锅炉汽包瞬态温度场在线监测

李 斌， 陈 丰， 史良宵

（华北电力大学 能源动力与机械工程学院，河北保定071003）

温度场计算是汽包应力计算和疲劳寿命分析的基础，根据瞬态温度场，可实时监测汽包疲劳寿命损耗，保证机组安全经济运行［1－5］.

传统的汽包温度场计算通常采用直接解法［6］，即在一定的初始条件和边界条件下求解导热微分方程，它需要已知求解区域的所有边界条件.由于汽包内有各种复杂的流动和换热过程，因而采用该方法时需要进行很多假设，从而带来较大的计算误差.另外，直接解法最致命的缺陷是需要知道汽包内壁与水和水蒸气的对流传热系数，而汽包内复杂的工况使得很难得到确切的传热系数，通常采用经验数据，但会导致误差［7－10］.此外还要已知热流密度和流体温度，尽管热流密度和流体温度都可以通过热流密度计和热电偶测得，然而该数据只适用于低压情况下，不适用于工程实际高温高压或者不便布置传感器等情况下瞬态传热的计算.

反演解法（Inverse Method）［11－12］采用控制容积法，在控制体上选择一些离散的测点，采用热电偶测温，通过测得的温度值实现整个温度场的求解，即导热反问题（Inverse Problem of Heat Conduction）求解［13－14］.采用反演解法求解大型高温高压设备温度场不需要进行直接解法中的假设，仅需要知道高温高压设备外壁的温度，即可快捷准确地计算出该设备的温度场.计算汽包的温度场时，汽包外壁的温度可以通过在汽包外壁布置热电偶进行测量，比较准确和方便.该方法克服了直接解法中的缺陷，不需要知道汽包内壁的传热系数便可求解汽包温度场.

1 反演解法

应用反演解法实现汽包瞬态温度场的在线监测，首先在汽包的外壁布置测点，安装热电偶测取外壁温度，然后根据测点对汽包横截面进行网格划分，最后利用控制容积法对导热微分方程进行离散，将局部热传导方程转化为常微分方程，由外向内逐次内推，求解得到汽包横截面的温度场.

汽包横截面示意图见图1，图中kw、ks分别表示汽包内饱和水.饱和水蒸气与汽包内壁的对流传热系数，W／（m2·K）.

图1 汽包横截面示意图Fig.1 Sectional view of the boiler drum

汽包模型相关参数如下：内径r1为0.400m，外径r5为0.426m；汽包材料的物性参数：导热系数λ为36.0W／（m·K），密度ρ为7 850kg／m3，比热容c为468J／（kg·K）.

1.1 网格划分

根据汽包横截面外壁热电偶的布置情况，结合反演解法的思想对其进行网格划分，见图2.图中，r1和r5分别为内径和外径，m；r2、r3和r4分别表示中间各层半径，m.由于汽包横截面的对称性，只绘制右半部分.沿着半径方向由内向外划分为三层，分别为内层、中间层和外层，沿圆周方向划分为13个节点.

图2 反演解法网格划分示意图Fig.2 Grid division of inverse method

1.2 反演解法

根据已知的外壁温度和热流密度，按照导热问题数值解法的思想，采用热平衡法，对汽包外壁每个节点所代表的控制容积用傅里叶定律列出能量守恒表达式.外层节点32、节点33和节点34的能量守恒表达式如下：

式中：Δφ为容积角度变化量，rad；Δr为容积径向变化量，m；Ti为节点i的温度，℃；qi为节点i处的热流密度，W／m2.

根据式（1）、式（2）和式（3）可以得出中间层节点19、节点20和节点21的温度表达式：

同理，根据图2，对中间层节点20列出能量守恒表达式，得到内层节点7的温度表达式：

将外层节点的温度逐次反演求解得到内层节点的温度，改变外层节点位置可求解得到整个内层节点的温度，从而得到整个汽包横截面的瞬态温度场.

2 Ansys分析验证

由于采用反演解法求解温度场需要知道求解区域外边界的温度及其边界换热条件，因此采用Ansys数值模拟计算结果进行相关验证.

根据导热微分方程、求解区域的边界条件以及初始条件，采用Ansys进行数值求解，将求解得到的外边界的温度作为反演解法的已知条件，通过反演解法求解得到内层节点温度，再与Ansys数值模拟得到的内层节点温度进行比较.

2.1 Ansys数值模拟数学模型的建立

导热微分方程：

边界条件：

初始温度：

汽包内饱和水（或饱和水蒸气）的温度随时间的变化关系如下：

式中：t为时间，s；T∞为流体温度，℃.

Ansys数值模拟的网格划分如图3所示.由于对象具有对称性，故取其右半部分进行分析.沿径向划分为5个节点，周向25个节点，共计96个单元，125个节点.

图3 Ansys数值模拟网格划分Fig.3 Grid division of Ansys numerical simulation

2.2 数值模拟

按照上述模型对汽包横截面瞬态温度场进行数值模拟，截面的初始温度为70℃，模拟的时间步长为10s，终止时刻为3 000s，得到汽包横截面的节点温度随时间的变化曲线如图4所示.图中节点29、节点33和节点37为外层节点，节点3和节点11为内层节点.

2.3 结果验证

将上述Ansys数值模拟得到的外层节点温度作为反演解法的已知条件（即测量外壁温度），时间步长同样取10s，逐渐内推得到相应的内层节点温度，并与Ansys模拟所得的内层节点温度进行比较验证.

图4 Ansys数值模拟所得汽包横截面的节点温度随时间的变化Fig.4 Variation of node temperature with time obtained by Ansys numerical simulation

由于采用反演解法求解温度值时需用到温度对时间的高阶导数，且测量的外边界的温度对时间的变化非常敏感，在逐层推进的过程中，温度值随时间的变化产生的误差被逐层放大，因此，为减小误差，在计算前对外壁绝热层的温度测量值及外壁温度对时间的导数进行修正，采用局部多项式方法对外层温度值进行空间和时间的平滑后再进行计算［8］.

Ansys模拟所得外层节点温度经过相应的平滑后，得到的外层节点温度随时间的变化曲线如图5所示.

图5 经空间、时间平滑后的外层节点29、节点33和节点37的温度Fig.5 Space－time averaged value of temperature at outer nodes 29，33and 37

根据图1，以汽包横截面几何中心线作为汽水分界面，下部饱和水（上部饱和水蒸气）与内壁发生对流传热，各点处温差相对比较小，因此在进行反演解法验证时，中心线以下、中心线处以及中心线以上三部分分别采用节点29、节点33和节点37的温度作为相应的外层温度.

反演解法的计算结果与Ansys数值模拟结果的对比如图6所示，2种算法计算所得的t＝1 000s和t＝2 000s时部分节点的温度值如表1和表2所示.由图6、表1和表2可以看出，反演解法的计算结果与Ansys数值模拟结果具有较高的吻合度，说明反演解法具有较高的精确度.

图6 反演解法和Ansys模拟所得内层节点3和节点11温度值的对比Fig.6 Comparison of temperature changes at internal nodes 3and 11by inverse method and Ansys numerical simulation

表1 t＝1 000s时2种算法计算所得部分节点的温度值Tab.1 Comparison of temperature changes at partial nodes obtained by two methods at t＝1 000s

表2 t＝2 000s时2种算法计算所得部分节点的温度值Tab.2 Comparison of temperature changes at partial nodes obtained by two methods at t＝2 000s

3 无限长圆柱理论解验证

3.1 模型的建立

半径为R的实心圆柱，其材料的导热系数为λ，热扩散率α为常数，初始温度为T0，将其放在温度为Tf并保持不变的流体中发生对流传热，流体与圆柱表面间的表面传热系数k为常数.其模型简化示意图如图7所示.

3.2 理论解析解

圆柱的无量纲过余温度解析解为

图7 理论解模型示意图Fig.7 Model for analytical solution

式中：FO为傅里叶数；η为半径比；τ为时间，s；r为计算半径，m；θ为过余温度，K；θ0为初始时刻过余温度，K；J0、J1分别为零阶和一阶的第一类贝塞尔（Bessel）函数；μn为超越方程的特征值；k为对流传热系数，W／（m2·K）；λ为材料导热系数，W／（m·K）；

选取初始温度为500℃，流体温度为20℃，计算时间步长为1s，对流传热系数为845W／（m2·K），根据理论解析解得到r＝r1和r＝r5处的温度值，如图8所示.

图8 无限长圆柱r＝r1和r＝r5处的理论解析解Fig.8 Analytical solution for infinitely long cylinder in the case of r＝r1and r＝r5

3.3 结果验证

将无限长圆柱r＝r5（即外层）处的理论解析解作为“测量外层温度”，用作反演解法的已知条件进行反演计算.将求解得到的结果与无限长圆柱r＝r1（即内层）处的理论解析解进行对比分析验证，如图9所示.

相对于实际热电偶测量的外层数据而言，理论解析解波动或者误差更小，因此为了验证应用测量数据进行反演解法计算的正确性，在外层理论解析解的基础上增加一个随机误差－0.5～0.5K，再作为“测量外层温度”进行反演解法计算.带扰动情况下内层理论解析解与反演解法结果的对比见图10.

由图9和图10可知，反演解法结果与理论解析解很接近，吻合度高，同时外加扰动情况下结果也比较接近，可以验证反演解法具有较高的精确度.

图9 内层反演解法结果与理论解析解的对比Fig.9 Comparison of temperature changes respectively obtained by inverse method and analytical solution

图10 带扰动情况下内层理论解析解与反演解法结果的对比Fig.1 0 Comparison of temperature changes respectively obtained by inverse method and analytical solution with disturbances

4 试验验证

为了验证反演解法的计算精度，在某小型锅炉外壁安装热电偶进行实际数据测量，将计算结果与测得的试验数据进行比较，对该方法进行验证.

选择锅炉的某一段工况，在锅炉外壁选择适当的位置布置相应的热电偶，布置的外壁测点编号为35～39，内壁测点编号为11，如图2所示.

由于实际现场测量数据存在较大的波动，首先将锅炉外壁的测量数据按照反演解法的步骤进行时间和空间的平滑处理，然后进行计算，即可求解得到内壁测点11的温度随时间的变化.再将计算结果与实际热电偶测量结果进行对比分析，如图11所示.

由于实际运行工况的复杂性，实际测量结果存在一定的波动，另外热电偶测量数据本身存在一定的误差，但是由图11可知，反演解法的计算结果与热电偶的实际测量结果吻合度较高.

图11 内壁测点11温度的计算结果与实测结果的对比分析Fig.1 1 Comparison of temperature changes at node 11obtained by inverse method and experimental test

5 结 论

（1）反演解法思路简单明了，求解过程不需要迭代，计算精度高，相对于直接解法而言，求解网格相对较少，但是精度却一致.另外计算不需要已知内壁换热条件，只需知道外层节点的温度即可反演得到内层节点温度，从而得到整个汽包的瞬态温度场分布.同时也避免了打孔对压力容器造成的设备寿命损耗，提高了设备的使用寿命.

（2）通过Ansys数值模拟，取汽包外壁边界换热条件为绝热，从二维瞬态的角度来验证反演解法的计算结果，证明该方法计算精度高，二者具有很好的吻合度.

（3）采用无限长圆柱理论解析解，圆柱体外壁边界与流体对流传热，从一维瞬态的角度验证反演解法的适用范围，同时验证了其计算的准确性；另外在附加扰动的情况下同样得到很好的吻合度.

（4）通过与实际试验数据的对比分析，验证了反演解法具有较高的计算精度，与实际试验数据具有很好的吻合度，满足工程应用的要求.

［1］管德清，莫江春，张学纶，等.电站锅炉汽包寿命在线监测系统［J］.动力工程，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］刘彤，徐钢，庞力平，等.锅炉炉内承压部件的蠕变分析及寿命计算［J］.动力工程，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］李立人，陈玮，盛建国，等.锅炉受压元件的高温蠕变－疲劳寿命设计计算方法［J］.动力工程，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］.热能动力工程，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］张健.超临界锅炉炉外承压部件的寿命分析及在线检测［D］.北京：华北电力大学，2008.

［6］赵铁成，沈月芬，朱国桢，等.电站锅炉锅筒温度场计算——三维非稳态变物性材料不均匀导热问题有限元分析［J］.中国电机工程学报，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.