张振国，李小畅，谭思超,*，田瑞峰，刘思超，张凤奎

(1.哈尔滨工程大学 核科学与技术学院，黑龙江 哈尔滨 150001；
2.哈尔滨工程大学 航天与建筑工程学院，黑龙江 哈尔滨 150001)

直流蒸汽发生器作为钠冷快堆关键设备，一方面承担钠、水两侧的热交换，另一方面是钠侧回路的压力边界，其热工水力特性直接影响快堆电站的安全性与经济性。快堆运行中，蒸汽发生器在流体激振力与热负荷的共同作用下，蒸汽发生器易诱发流致振动、流弹失稳、热应力、泥渣沉积及由此引起的微动磨损、疲劳失效、腐蚀等传热管失效甚至破裂问题。对换热管失效因素进行预测，需对蒸汽发生器开展热-流-固多物理场耦合分析，蒸汽发生器的三维温度场及流场则是必不可少的输入条件。因此，对蒸汽发生器开展三维精细化热工水力数值分析，不仅可为蒸汽发生器热-流-固耦合分析提供输入条件，也能为结构设计与力学校核奠定基础，同时也可指导实验研究。

蒸汽发生器结构复杂，包含数百甚至数千根细长换热管束，管束尺寸跨度大，且轴向存在数十个截面形状复杂的支撑结构，对全尺寸蒸汽发生器进行试验存在难度大、周期长、成本高且需投入巨大的人力物力等问题，诸多关键局部参数也难以通过试验获取。数值计算作为一种高效且相对低成本的分析手段，可很好弥补全尺寸试验的不足。由于蒸汽发生器具有几何结构和物理现象的双复杂特征，早期针对蒸汽发生器的数值分析主要采用一维程序，但该方法无法提供蒸汽发生器局部流场及温度场信息，分析结果较粗糙。近年来，随着计算流体力学、数值分析及多相流理论的发展，数值模拟越来越多应用于全流型管内沸腾换热数值模拟，如单流体多相流模型[1]、两流体多相流模型[2]及“两流体+三流场”模型[3-4]等CFD分析方法。尽管CFD方法很好弥补了一维程序分析精细度不够的问题，但存在一共同的问题：由于采用了复杂的三维多相流数值模型，对包含大规模换热管束且存在两侧耦合换热、多相流动及相变换热等复杂物理现象的蒸汽发生器进行三维精细化分析时，普遍存在计算量大、分析稳定性差及分析效率低等问题，难以实现全尺寸蒸汽发生器的三维精细化数值分析。因此，现有针对全尺寸蒸汽发生器的CFD分析大多采用多孔介质模型[5-10]。多孔介质模型忽略换热管及支撑结构的物理边界，将管侧以热源的形式予以考虑，支撑板及换热管本身对壳侧热工水力特性的影响，则用多孔介质模型予以考虑。多孔介质模型很好解决了全三维精细化CFD分析数值计算成本高、效率低、收敛性差等问题。然而，多孔介质模型也面临其他问题：无法很好模拟局部结构对流场及温度场的影响；
无法精细得到每根换热管内、外壁面温度，无法提供同截面壳、管两侧及支撑结构的温度梯度，因此其分析结果难以为热-流-固耦合分析提供精准的输入条件，难以进行堵管工况分析。

本文针对管内沸腾的钠冷快堆直流蒸汽发生器，建立壳侧三维CFD模型与管侧一维沸腾传热模型的跨尺度实时耦合分析模型，以实现每根换热管与壳侧的独立耦合传热分析，既精细获得壳侧三维精细化流场与温度场、每根换热管内外壁及管侧流体温度的三维精细化分布，又克服一维程序、全三维精细化数值模拟及多孔介质模型中的各种问题。

本文建立的跨尺度耦合方法主要针对管内沸腾的直流蒸汽发生器三维温度场的求解，核心思路为壳、管两侧传热的实时跨尺度耦合。跨尺度耦合指将壳侧三维精细化CFD流场与管侧每根换热管的一维沸腾传热模型进行耦合，实现三维温度场的精细化求解。这种方法的好处在于避免了直接采用复杂的三维多相流模型，将计算简化为单相流问题，大幅提高计算量；
另一方面可获得壳侧三维空间任意点的热工水力参数、水侧每根换热管任意高度处的热工水力参数，同时可获取任意换热管内外壁面及管壁固体域温度。

1.1 壳侧数学模型

壳侧流体为液态金属钠，其流动与传热遵循质量守恒、动量守恒及能量守恒，分别由连续方程、动量方程及能量方程来描述。将壳侧液态金属钠视为不可压单相流体，其稳态控制方程如下。

连续性方程：

(1)

动量方程：

(2)

能量方程：

(3)

对于实际工程问题，上述控制方程组无法直接求解，需对方程组进行时间平均处理。时均化后的控制方程组中，动量方程将出现包含脉动速度分量乘积的雷诺应力项，能量方程中则出现包含脉动速度与脉动温度乘积的湍流热通量项，因此需对雷诺应力项与湍流热通量项进行封闭。目前对于雷诺应力项的封闭主要采用涡黏性湍流模型，湍流热通量项则基于湍流普朗特数模型。本计算中，壳侧流体为液态金属钠，在快堆实际运行工况下，其普朗特数低至10-2量级，属于典型的低普朗特数流体，这意味流体流动与传热过程中，导热占主导地位。目前常见的两方程涡黏模型处理湍流热通量项时，大多将湍流热通量项合并至扩散项，并引入湍流黏性系数μt和湍流普朗特数Prt重新定义有效导热系数Keff。但在传统数值处理方法中，将湍流普朗特数Prt取常数0.85的方法已不再适用导热占主导地位的低普朗特数金属流体。已有研究表明，当湍流模型中的湍流普朗特数同时考虑流体分子黏性系数及湍流黏性系数的影响时，对低普朗特数流体具有相对更好的计算效果[11]。ANSYS FLUENT求解器中的RNGk-ε模型具备上述特点，因此本文选择RNGk-ε湍流模型结合壁面函数的方法模拟壳侧液态金属钠的湍流流动与传热。鉴于壳侧液态金属钠与管外壁的对流换热系数对整个蒸汽发生器三维温度场有显著影响，本文未使用CFD中的近壁面处理方法，如壁面函数或低雷诺数模型，而是将液态金属钠平行流过管束的换热公式嵌入三维CFD模型中求解液态钠与换热管外壁面的对流换热。对于液态金属平行流过换热管束的对流换热，文献[12]给出了8个推荐关系式，本文采用对液态金属钠具有更好适应性的Subbotin关系式[13]：

(4)

式中：Nu为努塞尔数；
E为换热管节距，m；
Pe为佩克莱数；
L为换热管长，m；
Dh为壳侧水力直径，m；
Do为换热管外径，m。

1.2 管侧数学模型

管侧沸腾传热计算模型用于求解水侧不同换热管不同高度处的热工水力参数，以实现与钠侧三维精细化CFD模型间的跨尺度耦合。本文对水侧流动沸腾的求解采用“四段式”，即将换热管内的流型分成4个区域，各区域采用不同的换热模型及热工水力参数确定模式。计算中，将全过程沸腾传热流型分成4个区域，分别为单相液换热区、泡核沸腾换热区、膜态沸腾换热区及过热蒸汽换热区。全过程沸腾传热及流型判别方法如图1所示。

图1 管侧传热模型与流型判别准则Fig.1 Heat transfer model and discrimination criteria of flow regime at tube side

如图1所示，各流型区判别准则如下：当水侧流体比焓h小于饱和液体比焓时，管侧处于单相液对流换热区，采用Sieder-Tate传热关系式；
当水侧流体比焓超过饱和蒸汽比焓时，仍采用Sieder-Tate传热关系式；
当水侧比焓大于饱和水比焓且小于饱和蒸汽比焓时，管侧处于两相换热区，此时水侧流体温度为饱和温度。两相换热区包括泡核沸腾区与膜态沸腾区，本文采用较成熟的临界质量含汽率来差别泡核沸腾与膜态沸腾的分界点，即干涸点处的质量含汽率。

两相段包括泡核沸腾区及膜态沸腾区，二者传热机理及传热模型均不相同，因此需采用合适的方法识别每根换热管中泡核沸腾与膜态沸腾的临界点。本文计算过程中，采用较成熟的临界质量含汽率判别泡核沸腾与膜态沸腾的分界点，即干涸点处的质量含汽率。干涸点处的质量含汽率采用古塔杰拉奇关系式[14]：

(5)

式中：xcrit为临界质量含汽率；
xcrit,D8为直径8 mm的换热管的临界质量含汽率；
D为换热管内径，mm；
G为质量流率，kg/(m2·s)；
ρ为密度，kg/m3；
μ为动力黏度，Pa·s；
下标l与v分别表示液相与气相。古塔杰拉奇关系式的适应范围：p=9.8～16.66 MPa；
qω=750～3 000 kg/(m2·s)。

管内过冷区单相液与过热区单相气的强迫对流换热，均采用Sieder-Tate关系式[15]。该关系式是Sieder与Tate两位学者于1936年在Dittus-Boelter关系式的基础上提出，能更好兼顾较大温差工况：

(6)

式中：Re与Pr分别为水侧雷诺数与普朗特数；
μb与μw分别为以水侧主流区及壁面为温度定性的黏性系数，Pa·s。Sieder-Tate关系式的适应范围：0.7≤Pr≤16,Re≥104,L/dw≥60。

泡核沸腾换热区水侧流体温度为饱和温度，此区域换热系数采用Chen关系式[16]，该式能预测欠热泡核沸腾及饱和核态沸腾，具有较好的工况适应性。该传热模型将壁面有效热流量表示成核态沸腾热流量与强迫对流热流量之和。主要基于两个假设，一是壁面蒸发过程中增加了液相流速，因此强迫对流换热能力相较于纯液相单相对流换热能力所有增强；
二是液相对流换热过程对泡核沸腾换热能力存在一定的抑制作用。Chen关系式具体形式如下：

hTP=hNB+hFC

(7)

其中：hTP为泡核沸腾区综合换热系数；
hFC与hNB分别为强迫对流换热系数与核态沸腾换热系数，如式(8)所示，单位均为W/(m2·K)。

(8)

其中，F为雷诺数因子：

(9)

Xtt为Lockhart-Martinelli参数：

(10)

S为泡核沸腾抑制因子：

(11)

其中，Retp的计算式如下：

(12)

式中：p为压强，Pa；
x为质量含汽率；
cp为比定压热容，J/(kg·K)；
k为导热系数，W/(m·K)；
下标s表示饱和状态。

膜态沸腾换热区水侧温度仍然为系统压力下的饱和温度，其传热模型采用Miropolskiy关系式[17]：

Nu=

(13)

该式基于经典单相对流Dittus-Boelter关系式，通过引入修正系数考虑了含汽率的影响。

基于以上数学模型，本文充分利用沸腾相变传热发生在细长型换热管中这一特点，将沸腾传热过程视为一维问题。根据沸腾传热关系式及流型判别方法建立每根换热管的一维沸腾传热求解模型与求解程序，通过Fluent二次开发接口将管侧一维沸腾传UDF热求子程序与壳侧三维精细化CFD模型建立实时耦合求解模型。

模型验证基于印度19管SGTF钠冷快堆模型实验[18]，取其中的30°区域作为分析对象，SGTF模型、换热管编号及计算区域如图2所示。由于本文主要给出跨尺度耦合方法的建模思路及求解可靠性，后续数值建模及计算中，未考虑支撑件的影响。该实验模型提供了15及16号换热管不同高度处近壁面液态钠介质温度测量数据。SGTF模型几何参数及运行工况参数列于表1。

图2 SGTF蒸汽发生器模型及计算区域Fig.2 Sectional view of SGTF model with indicating domain of computation

表1 SGTF模型参数Table 1 Parameter of SGTF model

对于19管SGTF模型，当径向基准网格尺度小于0.2 mm时，可避免网格尺寸对计算结果的影响[18]。考虑到液态金属钠的低普朗特数特性，本文选择ANSYS FLUENT中对低普朗特数具有更好计算效果的RNGk-ε模型，且保证换热管壁面y+处于30～50之间，此时壁面第1层网格高度为0.36 mm，轴向网格数取650。所建立的网格模型如图3所示，包括钠侧流体域及换热管固体域，两个侧边界为对称边界。由于固体域为纯导热，且管壁较薄，温度近似线性分布，因此建立一层网格即可保证精度。计算过程中，水侧所有参数储存到固体域单元，通过换热量实现钠、水两侧跨尺度耦合。由于本文直接对换热管固体域进行了三维建模，因此换热管固体域导热及热容对换热量与温度场的影响直接在耦合计算中体现。

图3 计算域网格Fig.3 Mesh of computational domain

数值求解时，首先采用定热流密度边界条件，迭代求解约20步后加载水侧沸腾传热计算程序进入耦合求解。图4为直流蒸汽发生器钠侧与水侧跨尺度耦合求解过程收敛性监测曲线，分别监测了4根换热管的水侧出口温度及壳侧液态金属钠出口温度。耦合过程中，采用了耦合传递参数的欠松弛方法，以保证计算收敛性。从图4可看出，耦合初期，水侧出口温度出现了大幅度振荡，随着迭代进行，振荡逐渐减小，约350个迭代步后钠、水两侧出口温度即达到稳定状态，耦合求解具有良好的收敛性与稳定性。从建模及求解过程来看，建立的跨尺度耦合方法的计算量与单相流动与传热问题相当，具有很高的求解效率。表2列出了钠侧及水侧出口温度的计算值与设计值的对比，钠侧与水侧出口温度计算结果与设计值的误差均小于1 K。

图4 耦合计算收敛性监测Fig.4 Convergence monitoring of coupling solution

表2 钠侧及水侧出口温度的计算值与设计值的对比Table 2 Comparison of sodium side and water side outlet temperatures between designed values and calculated values

SGTF钠冷快堆模型实验给出了部分换热管轴向不同位置的钠侧温度。图5对比了本文计算的换热管15及换热管16外管壁附近液态金属钠温度与文献[18]计算结果及实验数据的对比，图中实验数据误差棒为热电偶的测量误差(±1%)[18]。该实验热电偶安装在换热管外壁附近，但未给出热电偶的具体安装位置，因此本文计算值取对应换热管同一高度处外管壁第1层网格温度数据的体积加权平均值。根据本文计算结果，同一换热管同一截面处的温差小于10 K，因此基于该处理方法的对比分析仍能验证计算结果的可靠性。从图5看，本文计算结果与实验数据吻合良好，尤其在轴向位置9 m附近进入膜态沸腾时的温度小幅度跳跃现象也与实验较吻合。对于换热管15，本文计算结果与实验数据的最大相对误差约为2.6%，位于轴向位置7.46 m处；
对于换热管16，本文计算结果与实验数据的最大相对误差约为2.4%，位于轴向位置8.79 m处，可见本文所计算的轴向温度与实验数据之间的最大相对偏差小于3%。本文计算结果相比于文献[18]的计算结果，在轴向分布趋势及数值上基本吻合，但对于换热管16，在非实验数据区存在一定偏差，文献[18]的计算结果更早进入膜态沸腾区。

图5 换热管近壁面钠温与实验数据的对比Fig.5 Comparison of sodium temperatures near wall of tubes with experimental data

图6为4根换热管水侧总吸热量随高度的变化，随着高度的增加，吸热量逐渐增加。由图6可看出，对于本计算模型，换热管15、16、19由于换热面积相同，总吸热量基本相当，换热管16的吸热量略高于15与19。换热管10的换热面积为其他换热管的1/6，其总吸热量最低。

图6 换热管吸热量随高度的分布Fig.6 Heat transfer rate distribution along tube length

图7为换热管10、15及19水侧温度、钠侧温度、管内壁温度、管外壁温度及水侧换热系数(HTC)随高度的变化，其中图7d标注了各换热区域。从图7可看出，对于本计算模型及工况，约在6.5 m处进入两相区，约在9 m处由泡核沸腾进入膜态沸腾，约在13.5 m处进入单相汽区域。由泡核沸腾转为膜态沸腾时，水侧进入缺液区，管内、外壁及钠侧温度均出现小幅度跃升。这种温度跃升相较于典型的单管定热流密度流动沸腾明显更小，主要是由于两侧耦合换热时，一侧的温度始终受另一侧温度的限制，这也是两侧耦合换热的特征之一。各换热管在单相液及泡核沸腾区域由于换热系数相对较大，水侧温度与管内壁温度差更小，膜态沸腾及单相汽区域则由于换热系数较小，温差更大，与物理事实相符。从图7d看，由于本计算模型所取换热管规模较小，各换热管钠、水两侧温度及水侧换热系数均十分接近。计算结果表明，各换热管内壁平均换热系数约为12 500 W/(m2·K)，泡核沸腾区换热系数最大，平均换热系数达到约32 500 W/(m2·K)，单相汽区换热系数相对较小，至管出口时降至约5 000 W/(m2·K)。

图7 不同换管温度及水侧换热系数随高度的变化Fig.7 Temperatures and water side heat transfer coefficients along length of each tube

图8为不同高度处横截面钠侧温度分布，图中z为轴向坐标，ΔT为截面最大温差。从图8可看出，各截面温度分布较相似，表现为四周温度高、中间温度低，这主要是由于四周较大空隙中的液态钠没有足够的换热管带走热量，即四周单位质量钠介质均摊的换热管较少，换热量较小，从而表现出高温。从各截面最大温差来看，即使本文仅选取了4根换热管模型，最大温差可达54 K，最小温差约为8 K。同时，越靠近钠侧出口截面温差越大，越靠近钠侧入口截面温差越小。这主要是因为钠侧入口温度为均匀分布，接近钠侧入口的截面处液态金属钠来不及与水侧换热，换热量较小，因此截面温差也较小。这也说明，如果工程上特别关心截面温差引起的结构热应力，应重点关注靠近钠侧出口位置的温度不均匀性。

图8 钠侧不同轴向位置横截面温度分布Fig.8 Temperature distribution at various sections of sodium side

蒸汽发生器在运行过程中，换热管长期受微动磨损、冲刷、腐蚀等作用，易造成传热管失效，因此工程上常通过检修对有失效风险的换热管进行封堵，避免出现安全事故。本文所建立的直流蒸汽发生器两侧跨尺度耦合分析模型，最显著特点在于可获得任意换热管的热工水力参数，也易于实现任意换热管堵管工况的模拟。本文利用建立的跨尺度耦合分析程序，对上述4管模型进行了堵管工况分析。共计算了4组堵管工况，分别为封堵传热管10、15、19及同时封堵传热管15与19，分析其热工水力特性的变化。在处理堵管工况时，直接将被封堵的传热管内壁面热边界条件给定为零热流密度条件，其他换热管则为耦合换热条件。分析堵管工况时，其他换热管水侧单管流量与正常工况相同，钠侧总流量不变。

表3为4组堵管工况下宏观计算结果与正常工况的对比。由于本文计算对象为19管模型中的30°扇形角，因此堵管10工况对质量流量、换热量及钠水两侧的出口温度影响最小。随着堵管数量的增加，总换热量减小，水侧与钠侧出口温度均逐渐增大。单独封堵15或19号换热管时，堵管对换热管束的影响基本相同，钠、水两侧出口温度及换热量均相近，此时传热面积减小31.6%，热功率减小29.9%。当同时封堵15与19号管时，水侧出口温度达到793.2 K，钠侧达到732.4 K，换热量也最低。需要说明的是，本文所计算的堵管工况，主要验证跨尺度耦合模型对堵管分析的适应性，同时分析堵管对局部热工参数的影响规律，而在实际工程中，当堵管数量接近设计裕量上限时，两侧流量均需做相应调整。

表3 堵管工况与正常工况计算参数对比Table 3 Comparison of values between plugged and normal conditions

图9为4组堵管工况下钠、水两侧温度沿轴向高度的变化与正常工况的对比。由于被封堵换热管无水侧数据，因此图中各堵管工况按以下原则选取换热管作为分析对象：钠侧温度取被封堵换热管及其相邻换热管数据，水侧温度选取被封堵换热管的相邻换热管数据，以更好分析封堵对该换热管本身及其邻近换热管钠、水两侧温度的影响。从图9可看出，被封堵换热管附近钠温相比于正常工况，温度更高且变化更平缓，且由于热量损失大幅降低，进出口温差也减小。被封堵换热管的邻近换热管由于周围钠温更高，壳、管两侧温差也相应增大，吸热量增加，因此更进早入两相区，出口温度也明显升高。当某根换热管周围存在两根堵管时，其吸热量急剧增加，出口温度显著升高。计算结果表明，本文建立的跨尺度耦合方法可得到堵管工况下任意换热管钠水两侧及管壁固体域的详细温度场数据，作为堵管工况下蒸汽发生器结构热应力分析评估的输入条件。

图10为4组堵管工况下，轴向高度15 m处钠侧横截面温度分布与正常工况的对比。由图10可知，随着堵管数量的增加，同一截面处钠侧温度也更高，这主要是由于钠侧放热量减少所致。尽管从图中来看，随着堵管数量的增加，截面温差明显降低，但这并不能说明实际蒸汽发生器中的堵管处理也会导致截面温差减小。主要是因为，实际蒸汽发生器换热管数量众多，堵管数量较少，堵管仅对局部温度场有明显影响，对整个截面温度场分布影响有限。堵管位置及堵管数量对实际蒸汽发生器换热量及温度分布的影响，需开展全尺寸级蒸汽发生器的三维精细化分析。

图10 不同堵管工况的截面温度分布对比Fig.10 Comparison of temperature distribution under different plugged conditions

本文采用跨尺度耦合方法，建立了钠冷快堆蒸汽发生器换热管束钠、水两侧耦合换热的三维精细化数值求解模型，并通过实验数据进行了验证，同时开展了堵管工况分析，得到如下结论：

1) 采用管侧一维全流型沸腾传热模型与壳侧三维CFD模型的跨尺度耦合方法，可获得壳侧三维精细化热工水力参数及任意换热管管壁及水侧热工水力参数分布，显著提高了三维精细化数值求解效率，且易于实现任意换热管的堵管分析；

2) 与设计值及实验数据的对比分析表明，本文所建立的跨尺度耦合方法具有良好的可靠性，能够准确预测全流型沸腾换热的各传热特征点，出口宏观温度计算值与设计值的偏差小于1%，局部温度与实验数据的偏差小于3%；

3) 对于两侧逆向流动的耦合传热管束，截面最大温差更可能出现在钠侧出口附近，该区域的温度梯度引起的热应力等问题应重点关注。

本文重点介绍了管内沸腾直流蒸汽发生器的三维跨尺度耦合求解方法及其有效性，从以上分析结果来看，无论是与宏观设计值还是局部实验数据对比，该方法均表现了良好的可靠性。该方法与全三维多相流CFD方法相比，不仅大幅降低了求解方程的数量，同时也能显著提高求解效率与稳定性。对于全尺寸蒸汽发生器，壳侧支撑件等结构可考虑采用局部多孔介质方法，整个换热管束的计算可简化为裸管束通道，可实现全尺寸蒸汽发生器三维温度场的精细化数值求解。

猜你喜欢 热管蒸汽耦合 基于增强注意力的耦合协同过滤推荐方法辽宁工业大学学报(自然科学版)(2022年4期)2022-09-19擎动湾区制高点，耦合前海价值圈！房地产导刊(2021年12期)2021-12-31复杂线束在双BCI耦合下的终端响应机理北京航空航天大学学报(2021年9期)2021-11-02热管余热锅炉的应用及设计分析科学家(2021年24期)2021-04-25煤、铁、蒸汽与工业革命科学文化评论(2021年5期)2021-04-23用于大功率航天器的3D打印钛水热管设计及试验研究航天器工程(2020年4期)2020-08-14热管的讲究高效散热器不能忽略的因素电脑爱好者(2020年5期)2020-05-11车坛往事3：蒸汽时代步入正轨！人民交通(2020年2期)2020-04-16基于磁耦合的高效水下非接触式通信方法研究电子制作(2019年16期)2019-09-27蒸汽闪爆军事文摘·科学少年(2017年1期)2017-04-26