李 斌， 陳 豐， 史良宵

（華北電力大學(xué) 能源動(dòng)力與機(jī)械工程學(xué)院，河北保定071003）

溫度場(chǎng)計(jì)算是汽包應(yīng)力計(jì)算和疲勞壽命分析的基礎(chǔ)，根據(jù)瞬態(tài)溫度場(chǎng)，可實(shí)時(shí)監(jiān)測(cè)汽包疲勞壽命損耗，保證機(jī)組安全經(jīng)濟(jì)運(yùn)行［1－5］.

傳統(tǒng)的汽包溫度場(chǎng)計(jì)算通常采用直接解法［6］，即在一定的初始條件和邊界條件下求解導(dǎo)熱微分方程，它需要已知求解區(qū)域的所有邊界條件.由于汽包內(nèi)有各種復(fù)雜的流動(dòng)和換熱過(guò)程，因而采用該方法時(shí)需要進(jìn)行很多假設(shè)，從而帶來(lái)較大的計(jì)算誤差.另外，直接解法最致命的缺陷是需要知道汽包內(nèi)壁與水和水蒸氣的對(duì)流傳熱系數(shù)，而汽包內(nèi)復(fù)雜的工況使得很難得到確切的傳熱系數(shù)，通常采用經(jīng)驗(yàn)數(shù)據(jù)，但會(huì)導(dǎo)致誤差［7－10］.此外還要已知熱流密度和流體溫度，盡管熱流密度和流體溫度都可以通過(guò)熱流密度計(jì)和熱電偶測(cè)得，然而該數(shù)據(jù)只適用于低壓情況下，不適用于工程實(shí)際高溫高壓或者不便布置傳感器等情況下瞬態(tài)傳熱的計(jì)算.

反演解法（Inverse Method）［11－12］采用控制容積法，在控制體上選擇一些離散的測(cè)點(diǎn)，采用熱電偶測(cè)溫，通過(guò)測(cè)得的溫度值實(shí)現(xiàn)整個(gè)溫度場(chǎng)的求解，即導(dǎo)熱反問(wèn)題（Inverse Problem of Heat Conduction）求解［13－14］.采用反演解法求解大型高溫高壓設(shè)備溫度場(chǎng)不需要進(jìn)行直接解法中的假設(shè)，僅需要知道高溫高壓設(shè)備外壁的溫度，即可快捷準(zhǔn)確地計(jì)算出該設(shè)備的溫度場(chǎng).計(jì)算汽包的溫度場(chǎng)時(shí)，汽包外壁的溫度可以通過(guò)在汽包外壁布置熱電偶進(jìn)行測(cè)量，比較準(zhǔn)確和方便.該方法克服了直接解法中的缺陷，不需要知道汽包內(nèi)壁的傳熱系數(shù)便可求解汽包溫度場(chǎng).

1 反演解法

應(yīng)用反演解法實(shí)現(xiàn)汽包瞬態(tài)溫度場(chǎng)的在線監(jiān)測(cè)，首先在汽包的外壁布置測(cè)點(diǎn)，安裝熱電偶測(cè)取外壁溫度，然后根據(jù)測(cè)點(diǎn)對(duì)汽包橫截面進(jìn)行網(wǎng)格劃分，最后利用控制容積法對(duì)導(dǎo)熱微分方程進(jìn)行離散，將局部熱傳導(dǎo)方程轉(zhuǎn)化為常微分方程，由外向內(nèi)逐次內(nèi)推，求解得到汽包橫截面的溫度場(chǎng).

汽包橫截面示意圖見(jiàn)圖1，圖中kw、ks分別表示汽包內(nèi)飽和水.飽和水蒸氣與汽包內(nèi)壁的對(duì)流傳熱系數(shù)，W／（m2·K）.

圖1 汽包橫截面示意圖Fig.1 Sectional view of the boiler drum

汽包模型相關(guān)參數(shù)如下：內(nèi)徑r1為0.400m，外徑r5為0.426m；汽包材料的物性參數(shù)：導(dǎo)熱系數(shù)λ為36.0W／（m·K），密度ρ為7 850kg／m3，比熱容c為468J／（kg·K）.

1.1 網(wǎng)格劃分

根據(jù)汽包橫截面外壁熱電偶的布置情況，結(jié)合反演解法的思想對(duì)其進(jìn)行網(wǎng)格劃分，見(jiàn)圖2.圖中，r1和r5分別為內(nèi)徑和外徑，m；r2、r3和r4分別表示中間各層半徑，m.由于汽包橫截面的對(duì)稱性，只繪制右半部分.沿著半徑方向由內(nèi)向外劃分為三層，分別為內(nèi)層、中間層和外層，沿圓周方向劃分為13個(gè)節(jié)點(diǎn).

圖2 反演解法網(wǎng)格劃分示意圖Fig.2 Grid division of inverse method

1.2 反演解法

根據(jù)已知的外壁溫度和熱流密度，按照導(dǎo)熱問(wèn)題數(shù)值解法的思想，采用熱平衡法，對(duì)汽包外壁每個(gè)節(jié)點(diǎn)所代表的控制容積用傅里葉定律列出能量守恒表達(dá)式.外層節(jié)點(diǎn)32、節(jié)點(diǎn)33和節(jié)點(diǎn)34的能量守恒表達(dá)式如下：

式中：Δφ為容積角度變化量，rad；Δr為容積徑向變化量，m；Ti為節(jié)點(diǎn)i的溫度，℃；qi為節(jié)點(diǎn)i處的熱流密度，W／m2.

根據(jù)式（1）、式（2）和式（3）可以得出中間層節(jié)點(diǎn)19、節(jié)點(diǎn)20和節(jié)點(diǎn)21的溫度表達(dá)式：

同理，根據(jù)圖2，對(duì)中間層節(jié)點(diǎn)20列出能量守恒表達(dá)式，得到內(nèi)層節(jié)點(diǎn)7的溫度表達(dá)式：

將外層節(jié)點(diǎn)的溫度逐次反演求解得到內(nèi)層節(jié)點(diǎn)的溫度，改變外層節(jié)點(diǎn)位置可求解得到整個(gè)內(nèi)層節(jié)點(diǎn)的溫度，從而得到整個(gè)汽包橫截面的瞬態(tài)溫度場(chǎng).

2 Ansys分析驗(yàn)證

由于采用反演解法求解溫度場(chǎng)需要知道求解區(qū)域外邊界的溫度及其邊界換熱條件，因此采用Ansys數(shù)值模擬計(jì)算結(jié)果進(jìn)行相關(guān)驗(yàn)證.

根據(jù)導(dǎo)熱微分方程、求解區(qū)域的邊界條件以及初始條件，采用Ansys進(jìn)行數(shù)值求解，將求解得到的外邊界的溫度作為反演解法的已知條件，通過(guò)反演解法求解得到內(nèi)層節(jié)點(diǎn)溫度，再與Ansys數(shù)值模擬得到的內(nèi)層節(jié)點(diǎn)溫度進(jìn)行比較.

2.1 Ansys數(shù)值模擬數(shù)學(xué)模型的建立

導(dǎo)熱微分方程：

邊界條件：

初始溫度：

汽包內(nèi)飽和水（或飽和水蒸氣）的溫度隨時(shí)間的變化關(guān)系如下：

式中：t為時(shí)間，s；T∞為流體溫度，℃.

Ansys數(shù)值模擬的網(wǎng)格劃分如圖3所示.由于對(duì)象具有對(duì)稱性，故取其右半部分進(jìn)行分析.沿徑向劃分為5個(gè)節(jié)點(diǎn)，周向25個(gè)節(jié)點(diǎn)，共計(jì)96個(gè)單元，125個(gè)節(jié)點(diǎn).

圖3 Ansys數(shù)值模擬網(wǎng)格劃分Fig.3 Grid division of Ansys numerical simulation

2.2 數(shù)值模擬

按照上述模型對(duì)汽包橫截面瞬態(tài)溫度場(chǎng)進(jìn)行數(shù)值模擬，截面的初始溫度為70℃，模擬的時(shí)間步長(zhǎng)為10s，終止時(shí)刻為3 000s，得到汽包橫截面的節(jié)點(diǎn)溫度隨時(shí)間的變化曲線如圖4所示.圖中節(jié)點(diǎn)29、節(jié)點(diǎn)33和節(jié)點(diǎn)37為外層節(jié)點(diǎn)，節(jié)點(diǎn)3和節(jié)點(diǎn)11為內(nèi)層節(jié)點(diǎn).

2.3 結(jié)果驗(yàn)證

將上述Ansys數(shù)值模擬得到的外層節(jié)點(diǎn)溫度作為反演解法的已知條件（即測(cè)量外壁溫度），時(shí)間步長(zhǎng)同樣取10s，逐漸內(nèi)推得到相應(yīng)的內(nèi)層節(jié)點(diǎn)溫度，并與Ansys模擬所得的內(nèi)層節(jié)點(diǎn)溫度進(jìn)行比較驗(yàn)證.

圖4 Ansys數(shù)值模擬所得汽包橫截面的節(jié)點(diǎn)溫度隨時(shí)間的變化Fig.4 Variation of node temperature with time obtained by Ansys numerical simulation

由于采用反演解法求解溫度值時(shí)需用到溫度對(duì)時(shí)間的高階導(dǎo)數(shù)，且測(cè)量的外邊界的溫度對(duì)時(shí)間的變化非常敏感，在逐層推進(jìn)的過(guò)程中，溫度值隨時(shí)間的變化產(chǎn)生的誤差被逐層放大，因此，為減小誤差，在計(jì)算前對(duì)外壁絕熱層的溫度測(cè)量值及外壁溫度對(duì)時(shí)間的導(dǎo)數(shù)進(jìn)行修正，采用局部多項(xiàng)式方法對(duì)外層溫度值進(jìn)行空間和時(shí)間的平滑后再進(jìn)行計(jì)算［8］.

Ansys模擬所得外層節(jié)點(diǎn)溫度經(jīng)過(guò)相應(yīng)的平滑后，得到的外層節(jié)點(diǎn)溫度隨時(shí)間的變化曲線如圖5所示.

圖5 經(jīng)空間、時(shí)間平滑后的外層節(jié)點(diǎn)29、節(jié)點(diǎn)33和節(jié)點(diǎn)37的溫度Fig.5 Space－time averaged value of temperature at outer nodes 29，33and 37

根據(jù)圖1，以汽包橫截面幾何中心線作為汽水分界面，下部飽和水（上部飽和水蒸氣）與內(nèi)壁發(fā)生對(duì)流傳熱，各點(diǎn)處溫差相對(duì)比較小，因此在進(jìn)行反演解法驗(yàn)證時(shí)，中心線以下、中心線處以及中心線以上三部分分別采用節(jié)點(diǎn)29、節(jié)點(diǎn)33和節(jié)點(diǎn)37的溫度作為相應(yīng)的外層溫度.

反演解法的計(jì)算結(jié)果與Ansys數(shù)值模擬結(jié)果的對(duì)比如圖6所示，2種算法計(jì)算所得的t＝1 000s和t＝2 000s時(shí)部分節(jié)點(diǎn)的溫度值如表1和表2所示.由圖6、表1和表2可以看出，反演解法的計(jì)算結(jié)果與Ansys數(shù)值模擬結(jié)果具有較高的吻合度，說(shuō)明反演解法具有較高的精確度.

圖6 反演解法和Ansys模擬所得內(nèi)層節(jié)點(diǎn)3和節(jié)點(diǎn)11溫度值的對(duì)比Fig.6 Comparison of temperature changes at internal nodes 3and 11by inverse method and Ansys numerical simulation

表1 t＝1 000s時(shí)2種算法計(jì)算所得部分節(jié)點(diǎn)的溫度值Tab.1 Comparison of temperature changes at partial nodes obtained by two methods at t＝1 000s

表2 t＝2 000s時(shí)2種算法計(jì)算所得部分節(jié)點(diǎn)的溫度值Tab.2 Comparison of temperature changes at partial nodes obtained by two methods at t＝2 000s

3 無(wú)限長(zhǎng)圓柱理論解驗(yàn)證

3.1 模型的建立

半徑為R的實(shí)心圓柱，其材料的導(dǎo)熱系數(shù)為λ，熱擴(kuò)散率α為常數(shù)，初始溫度為T0，將其放在溫度為Tf并保持不變的流體中發(fā)生對(duì)流傳熱，流體與圓柱表面間的表面?zhèn)鳠嵯禂?shù)k為常數(shù).其模型簡(jiǎn)化示意圖如圖7所示.

3.2 理論解析解

圓柱的無(wú)量綱過(guò)余溫度解析解為

圖7 理論解模型示意圖Fig.7 Model for analytical solution

式中：FO為傅里葉數(shù)；η為半徑比；τ為時(shí)間，s；r為計(jì)算半徑，m；θ為過(guò)余溫度，K；θ0為初始時(shí)刻過(guò)余溫度，K；J0、J1分別為零階和一階的第一類貝塞爾（Bessel）函數(shù)；μn為超越方程的特征值；k為對(duì)流傳熱系數(shù)，W／（m2·K）；λ為材料導(dǎo)熱系數(shù)，W／（m·K）；

選取初始溫度為500℃，流體溫度為20℃，計(jì)算時(shí)間步長(zhǎng)為1s，對(duì)流傳熱系數(shù)為845W／（m2·K），根據(jù)理論解析解得到r＝r1和r＝r5處的溫度值，如圖8所示.

圖8 無(wú)限長(zhǎng)圓柱r＝r1和r＝r5處的理論解析解Fig.8 Analytical solution for infinitely long cylinder in the case of r＝r1and r＝r5

3.3 結(jié)果驗(yàn)證

將無(wú)限長(zhǎng)圓柱r＝r5（即外層）處的理論解析解作為“測(cè)量外層溫度”，用作反演解法的已知條件進(jìn)行反演計(jì)算.將求解得到的結(jié)果與無(wú)限長(zhǎng)圓柱r＝r1（即內(nèi)層）處的理論解析解進(jìn)行對(duì)比分析驗(yàn)證，如圖9所示.

相對(duì)于實(shí)際熱電偶測(cè)量的外層數(shù)據(jù)而言，理論解析解波動(dòng)或者誤差更小，因此為了驗(yàn)證應(yīng)用測(cè)量數(shù)據(jù)進(jìn)行反演解法計(jì)算的正確性，在外層理論解析解的基礎(chǔ)上增加一個(gè)隨機(jī)誤差－0.5～0.5K，再作為“測(cè)量外層溫度”進(jìn)行反演解法計(jì)算.帶擾動(dòng)情況下內(nèi)層理論解析解與反演解法結(jié)果的對(duì)比見(jiàn)圖10.

由圖9和圖10可知，反演解法結(jié)果與理論解析解很接近，吻合度高，同時(shí)外加擾動(dòng)情況下結(jié)果也比較接近，可以驗(yàn)證反演解法具有較高的精確度.

圖9 內(nèi)層反演解法結(jié)果與理論解析解的對(duì)比Fig.9 Comparison of temperature changes respectively obtained by inverse method and analytical solution

圖10 帶擾動(dòng)情況下內(nèi)層理論解析解與反演解法結(jié)果的對(duì)比Fig.1 0 Comparison of temperature changes respectively obtained by inverse method and analytical solution with disturbances

4 試驗(yàn)驗(yàn)證

為了驗(yàn)證反演解法的計(jì)算精度，在某小型鍋爐外壁安裝熱電偶進(jìn)行實(shí)際數(shù)據(jù)測(cè)量，將計(jì)算結(jié)果與測(cè)得的試驗(yàn)數(shù)據(jù)進(jìn)行比較，對(duì)該方法進(jìn)行驗(yàn)證.

選擇鍋爐的某一段工況，在鍋爐外壁選擇適當(dāng)?shù)奈恢貌贾孟鄳?yīng)的熱電偶，布置的外壁測(cè)點(diǎn)編號(hào)為35～39，內(nèi)壁測(cè)點(diǎn)編號(hào)為11，如圖2所示.

由于實(shí)際現(xiàn)場(chǎng)測(cè)量數(shù)據(jù)存在較大的波動(dòng)，首先將鍋爐外壁的測(cè)量數(shù)據(jù)按照反演解法的步驟進(jìn)行時(shí)間和空間的平滑處理，然后進(jìn)行計(jì)算，即可求解得到內(nèi)壁測(cè)點(diǎn)11的溫度隨時(shí)間的變化.再將計(jì)算結(jié)果與實(shí)際熱電偶測(cè)量結(jié)果進(jìn)行對(duì)比分析，如圖11所示.

由于實(shí)際運(yùn)行工況的復(fù)雜性，實(shí)際測(cè)量結(jié)果存在一定的波動(dòng)，另外熱電偶測(cè)量數(shù)據(jù)本身存在一定的誤差，但是由圖11可知，反演解法的計(jì)算結(jié)果與熱電偶的實(shí)際測(cè)量結(jié)果吻合度較高.

圖11 內(nèi)壁測(cè)點(diǎn)11溫度的計(jì)算結(jié)果與實(shí)測(cè)結(jié)果的對(duì)比分析Fig.1 1 Comparison of temperature changes at node 11obtained by inverse method and experimental test

5 結(jié) 論

（1）反演解法思路簡(jiǎn)單明了，求解過(guò)程不需要迭代，計(jì)算精度高，相對(duì)于直接解法而言，求解網(wǎng)格相對(duì)較少，但是精度卻一致.另外計(jì)算不需要已知內(nèi)壁換熱條件，只需知道外層節(jié)點(diǎn)的溫度即可反演得到內(nèi)層節(jié)點(diǎn)溫度，從而得到整個(gè)汽包的瞬態(tài)溫度場(chǎng)分布.同時(shí)也避免了打孔對(duì)壓力容器造成的設(shè)備壽命損耗，提高了設(shè)備的使用壽命.

（2）通過(guò)Ansys數(shù)值模擬，取汽包外壁邊界換熱條件為絕熱，從二維瞬態(tài)的角度來(lái)驗(yàn)證反演解法的計(jì)算結(jié)果，證明該方法計(jì)算精度高，二者具有很好的吻合度.

（3）采用無(wú)限長(zhǎng)圓柱理論解析解，圓柱體外壁邊界與流體對(duì)流傳熱，從一維瞬態(tài)的角度驗(yàn)證反演解法的適用范圍，同時(shí)驗(yàn)證了其計(jì)算的準(zhǔn)確性；另外在附加擾動(dòng)的情況下同樣得到很好的吻合度.

（4）通過(guò)與實(shí)際試驗(yàn)數(shù)據(jù)的對(duì)比分析，驗(yàn)證了反演解法具有較高的計(jì)算精度，與實(shí)際試驗(yàn)數(shù)據(jù)具有很好的吻合度，滿足工程應(yīng)用的要求.

［1］管德清，莫江春，張學(xué)綸，等.電站鍋爐汽包壽命在線監(jiān)測(cè)系統(tǒng)［J］.動(dòng)力工程，2002，22（6）：2044－2047.GUAN Deqing，MO Jiangchun，ZHANG Xuelun，et al.Life on－line monitoring system for boiler drum of power station［J］.Power Engineering，2002，22（6）：2044－2047.

［2］劉彤，徐鋼，龐力平，等.鍋爐爐內(nèi)承壓部件的蠕變分析及壽命計(jì)算［J］.動(dòng)力工程，2004，24（5）：631－635.LIU Tong，XU Gang，PANG Liping，et al.Creep analysis and life calculation of the pressure components inside boiler［J］.Power Engineering，2004，24（5）：631－635.

［3］李立人，陳瑋，盛建國(guó)，等.鍋爐受壓元件的高溫蠕變－疲勞壽命設(shè)計(jì)計(jì)算方法［J］.動(dòng)力工程，2009，29（5）：409－416.LI Liren，CHEN Wei，SHENG Jianguo，et al.Creepfatigue life design and calculation method for boiler pressure elements under elevated temperature［J］.Journal of Power Engineering，2009，29（5）：409－416.

［4］鄭心偉，孫瑜，王曉軍.增壓鍋爐汽包低周疲勞壽命計(jì)算方法研究［J］.熱能動(dòng)力工程，2010，25（2）：184－189.ZHENG Xinwei，SUN Yu，WANG Xiaojun.Study of the methods for calculating the low－cycle fatigue life of a supercharged boiler drum ［J］.Journal of Engineering for Thermal Energy and Power，2010，25（2）：184－189.

［5］張健.超臨界鍋爐爐外承壓部件的壽命分析及在線檢測(cè)［D］.北京：華北電力大學(xué)，2008.

［6］趙鐵成，沈月芬，朱國(guó)楨，等.電站鍋爐鍋筒溫度場(chǎng)計(jì)算——三維非穩(wěn)態(tài)變物性材料不均勻?qū)釂?wèn)題有限元分析［J］.中國(guó)電機(jī)工程學(xué)報(bào)，1997，17（4）：217－220.ZHAO Tiecheng，SHEN Yuefen，ZHU Guozhen，et al.Temperature field calculation of boiler drum—finite element analysis of 3－D unsteady variable character uneven material heat conduction problem［J］.Proceedings of the CSEE，1997，17（4）：217－220.

［7］HàO D N，THANH P X，LESNIC D.Determination of the heat transfer coefficients in transient heat conduction［J］.Inverse Problems，2013，29（9）：095020.

［8］TALER J.Determination of local heat transfer coefficient from the solution of the inverse heat conduction problem ［J］.Forschung im Ingenieurwesen，2007，71（2）：69－78.

［9］CEBULA A，TALER J.Determination of transient temperature and heat flux on the surface of a reactor control rod based on temperature measurements at the interior points ［J］.Applied Thermal Engineering，2014，63（1）：158－169.

［10］HU Wensen，LI Bin，CAO Zidong，et al.An inverse method for online stress monitoring and fatigue life analysis of boiler drums［J］.Journal of Chongqing University：English Edition，2009，8（2）：89－96.

［11］JAREMKIEWICZ M，TALER D，SOBOTA T.Measuring transient temperature of the medium in power engineering machines and installations［J］.Applied Thermal Engineering，2009，29（16）：3374－3379.

［12］TALER J.A new space marching method for solving inverse heat conditions problems ［J］.Forschung im Ingenieurwesen，1999，64（11）：296－306.

［13］TALER J，ZIMA W.Solution of inverse heat conduction problems using control volume approach［J］.International Journal of Heat and Mass Transfer，1999，42（6）：1123－1140.

［14］WANG Peng，LI Bin，WANG Songling.An On－line fatigue life monitoring system for boiler drums based on the inverse problem of heat conduction method［J］.Advanced Materials Research，2012，413：361－366.