間歇熱流下管式吸熱器溫度場及熱應力場分析

王富強，帥 永，談和平，林日億

(1.中國石油大學(華東)儲運與建筑工程學院，266580山東青島;2.哈爾濱工業(yè)大學能源科學與工程學院，150001哈爾濱)

太陽能熱利用是可再生能源發(fā)展的有效途徑之一.目前國際上太陽能熱發(fā)電技術(shù)主要有槽式、塔式和碟式［1］3種形式，其中塔式和碟式太陽能熱發(fā)電技術(shù)仍處于試驗和示范階段，而槽式太陽能熱發(fā)電技術(shù)已成為成熟的商業(yè)化技術(shù).槽式聚光器利用拋物線的幾何光學特性將低輻射熱流密度的太陽光匯聚后，形成高匯聚熱流.安置在槽式聚光器焦平面處的管式吸熱器吸收高匯聚熱流并通過導熱、對流等耦合換熱的方式將高匯聚熱流轉(zhuǎn)換為熱能.

由于太陽能一方面受到晝夜、季節(jié)和地理緯度等規(guī)律性因素的影響，另一方面又受到晴朗、陰云和雨雪天氣等隨機因素的制約，因而太陽能的熱利用具有顯著的間斷性和不均勻性，吸熱器容易承受高的熱應力，引發(fā)玻璃管罩的破裂、管式吸熱器的彎曲并最終失效［2］，如:墨西哥國立大學太陽能熱利用電站在運行和試驗中，多次發(fā)生不銹鋼管式吸熱器的大彎曲變形并引發(fā)玻璃罩破裂的事故［3～5］.

管式吸熱器熱應力場的研究可以為提高吸熱器的可靠性提供理論指導，Kim等［6］采用有限體積與有限元聯(lián)合計算法對管式吸熱器的熱變形進行了分析;王志峰等［7］采用商業(yè)軟件Ansys對管式吸熱器密封節(jié)的殘余熱應力進行了計算;王富強等［8］采用順序解耦的計算方法對管式吸熱器的熱應變進行分析，指出采用偏心管式吸熱器替代同心管式吸熱器可以有效地降低吸熱器熱應變.

在對吸熱器的溫度場及熱應力場分析時，目前大多數(shù)學者都假定入射太陽能熱流與吸熱器表面呈垂直關(guān)系［9～12］，沒有考慮到太陽能熱流入射方向與吸熱器表面之間夾角關(guān)系.本文采用命令流編程方式，在對管式吸熱器進行網(wǎng)格劃分之后，計算吸熱器外表面每一單元網(wǎng)格面的法線方向與太陽能熱流入射方向之間的余弦關(guān)系與匯聚熱流的乘積，并將其作為溫度場分析求解的邊界條件.

1 幾何模型及模擬方法

1.1 幾何模型

圖1 管式吸熱器網(wǎng)格劃分

1.2 計算方法

與 Submodeling方法相比，Coupled-field element方法在溫度場和應力場分析過程可以采用相同的網(wǎng)格并同時對溫度場和應力場方程矩陣進行計算求解.因此，Coupled-field element方法求解過程更簡便，更適于含有函數(shù)邊界條件下的命令流編程分析.在熱應力場分析過程中沒有采用任何的外部約束，因此，管式吸熱器的熱應力是由于吸熱器管壁上不均勻溫度分布和吸熱器自身幾何形狀的約束所產(chǎn)生［13］.

為了驗證計算方法的可靠性，在采用相同的幾何模型尺寸及邊界條件下，將本文計算得到的徑向應力及切向應力沿長度變化曲線與文獻［14］進行了對比驗證.驗證算例為陶瓷管式換熱器在管壁內(nèi)側(cè)施加均勻?qū)α鲹Q熱系數(shù)及非均勻?qū)α鲹Q熱系數(shù)兩種邊界條件下的熱應力場分析，其中管壁外側(cè)施加均勻熱流密度(800 W/m2)的邊界條件.由圖2可知，本文計算得到的徑向應力及切向應力沿長度方向變化曲線與文獻［14］吻合良好.

圖2 本文計算結(jié)果與文獻［14］的對比驗證

1.3 邊界條件

在數(shù)值模擬過程中，假設(shè)管壁材料物性不隨溫度變化，分析不銹鋼管式吸熱器的溫度場以及熱應力場的分布.

1.3.1 吸熱器內(nèi)壁面

由于在進行熱應力場分析過程中，管式吸熱器內(nèi)部流體不具有剛度矩陣，無法進行熱應力場分析，因此采用對管式吸熱器內(nèi)壁面施加一恒定對流換熱系數(shù)方法進行分析.管式吸熱器內(nèi)壁面施加邊界條件如下:

1.3.2 吸熱器外壁面

考慮到太陽能熱流入射方向與吸熱器表面之間夾角關(guān)系，每一個單元網(wǎng)格上的太陽能熱流入射方向與單元網(wǎng)格法線方向之間角度變化關(guān)系如圖3所示.

圖3 單元網(wǎng)格上熱流入射方向與網(wǎng)格法線之間夾角關(guān)系

由文獻［6］可知，吸熱器外壁面上每一個單元網(wǎng)格面受到熱流密度為

A：通過這幾年的實踐，我們認為，VOCs治理應該從源頭入手，對過程進行控制，并輔以末端治理。源頭，是我們治理的重點，目前我們使用的是無醇潤版液，因為不用酒精了，整個生產(chǎn)車間的空氣質(zhì)量大大改善，幾乎聞不到刺鼻的氣味。過程控制，主要是做好各項管理，使VOCs排放源處于管控之中。末端治理，即VOCs的收集系統(tǒng)，這也是非常關(guān)鍵的一環(huán)。

其中，θi為網(wǎng)格節(jié)點i與水平面之間的夾角，Esun為太陽輻照度(Esun=1000W/m2)，Es和Er分別為大氣中的散射輻射和反射輻射.由于大氣中散射輻射和反射輻射相對較小且影響因素很多［15］，本文計算中忽略不計.α為吸熱器表面吸收系數(shù).由于管式吸熱器的表面具有高吸收率涂層(可高達0.97)及玻璃罩抑制輻射損失的作用［15］，本文選擇理想性吸收表面進行計算，即α=1.0.

考慮到管式吸熱器受到間歇太陽能熱流的特點，對管式吸熱器外壁面下表面施加一移動熱源.移動熱源的大小也同樣考慮到每一個單元網(wǎng)格上的太陽能熱流入射方向與單元網(wǎng)格法線方向之間角度變化關(guān)系，移動熱源的大小沿管壁圓周方向隨著夾角變化而不斷變化.移動熱源位置沿吸熱器長度方向隨時間步長不斷變化來模擬間歇熱流的影響，如由陰天突變?yōu)榍缣?、日升日落?由文獻［6］可知:

移動熱源為

其中Cr為幾何匯聚比;

移動熱源位置為

其中:為t為時間步長，t=0.1s;S為時間步數(shù);T為總的計算時間，T=t×St;St為總的計算時間步數(shù).由文獻［13-14］可知，熱應力產(chǎn)生的主要原因是溫度場的突變，均勻、線性變化溫度場均不會引起熱應力.因此，本文選用總的計算時間為10s來模擬間歇熱流的作用，即St=100.

1.3.3 吸熱器端面

為了使溫度場變化對熱應力場的影響這一研究更清晰，在熱應力場分析過程中管式吸熱器兩端面沒有采用任何的外部約束.

2 計算結(jié)果與分析

2.1 不銹鋼管式吸熱器的溫度場分布

圖4為不同時間步下不銹鋼管式吸熱器外壁面在y=0位置處溫度場沿長度方向變化曲線.由于移動熱源的加載是隨時間沿管式吸熱器長度方向不斷前進的，因此在所有的時間步下，最高溫度出現(xiàn)在z=0位置處，沿著長度方向溫度不斷遞減.隨著時間的推移，吸熱器的最高溫度不斷上升，在時間步s=100時，不銹鋼管式吸熱器最高溫度已經(jīng)達到了572K.這是由于隨著時間的增長，受到高匯聚比太陽能熱流的面積不斷增大.由文獻［13］可知，熱應力場與溫度的二階導數(shù)呈函數(shù)關(guān)系，因此本文將溫度場沿長度方向進行二次求導來揭示熱應力場與溫度場的變化關(guān)系.

圖4 不同時間步下管式吸熱器溫度場沿長度變化

圖5所示為不同時間步下不銹鋼管式吸熱器外壁面在y=0位置處溫度場的二階導數(shù)沿長度方向變化曲線.

圖5 不同時間步下管式吸熱器溫度二階導數(shù)沿長度變化

由圖5可知，不同時間步下溫度二階導數(shù)曲線形狀基本一致，不同時間步下的溫度二階導數(shù)曲線的波峰與波谷的數(shù)值也相差很小.不同時間步下溫度二階導數(shù)曲線的波峰與波谷值出現(xiàn)的位置隨著時間步增長而沿著管式吸熱器長度方向推移，但是溫度二階導數(shù)的波峰與波谷值并不是出現(xiàn)在移動熱源加載位置邊界處，而是在移動熱源加載位置邊界稍靠前位置處.

2.2 不銹鋼管式吸熱器的應力場分布

圖6為不同時間步下不銹鋼管式吸熱器外壁面在y=0位置處軸向熱應力場沿長度方向變化曲線.由圖6可知，不同時間步下軸向應力沿管式吸熱器長度方向變化趨勢是一致的:軸向應力沿著管式吸熱器長度方向首先是不斷下降，一直達到曲線的波谷值;然后軸向應力沿著長度的方向增加呈現(xiàn)出線性增長到達軸向應力曲線的第二拐點處，此后軸向應力急速增長并趨近于零.對比軸向應力與溫度二階導數(shù)沿管式吸熱器長度方向變化曲線可以發(fā)現(xiàn)，在軸向應力曲線的拐點處均是溫度二階導數(shù)曲線拐點處.以時間步S=20為例:溫度二階導數(shù)曲線在L=0.18 m處出現(xiàn)了曲線的拐點，溫度二階導數(shù)數(shù)值由負值變?yōu)榱?，而軸向應力曲線在L=0.18 m處也出現(xiàn)了曲線的拐點，此拐點為軸向應力曲線波谷值;溫度二階導數(shù)曲線在L=0.25 m處出現(xiàn)了曲線的拐點，而軸向應力曲線在L=0.25 m處也出現(xiàn)了曲線的拐點，經(jīng)歷此拐點后，軸向應力曲線繼續(xù)上升并趨近于零;溫度二階導數(shù)曲線在L=0.42 m處出現(xiàn)了曲線的拐點，即為溫度二階導數(shù)曲線的峰值，而軸向應力曲線在L=0.42 m處也出現(xiàn)了曲線的第3個拐點，在此拐點處，軸向應力曲線數(shù)值變?yōu)?.

圖6 不同時間步下管式吸熱器軸向應力沿長度變化

對比圖5與圖6可知，雖然不同時間步下的溫度二階導數(shù)曲線的波峰波谷數(shù)值大小相差不大，但是不同時間步下的軸向應力的大小卻有著比較大的變化.這是由于不同時間步下，溫度二階導數(shù)的拐點出現(xiàn)位置不同而引起的.由彈性力學理論可知，在不考慮材料物性隨溫度場變化時，對管式吸熱器不同位置施加相同的載荷而引起管式吸熱器變形大小也不相同，因此產(chǎn)生的應力的數(shù)值大小也不相同，因此發(fā)生了不同時間步下相同位置的軸向應力值變化較大的現(xiàn)象.

圖7為不同時間步下不銹鋼管式吸熱器外壁面下表面在y=0位置處徑向熱應力場沿長度方向變化曲線.不同時間步下的徑向應力數(shù)值最大值均出現(xiàn)在z=0的端面處，徑向應力表現(xiàn)為壓應力.這是由于在z=0的端面處出現(xiàn)溫度峰值，且在此端面處溫度沿半徑方向變化最為劇烈.不同時間步下徑向應力沿管式吸熱器長度方向變化曲線形狀基本一致.除z=0的端面處外，不同時間步下的徑向應力變化曲線的波峰與波谷的數(shù)值也相差很小.同溫度二階導數(shù)沿長度方向變化曲線一樣，不同時間步下徑向應力變化曲線的波峰與波谷值出現(xiàn)的位置隨著時間步增長而沿著管式吸熱器長度方向推移.與軸向應力相比，相同時間步下的最大徑向應力僅為軸向應力的1/10左右.

圖7 不同時間步下管式吸熱器徑向應力沿長度變化

圖8為不同時間步下不銹鋼管式吸熱器外壁面在y=0位置處切向熱應力場沿長度方向變化曲線.同徑向應力沿長度方向變化曲線一樣，切向應力最大值均出現(xiàn)在z=0的端面處且均表現(xiàn)為壓應力.切向應力沿長度變化曲線的變化趨勢與徑向應力變化趨勢基本上是一致的.對比切向應力與溫度二階導數(shù)沿管式吸熱器長度方向變化曲線可以發(fā)現(xiàn)，切向應力曲線發(fā)生拐點處均是溫度二階導數(shù)曲線拐點發(fā)生處.

圖9為不同時間步下不銹鋼管式吸熱器外壁面下表面在y=0位置處等效應力Von-Mises熱應力場沿長度方向變化曲線.不同時間步下等效應力Von-Mises熱應力沿管式吸熱器長度方向變化趨勢是一致的.由等效應力Von-Mises應力理論可知［13-14］，軸向應力和切向應力對 Von-Mises熱應力場影響最大，徑向應力影響最小.等效應力Von-Mises熱應力沿長度方向變化曲線的拐點處均是溫度二階導數(shù)沿長度方向變化曲線的拐點發(fā)生處，這與文獻［13-14］研究得到結(jié)論相同.

圖8 不同時間步下管式吸熱器切向應力沿長度變化

圖9 不同時間步下管式吸熱器等效應力沿長度方向變化

3 結(jié)論

1)本文采用命令流編程方式，計算吸熱器外表面每一單元網(wǎng)格面的法線方向與太陽能熱流入射方向之間的余弦關(guān)系，并以此余弦關(guān)系作為溫度場分析求解的邊界條件.采用 Coupled-field element方法對管式吸熱器的溫度場和應力場進行了耦合求解計算.

2)研究了管式吸熱器在不均勻熱流密度以及間歇太陽能熱流邊界條件下的溫度場分布以及熱應力場分布.

3)研究了溫度場沿長度方向二階導數(shù)變化規(guī)律，發(fā)現(xiàn)熱應力曲線拐點處均是溫度二階導數(shù)曲線拐點發(fā)生處.

［1］方嘉賓，魏進家，董訓偉，等.腔式太陽能吸熱器熱性能的模擬計算［J］.工程熱物理學報，2009，30(3):428－432.

［2］LATA J M，RODRIGUEZ M A，de LARA M A.High flux central receivers of molten salts for the new generation of commercial stand-alone solar power plants［J］.ASME Journal of Solar Energy Engineering，2008，130(2):0210021－0210025.

［3］ ALMANZA R F，LENTZ A，JIMENEZ G.Receiver behavior in direct steam generation with parabolic toughs［J］.Solar Energy，1997，61:275 －278.

［4］ ALMANZA R F，F(xiàn)LORES V C，LENTZ A，et al.Compound wall receiver for DSG in parabolic troughs［C］//Proceedings of the 10th international symposium of solar thermal.Sydney:Solar PACES，2002:131 -135.

［5］FLORES V C，ALMANZA R F.Behavior of compound wall copper-steel receiver with stratified two-phase flow regimen in transient states when solar irradiance is arriving on one side of receiver［J］.Solar Energy，2004，76:195-198.

［6］ KIM D C，CHOI Y D.Analysis of conduction-natural convection conjugate heat transfer in the gap between concentric cylinders undersolarirradiation ［J］.International Journal of Thermal Science，2009，48(6):1247-1258.

［7］LEI Dongqiang，WANG Zhifeng，LI Jian.The analysis of residual stress in glass-to-metal seals for solar receiver tube［J］.Materials and Design，2010，31(4):1813 -1820.

［8］ WANG Fuqiang，SHUAI Yong，YUAN Yuan，et al.Thermal stress analysis of eccentric tube receiver using concentrated solar radiation［J］.Solar Energy，2010，85(10):1809-1815.

［9］劉斌，吳玉庭，楊征，等.螺旋盤管腔式熔鹽太陽能吸熱器的分析優(yōu)化和研制［J］.工程熱物理學報，2009，30(6):995－998.

［10］KUMAR N S，REDDY K S.Comparison of receivers for solar dish collector system ［J］.Energy Conversion and Management，2008，49:812 －819.

［11］YAPICI H S，ALBAYRAK B.Numerical solutions of conjugate heat transfer and thermal stresses in a circular pipe externally heated with concentrated solar irradiation heat flux ［J］.Energy Conversion and Management，2004，45:927 －937.

［12］KUMAR N S，REDDY K S.Thermal analysis of solar parabolic collector with porous disc receiver ［J］.Applied Energy，2009，86:1804 -1812.

［13］IFRAN M A，CHAPMAN W C.Thermal stresses in radiant tubes due to axial，circumferential and radial temperature distributions ［J］.Applied Thermal Engineering，2009，29:1913 -1920.

［14］ISLAMOGLU Y S.Finite element model for thermal analysis of ceramic heat exchanger tube under axial concentrated solar irradiation convective heat transfer coefficient［J］.Materials and Design，2004，25:479 －482.

［15］何梓年.太陽能熱利用［M］.合肥:中國科學技術(shù)大學出版社，2009:31－98.