基于風電消納電熱儲能裝置熱流固耦合研究

張雪平，齊鳳升，邢作霞，單建標，李寶寬

（1東北大學冶金學院，遼寧 沈陽 110819；2沈陽工業大學，遼寧 沈陽 110870）

風能由于其具有分布地區廣泛、不造成環境污染、資源蘊藏量豐富等優勢，現已成為可再生能源發展的主力軍。然而，截至目前風電的發展還有諸多問題有待解決，如風電并網、消納問題[1]依然嚴重，尤其是在三北地區，棄風現象較為普遍。我國電力系統負荷谷峰比例達0.4[2]以上，電力谷峰問題較為突出。在供暖地區，電蓄熱技術是解決風電并網問題，提高電網調峰能力[3]的有效途徑之一。系統利用低谷時期的電能或者無法并網消納的棄風電通過電加熱裝置將熱量儲存于蓄熱材料中，當需要時對熱用戶進行供熱。

國內外針對固體電蓄熱技術做了大量研究，尤其是在蓄熱體經濟性分析、蓄熱材料的選擇以及影響蓄熱性能因素等方面進行分析[4-5]。KHARE等[6]對多種高溫固體電蓄熱材料進行了物性參數和經濟性分析，發現鎂磚以及高鋁混凝土成本最低、熱容量較高。哈爾濱工業大學的廖晉[7]深入研究了固體電蓄熱技術理論及其裝置的熱工特性，提出了利用固體電蓄熱裝置消納過剩風電的方法。胡思科、周林林等[8]將蓄熱體模型簡化為二維模型，分析了放熱孔形及孔數對蓄放熱特性的影響，給定蓄熱體固體材料溫度作為放熱熱源，研究結果為蓄熱設備的涉及提供一定的參考依據。以往針對蓄熱裝置蓄放熱過程的研究運用的方法或者將蓄熱體簡化，選取其中包含單個換熱通孔進行二維模擬，或者通過運用傳熱學知識計算得到蓄熱材料與換熱介質邊界處（空氣、水等）熱流密度，將流固邊界設為熱流邊界條件。兩種方法只能簡單說明蓄熱體溫度場的變化規律，建立的傳熱模型與實際相差加大，不能全面的說明問題。裝置內傳熱過程較為復雜，涉及輻射、對流以及導熱傳熱，且涉及固體區域與流體區域的傳熱耦合，蓄熱體周期性地蓄熱和釋熱，溫度分布不均勻，企業投入運行的蓄熱裝置易出現開裂、坍塌等問題影響裝置工作。本文在上述研究的基礎上，對某企業運行中的固體電蓄熱裝置建立了流固耦合裝置內熱流固多場三維耦合傳熱數學模型，采用流固耦合傳熱的整場求解法[9]，詳細分析了裝置內蓄熱材料的溫度場及應變場，并在此基礎上對比了三種不同孔隙率、電熱絲排布方式對蓄熱裝置溫度分布均勻性的影響。

1 模 型

本文研究的固體電蓄熱裝置可用于棄風電消納與電力調峰系統中。裝置主要部件有蓄熱體、電熱絲、換熱管、循環風機等。蓄熱時，無法并網消納的風電通過加熱絲以輻射換熱的方式將電能轉換成熱量儲存于蓄熱體中；釋熱時，冷空氣流過通道與蓄熱材料表面直接接觸，蓄熱體蓄存的熱量通過對流換熱的方式將空氣加熱，被加熱的空氣通過熱交換器將熱量轉移給熱用戶。裝置工作原理示意圖如圖1所示。本文僅研究蓄熱過程。

圖1 固體電蓄熱裝置工作原理示意圖Fig.1 Schematic diagram of working principle of solid electric heat storage device

1.1 幾何及網格模型

整個蓄熱體組件由2028塊氧化鎂磚堆砌而成，磚型為異型磚，蓄熱體的總體積為8.355m3。蓄熱體內流體通道數為78個，其中36個通道內布置電熱絲，孔隙率為15%（保持蓄熱總體積不變，改變通道尺寸，通道總截面積與裝置截面積之比定義為孔隙率）。為了確定影響研究對象運行特性的主要因素，對研究模型進行合理簡化：①蓄熱材料的內部是均勻且連續填充的固態導熱介質，蓄熱單元看作一個長寬高固定的一個立體結構；②由于電熱絲的橫截面積相對于整個蓄熱單元而言非常小，故將其看作棒熱源處理。蓄熱體長2300mm，寬1500mm，高1874mm。為了使流體全部流入流體通道，蓄熱體上面以及兩個側面設有磚體擋板，用來阻擋繞流蓄熱體外部的流體。固體電蓄熱裝置被分為6個區域，每個區域均布有熱電偶監測溫度，幾何結構如圖2所示。

圖2 固體電蓄熱裝置結構圖Fig.2 Structural chart of solid electric heat storage device

圖3 固體電蓄熱裝置網格圖Fig.3 Grid of solid electric heat storage device

系統工作時，流體流過通道與蓄熱材料表面直接接觸，傳熱方式為流-固耦合對流傳熱，在結構化網格劃分的過程中，將流體區域和固體區域分離開來，流體區域與固體區域網格節點一一對應，數據通過流固耦合界面直接傳遞。考慮計算速度時間因素，確定網格數為106萬，如圖3所示。

1.2 控制方程

在整個裝置內，共存在三種不同形式的傳熱過程，即電熱絲與蓄熱磚的輻射換熱，蓄熱磚內部的導熱以及蓄熱磚與空氣的流固耦合對流換熱。經計算確定，空氣在通道內流動屬于湍流流動，近壁面處，采用壁面函數法[10]處理邊界層流動，湍流模型選用標準κ-ε湍流模型[11]。DO輻射模型[12]考慮了所有光學深度區間的輻射，且占用計算機內存適中，本文輻射模型選用DO輻射模型。流體區域、固體區域以及界面處控制方程如下。

（1）流體區域控制方程質量守恒方程

能量守恒方程

動量守恒方程

其中，ρf為流體密度，v為流體時均速度，v'為脈動速度，h為流體焓，p為流體的靜壓，μ為湍流黏度，SE、Si分別為能量源項以及動量源項，在本研究中為0。

（2）固體區域只有熱量傳遞過程，其導熱微分控制方程

其中，ρ為密度，τ為時間，λ為導熱系數，c為比熱容，Ts為固體溫度，Φ為固體內部體積熱源。

（3）流固交界面上滿足能量連續性條件，即溫度、熱流密度、位移相等。熱流密度控制方程式為

其中，qf、qs分別為流-固交界面上流體側和固體側的熱流密度，k為對流換熱系數，Tf為靠近交界面處流體溫度，Tm為交界面處溫度，n為流-固交界面法向量。

1.3 應變/應力理論

物質被加熱或冷卻時會產生熱脹冷縮，即發生“熱變形”，變形量和材料的膨脹系數有關。而由于存在約束，熱變形不能自由發生。這種因溫度變化而產生的應力稱之為熱應力。對蓄熱體作熱應力應變分析計算遵循材料力學的熱彈性理論。應變與節點位移矢量關系方程為

應力應變關系方程

式中，{σ}為應力矢量，[D]為彈性矩陣，{ε}為總應變矢量，{ε0}為熱應變矢量，[B]為節點應變-位移矩陣，{u}為節點位移矢量。

由式（6）～（7）可得：

熱應變矢量方程為

其中，T為單元體當前溫度（K），T0為初始溫度（K），αx、αy、αz分別為三個方向的熱膨脹系數（1/K）。考慮到蓄熱磚體在其溫度變化范圍內熱膨脹系數變化較小，在計算時將其設為一常值。式（9）可簡化為

則應變總矢量為

引起蓄熱體材料產生變形的原因主要有兩個，溫度分布不均以及流場對材料界面處的壓強。根據式（11）聯立平衡方程、幾何方程即可求解出對應的熱變形以及熱應力。

1.4 均勻性評價指標

基于不同工況下蓄熱裝置溫度場分布規律相似，采用相對標準偏差指標[13]衡量溫度場的均勻程度，其計算公式為

式中，CV為相對標準偏差，取0～1,越接近0表示溫度均勻性越好，n為數據點個數，Ti為某點溫度值，為數據平均值。

1.5 邊界條件和計算方法

邊界條件根據裝置實際運行過程中操作數據給出。蓄熱時，加熱絲總功率為270 kW，為了防止加熱絲表面過熱以及風機不斷啟停對設備影響，向裝置內通入少量空氣。入口為速度入口，空氣流速為0.01 m/s，溫度根據現場實際測試溫度編譯UDF導入確定。出口設為壓力出口，壁面設為絕熱無滑移邊界條件。流體與固體表面的交界處設為耦合壁面。初始蓄熱體內溫度均勻為482 K，總蓄熱時間為24000 s。氧化鎂磚熱物理性質如表1所示。

表1 氧化鎂磚的熱物理性質Table1 The thermal physical properties of magnesia brick

2 模擬結果與討論

2.1 模型驗證

為了驗證所建三維模型的正確性，提取6個熱電偶處的計算溫度與監測進行對比，如圖4所示。模擬初始條件根據企業提供的數據設置，6個熱電偶初始溫度有所不同，取裝置內平均溫度作為模擬初始溫度。誤差產生的來源主要是初始溫度的不同，熱電偶測得的溫度為局部區域溫度，而計算初始化取裝置內平均溫度為482 K。4號熱電偶的初始溫度為442 K，與裝置初始平均溫度最為接近，計算過程中最大為8%。各點的實驗溫度與計算溫度的變化趨勢基本一致，溫度上升速率也差別不大。

2.2 溫度場及應變場分析

圖5為蓄熱體在加熱過程中不同時間段（7200s、14400s、21600s以及24000s）溫度云圖（z=0.15m）。可以看出，隨著蓄熱過程的進行，布置加熱絲的通道表面由于直接接受電熱絲輻射，溫度最先升高，周圍磚體通過導熱傳熱，溫度逐漸降低，傳熱速率主要取決于蓄熱材料的導熱系數以及傳熱溫差。布置加熱絲的通道表面溫度最高可達932K，上部隔板處溫度上升較慢，無明顯升高。

圖4 計算值與現場數據對比Fig.4 Comparison between calculated values and measurement date

將模擬得到的固體溫度導入到Workbench熱應力分析模塊，將裝置底部以及四周設置為固定約束，利用等效應力以及總膨脹量得到變形、應力云圖分布如圖6、圖7所示。蓄熱體受熱向四周膨脹，加上底部約束，從下端到上端，熱膨脹量量呈拋物線型分布，中部中心區域熱膨脹量量最大。約束部位蓄熱體受熱無法自由膨脹，產生的熱應力最大，在這些區域磚體易出現斷裂現象。圖8給出了蓄熱體平均溫度以及膨脹量隨時間變化的曲線。蓄熱體平均溫度以及最大膨脹量均隨著時間呈線性上升。蓄熱結束后，蓄熱體平均溫度達760K，最大膨脹量為0.015m，與蓄熱體原寬相比，最大變形率最大為1%，影響可以忽略。

圖5 蓄熱體（z=0.15m截面）加熱過程中溫度分布Fig.5 Temperature distribution of heat storage (z=0.15m section) during heating process

圖6 蓄熱體變形量分布Fig.6 Distribution of deformation in heat storage

圖7 蓄熱體等效應力分布Fig.7 Distribution of equivalent stress

2.3 電熱絲個數對蓄熱溫度均勻性影響

圖8 蓄熱體最大變形量、溫度隨時間變化曲線Fig.8 Distribution of deformation and maximum stress during heating process

圖9 相對標準偏差隨時間變化曲線Fig.9 Curve of relative standard deviation during heating process

在蓄熱過程中，電熱絲的分布密度及單個功率對蓄熱溫度均勻性產生影響。根據傅里葉導熱定律，熱流量與溫差成正比，而與距離成反比。當單個電熱絲的功率不變時，電熱絲分布密度越大，溫度均勻性越好。當電熱絲分布密度不變時，電熱絲功率越大，則溫度均勻性較差。在保持電熱絲總加熱功率不變情況下，合理選擇電熱絲的個數可獲得較好蓄熱溫度均勻性。圖9為不同電熱絲個數下，相對標準偏差隨時間的變化曲線。三種情況下，在蓄熱初期，蓄熱體內部溫度梯度較大，相對標準偏差增大速度較快，后期電熱絲為18的工況出現緩慢下降的現象，最終趨于不變。電熱絲個數為72的工況下，蓄熱體均勻性程度最大，個數為18的工況均勻性較弱。

為進一步明確電熱絲個數對蓄熱溫度均勻性的影響，圖10給出了同一截面下，蓄熱體溫度分布隨電熱絲個數變化的云圖。電熱絲個數為18的工況下，高溫區與低溫區區別明顯，靠近電熱絲周圍的蓄熱磚溫度較高，而隔板處以及下部通道處，溫度呈波峰波谷狀分布，這是由于遠離電熱絲，蓄熱體內部通過導熱傳播熱量，而又受周圍左右兩個電熱絲的影響，距離越遠，低溫區越大。隨著電熱絲個數增多，上部低溫區面積變少，下部溫度趨于均勻，最高溫由1149K降為858K，最低溫由524K升為565K。電熱絲個數為72的工況下蓄熱溫度均勻性程度最好。

圖10 不同電熱絲個數在 同一截面處溫度分布云圖Fig.10 Temperature distribution at same section with number of electric heating wire (a) 18,(b) 36,(c) 72

2.4 孔隙率對蓄熱溫度均勻性影響

圖11為不同孔隙率下，相對標準偏差隨時間的變化。在蓄熱初期（1800s內），電熱絲溫度與蓄熱磚溫度相差較大，輻射熱流量較大，而蓄熱體內部導熱熱流傳播速率取決于材料的熱擴散系數，導熱需要一定的時間，蓄熱體內部溫度梯度較大，故相對標準偏差增長速度較快。大約加熱7000s之后，三種工況下的標準偏差基本不再發生變化，孔隙率為20%的工況相對標準偏差有所降低。根據此評價指標，孔隙率為20%的工況蓄熱溫度均勻性較好，而15%的工況蓄熱溫度均勻性相對較差。

圖11 相對標準偏差隨時間變化曲線Fig.11 D Curve of relative standard deviation during heating process

圖12所示為同一位置處，孔隙率為15%、20%以及25%的溫度分布云圖。孔隙率為25%的工況，其最高溫與最低溫相差最小，然而其低溫區面積最大，整體來看，其內部溫度梯度較大。孔隙率為20%的工況其最高溫僅為914K,而其溫度分布與15%的工況相似，溫度分布均勻性較其他兩個工況好。

3 結 論

文章以某企業實際運行的固體電蓄熱裝置為原型建立了三維流動與傳熱的數學模型，采用整場流固耦合法，將流固界面處難以確定的熱流邊界條件轉化為內部邊界，詳細分析了裝置蓄熱過程中內部溫度分布均勻性以及變形量，結論如下。

（1）將模擬結果與實驗數據對比，誤差產生的來源主要是初始溫度的不同，各點的實驗溫度與計算溫度的變化趨勢基本一致，溫度上升速率也差別不大。采用流固耦合傳熱模型較為精確。

（2）在蓄熱的過程中，蓄熱體平均溫度近似呈線性上升，最終達760K左右，蓄熱量為蓄熱體最大膨脹量為0.015m，最大應力出現在擋板處，與蓄熱體原寬相比，最大變形率最大為1%，影響可以忽略。

（3）運用相對標準偏差評價蓄熱溫度均勻性，并結合溫度云圖，得到電熱絲個數為72、孔隙率為20%的蓄熱體溫度分布均勻性較好，此結果為裝置結構優化提供參考。