CN101476715A - 一种电站锅炉水冷壁失效预警方法 - Google Patents

Abstract

本发明涉及一种水冷壁在线失效预警方法，包括以下步骤：水动力计算、金属壁温计算、应力解析、疲劳寿命消耗解析、蠕变寿命消耗解析、总寿命消耗解析。通过流动传热规律在线实时计算该管组所有管子的流量、进出口压力、沿长度各段管子的炉内壁温值、应力变化幅值等参数，结合在各辐射区具有流动传热代表性的管段布置少量的测点，和已有的炉外进出口测点作为计算、校核点，输入计算机计算处理，以判断是否出现水动力不稳定性、超温、疲劳蠕变破坏等不安全因素。计算结果实时显示、记录有关数据，作为判断水冷壁是否安全运行的依据，实现设备的优化运行和依据状态安排检验与维修管理，全面实施设备状态检修。

Description

一种电站锅炉水冷壁失效预警方法

技术领域

本发明涉及电站锅炉安全技术领域，涉及一种电站锅炉水冷壁失效预警方法。

背景技术

超(超)临界压力发电技术现在已经是成熟的技术，具有效率高、排放少、易于调峰、运行稳定的特点，但是实际运行中，水冷壁热偏差将导致流量偏差扩大，致使偏差管内工质热物理性质剧烈变化，进而产生流量偏差和传热恶化，使水冷壁壁温偏差增大，其危害比亚临界自然循环锅炉的程度严重得多。而传统的测量方法由于受到技术的约束，无法在炉膛设置大量的测点，测量整个炉膛工况的实时变化，且计算时依赖大量的经验参数，使得推算的结果误差值过大。

因此，应针对性地对水冷壁运行参数进行在线监测，利用管内流动传热规律及少量的测量参数对部件不同位置的实时状态(温度、流量、压力及残余寿命)进行评估，及时正确地将状态和寿命评估结果应用到设备管理决策中，可明显提高设备运行的安全性、可靠性，实现设备的优化运行和依据状态安排检验与维修管理，全面实施设备状态检修。

发明内容

本发明目的在于，提供一种水冷壁失效预警方法，该方法能够在线监测少量的电站锅炉水冷壁特征参数作为计算、校核点，通过流动传热计算在线判断水冷壁是否发生流量脉动、壁温超温，通过疲劳蠕变计算在线判断水冷壁剩余寿命。

为了实现上述任务，本发明采取如下技术解决方案，

首先，根据电站锅炉结构、运行数据划分水冷壁回路，结合运行时水冷壁总流量、进口压力数据，进行水动力计算，得到各回路流量、压力数据；

其次，在结构数据及各回路流量、压力数据基础上，采用数值方法，进行金属壁温计算，获得整个鳍片管的温度场；

再次，在已得到水冷壁各回路压力、流量、整个鳍片管壁温的基础上，根据公式σ_p＝K_PPD/2S_S得到切向薄膜应力σ_p，其中P为工作压力，MPa；D为管子外径与内径平均值，mm；S_s为管子实际壁厚，mm；K_p为综合减弱系数和应力集中系数；

根据公式 ${\mathrm{\σ}}_T=K_t\frac{E_a\·\mathrm{\ΔT}}{1-v}[1/\ln (R_2/R_1)-\frac{{2R}_2^2}{R_2^2-R_1^2}]$ 得到切向热应力σ_T。其中，E为当时温度下的弹性模数；α为当时温度下的微分线形膨胀系数；v为泊桑比；R₂为外径；R₁为内径；K_t为热应力集中系数；

最大应力值由每次循环中的最大切向薄膜应力σ_p和最大切向热应力σ_T相迭加而得，最小应力值则由最小切向薄膜应力σ_p和最小切向热应力σ_T相迭加而得，从而形成应力履历数据库，为疲劳损伤计算提供应力振幅，为蠕变损伤计算提供有效蠕变应力数据；

最后，在应力履历数据库的基础上，结合试验所建立的疲劳设计曲线S-N曲线，先统计出部件承受的各种循环应力的幅值，再根据S-N曲线找出对应的断裂周次，由公式 $D=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{n_i}{N_i}$ 得到部件的累积损伤程度，当D≥1时，判定水冷壁管因疲劳发生失效。其中，n_i为各种循环应力的循环次数；N_i为相应循环应力的允许循环次数；同时进行蠕变寿命损伤的评定，对于某应力，绝对温度T、蠕变破坏时间τ和应力σ采用公式T(C+lgτ)＝a₀+a₁lgσ+a₂lg²σ+a₃lg³σ计算，得到蠕变破坏时间τ。对于给定材料，可根据其实测对应不同温度和时间的材料持久强度数据求出系数a₀、a₁、a₂、a₃和常数C，再根据公式 ${\mathrm{\Φ}}_r=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\Δ\τ}}_i}{{\mathrm{\τ}}_i}$ 进行计算，当Φ_r≥1时，判定水冷壁管因蠕变发生失效，其中，φ_r为材料的蠕变寿命损耗量，⊿τ_i为在i参数下部件的运行时间，τ_i为在i参数下部件的蠕变断裂时间；若D≤1且Φ_r≤1，但φ＝D+Φ_r≥1时，判定水冷壁因疲劳、蠕变失效，实现水冷壁在线失效预警。

所述的水动力计算，通过逐步计算，获得个回路，包括节流元件、局部阻力元件、水冷壁等，不同流量下的压降，这样获得回路分散点的数据，采用四次方切比雪夫多项式拟合，在以获得的压降-流量数据基础上，得到回路的水动力曲线方程，在串联***中，根据工质流量相等，压降叠加，在并联***中，根据压降相等，流量叠加来获得串、并联回路的压降-流量或流量-压降数据，据此通过切比雪夫多项式拟合获得串并联***的水动力曲线方程，对于直流锅炉，入口总流量M_sum总是已知的，即可求得压降ΔP。

所述的金属壁温计算采用区域分离法建立了膜式水水冷壁稳态无内热源而且变物性的二维温度场计算模型，获得整个鳍片管的温度场，可用于计算分析膜式水冷壁界面温度场分布、水冷壁管壁温沿炉膛高度的分布和炉膛同一截面各水冷壁管的温度分布等问题。

本发明的水冷壁失效预警方法的技术效果是：

1、在水动力计算方法上，我国水动力计算常用的方法弦截法和拟牛顿法，用以求解复杂的多极串并联回路水动力问题时，难以保证方程求解的收敛性，并且不经济。针对以上问题，采用四次方切比雪夫多项式拟合法，具有足够的精度，但计算工作量只有上述方法的几十分之一，适用于超临界和超超临界锅炉的水冷壁；

在金属壁温计算方法上，针对以往计算方法只能计算水冷壁某些特殊点的温度，而不能计算整个鳍片管的温度场，本专利在有限差分方法的基础上，进一步完善膜式水冷壁温度场计算的理论模型和程序，用数值计算获得整个鳍片管的温度场，适用于亚临界、超临界和超超临界锅炉的膜式水冷壁；

2、能以少量的特征参数实时地在线计算各回路管段的流量、壁温、剩余寿命；

3、实现设备的优化运行和依据状态安排检验与维修管理，全面实施设备状态检修。

附图说明

图1是在线寿命监测***框图

图2是并联回路***图。

图3是并联回路流量叠加原理图。

图4是串联回路压降叠加原理图。

图5是水动力曲线方程计算流程图。

图6是单管区段压降计算子程序框图

图7是水冷壁水动力程序框图

图8是膜式水冷壁数理问题示意图。

图9是极坐标系下的网格***图。

图10是直角坐标系下的网格***图。

图11是e界面两侧的几何关系图。

图12是附加源项法图。

图13是区域扩充法示意图。

图14是区域分离法示意图。

图15是膜式水冷壁结构图。

图16是膜式壁向火面角系数计算示意图。

图17是角系数交叉线法示意图。

图18是焊缝上角系数计算示意图。

以下结合附图和发明人给出的具体实例对本发明作进一步的描述。

具体实施方式

参照图1所示，首先，根据电站锅炉结构、运行数据划分水冷壁回路，结合运行时水冷壁总流量、进口压力数据，进行水动力计算，得到各回路流量、压力数据；

其次，在结构数据及各回路流量、压力数据基础上，采用数值方法，进行金属壁温计算，获得整个鳍片管的温度场；

再次，在已得到水冷壁各回路压力、流量、整个鳍片管壁温的基础上，根据公式σ_p＝K_PPD/2S_S得到切向薄膜应力σ_p，其中P为工作压力，MPa；D为管子外径与内径平均值，mm；S_s为管子实际壁厚，mm；K_p为综合减弱系数和应力集中系数；

根据公式 ${\mathrm{\σ}}_T=K_t\frac{E_a\·\mathrm{\ΔT}}{1-v}[1/\ln (R_2/R_1)-\frac{{2R}_2^2}{R_2^2-R_1^2}]$ 得到切向热应力σ_T。其中，E为当时温度下的弹性模数；α为当时温度下的微分线形膨胀系数；v为泊桑比；R₂为外径；R₁为内径；K_t为热应力集中系数；

最大应力值由每次循环中的最大切向薄膜应力σ_p和最大切向热应力σ_T相迭加而得，最小应力值则由最小切向薄膜应力σ_p和最小切向热应力σ_T相迭加而得，从而形成应力履历数据库，为疲劳损伤计算提供应力振幅，为蠕变损伤计算提供有效蠕变应力数据；

最后，在应力履历数据库的基础上，结合试验所建立的疲劳设计曲线S-N曲线，先统计出部件承受的各种循环应力的幅值，再根据S-N曲线找出对应的断裂周次，由公式 $D=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{n_i}{N_i}$ 得到部件的累积损伤程度，当D≥1时，判定水冷壁管因疲劳发生失效。其中，n_i为各种循环应力的循环次数；N_i为相应循环应力的允许循环次数；同时进行蠕变寿命损伤的评定，对于某应力，绝对温度T、蠕变破坏时间τ和应力σ采用公式T(C+lgτ)＝a₀+a₁lgσ+a₂lg²σ+a₃lg³σ计算，得到蠕变破坏时间τ。对于给定材料，可根据其实测对应不同温度和时间的材料持久强度数据求出系数a₀、a₁、a₂、a₃和常数C，再根据公式 ${\mathrm{\Φ}}_r=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\Δ\τ}}_i}{{\mathrm{\τ}}_i}$ 进行计算，当Φ_r≥1时，判定水冷壁管因蠕变发生失效，其中，φ_r为材料的蠕变寿命损耗量，⊿τ_i为在i参数下部件的运行时间，τ_i为在i参数下部件的蠕变断裂时间；若D≤1且Φ_r≤1，但φ＝D+Φ_r≥1时，判定水冷壁因疲劳、蠕变失效，实现水冷壁在线失效预警。

一、水冷壁水动力研究

水动力计算的步骤为：回路划分——读入设计节流压降数据——计算并联回路(各辐射区)水动力曲线——计算串联回路(每面墙)水动力曲线——获得整个水冷壁总的水动力曲线——根据总压降，计算各面墙流量——根据各墙流量，代入各墙各辐射区水动力曲线方程，获得压降——根据压降，计算各回路流量——做出水动力不稳定性判断

测点布置：总流量测点、水冷壁各回路进口压力测点

1 水动力计算方法与模型

水动力计算的目的在于，在给定负荷率、节流方式下(节流压降或孔圈阻力系数)，计算各负荷率下的各回路流量，进而获得出口汽温、各种偏差及为壁温计算提供流量参数。

弦截法和拟牛顿法是我国水动力计算常用的方法，实践表明，诸如上述迭代方法用以求解复杂的多极串并联回路水动力问题，难以保证方程求解的收敛性，并且不经济。目前，复杂并联管路***的流量分配主要采用迭代法，因涉及的数组变量多，占用较大的内存，计算时需要多次迭代，CPU耗时大。

本文水动力计算时，是通过逐步计算，获得个回路(包括节流元件、局部阻力元件、水冷壁等)不同流量下的压降。这样获得回路分散点的数据，采用四次方切比雪夫多项式拟合，在以获得的压降-流量数据基础上，得到回路的水动力曲线方程。在串联***中，根据工质“流量相等，压降叠加”的原则，在并联***中，可以根据“压降相等，流量叠加”的原则来获得串、并联回路的压降-流量或流量-压降数据，据此通过切比雪夫多项式拟合获得串并联***的水动力曲线方程。对于直流锅炉，入口总流量M_sum总是已知的，即可求得压降ΔP。

2 回路***区段的划分

按结构尺寸和吸热相差不多的原则，将水冷壁划分为若干个计算回路，在高度方向划分若干个区段。用回路中的一根管的水动力特性代表所在回路的水动力特性，沿高度各个区段的压降之和即为总压降。水冷壁***构成复杂的串并联回路***，其中每一面墙的每一辐射区由若干回路构成并联***，自下而上各个辐射区构成串联***。四面墙组成一个大的并联***，下公共点可认为是水冷壁入口集箱，而上公共点是顶棚出口集箱。

将管组的各种结构数据输入数据库，包括管组管子数量、管长、倾角、内径、材料、集箱结构等。

3 复杂串并联回路水动力特性曲线方程

(1)单回路的水动力特性曲线方程

i回路j段在流量G_i下的出口焓值i_i，j为：

i_i，j＝i_i，j-1+Q_i，j/G_i        (1)

其中，Q_i，j——i回路j段的吸热量。

j段的摩擦压降ΔP_f(i，j，G_i)和重位压降ΔP_h(i，j，G_i)分别为：

ΔP_f(i，j，G_i)＝f₁(P，i_i，j-1，i_i，j，G_i，D_i，j，h_i，j)  (2)

ΔP_h(i，j，G_i)＝f₂(P，i_1，j-1，i_i，j，G_i，D_i，j，h_i，j，)  (3)

在某一级辐射区上累加，可得：

${\mathrm{\ΔP}}_f(i,G_i)=\underset{j}{\mathrm{\Σ}}{\mathrm{\ΔP}}_f(i,j,G_i)---\left(4\right)$

${\mathrm{\ΔP}}_h(i,G_i)=\underset{j}{\mathrm{\Σ}}{\mathrm{\ΔP}}_h(i,j,G_i)---\left(5\right)$

ΔP(i，G_i)＝ΔP_f(i，G_i)+ΔP_h(i，G_i)   (6)

对多个流量重复以上过程的计算，就可得到摩擦压降和重位压降之和与流量的对应关系，通过切比雪夫(chebyshev)曲线拟合数学处理，即可获得i回路的水动力曲线方程。同时得到ΔP＝F_i，1(G)和G＝F_i，2(ΔP)，以备后面计算方便。经验表明，四次chebyshev拟合以足够满足精度要求。

当计算下辐射区域时，必须考虑进口节流压降和三叉管局部阻力压降。其压降为：

ΔP(i，G_i)＝ΔP_f(i，G_i)+ΔP_h(i，G_i)+ΔP_JL(i，G_i)+ΔPSC(i，G_i)   (7)

(2)并联回路的水动力特性曲线方程

参照图2所示，I_I表示进口集箱内工质焓值，I₂表示出口集箱内工质焓值，M_sum表示工质总流量，P_in表示进口集箱内工质压力，P_out表示出口集箱内工质压力，N表示并联的第N根水冷壁管。设并联回路由N个单回路组成，已知入口工质的总流量G_sum，入口焓i_in和入口压力P_in。在水冷壁结构已知的情况下，可以认为有：

ΔP＝f(G_sum，i_in，P_in)     (8)

ΔP——并联回路进、出口公共点的总压差。

按单回路水动力特性曲线方程的计算方法，分别求出并联回路中的每个回路在i_in和P_in条件下的水动力曲线方程。

参照图3所示，(ρω)₁表示管路1的质量流量，(ρω)₂表示管路2的质量流量，(ρω)₁+(ρω)₂表示并联管路1和管路2流量叠加，△P表示并联管路共同的压力降。由于各回路并联，进、出口公共点之间在同一压降下工作，并联回路水动力特性曲线可以根据“压降相等，流量叠加”原则得到。

对回路i，流量为G_i时，已知其压降ΔP_i，则并联回路中其余回路的流量分别为：

G₁＝F_1，2(ΔP_i)G₂＝F_2，2(ΔP_i)G_N＝F_N，2(ΔP_i)  (9)

并联回路的总流量为： $G=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{G}_i---\left(10\right)$

对其他的G_i计算就可获得整个并联回路的ΔP—G数据，同样采用chebyshev拟合，即可获得并联回路的水动力曲线方程。

ΔP＝F₃(G)          (11)

获得并联回路的压降。

已知的G_sum，由式(10)即可获得ΔP_sum，同时也可得到出口处的i_out和P_out。

(3)串联回路的水动力特性曲线方程

参照图4所示，⊿P₁表示管路1的压力降，△P₂表示管路1的压力降，⊿P₁+⊿P₂表示串联管路1和管路2的压力降叠加，ρω表示串联管路共同的质量流量。串联***中，总水动力特性曲线可以根据工质“流量相等，压降叠加”的原则得到。

参照图5所示，方框I、II表示两个相接的并联回路或者单回路或者单管，i₁表示方框I前的工质焓值，i₂表示方框I后、方框II前的工质焓值，i₃表示方框II后的工质焓值；P₁表示方框I前的工质压力，P₂表示方框I后、方框II前的工质压力，P₃表示方框II后的工质压力；G_sum表示方框I、II前后共同的工质流量。不论方框代表的是并联回路、单回路，还是单管，已知方框I前的i₁，G_sum，P₁，依据单回路的水动力特性曲线方程和并联回路的水动力特性曲线方程的方法，既可求出方框I后的i₂，G_sum，P₂和方框II后的i₃，G_sum，P₃，依次类推，便可得到G_sum与多级串联回路的上下公共点之间的压降的对应关系。多对G_sum—ΔP数据，即可获得串联回路的水动力曲线方程。

4 压降计算补充方程

锅炉水冷壁管中加速压降可忽略，只计算摩擦压降和重位压降。水冷壁有光管和内螺纹管组成。

(1)单相压降计算

对光管，处于旺盛紊流区单相摩擦系数为：

$f=\frac{1}{4{\left(\lg (3700D_{\mathrm{in}}/k)\right)}^2}---\left(12\right)$

其中，k为粗糙度，对于SA213T12，k＝0.008。

ΔP_f＝f·(L/D_in)·(G²/2)·v    (13)

对于四头内螺纹管，CE公司推荐将光管的摩擦阻力乘以1.66得到内螺纹管的摩擦阻力：

${\mathrm{\ΔP}}_f^{*}=1.66*{\mathrm{\ΔP}}_f---\left(14\right)$

Kohler给出的内螺纹管摩擦阻力计算公式为：

${\mathrm{\ΔP}}_f^{*}=f^{*}\·(L/D_{\mathrm{in}})\·(G^2/2)\·v---\left(15\right)$

f^*＝0.0213+1.01×10⁴/Re^1.2   (16)

通过实验进行四头内螺纹管的阻力特性研究给出的摩擦压降计算公式为：

$f=\frac{2.6082}{{\mathrm{Re}}^{0.5546}}+0.03716---\left(17\right)$

在亚临界参数的广阔范围内，利用式(13)和(14)计算结果相差小于14％；式(13)和(17)计算结果相差约为11％。

对于超临界压力，单相摩擦压降计算公式中的平均比容按各区段进、出口平均焓值选取，即：i＝(i_in+i_out)/2

若i处在(i＝1670kJ/kg～2720kJ/kg)，且在区段内Δi>210kJ/kg时，平均比容按下式选取：

v＝(v_out·i_out-v_in·i_in)/(i_in-i_out)    (18)

单相重位压降计算公式为：

ΔP_h＝ρ·g·h＝g·h/v     (19)

超临界压力时比容v的选取同上。

(2)汽液两相压降计算

文献中已发表了许多垂直光管的两相流摩擦压降，普遍认为Chisholm方法较为准确。

引进两相压降倍律

：

ΦL²＝ΔP_f/ΔP_L＝1+c/x+1/x²    (20)

x²＝ΔP_L/ΔP_g    (21)

△P_L，△P_g分别为两相流中液相单独流过管道时的阻力及汽相单独流过管道时的阻力。

ΔP_L＝λ_L(h/D_i){[G(1-x)]²/2}v_L   (22)

ΔP_g＝λ_g(h/D_i){(Gx)²/2}v_g      (23)

进入旺盛紊流区后

λ_L＝λ_g＝1/{4[lg(3700D_i/k)]}²      (24)

系数C值按下式计算

C＝{λ+(C₁-λ)[(v_g-v_L)/v_g]^0.5}[(v_g/v_L)^0.5+(v_L/v_g)^0.5]   (25)

1.G≤G^*＝1500kg/m²·s(粗糙管)

取λ＝1.0，C₁＝G^*/G，计算C值。

2.G>G^*＝1500kg/m²·s

Φ_L ²＝(1+c₂/x+1/x²)Ψ    (26)

Ψ＝(1+c/T+1/T²)/(1+c₁/T+1/T²)  (27)

C₂＝(v_g/v_L)^0.5+(v_L/v_g)^0.5   (28)

T＝[x/(1-x)]^(2-n)/2(μ_L/μ_g)^n/2(v_L/v_g)^1/2   (29)

求得后，ΔP_f＝ΔP_L·Ψ_L ²，从以上计算过程可以看出，Chisholm方法综合考虑了物性、质量流速的影响。

CE推荐内螺纹管中摩擦阻力为光管中的1.66倍，即

ΔP_f ^*＝1.66ΔP_f    (30)

Kohler给出的内螺纹管汽液两相摩擦阻力计算方法如下：

R＝(1-x)²+x²(v_L/v_g)(λ_g/λ_L)+6x^1.2(1-x)^0.41(v_g/v_L)

                                                (31)

＝(μ_g/μ_L)^0.4[1-(μ_g/μ_L)]F_rL ^-0.05W_oL ^-0.033

这里F_rL＝(G²v_L ²)/(gD_i)，W_oL＝(G²D_iv_L)/σ

λ_g，λ_L由以下两式确定：

λ_g＝0.0213+(1.01×10⁴/Re_g ^1.2)     (32)

λ_L＝0.0213+(1.01×10⁴/ReL^1.2)     (33)

Re_g＝GD_i/μ_g，Re_L＝GD_i/μ_L    (34)

ΔP_f ^*＝RΔP_L＝Rλ(h/D_i)(G²/2)v_L   (35)

过实验进行四头内螺纹管的阻力特性研究给出的汽液两相摩擦压降计算公式：

当P＝9～19MPa时：

C(x)＝1.2803·x^0.7063·(1-x)^0.2515    (36)

当P＝19～21.5MPa时：

C(x)＝0.9885·x^0.6198·(1-x)^0.1976    (37)

两相区重位压降的计算方法如下：

ΔP_h＝ρ_mgh    (38)

混合物的密度为：

ρ_m＝α/v_g+(1-α)/v_L    (39)

α为界面平均含气率。

截面含气率和体积含汽率从理论上存在以下关系：

α＝1/{1+S[(1-β)/β]}            (40)

由体积含汽率β和质量含汽率x的关系

β＝1/{1+(v_L/v_g)[(1-x)/x]}    (41)

可得：

α＝1/{1+S(v_L/v_g)[(1-x)/x]}   (42)

其中S为汽速和水速之比。

滑速比S为多种计算方法，现采用电站锅炉水动力计算方法：

S＝1+(0.4+β²)/(G·v_L)0.5·(1-P/221.5)   (43)

具体计算时，首先确定体积含汽率，进而确定滑速比S，求得截面含汽率α，然后求得ρ_m，即可确定重位压降ρ_m。

参照图6所示，是单管区段压降计算子程序框图，它提供了单管中不同流态的流体的压降计算。首先输入管长、倾角、内径、流量、压力等数据，判断是否为超临界，如果是，直接运用超临界压降计算方法计算压降；如果不是，则进入下一级程序，依次判断进出口均为两相，进口两相出口过热，出口欠热进口过热，进口欠热，出口两相等情况，依次运用相关方法计算压降。

参照图7所示，是水冷壁水动力程序框图，首先读入设计节流压降数据，然后计算各个回路下包括节流在内的分散点数据，然后拟合并联回路(各面墙)的水动力曲线，接着拟合各串联回路(各辐射区)的水动力曲线，接着通过各串并联回路获得整个水冷壁总的水动力曲线，然后根据水冷壁总流量，代到水冷壁水动力曲线方程，获得水冷壁总压降，然后根据总压降，计算各面墙流量，然后根据各墙的流量，代入各墙各辐射区的水动力曲线方程，获得压降，接着根据压降，计算各回路流量、出口焓值、出口压力，通过得到各回路流量、出口焓值、出口压力，可以得到水冷壁壁温，最后计算完毕，进入下一负荷的计算准备。

二、膜式水冷壁温度场计算

本发明在有限差分方法的基础上通过数值计算可以获得整个鳍片管的温度场。

对于短期超温，由于其时间相对于水冷壁寿命而言极其短暂，所以剩余寿命计算的实际意义较小，这里以金属监督及检查为主，达到预防和调整运行的目的。

金属壁温计算的步骤为：数学模型的建立——外部边界条件计算——内部边界条件的计算——膜式水冷壁温度场迭代求解——做出超温判断。

测点布置：延炉膛高度、宽度方向背火侧设置少量测温点，利用数值计算得到的温度场间接得到向火侧壁温，避免了高温和磨损的问题引起测量装置的脱落或失效，适合长期监测；在热负荷较高处(如燃烧器附近)，在炉膛高度、宽度方向有代表性的管束处(如各回路中间管束)布置少量向火侧测点，用以预警短期超温，同时对数值计算得到的温度场分布进行校核。

1 数学物理问题的分析

在锅炉运行中，膜式水冷壁外部一面受到炉膛火焰的强烈辐射，另一面则敷设绝热材料，而管内则是有一定质量流速的水及蒸汽流过。要得到膜式水冷壁管壁稳定的分布规律就是求解在稳定状态下满足给定边界条件的调和方程ΔT＝0的解。该调和方程是一个线性偏微分方程，其中无关变量是空间坐标，而相关变量则是温度，为简化求解，引入如下假设：

(1)炉膛热负荷在一定局部范围内是均匀的；

(2)在锅炉机组水冷壁中，沿管件高度方向的热流可以忽略，因而问题的解域可看作二维；

(3)管壁与管内介质的传热系数沿管子内周界取为常数；

(4)在炉膛内，管子的对流吸热量相对辐射吸热量来说很小，可以忽略不计；

(5)经炉墙绝热材料所散失的热量忽略不计。

为了估计以上简化假设对计算结果可靠性的影响程度，做如下讨论：

(1)虽然炉膛热负荷随炉膛空间位置而变化，但对于所计算的某段位置上的膜式水冷壁而言，就其固定位置上而言热负荷变化不大，所以可以认为在所确定的局部范围内热负荷是均匀的；

(2)相对沿膜式水冷壁管截面径向及周向的传热而言，沿其长度方向的导热较小，而且膜式水冷壁长度方向尺寸相对较大，这样把沿膜式水冷壁长度方向的导热忽略掉，对其计算结果影响很小；

(3)炉膛中辐射换热强度远远超过膜式水冷壁管外壁与炉膛烟气之间的对流换热。因此，在计算中不考虑膜式水冷壁管外壁的对流换热，这对计算结果的影响不大。

(4)在超临界直流锅炉实际运行过程中，通过炉墙绝热材料的散热量很小，将其忽略后的计算结果仍然可以很好的与实际工况相吻合。

综上所述，根据安全性与实际可行性两方面的考虑，可以把超临界直流螺旋管圈膜式水冷壁的计算区域简化为二维温度场计算模型来进行处理。

2 膜式水冷壁温度场数理模型及算法

(1)物理问题和控制方程及其离散

参照图8所示，q表示截面外壁辐射均匀热负荷，1、1’、2、2’、3、4、5、6为所设边界点，A为边界1-2-3-4-1的向火面，B为边界1-2-3-4-1的背火面。

根据上节中的合理假设，可以认为某根管截面外壁辐射热负荷q是均匀的，如果膜式水冷壁相邻管的几何特征和传热特性相同，即：相邻管的管内换热系数相等、管型和结构尺寸相同，认为膜式水冷壁管截面和鳍片的对称面是绝热面，则可选取参照图8(a)所示的斜线填充区为计算对象，这样便可大大减少计算复杂程度及计算时间。我们知道，在超临界压力及亚临界压力核态沸腾区是基本满足内壁放热系数α左右相同这一条件的，但在亚临界传热恶化区，则差别较大。

参照图8(a)所示，斜线填充区所示的计算区域中，可以将垂直管圈膜式水冷壁温度场分布问题简化为二维、稳态无内热源而且变物性(管壁金属导热系数随温度变化)的导热问题进行的求解。由于计算区域是非规则区域，采用区域分离法将计算区域进行分解，一部分为半圆环形，如图8(b)左侧所示；另一部分经扩充后为长方形，如图8(b)右侧所示。计算区域所分成的两部分有一重叠区域，在二者的重叠区实现半圆环区域和长方形区域的耦合迭代，从而求解整个膜式水冷壁温度场分布。下面分别讨论半圆环区域、长方形区域和重叠区域的数理控制方程及其方程的离散。

1)半圆环区域的控制方程及其离散

$\frac{{\∂}^{2}t}{{\∂r}^2}+\frac{1}{r}\·\frac{\∂t}{\∂r}+\frac{1}{r^2}\·\frac{{\∂}^{2}t}{{\∂\mathrm{\θ}}^2}=0---\left(44\right)$

参照图9所示，在极坐标系中，采用内节点法对半圆环形区域进行离散化。O为极坐标原点，n、e、s、w为控制体的4个边；N、E、S、W为所选4个节点；△r为控制体r方向范围，△θ为控制体θ方向范围；(δr)_n为N点r方向坐标，(δr)_s为S点r方向坐标，(δθ)_w为W点θ方向坐标，(δθ)_e为E点θ方向坐标。

在离散化区域中，式(44)的离散化形式为：

a_PT_P＝a_ET_E+a_wT_w+a_NT_N+a_ST_S+b    (45)

式中：

$a_E=\frac{\mathrm{\Δr}}{r\·{\left(\mathrm{\δ\θ}\right)}_e/k_e}$ $a_w=\frac{\mathrm{\Δr}}{r_p\·{\left(\mathrm{\δ\θ}\right)}_w/k_w}$

$a_N=\frac{r_n\·\mathrm{\Δ\θ}}{{\left(\mathrm{\δr}\right)}_n/k_n}$ $a_s=\frac{r_s\·\mathrm{\Δ\θ}}{{\left(\mathrm{\δr}\right)}_s/k_s}$

a_p＝a_E+a_W+a_N+a_S-S_p·△_V

b＝S_c·△V

△V＝0.5(r_n+r_s)·△_r·△θ

其中，k_e、k_w、k_n及k_s为控制容积界面上的导热系数，可按调和平均方法计算，以k_e为例， $k_e={\left(\mathrm{\δ\θ}\right)}_e/[\frac{{\left(\mathrm{\δ\θ}\right)}_{e^{-}}}{k_P}+\frac{{\left(\mathrm{\δ\θ}\right)}_{e^{+}}}{k_E}];$ 源项S在任一控制容积中之值可以表示为温度的线性函数，即S＝S_C+S_P·T_P，S_C为常数；S_P为曲线S＝f(T)在P点的斜率。

2)矩形区域的控制方程及其离散

$\frac{{\∂}^{2}t}{{\∂X}^2}+\frac{{\∂}^{2}t}{{\∂Y}^2}=0---\left(46\right)$

参照图10所示，在直角坐标系中，采用内节点法对矩形区域进行离散化。O为直角坐标原点，n、e、s、w为控制体的4个边；N、E、S、W为所选4个节点；△x为控制体在X方向长度，△Y为控制体在Y方向长度；(δy)_n为N点Y方向坐标，(δy)_s为S点Y方向坐标，(δx)_w为W点X方向坐标，(δx)_e为E点X方向坐标。

在离散化区域中，式(46)的离散化形式为：

式中：

$a_E=\frac{\mathrm{\ΔY}}{{\left(\mathrm{\δX}\right)}_e/k_e}$ $a_w=\frac{\mathrm{\ΔY}}{{\left(\mathrm{\δX}\right)}_w/k_w}$

$a_N=\frac{\mathrm{\ΔX}}{{\left(\mathrm{\δY}\right)}_n/k_n}$ $a_s=\frac{\mathrm{\ΔX}}{{\left(\mathrm{\δY}\right)}_s/k_s}$

a_p＝a_E+a_W+a_N+a_S-S_p·△V

b＝S_c·△V

△V＝△X·△Y

其中，k_e、k_w、k_n及k_s为控制容积界面上的导热系数，可按调和平均方法计算。

参照图11所示，P、E为两节点，(δx)_P ^-为P点横坐标，(δx)_E ⁺为E点横坐标，δx为PE两点在X方向长度。以k_e为例， $k_e={\left(\mathrm{\δx}\right)}_e/[\frac{{\left(\mathrm{\δx}\right)}_{e^{-}}}{k_P}+\frac{{\left(\mathrm{\δ\θx}\right)}_{e^{+}}}{k_E}].$

3)边界条件的处理方法

对于半圆形区域，如图8(b)，边界1-6、3-4、4-5均为绝热条件，边界1-2为给定热负荷条件，边界2-3为给定温度条件，边界5-6为第三类边界条件，即对流换热边界条件。

对于矩形左侧区域，如图8(b)，边界1-2、2-3、3-4均为绝热边界条件，边界1-4为给定温度条件。

边界条件处理采用附加源项法，把由第二类和第三类边界条件所规定的进入或导出计算区域的热量作为与边界相邻的控制容积的当量源项。

对第二类边界条件，参照图12所示，P表示控制容积，N、E、S为所选节点，△y为控制容积在Y方向长度，q_B为边界热流密度。与边界相邻的控制容积P的附加常数源项为：

S_c，ad＝(A/△V)·q_B   (47)

且同时令a_w＝0。其中A为所研究控制容积在边界上的换热面积，△V为控制容积P的体积。

对第三类边界条件，参照图12所示，边界控制容积的附加源项为：

$S_{c,\mathrm{ad}}=-\frac{A}{\mathrm{\ΔV}}\·\frac{T_f}{1/\mathrm{\α}+{\left(\mathrm{\δX}\right)}_w/k_B}---\left(48\right)$

$S_{P,\mathrm{ad}}=-\frac{A}{\mathrm{\ΔV}}\·\frac{1}{1/\mathrm{\α}+{\left(\mathrm{\δX}\right)}_w/k_B}---\left(49\right)$

且同时令a_w＝0。

由于令a_w＝0，就把未知的边界节点温度排除在外，在获得内部节点的解以后，再按边界条件来确定边界节点温度。

4)扩充区域的处理

由于图8(b)中的阴影部分区域不够规则，采用区域扩充法把计算区域扩充到直角坐标系易于描述的矩形区域。真实边界上作用有一定的热流，把它以附加源项的形式送到与真实边界相邻的控制容积中去，而扩充区域则处于绝热状态。参照图13所示，域扩充法，c、d、e、f为控制容积P的四个边界点。

在图中控制容积P与真实边界相交。于是P控制容积的附加源项为：

S_ad＝q·L_ef/△V_cdef   (50)

其中：L_ef表示真实边界与控制容积P的两条边界相交部分的长度，△V_cdef是控制容积P的体积，q是给定的边界热流密度。在引入附加热源的同时再令扩充区中的导热系数为零，于是控制容积P中的热量就无法向扩充区域中传递，这就实现了真实边界上的均匀加热条件。

5)分离区域的耦合

区域分离法参照图14所示，1、2、3、4、5、6、7、8表示所选的8个节点。

在图中假设

上各节点温度，算得其内部各节点(环形区域)温度分布后，ab上各节点温度通过插值进行计算。

$T^{*}=\frac{X_{0-8}}{X_{6-8}}\·(T\′-T\′\′)+T\′\′---\left(51\right)$

其中： $T\′=\frac{X_{0-5}}{X_{6-5}}\·(T_6-T_5)+T_5$

$T\′\′=\frac{X_{0-7}}{X_{7-8}}\·(T_8-T_7)+T_7$

如此便求得了ab上各节点(矩形区域)的温度，把这些温度作为已知条件，算得矩形区域内各节点温度，随之可通过插值求得

上各节点的温度。

$T_{*}=\frac{Y_{0-3}}{Y_{1-3}}\·(T_a-T_b)+T_b---\left(52\right)$

其中， $T_a=\frac{X_{0-1}}{X_{2-1}}\·(T_2-T_1)+T_1;$

$T_b=\frac{X_{0-3}}{X_{4-3}}\·(T_4-T_3)+T_3.$

上式中，X_0-5为被求节点与第5点在X方向上的距离；Y_0-3为被求节点与第3点在Y方向上的距离。比较前后两次算得

上各节点温度，如误差小于规定精度，则终止迭代。

3 膜式水冷壁边界条件的确定

根据第1节中的论述可知，膜式水冷壁的边界条件包括管内壁对流边界条件、管外壁边界条件包括膜式壁向火面的辐射边界和背火面及鳍片两端面的绝热边界，下面分别讨论。

(1)外部边界条件

1)膜式壁向火面辐射边界

根据设计厂家提供的沿炉膛高度方向上热负荷分布曲线，将此曲线拟合成多项式。将沿高度方向的平均热负荷q乘以沿宽度方向上的热负荷分配系数η，即可得到炉膛内任一处的壁面热负荷，即：

q＝q·η         (53)

2)膜式水冷壁表面热流密度分布的精确解

膜式水冷壁外部在炉膛中主要受辐射换热的影响。膜式水冷壁辐射边界条件少有精确的计算方法。为了提高其壁温计算的精确性，本文采用如下膜式水冷壁角系数精确解的计算方法。

A 结构参数

根据第1节中的假设，不考虑其沿管长方向的导热(假定在一局部范围内，炉膛热负荷是均匀的)，可将垂直管圈膜式水冷壁的温度场简化为二维稳态导热问题。参照图15所示，膜式水冷壁主要由管子、鳍片以及焊缝三部分组成，所以需要分别确定以上三个部分的热流密度分布。图中P为焊缝位置参数；s为管子节距；2δ为鳍片厚度；θ为焊缝交角；r为管子半径。

B 辐射边界条件分析

根据模型假设，锅炉炉膛的辐射面积热负荷q_L在一定工况下是不变的，因此图15中

面上所吸收的热量为：

Q_L＝q_L·s                 (54)

在管壁、焊缝及鳍片上任一微元长度dx上所吸收的热量为：

Q_dx＝q′·dx         (55)

q′为该微元长度dx上的热流密度。

ab面对dx的角系数为

因此，微元长度上的热流密度为

从上式可知，求解微元长度上的热流密度之前，首先需要计算角系数

C 膜式水冷壁热流密度分布函数

a 鳍片上热流密度分布函数

参照图16所示，在鳍片距O₁点x处取微元长度dX，过c点分别作圆O₁、O₂的切线交圆O₁于M点，交圆O₂于N点；过d点分别作圆O₁、O₂的切线交圆O₁于M’点，交圆O₂于N’点。则ab、

cd、构成了由四个非凹表面组成的封闭***。由角系数的相对性与完整性可得出确定角系数的交叉线法如下：

上式中，

为表面1对表面2的角系数。参照图17所示，角系数

采用下式计算。

根据角系数交叉线法，可以确定图16中ab(火焰)对cd(dx)的角系数为：

据几何关系可以得：

ab＝S(60)

将式(60)～(64)代入(59)式，根据f(u+du)-f(u)＝f′(u)du整理可得：

g_g，dx＝(1/2S)·[(r·ctgα-δ·cscα)·ctgα·dα+(r·ctgβ-δ·cscβ)·ctgβ·dβ](65)

由几何条件可以得α和X的关系：

tgα·(X-r·sinα)＝(r·cosα-δ)    (66)

由式(66)可得：

$\mathrm{\α}=\arcsin \frac{r\·X-\mathrm{\δ}\·{(X^2+{\mathrm{\δ}}^2-r^2)}^{1/2}}{X^2+{\mathrm{\δ}}^2}---\left(67\right)$

将等式(67)两边对x求导，可得：

${\sec }^2\mathrm{\α}\·\frac{\mathrm{d\α}}{\mathrm{dX}}\·(X-r\·\sin \mathrm{\α})+\mathrm{tg\α}\·(1-r\·\cos \mathrm{\α}\·\frac{\mathrm{d\α}}{\mathrm{dX}})$

$=-r\·\sin \mathrm{\α}\·\frac{\mathrm{d\α}}{\mathrm{dX}}$

$\frac{\mathrm{d\α}}{\mathrm{dX}}=\frac{-\sin \mathrm{\α}}{r\·\mathrm{tg\α}-X\·\sec \mathrm{\α}}---\left(68\right)$

同理可以推导出β和X的关系如下：

tgβ·[(S-X)-r·sin β]＝(r·cos β-δ)          (69)

$\mathrm{\β}=\arcsin \left[\frac{r\·(S-X)-\mathrm{\δ}\·{[{(S-X)}^2+{\mathrm{\δ}}^2-r^2]}^{1/2}}{{(S-X)}^2+{\mathrm{\δ}}^2}\right]---\left(70\right)$

将式(70)两边对X求导，可得：

$\frac{\mathrm{d\β}}{\mathrm{dX}}=\frac{\sin \mathrm{\β}}{[(S-X)\·\sec \mathrm{\β}-r\·\mathrm{tg\β}]}---\left(71\right)$

将α、dα、β及dβ的表达式代入式(65)即可以得出的表达式。则膜式水冷壁焊缝上的热流密度分布函数为：

b膜式水冷壁焊缝上的热流密度分布函数

参照18所示，在焊缝上取微元cd建立以点O₁为原点，

为X轴，

为Y轴的直角坐标系。c点X轴坐标为X，d点为X+dX，则cd长度cd＝dX/cosθ。分别过c、d作圆O₁的切线交圆O₁于M、M’，作圆O₂的切线交圆O₂于N、N’。

则ab、

cd、

构成了由四个非凹表面组成的封闭***。根据角系数交叉线法，ab(火焰)对cd(dX/cosθ)的角系数为：

根据几何关系可得：

将式(74)～(77)代入(73)式，根据

f(u+du)-f(u)＝f′(u)du整理可得：

$\frac{{\mathrm{\φ}}_{g,\mathrm{dX}}}{\cos \mathrm{\θ}}=(1/2S)\·\{[r\·\mathrm{rg\β}-(S-X)\·\sec \mathrm{\β}]\·\mathrm{tg\β}\·\mathrm{d\β}+(\sec \mathrm{\α}+\sec \mathrm{\β})\·\mathrm{dX}---\left(78\right)$

$+(r\·\mathrm{tg\α}-X\·\sec \mathrm{\α})\·\mathrm{tg\α}\·\mathrm{d\α}\}$

由几何条件可以得α和X的关系：

tgα·(r·cscα-X)＝(P-X)·tgθ                (79)

则 $\mathrm{\α}=\arcsin \frac{r\·X-W\·{(X^2+W^2-r^2)}^{1/2}}{X^2+W^2}---\left(80\right)$

其中：

W＝(P-X)·tgθ                                 (81)

将式(80)两边对X求导，进而求得dα/dX

$\frac{\mathrm{d\α}}{\mathrm{dX}}=\frac{\mathrm{tg\θ}-\mathrm{tg\α}}{\sec \mathrm{\α}\·(r\·\mathrm{tg\α}-X\·\sec \mathrm{\α})}---\left(82\right)$

同理可得β和X的关系：

tgβ·[r·cscβ-(S-X)]＝(P-X)·tgθ            (83)

解得：

$\mathrm{\β}=\arcsin \left[\frac{r\·(S-X)-W\·{[{(S-X)}^2+W^2-r^2]}^{1/2}}{{(S-X)}^2+W^2}\right]---\left(84\right)$

其中：W＝(P-X)·tgθ

将式(84)两边对X求导，进而求得dβ/dX

$\frac{\mathrm{d\β}}{\mathrm{dX}}=\frac{\mathrm{tg}\mathrm{\β}+\mathrm{tg\θ}}{\sec \mathrm{\β}\·[(S-X)\·\sec \mathrm{\β}-r\·\mathrm{tg\β}]}---\left(85\right)$

将α、dα、β、dβ的表达式代入(78)式，即可得到由r、S、P、θ及X和dX表达的φ_g，dX/cosθ的显示关系式。

焊缝上的热流密度分布函数为：

$q\′\left(X\right)=\frac{{\mathrm{\φ}}_{g,\mathrm{dX}.\cos \mathrm{\θ}}}{\mathrm{dX}/\cos \mathrm{\θ}}\·q_i\·S---\left(86\right)$

c 管壁上的热流密度分布函数

如图18所示，在圆O₁上取微元

分别从c、d作圆O₂的切线交圆O₂于c′、d′，则

abab构成了由四个非凹表面组成的封闭***。

根据角系数交叉线法，ab(火焰)对cd(r·dα)的角系数为：

据几何关系：

r·cos(α+dα)]·csc(β+dβ)

将式(88)～(90)代入(2-91)式，根据f(u+du)-f(u)＝f′(u)du得：

φ_g，r·dα＝(1/2S)·[r·(1-sinα·cscβ)·dα]·dα+   (91)

r·ctgβ·cscβ·(cosβ+cosα)·dβ

下面确定β和α的关系：

r·cos β-r·cos α＝tg β·(S-r·sin α-r·sin β)         (92)

则有：

$\mathrm{\β}=\arcsin \frac{-r(r\·\sin \mathrm{\α}-S)-{[r^2\·S\·{\cos }^2\mathrm{\α}\·(S-2r\·\sin \mathrm{\α})]}^{1/2}}{r^2+S^2-2r\·S\·\sin \mathrm{\α}}---\left(93\right)$

将(93)式两边对α求导的dβ/dα

$\frac{\mathrm{d\β}}{\mathrm{d\α}}=\frac{\sin (\mathrm{\α}+\mathrm{\β})}{-\sin (\mathrm{\α}+\mathrm{\β})+(S/r)\·\cos \mathrm{\β}}---\left(94\right)$

将β、dβ的表达式代入(91)式，即可得出由r、s及α、dα表达的φ_g，r·dα显示表达式。

管壁上的热流密度分布函数为：

q′(α)＝[φ_g，r·dα/(r·dα)]·q_i·S      (95)

(2)膜式水冷壁内部边界条件

膜式水冷壁的内部边界条件是第三类对流边界，西安交通大学在广泛的压力范围内研究了内螺纹管与介质水之间的传热特性和管壁温度的特性，获得了管内壁与水的对流换热系数，从而确定膜式水冷壁的内部边界条件。本节重点论述在亚临界及超临界及超超临界压力区内，内螺纹管内壁与水之间传热系数的确定。

1)垂直管圈水冷壁亚临界压力区

研究结果表明，在亚临界压力区，内螺纹管的传热强化作用明显，管内工质与管壁之间的换热系数主要分单相对流区、核态沸腾区和干涸后传热区分别计算。

A 单相对流区

水冷壁在锅炉低负荷下运行时，处于亚临界压力区。由于水冷壁入口有一定的过冷度，所以处于锅炉下炉膛的垂直上升内螺纹管水冷壁存在一段单相对流区域。在此区域，可采用西安交通大学根据试验数据拟合得到的垂直上升内螺纹管内单相流体对流换热计算关联式

Nu＝0.05298·Re^0.7417·Pr^0.4        (96)

上式的适用范围为：压力P＝9.0MPa～19MPa。

B 核态沸腾区

锅炉在亚临界压力区运行时，垂直管圈内螺纹管水冷壁在相当长范围内处于核态沸腾状态两相沸腾区。在此区域，可以采用垂直上升内螺纹管内两相流体沸腾传热系数关联式

$\frac{{\mathrm{\α}}_{\mathrm{tp}}}{{\mathrm{\α}}_l}=0.8936\·\frac{1}{{X_{\mathrm{tt}}}^{0.2627}}\·{\left(\frac{P}{P_{\mathrm{cr}}}\right)}^{-1.3889}\·{\left(\frac{G}{2000}\right)}^{-0.6731}---\left(97\right)$

上式的适用范围：压力P＝9.0～19.0MPa。

当汽液两相均为紊流时，参数X_tt为：

$X_{\mathrm{tt}}={\left(\frac{1-x}{x}\right)}^{0.9}\·{\left(\frac{{\mathrm{\ρ}}_g}{{\mathrm{\ρ}}_l}\right)}^{0.5}\·{\left(\frac{{\mathrm{\μ}}_l}{{\mathrm{\μ}}_g}\right)}^{0.1}---\left(98\right)$

C 干涸后传热区

垂直上升内螺纹管在亚临界压力的单相区和两相区都有显著得的强化传热，推迟传热恶化和降低恶化的后果的作用。西安交通大学进行实验研究发现，垂直上升内螺纹管中发生的传热恶化是第二类传热恶化，也即干涸。当处于干涸后的传热区，可以采用西安交通大学所给出的垂直上升内螺纹管干涸后传热系数计算关联式

${\mathrm{Nu}}_g=0.012\·{\{{\mathrm{Re}}_g\·[x+\frac{{\mathrm{\ρ}}_g}{{\mathrm{\ρ}}_l}\·(1-x)]\}}^{0.86}\·{{\mathrm{Pr}}_{g,w}}^{2.1}\·q^{0.95}\·{\left(\frac{k_g}{k_{\mathrm{cr}}}\right)}^{-1.3}\·{\left(\frac{P}{P_{\mathrm{cr}}}\right)}^{-0.98}$

上式的适用范围：压力P＝9.0～19.0MPa        (99)

2)垂直管圈水冷壁超临界及超超临界压力区

在超临界及超超临界压力区，内螺纹管的传热特性在不同焓值区域不同，在远离临界焓值的低焓值区和高焓值区，呈现相似的规律：管壁温度随焓值增大而增大，但处于低焓值区的超超临界水的换热强于处于高焓值区的超超临界汽，相同压力、质量流速和壁面热负荷条件下，超超临界水的传热系数高于超超临界汽的传热系数。内螺纹管内壁与流体之间的换热系数可根据管内流体的焓值，分别选用关联式

在拟临界温度前的低焓值区：

${\mathrm{Nu}}_w=0.7469\·{(\frac{h_w-h_f}{t_w-t_f}\·{\mathrm{Pr}}_w)}^{1.1109}\·{{\mathrm{Re}}_w}^{0.459}\·{\left(\frac{v_f}{v_w}\right)}^{0.468}---\left(100\right)$

在拟临界温度后的高焓值区域：

${\mathrm{Nu}}_w=47.39\·{(\frac{h_w-h_f}{t_w-t_f}\·{\mathrm{Pr}}_w)}^{1.183}\·{{\mathrm{Re}}_w}^{0.203}\·{\left(\frac{v_f}{v_w}\right)}^{1.696}---\left(101\right)$

上2式中的下标w和f分别表示以壁面温度和工质温度为定性温度；其中Re数的特征尺寸取内螺纹管的当量内径。式(100)和式(101)的平均相对误差分别为16.96％和19.34％。其适用范围：压力P＝24.5～32.0MPa，质量流速G＝400～1500kg/m²·s，热负荷q＝200～600kW/m²。

3)螺旋管圈水冷壁亚临界压力区

A.倾斜上升内螺纹管内单相流体对流换热系数

采用Dittus-Boelter公式的的形式，对试验数据进行拟合，得到亚临界压力区倾斜上升内螺纹管顶部内壁与单相流体之间的传热系数关联式如下

${\mathrm{Nu}}_f=0.0218{\mathrm{Re}}_f^{0.792}{\mathrm{Pr}}_f^{0.561}---\left(102\right)$

式中下标f表示流体参数。关联式与试验数据之间的均方根误差为12.6％。关联式的适用范围为：压力p＝9～22MPa，质量流速G＝600～1200kg/(m²s)，热负荷为q＝200～500kW/m²。

B.倾斜上升内螺纹管内两相流体传热系数

根据实验数据拟合得到倾斜上升内螺纹管顶点内壁与两相流体之间的传热系数关联式如下：

$\frac{h_{\mathrm{tp}}}{h_l}=96.8{[B_0+0.00015{\left(\frac{1}{X}\right)}^{0.67}]}^{-0.0168}{q}^{-0.196}{G}^{-0.534}{\left(\frac{{\mathrm{\ρ}}_i}{{\mathrm{\ρ}}_g}\right)}^{0.679},---\left(103\right)$

上式与试验数据之间的均方根误差为14.6％。上式的适用范围为：压力p＝9～22MPa，质量流速G＝600～1200kg/(m²s)，热负荷为q＝200～500kW/m²。

式中： $X={\left(\frac{1-x}{x}\right)}^{0.9}{\left(\frac{{\mathrm{\ρ}}_g}{{\mathrm{\ρ}}_l}\right)}^{0.5}{\left(\frac{{\mathrm{\μ}}_l}{{\mathrm{\μ}}_g}\right)}^{0.1}---\left(104\right)$

C.倾斜上升内螺纹管传热恶化后传热系数

用于计算传热恶化后倾斜上升内螺纹管顶点内壁与工质之间传热系数的试验关联式如下：

${\mathrm{Nu}}_g=0.0498\·{\left\{{\mathrm{Re}}_g[x+\frac{{\mathrm{\ρ}}_g}{{\mathrm{\ρ}}_l}\·(1-x)]\right\}}^{0.0215}\·{\mathrm{Pr}}_{g,w}^{4.18}\·q^{0.883}\·{\left(\frac{{\mathrm{\λ}}_{\mathrm{cr}}}{{\mathrm{\λ}}_g}\right)}^{1.22}---\left(105\right)$

式中：Pr_g，w表示以内壁温度为定性温度得到汽相Pr数；

λ为蒸汽导热系数；

λ_cr＝0.914kW/(m.K)为热力学临界点的导热系数。

上式与试验数据之间的均方根误差为12.6％。上式的适用范围为：压力p＝9～22MPa，质量流速G＝600～1200kg/(m²s)，热负荷为q＝200～500kW/m²。

4)螺旋管圈水冷壁超临界及超超临界压力区

拟合试验数据，得到如下超临界压力区倾斜上升内螺纹管顶部内壁与工质之间传热系数的关联式。

${\mathrm{Nu}}_b=0.0198{{\mathrm{Re}}_b}^{0.783}{{\mathrm{Pr}}_b}^{0.022}{\left(\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_b}\right)}^{-0.106}{\left(\frac{{\mathrm{\λ}}_w}{{\mathrm{\λ}}_b}\right)}^{0.918}---\left(106\right)$

式中：下标b为流体参数；

下标w为管壁参数；

上式与试验数据之间的均方根误差为10.4％。上式的适用范围为：压力p＝23～28MPa，质量流速G＝600～1200kg/(m²s)，热负荷为q＝200～500kW/m²。

4 膜式水冷壁温度场二维计算模型的程序步骤

对管壁温度进行计算时，由于内壁放热系数不仅与“当地”流体温度有关，还与内壁温有关，因此是一迭代过程，先假设一个内壁温度，由此计算放热系数，再求内壁温度，当前次内壁温度与当次内壁温度相差很小时，即终止迭代。膜式水冷壁温度场分布求解步骤如下：

(1)假设一个温度场分布；

(2)根据给定的温度场分布，得到内壁换热系数以及计算区域内各点上的导热系数；

(3)进行全场迭代求解，求出新的温度场分布；

(4)将算得的新温度场代替旧温度场，重复步骤(2)和(3)。直至相邻两次迭代的节点最大温度差小于规定精度，就可以认为迭代收敛，也即所得温度场为真实温度场。

(5)将同时求得的环形区域与矩形区域耦合求出膜式水冷壁全场壁温分布。

三 受热面应力在线监测

应力在线监测只需监测实时压力，并通过下式计算求得。应力误差对寿命的影响一般远小于温度对寿命的影响。

σ_p＝K_PPD/2S_S      (107)

式中，P为工作压力，MPa；D为管子原始中径(公称外径与内径平均值)，mm；S_s为管子实际壁厚，mm；K_p为综合减弱系数和应力集中系数。

${\mathrm{\σ}}_T=K_t\frac{E_a\·\mathrm{\ΔT}}{1-v}[1/\ln (R_2/R_1)-\frac{{2R}_2^2}{R_2^2-R_1^2}]---\left(108\right)$

式中：E——当时温度下的弹性模数；

α——当时温度下的微分线形膨胀系数；

v——泊桑比；

R₂——外径；

R₁——内径；

K_t——热应力集中系数。

实际运行中，由于腐蚀和磨损，不可避免地存在管子壁厚减薄的现象，为此：

S_s＝S_b-S_w-S_n-S_m    (109)

式中，S_b为原始壁厚，mm；S_w为外壁腐蚀减薄厚度，mm；S_n为内壁腐蚀减薄厚度，mm；S_m为飞灰磨损减薄厚度，mm。在线监测***中，S_w、S_n、S_m采用水管锅炉强度计算标准推荐的最大值，并根据锅炉检验实测壁厚定期校核修正。

四 疲劳损伤的在线寿命评定

锅炉承压部件在运行期间受到的反复交变应力，主要来源于机组的启停以及变负荷运行时压力和温度的变化和波动。目前常用按疲劳试验破坏前材料所经历的载荷循环周次的多少对疲劳进行分类。通常，把破坏周次大于10⁴～10⁵次的称为高周疲劳，低于此界限的称为低周疲劳。锅炉设备在启停和变负荷运行中产生的热疲劳多属于低周疲劳。

疲劳寿命具有统计特性，因此疲劳寿命估算多是在通过大量试验所建立的疲劳设计曲线的基础上进行的。传统的试验方法是在一系列的循环载荷下，测得无裂纹光滑试件的相应的断裂次数，由此获得应力与寿命的关系(σ-N)曲线或应力与寿命(ε-N)曲线，统称为S-N曲线。对部件进行疲劳寿命评价时，首先统计出部件承受的各种循环应力的幅值，再根据S-N曲线找出对应的断裂周次，由下式即得到部件的累积损伤程度：

$D=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{n_i}{N_i}---\left(110\right)$

式中：n_i为各种循环应力的循环次数；N_i为相应循环应力的允许循环次数，当D＝1时，即认为发生了疲劳破坏。

五 蠕变寿命损伤的评定

在高温蠕变条件下，应力越大，蠕变进行的越快，破坏得越早。试验表明：在一定温度条件下，应力σ与蠕变破坏时间τ存在如下关系：

τ＝Aσ^-B               (111)

式中，A、B为与钢种及温度有关的常数，可通过恒温条件下，采用较高的应力进行短期的试验获得。

采用经典的蠕变计算公式(L-M公式)，对于某应力，绝对温度T、蠕变破坏时间τ和应力σ有如下关系：

T(C+lgτ)＝a₀+a₁lgσ+a₂lg²σ+a₃lg³σ    (112)

对于给定材料，可根据其实测对应不同温度和时间的材料持久强度数据求出系数a₀、a₁、a₂、a₃和常数C。

${\mathrm{\Φ}}_r=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\Δ\τ}}_i}{{\mathrm{\τ}}_i}---\left(113\right)$

φ_r——材料的蠕变寿命损耗量

⊿τ_i——在i参数下部件的运行时间

τ_i——在i参数下部件的蠕变断裂时间

根据以上分析，由材料的工作温度和工作应力，可以确定部件在该工作温度和工作应力下的使用寿命。在线监测中，通过实时计算管壁温度和应力，可得到寿命损耗率。对于蠕变损伤，满足线性累加原则，即可将寿命损耗率累加，得到总的寿命损耗。

六 蠕变和疲劳的交互作用

在电厂锅炉承压部件中，材料都会受到蠕变和疲劳的交互作用，材料的软化会降低材料的蠕变性能和低周疲劳性能，脆化会降低材料的冲击韧性值和断裂韧性，从而严重降低锅炉承压部件的致裂寿命和剩余寿命。而蠕变-疲劳交互作用，又会进一步缩短锅炉承压部件的致裂寿命和剩余寿命。从上面的介绍已经知道了材料蠕变和疲劳的寿命损耗计算方法，它们都是在单独作用情况下得出的。

在评价锅炉承压部件蠕变一疲劳交互作用下损伤时，通常采用线性累计损伤法，即应用蠕变损伤的Robbinson法则和疲劳损伤的Miner准则，并认为蠕变与疲劳独立发生作用，其总损伤Φ为：

$\mathrm{\Φ}={\mathrm{\Φ}}_r+D=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\Δ\τ}}_i}{{\mathrm{\τ}}_i}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{n_i}{N_i}---\left(114\right)$

式中Φr和D分别为为蠕变损伤和疲劳损伤；Δτ_i为第i种工作条件下累积运行时间；τ_i为第i种工作条件下金属部件蠕变断裂时间。n_i为在应力σi下经过n_i反复循环作用；Ni为在应力σ经过Ni次反复循环作用材料发生裂纹。

Claims (3)

1、一种水冷壁在线失效预警方法，其特征在于，首先，根据电站锅炉结构、运行数据划分水冷壁回路，结合运行时水冷壁总流量、进口压力数据，进行水动力计算，得到各回路流量、压力数据；

其次，在结构数据及各回路流量、压力数据基础上，采用数值方法，进行金属壁温计算，获得整个鳍片管的温度场；

再次，在已得到水冷壁各回路压力、流量、整个鳍片管壁温的基础上，根据公式σ_p＝K_PPD/2S_S得到切向薄膜应力σ_p，其中P为工作压力，MPa；D为管子外径与内径平均值，mm；S_s为管子实际壁厚，mm；K_p为综合减弱系数和应力集中系数；

根据公式 ${\mathrm{\σ}}_T=K_t\frac{E_{\mathrm{\α}}\·\mathrm{\ΔT}}{1-v}[1/\ln (R_2/R_1)-\frac{2R_2^2}{R_2^2-R_1^2}]$ 得到切向热应力σ_T。其中，E为当时温度下的弹性模数；α为当时温度下的微分线形膨胀系数；v为泊桑比；R₂为外径；R₁为内径；K_t为热应力集中系数；

最大应力值由每次循环中的最大切向薄膜应力σ_p和最大切向热应力σ_T相迭加而得，最小应力值则由最小切向薄膜应力σ_p和最小切向热应力σ_T相迭加而得，从而形成应力履历数据库，为疲劳损伤计算提供应力振幅，为蠕变损伤计算提供有效蠕变应力数据；

最后，在应力履历数据库的基础上，结合试验所建立的疲劳设计曲线S-N曲线，先统计出部件承受的各种循环应力的幅值，再根据S-N曲线找出对应的断裂周次，由公式 $D=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{n_i}{N_i}$ 得到部件的累积损伤程度，当D≥1时，判定水冷壁管因疲劳发生失效。其中，n_i为各种循环应力的循环次数；N_i为相应循环应力的允许循环次数；同时进行蠕变寿命损伤的评定，对于某应力，绝对温度T、蠕变破坏时间τ和应力σ采用公式T(C+lgτ)＝a₀+a₁lgσ+a₂lg²σ+a₃lg³σ计算，得到蠕变破坏时间τ。对于给定材料，可根据其实测对应不同温度和时间的材料持久强度数据求出系数a₀、a₁、a₂、a₃和常数C，再根据公式 ${\mathrm{\Φ}}_r=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{\mathrm{\τ}}_i}{{\mathrm{\τ}}_i}$ 进行计算，当Φ_r≥1时，判定水冷壁管因蠕变发生失效，其中，φ_r为材料的蠕变寿命损耗量，为在i参数下部件的运行时间，τ_i为在i参数下部件的蠕变断裂时间；若D≤1且Φ_r≤1，但φ＝D+Φ_r≥1时，判定水冷壁因疲劳、蠕变失效，实现水冷壁在线失效预警。

2、根据权利要求1所述的一种水冷壁在线失效预警方法，其特征在于，水动力计算，通过逐步计算，获得个回路，包括节流元件、局部阻力元件、水冷壁，不同流量下的压降，这样获得回路分散点的数据，采用四次方切比雪夫多项式拟合，在以获得的压降-流量数据基础上，得到回路的水动力曲线方程，在串联***中，根据工质流量相等，压降叠加，在并联***中，根据压降相等，流量叠加来获得串、并联回路的压降-流量或流量-压降数据，据此通过切比雪夫多项式拟合获得串并联***的水动力曲线方程，对于直流锅炉，入口总流量M_sum总是已知的，即可求得压降ΔP。

3、根据权利要求1所述的一种水冷壁在线失效预警方法，其特征在于，所述的金属壁温计算采用区域分离法建立了膜式水水冷壁稳态无内热源而且变物性的二维温度场计算模型，获得整个鳍片管的温度场，用于计算分析膜式水冷壁界面温度场分布、水冷壁管壁温沿炉膛高度的分布和炉膛同一截面各水冷壁管的温度分布。

Images