基于風(fēng)電消納電熱儲能裝置熱流固耦合研究

張雪平，齊鳳升，邢作霞，單建標(biāo)，李寶寬

（1東北大學(xué)冶金學(xué)院，遼寧 沈陽 110819；2沈陽工業(yè)大學(xué)，遼寧 沈陽 110870）

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

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

1 模 型

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

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

1.1 幾何及網(wǎng)格模型

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

圖2 固體電蓄熱裝置結(jié)構(gòu)圖Fig.2 Structural chart of solid electric heat storage device

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

系統(tǒng)工作時，流體流過通道與蓄熱材料表面直接接觸，傳熱方式為流-固耦合對流傳熱，在結(jié)構(gòu)化網(wǎng)格劃分的過程中，將流體區(qū)域和固體區(qū)域分離開來，流體區(qū)域與固體區(qū)域網(wǎng)格節(jié)點一一對應(yīng)，數(shù)據(jù)通過流固耦合界面直接傳遞?？紤]計算速度時間因素，確定網(wǎng)格數(shù)為106萬，如圖3所示。

1.2 控制方程

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

（1）流體區(qū)域控制方程質(zhì)量守恒方程

能量守恒方程

動量守恒方程

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

（2）固體區(qū)域只有熱量傳遞過程，其導(dǎo)熱微分控制方程

其中，ρ為密度，τ為時間，λ為導(dǎo)熱系數(shù)，c為比熱容，Ts為固體溫度，Φ為固體內(nèi)部體積熱源。

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

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

1.3 應(yīng)變/應(yīng)力理論

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

應(yīng)力應(yīng)變關(guān)系方程

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

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

熱應(yīng)變矢量方程為

其中，T為單元體當(dāng)前溫度（K），T0為初始溫度（K），αx、αy、αz分別為三個方向的熱膨脹系數(shù)（1/K）?？紤]到蓄熱磚體在其溫度變化范圍內(nèi)熱膨脹系數(shù)變化較小，在計算時將其設(shè)為一常值。式（9）可簡化為

則應(yīng)變總矢量為

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

1.4 均勻性評價指標(biāo)

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

式中，CV為相對標(biāo)準(zhǔn)偏差，取0～1,越接近0表示溫度均勻性越好，n為數(shù)據(jù)點個數(shù)，Ti為某點溫度值，為數(shù)據(jù)平均值。

1.5 邊界條件和計算方法

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

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

2 模擬結(jié)果與討論

2.1 模型驗證

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

2.2 溫度場及應(yīng)變場分析

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

圖4 計算值與現(xiàn)場數(shù)據(jù)對比Fig.4 Comparison between calculated values and measurement date

將模擬得到的固體溫度導(dǎo)入到Workbench熱應(yīng)力分析模塊，將裝置底部以及四周設(shè)置為固定約束，利用等效應(yīng)力以及總膨脹量得到變形、應(yīng)力云圖分布如圖6、圖7所示。蓄熱體受熱向四周膨脹，加上底部約束，從下端到上端，熱膨脹量量呈拋物線型分布，中部中心區(qū)域熱膨脹量量最大。約束部位蓄熱體受熱無法自由膨脹，產(chǎn)生的熱應(yīng)力最大，在這些區(qū)域磚體易出現(xiàn)斷裂現(xiàn)象。圖8給出了蓄熱體平均溫度以及膨脹量隨時間變化的曲線。蓄熱體平均溫度以及最大膨脹量均隨著時間呈線性上升。蓄熱結(jié)束后，蓄熱體平均溫度達(dá)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 蓄熱體等效應(yīng)力分布Fig.7 Distribution of equivalent stress

2.3 電熱絲個數(shù)對蓄熱溫度均勻性影響

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

圖9 相對標(biāo)準(zhǔn)偏差隨時間變化曲線Fig.9 Curve of relative standard deviation during heating process

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

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

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

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

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

圖11 相對標(biāo)準(zhǔn)偏差隨時間變化曲線Fig.11 D Curve of relative standard deviation during heating process

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

3 結(jié) 論

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

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

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

（3）運用相對標(biāo)準(zhǔn)偏差評價蓄熱溫度均勻性，并結(jié)合溫度云圖，得到電熱絲個數(shù)為72、孔隙率為20%的蓄熱體溫度分布均勻性較好，此結(jié)果為裝置結(jié)構(gòu)優(yōu)化提供參考。