CN116562153A - 一种液态金属热分层特性的计算方法 - Google Patents

Abstract

本发明公开了一种液态金属热分层特性的计算方法，属于液态金属冷却核反应堆热工水力领域，解决了传统雷诺比拟假设方法计算液态金属热分层特性的数值精度不佳的问题。本发明方法通过引入热膨胀系数和浮力关联项有效地将一阶两方程湍流模型和一阶两方程换热模型相结合，从而建立起适配液态金属低普朗特数浮力主导流动传热计算的高保真三维流动传热分析平台。该方法使用先进的一阶两方程换热模型修正液态金属的湍流热扩散计算，同时将液态金属湍流热扩散的输运特性通过浮升力反馈至一阶两方程湍流模型，从而获得液态金属在某些浮力主导流动状态下的热分层传热特征，尤其是在自然对流和混合对流工况下。

Description

一种液态金属热分层特性的计算方法

技术领域

本发明属于液态金属冷却核反应堆热工水力领域，具体涉及一种液态金属热分层特性的计算方法。

背景技术

液态金属(如铅、铅铋合金、钠、钠钾合金等冷却剂)中子吸收截面小、导热系数高、自然循环潜能大，是第四代反应堆中的铅冷快堆和钠冷快堆以及先进加速器驱动次临界反应堆的主选冷却剂，可进一步提升核反应堆的可持续性、安全性及经济性。然而液态金属冷却核反应堆往往伴随着由高温度梯度和大密度差引起的热分层现象，并伴随热振荡的产生，这使得全工况下反应堆堆腔或管道内面临着由温度波动引起的设备热疲劳损坏、冷却剂流量振荡、阻碍事故工况下堆芯自然循环重新建立等潜在威胁。因此，准确掌握液态金属在反应堆内的热分层特性对科学评估液态金属冷却核反应堆的安全性具有重要意义。

然而如图1和图2所示，液态金属具有强烈的低普朗特数(Pr)流动传热特征，其热边界层厚度(δ_h)远大于其速度边界层厚度(δ_v)。液态金属具有的强分子导热作用在多数流动工况下可能影响至湍流核心区域，导致适用于常规流体热分层特性研究的雷诺比拟假设在计算液态金属热分层特性时适用性较差。当前，国内外针对液态金属的热分层特性和对流传热机理的数值研究多处于基于一阶两方程k-ε湍流模型与雷诺比拟假设(湍流普朗特数Pr_t为常数的模型)相结合的模拟水平，难以更进一步获得热分层现象中液态金属的温度脉动及温度振荡特性。因此，面向我国先进液态金属反应堆研发与设计的迫切需求，亟待系统地建立适配液态金属低普朗特数热分层特性和传热机制研究的有效湍流换热计算方法，并准确厘清液态金属的热分层特性和传热机理。

如图3所示，为克服传统雷诺比拟假设计算液态金属热分层特性的局限性，国内外学者对湍流热通量项直接建立微分输运方程或代数输运方程，形成二阶微分热通量模型(DHFM)或一阶半代数热通量模型(AHFM)，研究了低普朗特数流体在简单几何下的自然对流和混合对流特性。尽管更高精度的DHFM模型和AHFM模型可以克服的计算假设更多，但因为输运方程形式复杂和模化系数敏感性强导致数值求解难度大，尚未在复杂几何结构内的热分层研究中得到广泛地应用。为此，许多学者类比一阶两方程k-ε湍流模型,对湍流热通量建立一阶两方程k_θ-ε_θ换热模型，形成一阶四方程k-ε-k_θ-ε_θ湍流换热模型。如图4所示，一阶四方程k-ε-k_θ-ε_θ模型具有理论精度高和数值效率佳的求解特性，并可克服雷诺比拟假设的局限性。

四方程k-ε-k_θ-ε_θ模型通过动力时间尺度和热时间尺度联合输运湍流运动粘度和湍流热扩散系数，从而获得随雷诺数、Pr和几何结构变化的Pr_t分布，进而优化液态金属的动量和能量的输运过程。近年来的研究发现相比于雷诺比拟假设，基于四方程模型可以有效预测液态金属在反应堆复杂几何结构下的近壁湍流换热行为和具有强分子导热特性的热通量输运过程。然而，当前的四方程模型多应用于液体金属强迫对流换热研究，引入浮力效应的四方程模型用于研究液态金属热分层特性及浮力对流换热过程的报道较少。面向制约先进液态金属冷却核反应堆数值传热技术研发和堆芯安全分析的科学问题之一，亟待深入研究适配液态金属低普朗特数热分层特性研究的浮力四方程模型计算方法。

为此，本发明基于开源计算流体力学程序OpenFOAM，结合一阶两方程k-ε湍流模型和一阶两方程k_θ-ε_θ换热模型，引入浮力影响机制，开发适用于液态金属低普朗特数浮力主导热分层特性研究的高保真数值模拟计算方法，可为建立液态金属低普朗特数流体湍流换热模型、研究液态金属热分层特性和浮力主导对流传热机理、丰富液态金属反应堆数值传热技术和堆芯安全分析提供恰当的计算方法。

发明内容

本发明的目的是提供一种液态金属热分层特性的计算方法，以解决传统雷诺比拟假设方法计算液态金属热分层特性的数值精度不佳的问题。

本发明的技术方案是：一种液态金属热分层特性的计算方法，包括以下步骤：

步骤1：基于OpenFOAM内嵌的适用于不可压缩和可压缩流体浮力湍流的单相流动传热求解器，建立用户自定义单相流动传热求解器；

步骤2：在用户自定义单相流动传热求解器中的场创建文件中，定义一阶两方程湍流模型、一阶两方程换热模型及浮力项相关变量，包括：速度u_i、湍流运动黏度ν_t、湍动能k及其耗散率ε、湍流热扩散系数α_t、温度脉动k_θ及其耗散率ε_θ、流体热膨胀系数β；调用热物性函数库中的相关变量，包括：温度T、压力P、流体密度ρ、流体比热容C_p、流体导热系数λ、流体分子运动黏度ν；

步骤3：在求解动量方程文件中，不采用内嵌湍流模型调用函数，引入速度u_i、分子运动黏度ν和湍流运动黏度ν_t、流体密度ρ相关联的拉普拉斯项，为调用自定义一阶两方程湍流模型提供输入接口，求解的动量方程如下：

其中：

t为流动时间，

u_i和u_j为速度分量，

x_i和x_j为三维坐标分量，

ρ为流体密度，

ν为流体分子运动黏度，

ν_t为流体湍流运动黏度，

P为流体压力，

g_i为重力加速度分量；

步骤4：在求解能量方程文件中，不采用内嵌有效热扩散系数调用函数，引入比焓h、分子热扩散系数α和湍流热扩散系数α_t、流体密度ρ相关联的拉普拉斯项，为调用自定义一阶两方程换热模型提供输入接口，求解的能量方程如下：

其中：

α为流体分子热扩散系数，

α_t为流体湍流热扩散系数，

h为比焓，

K为比动能；

步骤5：保留使用基于SIMPLE算法的压力泊松方程求解文件，求解流体的连续性方程如下：

步骤6：不调用OpenFOAM内嵌湍流模型和壁面函数，定义求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件，然后添加更新湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界值的程序段，该程序段实现如下关系：

对于无滑移壁面，湍动能k及其耗散率ε在壁面上满足：

对于恒温壁面，温度脉动k_θ及其耗散率ε_θ在壁面上满足：

其中：

下标w代表该变量在壁面上的值，

δ为离开壁面的距离；

具体步骤如下：

步骤6-1：分别建立等式(4)和等式(5)等号右边的体积标量场；

步骤6-2：根据输入计算网格，定义存放壁面网格中心坐标值的变量；

步骤6-3：分别获取步骤6-1建立的体积标量场在壁面上的值，然后获取湍动能耗散率ε和温度脉动耗散率ε_θ在壁面上的值；

步骤6-4：遍历步骤6-2获取的壁面网格中心坐标变量，依次将步骤6-3获取的体积标量场在壁面上的值赋予给步骤6-3获取的湍动能耗散率ε和温度脉动耗散率ε_θ在壁面上的值；

步骤7：建立求解一阶两方程湍流模型的程序文件，包括如下考虑浮升力效应的一阶两方程湍流模型：

ν_t＝C_uf_ukτ_u

其中：

C_u、σ_k、σ_ε、C_ε1、C_ε2为模型系数；

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

R_t＝k²/(νε)为湍流雷诺数；

R_ε＝δu_ε/v为特征雷诺数，u_ε＝(vε)^0.25为Kolmogorov速度尺度；

为湍动能产生项，/>为雷诺应力项；

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

为浮升力产生项，β为热膨胀系数；

为沿重力方向上的速度分量值，/>为重力矢量，/>为速度矢量；

为沿非重力方向上的速度分量值；

步骤8：建立求解一阶两方程换热模型的程序文件，包括如下适用于液态金属低普朗特数流动传热的一阶两方程换热模型：

α_t＝C_λf_λkτ_u

其中：

C_λ、C_m、Pr_t∞、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)²]为修正的湍动能耗散阻尼函数；

为温度脉动产生项，/>为雷诺热通量项；

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

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

步骤9：添加用于计算一阶两方程湍流模型及一阶两方程换热模型的相关OpenFOAM壁面处理函数头文件至步骤1所建立的用户自定义单相流动传热求解器的主程序文件中；

步骤10：建立存放用户算例的计算文件夹，包括初始文件夹、固定文件夹、系统文件夹；

步骤11：使用网格划分软件对计算模型划分网格后，通过OpenFOAM内置网格转换工具转化成OpenFOAM可识别网格文件至用户算例固定文件夹下，并记录无滑移壁面、恒温壁面、恒热流加热壁面或绝热壁面的壁面边界信息；

步骤12：根据用户实际计算问题，在用户算例初始文件夹下和步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件内施加恰当的数值边界条件，具体步骤如下：

步骤12-1：在用户算例初始文件夹下对速度、温度、压力、湍流运动黏度、湍流热扩散系数施加标准的OpenFOAM边界条件；

步骤12-2：在用户算例初始文件夹下对湍动能k及其耗散率ε施加固定值边界条件，并调用在步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件来更新湍动能耗散率ε的在无滑移壁面上的边界值；

步骤12-3：在用户算例初始文件夹下对恒温壁面上的温度脉动k_θ及其耗散率ε_θ施加固定值边界条件，并调用在步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件来更新温度脉动耗散率ε_θ的在恒温壁面上的边界值；

步骤12-4：在用户算例初始文件夹下对恒热流加热或绝热壁面上的温度脉动k_θ及其耗散率ε_θ施加零梯度值边界条件；

步骤13：在用户算例固定文件夹下的物性字典文件中对流体密度、比热、导热系数、动力黏度等物性变量施加标准的OpenFOAM物性条件；在用户算例固定文件夹下的湍流模型字典文件中不调用OpenFOAM内嵌湍流模型；在用户算例固定文件夹下的重力字典文件中指定重力方向和大小；在用户自定义单相流动传热求解器中更新步骤2定义的流体热膨胀系数；

步骤14：在用户算例系统文件夹下的计算控制文件、离散格式文件、求解控制文件中分别对用户计算工况的计算步长及计算数据写入写出格式、数值离散格式、数值求解算法进行指定；

步骤15：用户自定义单相流动传热求解器开始计算后，读取步骤11转化的OpenFOAM可识别网格文件，读取步骤12建立的速度、温度、压力、湍流运动黏度、湍流热扩散系数、湍动能及其耗散率、温度脉动及其耗散率的边界条件，读取步骤13建立的流体物性条件、重力条件方向和大小，读取步骤14建立的计算步长及计算数据写入写出格式、数值离散格式、数值求解算法；

步骤16：求解步骤3建立的动量方程程序文件，得到预测速度场；

步骤17：求解步骤4建立的能量方程程序文件，更新流体温度场及热物性参数；

步骤18：求解步骤5建立的基于SIMPLE算法的压力泊松方程求解文件，更新速度场与压力场；

步骤19：求解步骤6建立的求解湍动能耗散率和温度脉动耗散率的壁面边界条件文件，更新湍动能耗散率和温度脉动耗散率的壁面边界值；

步骤20：求解步骤7建立的考虑浮升力效应的一阶两方程湍流模型，迭代求解湍动能及其耗散率微分输运方程，并更新液态金属湍流运动黏度；

步骤21：求解步骤8建立的适用于液态金属低普朗特数流动传热的一阶两方程换热模型，迭代求解温度脉动及其耗散率微分输运方程，并更新液态金属湍流热扩散系数；

步骤22：判断经步骤16-21迭代求解的速度、比焓、压力、湍动能及其耗散率、温度脉动及其耗散率的输运方程的计算残差是否达到设定的残差阈值，若达到，则判定迭代计算结束，若未达到，则重复步骤15-22，直至达到迭代设定的残差阈值。

最终液态金属的速度、比焓、压力、湍动能及其耗散率、温度脉动及其耗散率等物理量达到设定的迭代收敛条件后，能够通过后处理软件获得基于适配液态金属低普朗特数流动传热的一阶两方程换热模型与考虑浮升力效应修正的一阶两方程湍流模型的液态金属浮力主导热对流分层现象的高保真三维流动传热特性，并可基于本计算方法进行液态金属反应堆堆腔内热分层特性、浮力主导对流传热机理、低普朗特数流体湍流换热模型及规律等研究。

本发明的有益效果是：

1.本发明在液态金属浮力主导单相对流传热的热分层计算中，引入了相比于传统雷诺比拟假设模型理论精度更高的液态金属一阶两方程换热模型，基于该模型，可以在计算液态金属热分层特性问题的同时，获得适配液态金属低普朗特数浮力主导的三维高保真湍流换热特性和现象。

2.本发明在液态金属浮力主导单相对流传热的热分层计算中，引入了考虑浮升力效应修正的一阶两方程湍流模型，通过浮升力项的引入有效地将基于一阶两方程换热模型输运后的液态金属湍流热扩散系数传递至一阶两方程湍流模型输运中，可大大地改善液态金属浮力主导湍流换热的数值计算精度，从而可获得适配液态金属低普朗特数传热的热分层特性及现象。

3.本发明方法将适配液态金属低普朗特数传热的一阶两方程换热模型和考虑浮升力效应修正的一阶两方程湍流模型同时置于开源计算流体力学程序OpenFOAM的用户自定义求解器内进行开发和编译，可以便捷地将液态金属湍流热扩散系数传递至液态金属湍流模型输运中，无须添加或修正OpenFOAM湍流模型底层代码，使用户在求解器内可以快速地自定义浮升力、一阶两方程换热模型、一阶两方程湍流模型，可为研究液态金属低普朗特数湍流换热模型提供高保真三维单相对流传热计算分析平台。

4.本发明方法是将适配液态金属低普朗特数传热的一阶两方程换热模型和考虑浮升力效应修正的一阶两方程湍流模型共同嵌入开源计算流体力学程序OpenFOAM中的内置适用于不可压缩和可压缩流体浮力湍流的单相流动传热求解器的二次开发和应用，克服了传统闭源商业计算流体力学程序基于先进液态金属湍流换热模型和浮力传递机制开发液态金属对流传热求解器的困难，可为开发适用于液态金属低普朗特数热分层特性研究的高保真数值方法提供思路。

附图说明

图1为液态铅铋与液态钠的分子普朗特数随流体温度变化示意图；

图2为不同流体的速度边界层与温度边界层分布示意图；

图3为低普朗特数流体湍流换热模型分类示意图；

图4为不同阶次湍流换热模型求解复杂流动传热的理论精度及数值特性示意图；

图5为本发明的计算流程示意图；

图6为液态金属强迫对流换热输运机制示意图；

图7为本发明液态金属浮力主导对流换热输运机制示意图；

图8为本发明实施例1的几何模型示意图；

图9为本发明实施例1的网格模型示意图；

图10为本发明方法Pr＝0.025时无量纲温度模拟结果对比示意图；

图11为本发明方法Pr＝0.025时无量纲雷诺热通量模拟结果对比示意图；

图12为本发明方法Pr＝0.025时无量纲温度波动模拟结果对比示意图；

图13为本发明方法Pr＝0.006时无量纲温度模拟结果对比示意图；

图14为本发明方法Pr＝0.006时无量纲温度波动模拟结果对比示意图；

图15为本发明方法Pr＝0.025时剖面速度分布示意图；

图16为本发明方法Pr＝0.025时剖面温度分布示意图；

图17为本发明方法Pr＝0.025时剖面湍动能分布示意图；

图18为本发明方法Pr＝0.025时剖面温度波动分布示意图；

图19为本发明方法Pr＝0.025时剖面湍流运动黏度分布示意图；

图20为本发明方法Pr＝0.025时剖面湍流热扩散系数分布示意图。

具体实施方式

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

基于一阶两方程湍流模型和一阶两方程换热模型耦合输运的液态金属强迫对流换热输运机制如图6所示；本发明提供的一种液态金属浮力主导对流换热输运机制如图7所示，与图6形成对比，表明本发明方法在引入适配液态金属低普朗特数传热特性的一阶两方程换热模型的基础上，同时引入了考虑浮升力效应修正的一阶两方程湍流模型，通过浮升力可有效地将两方程换热模型输运后的液态金属湍流热扩散系数传递至液态金属湍动能及其耗散率输运过程中。上述过程可有效地考虑液态金属在浮力主导对流换热过程中的低普朗特数流动传热特性，进而有效地获得液态金属热分层特性及现象。为了方便理解本发明提出的一种液态金属热分层特性的计算方法的计算效果，接下来在实施例1中基于本发明方法对液态金属典型的浮力主导自然对流热分层过程进行数值模拟和计算分析。

实施例1、

结合如图8所示的液态金属在密闭方腔内的自然对流热分层计算几何模型，本实施例的具体条件如下：

几何参数：高为0.8m、宽为0.1m；

边界条件：左侧为恒温300K的壁面、右侧为恒温301K的壁面、上下侧为绝热壁面；

流动工况：重力大小为9.81m/s²且方向竖直向下；

流体物性：密度假定为1kg/m³、导热系数假定为0.01W/m/K、分子动力黏度假定为0.00001Pa·s；本实施案例分别对两种液态金属典型常物性进行计算，当Pr＝0.025时，比热为25J/kg/K；当Pr＝0.006时，比热为6J/kg/K；

计算工况：当Pr＝0.025时，取瑞利数为100000，此时对应的热膨胀系数为0.040775(1/K)；当Pr＝0.006时，取瑞利数为24000，此时对应的热膨胀系数为0.040775(1/K)；

计算时对密度引入自然对流Boussinesq假设，将密度写为1-β(T-T_ref)形式，T_ref取300K。

根据上述实施条件，按照以下步骤进行计算：

步骤1：基于OpenFOAM内嵌的适用于不可压缩和可压缩流体浮力湍流的单相流动传热求解器，建立用户自定义单相流动传热求解器。

步骤2：在用户自定义单相流动传热求解器中的场创建文件中，定义一阶两方程湍流模型、一阶两方程换热模型及浮力项相关变量，包括：速度u_i、湍流运动黏度ν_t、湍动能k及其耗散率ε、湍流热扩散系数α_t、温度脉动k_θ及其耗散率ε_θ、流体热膨胀系数β；调用热物性函数库中的相关变量，包括：温度T、压力P、流体密度ρ、流体比热容C_p、流体导热系数λ、流体分子运动黏度ν。

步骤3：在求解动量方程文件中，不采用内嵌湍流模型调用函数，引入速度u_i、分子运动黏度ν和湍流运动黏度ν_t、流体密度ρ相关联的拉普拉斯项，为调用自定义一阶两方程湍流模型提供输入接口，求解的动量方程如下：

其中：

t、u_i和u_j、x_i和x_j分别为流动时间、速度分量、三维坐标分量；

ρ、ν、ν_t分别为流体密度、流体分子运动黏度、流体湍流运动黏度；

P、g_i分别为流体压力、重力加速度分量。

步骤4：在求解能量方程文件中，不采用内嵌有效热扩散系数调用函数，引入比焓h、分子热扩散系数α和湍流热扩散系数α_t、流体密度ρ相关联的拉普拉斯项，为调用自定义一阶两方程换热模型提供输入接口，求解的能量方程如下：

其中：

α、α_t分别为流体分子热扩散系数、流体湍流热扩散系数；

h、K分别为比焓、比动能。

步骤5：保留使用基于SIMPLE算法的压力泊松方程求解文件，求解流体的连续性方程如下：

步骤6：不调用OpenFOAM内嵌湍流模型和壁面函数，定义求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件，然后添加更新湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界值的程序段，该程序段实现如下关系：

对于无滑移壁面，湍动能k及其耗散率ε在壁面上满足：

对于恒温壁面，温度脉动k_θ及其耗散率ε_θ在壁面上满足：

其中：

下标w代表该变量在壁面上的值，

δ为离开壁面的距离；

具体步骤如下：

步骤6-1：分别建立等式(4)和等式(5)等号右边的体积标量场；

步骤6-2：根据输入计算网格，定义存放壁面网格中心坐标值的变量；

步骤6-3：分别获取步骤6-1建立的体积标量场在壁面上的值，然后获取湍动能耗散率ε和温度脉动耗散率ε_θ在壁面上的值；

步骤6-4：遍历步骤6-2获取的壁面网格中心坐标变量，依次将步骤6-3获取的体积标量场在壁面上的值赋予给步骤6-3获取的湍动能耗散率ε和温度脉动耗散率ε_θ在壁面上的值。

步骤7：建立求解一阶两方程湍流模型的程序文件，包括如下考虑浮升力效应的一阶两方程湍流模型：

ν_t＝C_uf_ukτ_u

其中：

C_u、σ_k、σ_ε、C_ε1、C_ε2为模型系数；

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

R_t＝k²/(νε)为湍流雷诺数；

R_ε＝δu_ε/v为特征雷诺数，u_ε＝(vε)^0.25为Kolmogorov速度尺度；

为湍动能产生项，/>为雷诺应力项；

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

为浮升力产生项，β为热膨胀系数；

为沿重力方向上的速度分量值，/>为重力矢量，/>为速度矢量；

为沿非重力方向上的速度分量值。

步骤8：建立求解一阶两方程换热模型的程序文件，包括如下适用于液态金属低普朗特数流动传热的一阶两方程换热模型：

α_t＝C_λf_λkτ_u

其中：

C_λ、C_m、Pr_t∞、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)²]为修正的湍动能耗散阻尼函数；

为温度脉动产生项，/>为雷诺热通量项；

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

Pr—流体的分子普朗特数。

步骤9：添加用于计算一阶两方程湍流模型及一阶两方程换热模型的相关OpenFOAM壁面处理函数头文件至步骤1所建立的用户自定义单相流动传热求解器的主程序文件中。

上述步骤1-9完成了将考虑浮力效应的一阶两方程湍流模型和适用于液态金属低普朗特数流动传热的一阶两方程换热模型进行耦合输运计算的可用于液态金属热分层研究的用户自定义单相流动传热求解器的二次开发。

步骤10：根据步骤1至步骤9进行的用户自定义单相流动传热求解器的开发过程，建立存放用户算例的计算文件夹，包括初始文件夹、固定文件夹、系统文件夹，如图5所示。

步骤11：使用第三方网格划分软件GAMBIT对图8所示的计算几何模型划分网格后，通过OpenFOAM内置网格转换工具fluentMeshToFoam转化成OpenFOAM可识别网格文件polyMesh至用户算例固定文件夹下，并记录无滑移壁面、恒温壁面、绝热壁面的壁面边界信息，划分后的网格模型如图9所示。

步骤12：根据用户实际计算问题，在用户算例初始文件夹下和步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件内施加恰当的数值边界条件，如表1所示，具体步骤如下：

步骤12-1：在用户算例初始文件夹下对速度、温度、压力、湍流运动黏度、湍流热扩散系数施加标准的OpenFOAM边界条件；

步骤12-2：在用户算例初始文件夹下对湍动能k及其耗散率ε施加固定值边界条件，并调用在步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件来更新湍动能耗散率ε的在无滑移壁面上的边界值；

步骤12-3：在用户算例初始文件夹下对恒温壁面上的温度脉动k_θ及其耗散率ε_θ施加固定值边界条件，并调用在步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件来更新温度脉动耗散率ε_θ的在恒温壁面上的边界值；

步骤12-4：在用户算例初始文件夹下对绝热壁面上的温度脉动k_θ及其耗散率ε_θ施加零梯度值边界条件。

表1

参数	左恒温壁面	右恒温壁面	上、下绝热壁面
				ui	无滑移条件	无滑移条件	无滑移条件
T	固定值条件	固定值条件	零梯度条件
				P	零梯度条件	零梯度条件	零梯度条件
νt	自动计算条件	自动计算条件	自动计算条件
				k	固定值条件	固定值条件	固定值条件
ε	式(4)	式(4)	式(4)
				αt	自动计算条件	自动计算条件	自动计算条件
kθ	固定值条件	固定值条件	零梯度条件
				εθ	式(5)	式(5)	零梯度条件

步骤13：在用户算例固定文件夹下的物性字典文件中对本实施例流体密度、比热、导热系数、动力黏度等物性变量施加标准的OpenFOAM多项式函数物性条件；在用户算例固定文件夹下的湍流模型字典文件中不调用OpenFOAM内嵌湍流模型，在本实施例中指定为laminar；在用户算例固定文件夹下的重力字典文件中指定重力方向和大小；在用户自定义单相流动传热求解器中更新步骤2定义的流体热膨胀系数，本实施例为0.040775(1/K)。

步骤14：在用户算例系统文件夹下的计算控制文件、离散格式文件、求解控制文件中分别对用户计算工况的计算步长及计算数据写入写出格式、数值离散格式、数值求解算法进行指定，具体地，在本实施例中：时间项采用稳态格式、梯度项采用高斯线性格式、对流项采用有界高斯迎风格式、拉普拉斯项采用高斯线性修正格式、压力采用代数多重网格法(GAMG)求解、其余各物理量采用稳定预条件共轭法(PBiCG)求解、各物理量残差阈值设定为10^-9、压力松弛因子为0.7、速度松弛因子为0.3。

上述10-14步骤完成使用用户自定义单相流动传热求解器进行液态金属热分层计算所需要的主要前处理流程。

在完成上述步骤之后，采用经步骤13更新热膨胀系数后的用户自定义单相流动传热求解器，对如图8所示的实施例进行液态金属热分层特性计算，即实施如图5所示计算流程中的单相流动传热求解计算过程，即步骤15-22。

步骤15：用户自定义单相流动传热求解器开始计算后，读取步骤11转化的OpenFOAM可识别网格文件，读取步骤12建立的速度、温度、压力、湍流运动黏度、湍流热扩散系数、湍动能及其耗散率、温度脉动及其耗散率的边界条件，读取步骤13建立的流体物性条件、重力条件方向和大小，读取步骤14建立的计算步长及计算数据写入写出格式、数值离散格式、数值求解算法。

步骤16：求解步骤3建立的动量方程程序文件，得到预测速度场。

步骤17：求解步骤4建立的能量方程程序文件，更新流体温度场及热物性参数。

步骤18：求解步骤5建立的基于SIMPLE算法的压力泊松方程求解文件，更新速度场与压力场。

步骤19：求解步骤6建立的求解湍动能耗散率和温度脉动耗散率的壁面边界条件文件，更新湍动能耗散率和温度脉动耗散率的壁面边界值。

步骤20：求解步骤7建立的考虑浮升力效应的一阶两方程湍流模型，迭代求解湍动能及其耗散率微分输运方程，并更新液态金属湍流运动黏度。

步骤21：求解步骤8建立的适用于液态金属低普朗特数流动传热的一阶两方程换热模型，迭代求解温度脉动及其耗散率微分输运方程，并更新液态金属湍流热扩散系数。

步骤22：判断经步骤16-21迭代求解的速度、比焓、压力、湍动能及其耗散率、温度脉动及其耗散率的输运方程的计算残差是否达到设定的残差阈值，若达到，则判定迭代计算结束，若未达到，则重复步骤15-22，直至达到迭代设定的残差阈值。

当最终液态金属的速度、比焓、压力、湍动能及其耗散率、温度脉动及其耗散率等物理量达到设定的迭代收敛条件后，进入到如图5所示的计算流程中的数据处理部分。对当前实施例，采用开源后处理软件ParaView中进行数据处理。

为了比对本发明液态金属热分层输运计算方法与参考文献中传统雷诺比拟假设输运计算方法(K-Epsilon Low Reynolds)和高阶代数热通量模型计算方法(AHFM-2005)的计算精度，以高精度直接数值模拟(DNS)数据为基准，选取本实施例中半高处沿水平方向无量纲距离变化的无量纲温度、无量纲雷诺热通量、无量纲温度波动做数据分析，计算结果分别示于图10-14。从图10和图13可以看出，本发明方法相较于使用传统雷诺比拟假设模型的计算方法(K-Epsilon Low Reynolds)获得的无量纲温度更接近于高精度直接数值模拟(DNS)数据，有效地解决了传统雷诺比拟假设计算方法处理液态金属浮力主导热分层问题精度不佳的难题。综合图10-14中可以看出，本发明方法在液态金属自然对流无量纲温度、无量纲雷诺热通量和无量纲温度波动的模拟结果上较优于高阶代数热通量模型计算方法(AHFM-2005)的模拟结果，表明本发明方法可为进一步丰富液态金属热分层特性的计算方法的发展提供思路和可能。

最后，基于本发明方法，可在本实施例计算的液态金属浮力主导热分层问题中，获得如图15-20所示的适配液态金属低普朗特数流动换热的热分层特性及现象。其中，图15-20分别为Pr＝0.025时剖面速度、温度、湍动能、温度脉动、湍流运动黏度、湍流热扩散系数分布示意图。

本发明方法通过引入热膨胀系数和浮力关联项有效地将一阶两方程湍流模型和一阶两方程换热模型相结合，从而建立起适配液态金属低普朗特数浮力主导流动传热计算的高保真三维流动传热分析平台。该方法使用先进的一阶两方程换热模型修正液态金属的湍流热扩散计算，同时将液态金属湍流热扩散的输运特性通过浮升力反馈至一阶两方程湍流模型，从而获得液态金属在某些浮力主导流动状态下的热分层传热特征，尤其是在自然对流和混合对流工况下。本发明方法可为研究液态金属反应堆内热分层特性、浮力主导对流传热机理、低普朗特数流体湍流换热模型等问题提供高保真三维流动传热分析计算方法，属于先进液态金属冷却核反应堆热工水力技术和低普朗特数流体数值传热技术在开源计算流体力学程序中的二次开发和应用。

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

Claims (1)

1.一种液态金属热分层特性的计算方法，其特征在于包括以下步骤：

步骤1：基于OpenFOAM内嵌的适用于不可压缩和可压缩流体浮力湍流的单相流动传热求解器，建立用户自定义单相流动传热求解器；

步骤2：在用户自定义单相流动传热求解器中的场创建文件中，定义一阶两方程湍流模型、一阶两方程换热模型及浮力项相关变量，包括：速度u_i、湍流运动黏度ν_t、湍动能k及其耗散率ε、湍流热扩散系数α_t、温度脉动k_θ及其耗散率ε_θ、流体热膨胀系数β；调用热物性函数库中的相关变量，包括：温度T、压力P、流体密度ρ、流体比热容C_p、流体导热系数λ、流体分子运动黏度ν；

步骤3：在求解动量方程文件中，不采用内嵌湍流模型调用函数，引入速度u_i、分子运动黏度ν和湍流运动黏度ν_t、流体密度ρ相关联的拉普拉斯项，为调用自定义一阶两方程湍流模型提供输入接口，求解的动量方程如下：

其中：

t为流动时间，

u_i和u_j为速度分量，

x_i和x_j为三维坐标分量，

ρ为流体密度，

ν为流体分子运动黏度，

ν_t为流体湍流运动黏度，

P为流体压力，

g_i为重力加速度分量；

步骤4：在求解能量方程文件中，不采用内嵌有效热扩散系数调用函数，引入比焓h、分子热扩散系数α和湍流热扩散系数α_t、流体密度ρ相关联的拉普拉斯项，为调用自定义一阶两方程换热模型提供输入接口，求解的能量方程如下：

其中：

α为流体分子热扩散系数，

α_t为流体湍流热扩散系数，

h为比焓，

K为比动能；

步骤5：保留使用基于SIMPLE算法的压力泊松方程求解文件，求解流体的连续性方程如下：

步骤6：不调用OpenFOAM内嵌湍流模型和壁面函数，定义求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件，然后添加更新湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界值的程序段，该程序段实现如下关系：

对于无滑移壁面，湍动能k及其耗散率ε在壁面上满足：

对于恒温壁面，温度脉动k_θ及其耗散率ε_θ在壁面上满足：

其中：

下标w代表该变量在壁面上的值，

δ为离开壁面的距离；

具体步骤如下：

步骤6-1：分别建立等式(4)和等式(5)等号右边的体积标量场；

步骤6-2：根据输入计算网格，定义存放壁面网格中心坐标值的变量；

步骤6-3：分别获取步骤6-1建立的体积标量场在壁面上的值，然后获取湍动能耗散率ε和温度脉动耗散率ε_θ在壁面上的值；

步骤6-4：遍历步骤6-2获取的壁面网格中心坐标变量，依次将步骤6-3获取的体积标量场在壁面上的值赋予给步骤6-3获取的湍动能耗散率ε和温度脉动耗散率ε_θ在壁面上的值；

步骤7：建立求解一阶两方程湍流模型的程序文件，包括如下考虑浮升力效应的一阶两方程湍流模型：

ν_t＝C_uf_ukτ_u

其中：

C_u、σ_k、σ_ε、C_ε1、C_ε2为模型系数；

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

R_t＝k²/(νε)为湍流雷诺数；

R_ε＝δu_ε/v为特征雷诺数，u_ε＝(vε)^0.25为Kolmogorov速度尺度；

为湍动能产生项，/>为雷诺应力项；

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

为浮升力产生项，β为热膨胀系数；

为沿重力方向上的速度分量值，/>为重力矢量，/>为速度矢量；

为沿非重力方向上的速度分量值；

步骤8：建立求解一阶两方程换热模型的程序文件，包括如下适用于液态金属低普朗特数流动传热的一阶两方程换热模型：

α_t＝C_λf_λkτ_u

其中：

C_λ、C_m、Pr_t∞、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)²]为修正的湍动能耗散阻尼函数；

为温度脉动产生项，/>为雷诺热通量项；

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

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

步骤9：添加用于计算一阶两方程湍流模型及一阶两方程换热模型的相关OpenFOAM壁面处理函数头文件至步骤1所建立的用户自定义单相流动传热求解器的主程序文件中；

步骤10：建立存放用户算例的计算文件夹，包括初始文件夹、固定文件夹、系统文件夹；

步骤11：使用网格划分软件对计算模型划分网格后，通过OpenFOAM内置网格转换工具转化成OpenFOAM可识别网格文件至用户算例固定文件夹下，并记录无滑移壁面、恒温壁面、恒热流加热壁面或绝热壁面的壁面边界信息；

步骤12：根据用户实际计算问题，在用户算例初始文件夹下和步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件内施加恰当的数值边界条件，具体步骤如下：

步骤12-1：在用户算例初始文件夹下对速度、温度、压力、湍流运动黏度、湍流热扩散系数施加标准的OpenFOAM边界条件；

步骤12-2：在用户算例初始文件夹下对湍动能k及其耗散率ε施加固定值边界条件，并调用在步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件来更新湍动能耗散率ε的在无滑移壁面上的边界值；

步骤12-3：在用户算例初始文件夹下对恒温壁面上的温度脉动k_θ及其耗散率ε_θ施加固定值边界条件，并调用在步骤6定义的求解湍动能耗散率ε和温度脉动耗散率ε_θ的壁面边界条件文件来更新温度脉动耗散率ε_θ的在恒温壁面上的边界值；

步骤12-4：在用户算例初始文件夹下对恒热流加热或绝热壁面上的温度脉动k_θ及其耗散率ε_θ施加零梯度值边界条件；

步骤13：在用户算例固定文件夹下的物性字典文件中对流体密度、比热、导热系数、动力黏度等物性变量施加标准的OpenFOAM物性条件；在用户算例固定文件夹下的湍流模型字典文件中不调用OpenFOAM内嵌湍流模型；在用户算例固定文件夹下的重力字典文件中指定重力方向和大小；在用户自定义单相流动传热求解器中更新步骤2定义的流体热膨胀系数；

步骤14：在用户算例系统文件夹下的计算控制文件、离散格式文件、求解控制文件中分别对用户计算工况的计算步长及计算数据写入写出格式、数值离散格式、数值求解算法进行指定；

步骤15：用户自定义单相流动传热求解器开始计算后，读取步骤11转化的OpenFOAM可识别网格文件，读取步骤12建立的速度、温度、压力、湍流运动黏度、湍流热扩散系数、湍动能及其耗散率、温度脉动及其耗散率的边界条件，读取步骤13建立的流体物性条件、重力条件方向和大小，读取步骤14建立的计算步长及计算数据写入写出格式、数值离散格式、数值求解算法；

步骤16：求解步骤3建立的动量方程程序文件，得到预测速度场；

步骤17：求解步骤4建立的能量方程程序文件，更新流体温度场及热物性参数；

步骤18：求解步骤5建立的基于SIMPLE算法的压力泊松方程求解文件，更新速度场与压力场；

步骤19：求解步骤6建立的求解湍动能耗散率和温度脉动耗散率的壁面边界条件文件，更新湍动能耗散率和温度脉动耗散率的壁面边界值；

步骤20：求解步骤7建立的考虑浮升力效应的一阶两方程湍流模型，迭代求解湍动能及其耗散率微分输运方程，并更新液态金属湍流运动黏度；

步骤21：求解步骤8建立的适用于液态金属低普朗特数流动传热的一阶两方程换热模型，迭代求解温度脉动及其耗散率微分输运方程，并更新液态金属湍流热扩散系数；

步骤22：判断经步骤16-21迭代求解的速度、比焓、压力、湍动能及其耗散率、温度脉动及其耗散率的输运方程的计算残差是否达到设定的残差阈值，若达到，则判定迭代计算结束，若未达到，则重复步骤15-22，直至达到迭代设定的残差阈值。