CN115862779A - 一种基于两方程热湍流模型的液态金属共轭传热计算方法 - Google Patents

Abstract

本发明公开了一种基于两方程热湍流模型的液态金属共轭传热计算方法，属于液态金属冷却反应堆热工水力计算领域，解决了现有计算方法精度不佳的问题。本发明在液态金属共轭传热计算中，引入了理论性较强、理论精度较高、理论边界较简单的液态金属两方程k_θ‑ε_θ热湍流模型，基于该模型，可以在计算液态金属共轭传热问题的同时，获得更为丰富的液态金属低普朗特数湍流换热行为和固体温度行为。本发明方法基于开源计算流体力学程序OpenFOAM保留了调用OpenFOAM内置湍流模型的功能函数，可调用不同类型的湍流模型来获得液态金属的湍流流动行为，可利用不同的湍流模型并结合液态金属两方程热湍流输运模型来改进液态金属的共轭传热计算过程。

Description

一种基于两方程热湍流模型的液态金属共轭传热计算方法

技术领域

本发明属于液态金属冷却反应堆热工水力计算领域，具体涉及一种基于两方程热湍流模型的液态金属共轭传热计算方法

背景技术

液态金属冷却反应堆(铅冷快堆和钠冷快堆)是第四代先进核能系统的重要堆型，因其冷却剂具有低压高沸点、导热性能好和中子学性能优异等特点而备受关注。在液态金属冷却反应堆中，燃料棒与液态金属冷却剂之间、换热器管壁与液态金属冷却剂之间的传热过程，通常属于大温差和变热流等严峻环境下的流体传热与固体导热的共轭传热过程，因此对上述过程需要进行精确的流-固共轭传热计算。在传统的雷诺时均(RANS)共轭传热计算方法中，通常采用简单梯度扩散假设(SGDH)将流体能量方程中的湍流热通量项转化为湍流热扩散系数项计算，然后引入雷诺比拟假设进一步将湍流热扩散系数项转化为湍流普朗特数项，进而封闭流体能量方程。然而，液态金属的导热系数较大且分子普朗特数较低，其速度边界层与温度边界层的剖面存在较大差异，导致传统的雷诺比拟假设不再适用于液态金属的共轭传热计算。

为提高液态金属RANS传热计算的精度，已有的一些改进计算方法是：1)在SGDH和雷诺比拟假设下，将基于全局流动参数或局部湍流参数拟合得到的湍流普朗特数关系式，代替恒定的湍流普朗特数，该方法能够在一定程度上提高液态金属的CFD传热计算精度，但仍受限于雷诺比拟假设的限制，难以捕捉到精细的液态金属近壁热湍流信息。2)在湍流热通量项的求解方面，不采用SGDH和雷诺比拟假设，而直接对湍流热通量项建立微分热通量模型(DHFM)方程或代数热通量模型(AHFM)方程来输运计算液态金属的湍流热通量，该方法理论性强、理论精度高，但由于该模型的输运方程数目较多、方程模化形式多样、参数敏感度较高等原因，目前在复杂几何中的应用尚待推广。3)在SGDH假设下，将湍流热通量项转化为湍流热扩散系数项之后，类比两方程k-ε湍流模型的理论方法构建两方程k_θ-ε_θ热湍流模型来输运计算湍流热扩散系数项，从而获得液态金属丰富的近壁热湍流行为信息，该方法理论性较强、实际精度较高，已在复杂的液态金属冷却反应堆燃料组件热工水力计算中得到应用。但这类两方程热湍流模型的输运形式和边界条件等开发难度较大，较难在闭源商业CFD软件中实现。

发明内容

本发明的目的是提供一种基于两方程热湍流模型的液态金属共轭传热计算方法，以解决现有计算方法精度不佳的问题。

本发明的技术方案是：一种基于两方程热湍流模型的液态金属共轭传热计算方法，包括以下步骤：

步骤1：在OpenFOAM中，基于内置流固共轭传热求解器，建立用户自定义共轭传热求解器；

步骤2：在步骤1建立的用户自定义共轭传热求解器中，建立湍流热扩散系数变量α_t和输运计算所述湍流热扩散系数的温度脉动变量k_θ及其耗散率变量ε_θ；调用热物性函数库中的温度-压力及物性变量，包括温度T、压力P、流体密度ρf、流体比热容C_pf、流体导热系数λ_f、流体分子运动黏度v_f、固体密度ρ_s、固体比热容C_ps和固体导热系数λ_s；调用湍流函数库中的速度及湍流变量，包括速度u_i、湍流运动黏度v_t和输运计算湍流运动黏度的湍动能k及其耗散率ε；

步骤3：添加液态金属两方程k_θ-ε_θ热湍流模型程序文件及相关函数头文件；

步骤4：根据实际计算问题，在用户算例文件夹下的流体物性字典文件中对步骤2中调用的流体及固体热物性变量施加合适的OpenFOAM标准物性条件；

步骤5：在用户算例初始文件夹下对步骤2中调用和建立的各物理变量施加合适的数值模拟边界条件，具体步骤如下：

步骤5-1：对以下物理变量：温度T、压力P、速度u_i、湍流运动黏度v_t、湍动能k及其耗散率ε或其比耗散率ω、湍流热扩散系数α_t调用标准的OpenFOAM壁面边界条件；

步骤5-2：对处于流-固耦合面上的温度脉动k_θ及其耗散率ε_θ施加OpenFOAM标准的零固定梯度值边界条件；

步骤6：采用OpenFOAM内置PIMPLE算法求解液态金属流体质量守恒方程(1)、流体动量守恒方程(2)和流体能量守恒方程(3)，

流体质量守恒方程(1)为：

流体动量守恒方程(2)为：

流体能量守恒方程(3)为：

其中：

h-液态金属的比焓；

K-液态金属的比动能；

gi-重力加速度；

具体步骤如下：

步骤6-1：求解流体动量守恒方程(2)，并更新速度场；

步骤6-2：利用步骤6-1更新的速度场，迭代求解流体能量守恒方程(3)，并更新温度场和流体热物性；

步骤6-3：基于步骤6-1和步骤6-2更新的数据，对PIMPLE算法下的流体质量守恒方程(1)与流体动量守恒方程(2)联立导出的压力泊松方程进行迭代求解，判断是否达到设定的PIMPLE迭代次数，如达到，则判定当前迭代步PIMPLE压力-速度耦合迭代求解结束，更新速度场和压力场后进行下一步骤；

步骤7：在用户算例文件夹下的湍流模型字典文件turbulenceProperties中，调用OpenFOAM内置湍流模型，包括但不限于：标准k-ε湍流模型、SSTk-ω湍流模型或RNGk-ε湍流模型等，然后迭代求解湍动能k及其耗散率ε或其比耗散率ω输运方程至设定的计算残差，并更新湍流运动黏度v_t；

步骤8：求解步骤3在用户自定义共轭传热求解器中建立的液态金属两方程k_θ-ε_θ热湍流输运模型至设定的残差，并更新湍流热扩散系数a_t，具体步骤如下：

步骤8-1：求解步骤3在用户自定义共轭传热求解器中建立的温度脉动k_θ及其耗散率ε_θ输运方程至设定的计算残差：

其中：

σ_kθ、σ_εθ、C_p1、C_p2、C_d1、C_d2-模型常数；

f_d2＝1/C_d2(C_ε2f_ε-1)[1-exp(-R_ε/5.7)]²-温度脉动耗散阻尼函数；

f_ε＝1-0.3exp[-(R_t/6.5)²]-湍动能耗散阻尼函数；

-温度脉动产生项；

-湍动能产生项；

步骤8-2：由步骤8-1更新后的温度脉动k_θ及其耗散率ε_θ，更新湍流热扩散系数α_t：

α_t＝C_λf_λkτ_u

其中：

C_λ、C_m、Pr_t∞--模型常数；

τ_u＝k/ε-动力时间尺度；

R＝τ_θ/τ_u-混合时间尺度，τ_θ为k_θ/ε_θ--热湍流时间尺度；

R_t为湍流雷诺数，k²/(vε)；R_ε为特征雷诺数，δu_ε/v，δ为离开壁面的距离；

Pr--流体的分子普朗特数；

步骤9：基于当前迭代步更新的流体温度场信息及流固耦合面温度信息，迭代求解固体能量守恒方程(4)至计算残差，

固体能量守恒方程(4)为：

步骤10：基于当前迭代步更新的固体温度场更新用户算例文件夹下的固体物性字典文件中的固体物性值；

步骤11：判断经步骤10更新后的流体质量守恒方程(1)、流体动量守恒方程(2)和流体能量守恒方程(3)、OpenFOAM内置湍流模型湍动能k及其耗散率ε或其比耗散率ω输运方程、液态金属两方程热湍流模型温度脉动k_θ及其耗散率ε_θ输运方程、固体能量守恒方程(4)的残差是否达到设定的外残差阈值，若达到，则判定外迭代收敛；若未达到，则按照步骤5设定的各物理量的边界条件更新壁面边界网格上的u_i、P、T、v_t、k、ε或ω、α_t、k_θ、ε_θ数值，并重复步骤6-11，直至达到外迭代设定的残差阈值。

最终流体域及固体域各物理量达到设定的外迭代收敛条件后，能够获得基于先进的液态金属两方程k_θ-ε_θ热湍流模型的液态金属冷却剂与固体材料之间的三维高保真低普朗特数流动传热特性及固体域热扩散行为，并可基于此数值计算方法获得液态金属冷却反应堆燃料组件、换热器等换热结构内部液态金属冷却剂低普朗特数传热情况及固体域内部热分布情况。

本发明基于开源计算流体力学程序OpenFOAM以及液态金属两方程k_θ-ε_θ热湍流模型，开发了适用于液态金属冷却剂与固体材料之间共轭传热问题的计算方法，实现了针对液态金属冷却反应堆内液态金属与固体材料之间的三维高保真共轭传热数值模拟计算，为液态金属冷却反应堆热工水力设计和安全分析提供精确的数值计算方法。该方法克服了雷诺比拟假设的局限性，避开使用经验性较强的湍流普朗特数关系式，将理论先进的液态金属两方程热湍流模型应用到液态金属流-固共轭传热过程中，为研究反应堆内处于大温差和变热流环境下的液态金属冷却剂与固体材料之间热工水力特性提供高保真数值计算方法，属于先进液态金属热湍流模型在开源CFD平台中的二次开发和应用。

与现有技术相比，本发明具有以下优点：

1.本发明在液态金属共轭传热计算中，引入了理论性较强、理论精度较高、理论边界较简单的液态金属两方程k_θ-ε_θ热湍流模型，基于该模型，可以在计算液态金属共轭传热问题的同时，获得更为丰富的液态金属低普朗特数湍流换热行为和固体温度行为。

2.本发明方法基于开源计算流体力学程序OpenFOAM保留了调用OpenFOAM内置湍流模型的功能函数，可以根据实际计算流动几何的复杂程度和计算资源等要求而调用不同类型的湍流模型来获得液态金属的湍流流动行为，可以便捷地综合利用不同的湍流模型并结合液态金属两方程热湍流输运模型来改进液态金属的共轭传热计算过程。

3.本发明方法是将先进的液态金属热湍流模型嵌入开源CFD软件中的共轭传热求解器的二次开发和应用，可为开发适用于液态金属冷却反应堆流-固共轭传热特性研究的高保真数值方法提供思路。

附图说明

图1为本发明方法共轭传热输运机制的示意图；

图2为传统共轭传热输运机制的示意图；

图3为本发明方法的计算流程图；

图4为本发明具体实施方式中的液态金属圆管流动示意图；

图5为本发明具体实施方式中的液态金属圆管网格示意图；

图6为本发明计算方法、传统计算方法与经验关系式的传热计算结果对比图；

图7为本发明具体实施方式中的液态金属圆管共轭传热结果示意图。

具体实施方式

以下结合附图和具体实施方式对本发明进行进一步详细说明。

本发明方法共轭传热输运机制如图1所示，传统共轭传热输运机制如图2所示，本发明提供的液态金属共轭传热输运机制，引入了液态金属两方程热湍流输运模型，通过输运计算温度脉动k_θ及其耗散率ε_θ，获得热湍流时间尺度，然后结合OpenFOAM内置湍流模型输运获得的动力湍流时间尺度，综合计算得到湍流热扩散系数α_t，从而改进输运液态金属的温度场。上述过程代替了液态金属传统共轭传热方法中所使用的理论性较差的湍流普朗特数模型。为了方便比较本发明液态金属共轭传热输运机制与传统液态金属共轭传热输运机制的计算效果，接下来在实施例中进行两种输运机制计算结果与液态金属圆管经验关系式进行对比。

实施例1、

液态金属圆管流动几何结构如图4所示，其外直径为19.05mm、内直径为16.56mm、流动长度为1219.2mm，流体为液态铅铋，固体为低碳钢。图5是液态金属圆管网格示意图。

一种基于两方程热湍流模型的液态金属共轭传热计算方法，计算流程如图3所示，具体步骤如下：

步骤1：在OpenFOAM中，基于内置流固共轭传热求解器，建立用户自定义共轭传热求解器。

步骤2：在步骤1建立的用户自定义共轭传热求解器中，建立湍流热扩散系数变量α_t和输运计算所述湍流热扩散系数的温度脉动变量k_θ及其耗散率变量ε_θ；调用热物性函数库中的温度-压力及物性变量，包括温度T、压力P、流体密度ρ_f、流体比热容C_pf、流体导热系数λ_f、流体分子运动黏度v_f、固体密度ρ_s、固体比热容C_ps和固体导热系数λ_s；调用湍流函数库中的速度及湍流变量，包括速度u_i、湍流运动黏度v_t和输运计算所述湍流运动黏度的湍动能k及其耗散率ε。

步骤3：在步骤1建立的用户自定义共轭传热求解器的流体文件夹下添加液态金属两方程k_θ-ε_θ热湍流模型程序文件，然后将OpenFOAM内置与壁面距离相关的函数头文件添加至步骤1建立的用户自定义共轭传热求解器的主程序文件中，并将液态金属两方程k_θ-ε_θ热湍流模型程序文件名添加至控制求解流体域的程序文件中，进一步检查、保存并编译当前步骤修改的用户自定义共轭传热求解器。

步骤1-3完成本计算方法中基于两方程热湍流模型的液态金属共轭传热求解器的用户自定义开发。

步骤4：根据实际计算问题，在用户算例文件夹下的流体物性字典文件中对步骤2中调用的流体及固体热物性变量施加合适的OpenFOAM标准物性条件。

步骤5：在用户算例初始文件夹下对步骤2中调用和建立的各物理变量施加合适的数值模拟边界条件；本实施例中，液态铅铋进口为流量入口条件，出口为大气压常压，在外壁面上施加恒定热流密度，流体各物理量在外壁面和内壁面上的边界条件类型如表1所示。

表1

参数	外壁面	内壁面
			k	-	基于湍动能的低雷诺数壁面函数
ω	-	标准的湍动能比耗散壁面函数
			v<sub>*</sub>	-	标准的低雷诺数湍流运动黏度壁面函数
k<sub>θ</sub>	-	零梯度条件
			ε<sub>θ</sub>	-	零梯度条件
α<sub>t</sub>	-	自动计算条件
			u<sub>i</sub>	-	无滑移条件
T	-	可压缩：湍流-温度流固耦合条件
			P	-	零梯度条件

在进行步骤6计算开始之前，完成如图3所示的计算流程中前处理的其余部分设置内容，包括：1)采用网格划分软件GAMBIT划分高质量mesh网格，采用网格转化工具将mesh网格生成OpenFOAM可识别网格文件polyMesh，生成的网格示意图，如图5所示；2)在用户算例文件夹下的流体湍流模型字典文件中指定计算所调用的湍流模型，本实施例使用SST k-ω湍流模型；3)在用户算例文件夹下的流体重力字典文件中指定重力大小及方向；4)在用户算例文件夹下的区域参数字典文件中指定流固类型；5)在用户算例文件夹下的数值离散格式字典文件中指定流体域和固体域时间项为稳态格式，梯度项为高斯线性格式，对流项为有界高斯迎风格式，拉普拉斯项为高斯线性修正格式；6)在用户算例文件夹下的流体数值求解格式字典文件中指定流体域的压力采用代数多重网格法(GAMG)求解，其余各物理量采用稳定预条件共轭法(PBiCG)求解；7)在用户算例文件夹下的固体数值求解格式字典文件中指定固体比焓采用对称共轭法(PCG)求解。

完成步骤4-5之后，采用步骤3编译完毕的用户自定义共轭传热求解器，进行液态金属圆管共轭传热计算，即实施如图3所示计算流程中的共轭传热求解计算过程，即步骤6-11。

步骤6：采用OpenFOAM内置PIMPLE算法求解液态金属流体质量守恒方程(1)、流体动量守恒方程(2)和流体能量守恒方程(3)，

流体质量守恒方程(1)为：

流体动量守恒方程(2)为：

流体能量守恒方程(3)为：

其中：

h-液态金属的比焓；

K-液态金属的比动能；

g_i-重力加速度；

具体步骤如下：

步骤6-1：求解流体动量守恒方程(2)，并更新速度场；

步骤6-2：利用步骤6-1更新的速度场，迭代求解流体能量守恒方程(3)，并更新温度场和流体热物性；

步骤6-3：基于步骤6-1和步骤6-2更新的数据，对PIMPLE算法下的流体质量守恒方程(1)与流体动量守恒方程(2)联立导出的压力泊松方程进行迭代求解，判断是否达到设定的PIMPLE迭代次数，如达到，则判定当前迭代步PIMPLE压力-速度耦合迭代求解结束，更新速度场和压力场后进行下一步骤。

步骤7：调用用户算例文件夹下的湍流模型字典文件中指定的SSTk-ω湍流模型，然后迭代求解湍动能k及其比耗散率ω输运方程至设定的计算残差，并更新湍流运动黏度v_t。

步骤8：求解步骤3在用户自定义共轭传热求解器中建立的液态金属两方程k_θ-ε_θ热湍流输运模型至设定的残差，并更新湍流热扩散系数α_t，具体步骤如下：

步骤8-1：求解步骤3在用户自定义共轭传热求解器中建立的温度脉动k_θ及其耗散率ε_θ输运方程至设定的计算残差：

其中：

σ_kθ、σ_εθ、C_p1、C_p2、C_d1、C_d2-模型常数；

f_d2＝1/C_d2(C_ε2f_ε-1)[1-exp(-R_ε/5.7)]²-温度脉动耗散阻尼函数；

f_ε＝1-0.3exp[-(R_t/6.5)²]-湍动能耗散阻尼函数；

-温度脉动产生项；

-湍动能产生项；

步骤8-2：由步骤8-1更新后的温度脉动k_θ及其耗散率ε_θ，更新湍流热扩散系数α_t：

α_t＝C_λf_λkτ_u

其中：

C_λ、C_m、Pr_t∞--模型常数；

τ_u＝k/ε-动力时间尺度；

R＝τ_θ/τ_u-混合时间尺度，τ_θ为k_θ/ε_θ--热湍流时间尺度；

R_t为湍流雷诺数，k²/(vε)；R_ε为特征雷诺数，δu_ε/v，δ为离开壁面的距离；

Pr-流体的分子普朗特数；

步骤9：基于当前迭代步更新的流体温度场信息及流固耦合面温度信息，迭代求解固体能量守恒方程(4)至计算残差，

固体能量守恒方程(4)为：

步骤10：基于当前迭代步更新的固体温度场更新用户算例文件夹下的固体物性字典文件中的固体物性值。

步骤11：判断经步骤10更新后的流体质量守恒方程(1)、流体动量守恒方程(2)和流体能量守恒方程(3)、OpenFOAM内置湍流模型湍动能k及其比耗散率ω输运方程、液态金属两方程热湍流模型温度脉动k_θ及其耗散率ε_θ输运方程、固体能量守恒方程(4)的残差是否达到设定的外残差阈值，若达到，则判定外迭代收敛；若未达到，则按照步骤5设定的各物理量的边界条件更新壁面边界网格上的u_i、P、T、v_t、k、ω、α_t、k_θ、ε_θ数值，并重复步骤6-11，直至达到外迭代设定的残差阈值。

当计算达到外迭代收敛后，进入到如图3所示的计算流程中的数据处理部分。本实施例中，采用开源后处理软件Parayiew进行数据处理。

为了比对本发明液态金属共轭传热输运机制与传统液态金属共轭传热输运机制的计算精度，选取本实施例中达到热充分发展段的距离进口1m位置处的横截面做传热数据处理，计算结果与经验关系式对比如6所示。图6中，分别将本发明计算方法、传统计算方法(Prt＝0.8和Prt＝Kays)与Cheng等和Kays等推荐的液态金属圆管传热经验关系式进行了对比，可以看出，本计算方法相较于传统计算方法(Prt＝0.8和Prt＝Kays)计算的传热努塞尔数更接近于Cheng等和Kays等推荐的液态金属圆管传热经验关系式，有效地解决了液态金属传统计算方法精度不佳的问题。

最终，基于本发明方法，可在液态金属共轭传热计算问题中，额外获得如图7所示的丰富的液态金属低普朗特数湍流换热特性(包括温度脉动分布、湍流热扩散系数分布及湍流普朗特数分布)以及固体域温度分布。其中图7为达到热充分发展段的距离进口1m位置处的横截面上的中心线的速度、温度、湍动能、温度脉动、湍流热扩散系数、湍流普朗特数分布示意图。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明未详细说明部分属本领域技术人员公知常识。

Claims (1)

1.一种基于两方程热湍流模型的液态金属共轭传热计算方法，其特征在于包括以下步骤：

步骤1：在OpenFOAM中，基于内置流固共轭传热求解器，建立用户自定义共轭传热求解器；

步骤2：在步骤1建立的用户自定义共轭传热求解器中，建立湍流热扩散系数变量α_t和输运计算所述湍流热扩散系数的温度脉动变量k_θ及其耗散率变量ε_θ；调用热物性函数库中的温度-压力及物性变量，包括温度T、压力P、流体密度ρ_f、流体比热容C_pf、流体导热系数λ_f、流体分子运动黏度v_f、固体密度ρ_s、固体比热容C_ps和固体导热系数λ_s；调用湍流函数库中的速度及湍流变量，包括速度u_i、湍流运动黏度v_t和输运计算湍流运动黏度的湍动能k及其耗散率ε；

步骤3：添加液态金属两方程k_θ-ε_θ热湍流模型程序文件及相关函数头文件；

步骤4：根据实际计算问题，在用户算例文件夹下的流体物性字典文件中对步骤2中调用的流体及固体热物性变量施加合适的OpenFOAM标准物性条件；

步骤5：在用户算例初始文件夹下对步骤2中调用和建立的各物理变量施加合适的数值模拟边界条件，具体步骤如下：

步骤5-1：对以下物理变量：温度T、压力P、速度u_i、湍流运动黏度v_t、湍动能k及其耗散率ε或其比耗散率ω、湍流热扩散系数α_t调用标准的OpenFOAM壁面边界条件；

步骤5-2：对处于流-固耦合面上的温度脉动k_θ及其耗散率ε_θ施加OpenFOAM标准的零固定梯度值边界条件；

步骤6：采用OpenFOAM内置PIMPLE算法求解液态金属流体质量守恒方程(1)、流体动量守恒方程(2)和流体能量守恒方程(3)，

流体质量守恒方程(1)为：

流体动量守恒方程(2)为：

流体能量守恒方程(3)为：

其中：

h-液态金属的比焓；

K-液态金属的比动能；

g_i-重力加速度；

具体步骤如下：

步骤6-1：求解流体动量守恒方程(2)，并更新速度场；

步骤6-2：利用步骤6-1更新的速度场，迭代求解流体能量守恒方程(3)，并更新温度场和流体热物性；

步骤6-3：基于步骤6-1和步骤6-2更新的数据，对PIMPLE算法下的流体质量守恒方程(1)与流体动量守恒方程(2)联立导出的压力泊松方程进行迭代求解，判断是否达到设定的PIMPLE迭代次数，如达到，则判定当前迭代步PIMPLE压力-速度耦合迭代求解结束，更新速度场和压力场后进行下一步骤；

步骤7：在用户算例文件夹下的湍流模型字典文件turbulenceProperties中，调用OpenFOAM内置湍流模型，包括但不限于：标准k-ε湍流模型、SST k-ω湍流模型或RNGk-ε湍流模型等，然后迭代求解湍动能k及其耗散率ε或其比耗散率ω输运方程至设定的计算残差，并更新湍流运动黏度v_t；

步骤8：求解步骤3在用户自定义共轭传热求解器中建立的液态金属两方程k_θ-ε_θ热湍流输运模型至设定的残差，并更新湍流热扩散系数α_t，具体步骤如下：

步骤8-1：求解步骤3在用户自定义共轭传热求解器中建立的温度脉动k_θ及其耗散率ε_θ输运方程至设定的计算残差：

其中：

σ_kθ、σ_εθ、C_p1、C_p2、C_d1、C_d2-模型常数；

f_d2＝1/C_d2(C_ε2f_ε-1)[1-exp(-R_ε/5.7)]²-温度脉动耗散阻尼函数；

f_ε＝1-0.3exp[-(R_t/6.5)²]-湍动能耗散阻尼函数；

-温度脉动产生项；

-湍动能产生项；

步骤8-2：由步骤8-1更新后的温度脉动k_θ及其耗散率ε_θ，更新湍流热扩散系数α_t：

α_t＝C_λf_λkτ_u

其中：

C_λ、C_m、

--模型常数；

τ_u＝k/ε-动力时间尺度；

R＝τ_θ/τ_u-混合时间尺度，τ_θ为k_θ/ε_θ--热湍流时间尺度；

R_t为湍流雷诺数，k²/(vε)；R_ε为特征雷诺数，δu_ε/v，δ为离开壁面的距离；

Pr-流体的分子普朗特数；

步骤9：基于当前迭代步更新的流体温度场信息及流固耦合面温度信息，迭代求解固体能量守恒方程(4)至计算残差，

固体能量守恒方程(4)为：

/>

步骤10：基于当前迭代步更新的固体温度场更新用户算例文件夹下的固体物性字典文件中的固体物性值；

步骤11：判断经步骤10更新后的流体质量守恒方程(1)、流体动量守恒方程(2)和流体能量守恒方程(3)、OpenFOAM内置湍流模型湍动能k及其耗散率ε或其比耗散率ω输运方程、液态金属两方程热湍流模型温度脉动k_θ及其耗散率ε_θ输运方程、固体能量守恒方程(4)的残差是否达到设定的外残差阈值，若达到，则判定外迭代收敛；若未达到，则按照步骤5设定的各物理量的边界条件更新壁面边界网格上的u_i、P、T、v_t、k、ε或ω、α_t、k_θ、ε_θ数值，并重复步骤6-11，直至达到外迭代设定的残差阈值。

Images