[发明专利]一种液态金属流动传热计算方法

权利要求书

1.一种液态金属流动传热计算方法，其特征在于，包括如下步骤：

步骤1：在OpenFOAM中建立求解过程中需要的各物理场，包括标量场：温度T、压力p、随温度变化的热物性场、湍动能k及其单位耗散率Ω、脉动温度方均值k_θ及其单位耗散率Ω_θ；以及矢量场：液态金属的速度矢量所述随温度变化的热物性场包括液态金属密度ρ、比热容C_p、动力黏度μ和热导率λ；

步骤2：对步骤1中建立的各物理场施加数值模拟边界条件，具体步骤如下：

步骤2-1：对以下物理场：液态金属的速度矢量温度T、压力p、湍动能k、湍动能单位耗散率Ω施加OpenFOAM标准边界条件；

步骤2-2：由脉动温度方均值k_θ及其单位耗散率Ω_θ两个变量对边界条件要求严格，标准边界条件不能满足需要，故利用OpenFOAM内置的网格遍历功能对所有壁面边界网格进行识别和标记，以便在步骤7中对所标记壁面边界网格上的脉动温度方均值k_θ及其单位耗散率Ω_θ进行赋值更新；

步骤3：采用OpenFOAM内置SIMPLE算法对液态金属的质量、动量守恒方程(1)进行压力—速度耦合迭代求解；利用获得的速度、压力场求解经不可压缩、忽略粘性力做功简化后的液态金属比能方程(2)：

其中：

ρ——液态金属密度，kg·m^-3；

e——液态金属的比内能，m²·s^-2·kg^-1；

K——液态金属的比机械能，m²·s^-2·kg^-1；

——液态金属的速度矢量，m·s^-1；

α_eff——有效热扩散率，m²·s^-1；

p——压力，kg·m^-1·s^-2；

——重力加速度，m·s^-2；

μ_eff——有效动力黏度，kg·m^-1·s^-1；

步骤4：比较方程(1)求解后的残差是否小于设定的残差精度10^-4，即判断压力—速度耦合求解是否迭代收敛：若未收敛则继续对方程(1)进行压力—速度耦合迭代求解；若已收敛则结束迭代，进行下一步骤；

步骤5：基于四因子模型，即k-Ω-k_θ-Ω_θ湍流换热模型，求解动量守恒与能量守恒方程中的湍流粘度和湍流热扩散率；具体为：调用k-ΩSST湍流模型，求解并更新有效粘度ν_eff；求解k_θ-Ω_θ热湍流方程(3)，更新有效热扩散率α_eff；

其中：

k——湍动能，m²·s^-2；

Ω——湍动能单位耗散率，s^-1；

k_θ——脉动温度方均值，K²；

Ω_θ——脉动温度方均值单位耗散率，s^-1；

α——层流热扩散率，m²·s^-1；

α_t——湍流热扩散率，m²·s^-1；

ν_t——湍流运动黏度，m²·s^-1；

C_p1、C_p2、C_μ、C_d1、C_d2、σ_θ——模型常数；

其中v为层流运动黏度，m²·s^-1；

P_kθ——随温度变化的可变源项，

P_k——随速度变化的可变源项，

解出k_θ与Ω_θ值后，再由下式计算有效热扩散率：

其中：

C_θ——模型常数；

k——湍动能，m²·s^-2；

τ_lθ——局部热特征时间尺度，s，按式(5)-(7)计算：

τ_lθ＝f_1θB_1θ+f_2θB_2θ (5)

其中：

Pr——液态金属的普朗特数，Pr＝C_pμ/λ；

R_t——湍流雷诺数，R_t＝k/(C_μe^Ων)；

R_d——无量纲壁面距离，R_d＝δ(C_μe^Ωk)^0.25/ν^0.75，其中δ为壁面到第一层网格中心的距离，m；

R——热特征时间与动量时间比值，

f_1θ、f_2θ、B_1θ、B_2θ——中间变量，无物理意义；

Pr_t,∞、C_γ——模型常数；

步骤6：利用步骤3求解得到的压力、速度、温度场更新方程(1)-(3)中的热物性与可变源项，具体步骤如下：

步骤6-1：在步骤1中已将各热物性场：动力黏度μ、密度ρ、比热容C_p和热导率λ建立为与温度相关的物理场，此时根据当前迭代步的最新温度T更新质量及动量守恒方程(1)、比能方程(2)、k_θ-Ω_θ方程(3)中的热物性值；

步骤6-2：根据当前迭代步的最新速度场和温度T计算并更新k_θ-Ω_θ方程(3)中的可变源项：P_k与P_kθ；

步骤7：比较经步骤6更新后的方程(1)–(3)的残差是否小于设定的残差精度10^-4，即判断是否达到外迭代收敛；若未收敛，则按照式(8)更新步骤2中所标记的壁面边界网格上的k_θ与Ω_θ值，并重复步骤3至步骤7，直至达到外迭代收敛；

其中：

|_wall——表示该变量为壁面上的数值；

最终收敛后，能够获得基于k-Ω-k_θ-Ω_θ湍流换热模型的低普朗特数液态金属高保真三维关键热工水力参数分布，并依此评估液态金属反应堆燃料组件内部冷却剂的载热能力以及液态金属反应堆燃料棒壁面的冷却情况。