Welander自然循环问题高精度数值模拟

巢 飞，杨 文，邰 云，邱金荣，单建强

(1.武汉第二船舶设计研究所，湖北 武汉 430064；2.西安交通大学 核科学与技术学院，陕西 西安 710049)

0 引 言

自然循环是指依靠密度差和高度差形成的驱动力驱动流体循环流动的现象。自然循环在核动力系统中非常常见，其对反应堆安全有重要影响，因此准确模拟自然循环现象是当前反应堆热工水力分析中最重要研究内容之一。

目前主流热工水力系统分析程序如RELAP5，CATHARE，TRACE 等普遍采用低阶差分法来求解基本守恒方程，其在模拟自然循环不稳定性问题时存在预测误差。为提高自然循环不稳定性问题的计算精度，本文采用基于两流体双压力两相流模型的高精度数值解法对Welander 自然循环进行计算分析。

1 两流体双压力模型

两流体双压力模型基本守恒方程：

体积份额输运方程

式中：αg表示汽相体积份额；Pg和Pf分别代表汽相和液相压力，Pa；vint为相间界面速度，m·s−1；µ为压力松弛因子，Pa−1·s−1，表示汽液两相压力达到平衡的松弛速率；A为管道截面积，m2；ρint为相间界面密度，k g·m−3；Γg为单位体积液-汽相间质量交换率，kg·m−3·s−1。汽、液两相体积份额关系应满足：

其中 τP为压力松弛时间，s。

质量守恒方程

式中：Γf=−Γg为单位体积汽-液相间质量交换率，kg·m−3·s−1；ρ为相密度，kg·m−3；v为相速度，m·s−1；f为液相标识；g为汽相标识。

动量守恒方程：

式中：Pint为相间界面压力，Pa；fwall,k为阻力系数;ρc为连续相密度，kg·m−3；CiD为相间阻力系数；Aint为单位体积相间界面面积，m−1；gx为重力加速度在流动方向的投影，m·s−2。

能量守恒方程

式中：ek为k相比内能/J·kg−1；Qwg为汽相单位体积壁面换热量，W·m−3；Qwf为液相单位体积壁面换热量，W·m−3；Hig为单位体积相间界面到汽相的换热系数，W·m−3·K−1；Hif为单位体积相间界面到液相换热系数，W·m−3·K−1；Ts为相间界面温度，Tk为k相温度，K；为相间界面汽相比焓，J·kg−1；为相间界面液相比焓，J·kg−1。

2 高精度数值算法

2.1 半隐式数值解法

采用参考文献[1]作者提出的半隐数值算法求解两流体双压力模型。为达到快速计算的效果，半隐算法只对时间常数小、传播速度快的量隐式处理，如两相速度、压力梯度、相间界面质量和动量交换项，其他量显式处理。为了解决两相和单相之间转变带来的数值不稳定性问题，动量守恒方程和质量守恒方程处理成和与差的微分方程，所有守恒方程瞬态项采用展开的形式，再进行数值差分。采用交错网格划分控制体，如图1 所示。K，L等表示求解质量、能量守恒方程的控制体，存储流场中的标量如压力、比内能、空泡份额等；j代表动量控制体（虚线包围），用于求解动量守恒方程，存储流场中速度矢量。

图1 交错网格示意图Fig.1 Schematic drawing of staggered mesh

2.2 瞬态项二阶差分格式

目前现有的系统程序普遍采用一阶时间差分（BDF1）离散守恒方程瞬态项。BDF1 离散误差大，给守恒方程的差分方程引入了“人工”扩散，宜采用精度更高的二阶时间差分（BDF2）来降低瞬态项的离散误差。二阶时间差分根据Taylor 级数展开法得到[2]：

其中：∂φ/∂t表示守恒方程中的瞬态项，Δt=tn+1−tn，Δtold=tn−tn−1。该式考虑了时间步长的变化，具有一般性，在系统程序中应用更为方便。BDF1 和BDF2 可处理成统一形式：

2.3 对流项三阶TVD 差分格式

目前主流热工水力程序采用一阶施主元法离散对流项。一阶施主元法本质上是一阶迎风，计算精度低，因此宜采用稳定高阶施主元法，以提高对流项的计算精度。双压力模型守恒方程中与速度相关的项均采用施主元法计算。以能量方程对流项为例，其离散形式为：

由于速度本身定义在接管上，因此对流项的计算精度很大程度上取决于接管上的标量（如）的计算方案。接管上标量数值计算的统一形式可以写成：

其中：Δx是控制体长度；∂φ/∂x是控制体边界上的梯度；ψ(r)为通量限制函数；r为梯度比，由下式给出：

空间离散方案统一形式（13）由两部分组成，第一部分为稳定的一阶迎风方案，第二部分为反扩散项。因此，物理量 φ˙计算精度由反扩散项中的通量限制函数ψ(r)决定。

稳定三阶TVD 差分格式（命名为TOU_TVD_CFL）为[2]：

式 中C为Courant 数。

3 计算分析

3.1 Welander 自然循环问题

Welander 提出了竖直回路中单相自然循环问题来研究流体的不稳定性行为[3]。Welander 自然循环回路由2 个平行的竖直绝热管，顶部冷源和底部热源组成，如图2 所示，回路总长为L，截面积为A，热阱和热源长度为 Δs。通过顶部热阱的冷却和底部热源加热，流体在回路顶部和底部产生密度差，从而驱动流体发生自然循环流动。通过选择合适的浮力与摩擦阻力的比值，流动可能会出现稳定、弱不稳定和强不稳定行为。在文献[3]中，Welander 推导了实验回路层流下流动不稳定性理论边界，Ambrosini 和Ferreri[4]将不稳定性边界扩展湍流情况，如图3 所示。Welander 自然循环稳定性边界绘制成无量纲重力系数 α和无量纲摩擦系数 ε的关系式图，无量纲量具体表达式为：

图2 Welander 自然循环问题示意图Fig.2 Schematic diagram of Welander natural circulation

图3 Welander 自然循环问题流动不稳定性边界Fig.3 Flow instability map of Welander natural circulation

式中：Welander 定义的层流摩擦系数为R=64/Ref·Welander 定义的换热系数为温差ΔT=当a=0.316 4/4,b=0.25时，表示湍流条件；当a=16,b=1时，表示层流条件。

对于Welander 单相振荡自然循环问题，选取Ambrosinia 和Ferreri[5]设置的几何参数和边界条件。管内径为0.1 m，热源和冷源长度为0.1 m。管道充满压力为0.1 MPa 的过冷水。热源温度为30℃，与流体换热系数为20 000 W·m−2·K−1，冷源温度维持20℃，与流体换热系数为20 000 W·m−2·K−1。本文研究Welander 单相振荡自然循环问题3 个不同的实验工况：强不稳定模式、弱不稳定模式和稳定模式，见表1。图4 给出了Welander 自然循环回路数值模拟节点图。在进行模拟计算时，时间步长为0.5 s，热源和冷源分别采用一个节点模拟，两垂直绝热管都均分为N个控制体，因此环路控制体总数为2N+2。

表1 Welander 自然循环问题实验工况Tab.1 Test conditions of Welander natural circulation

图4 Welander 自然循环节点图Fig.4 Nodalization of Welander natural circulation

3.2 数值结果

对于强不稳定模式工况，环路长度取L=20.2 m，此时 (α,ε)为 (343.79,2.346 84)。强不稳定模式时，实验回路流量出现流动反转现象，并作平均速度为零的等幅振荡。为便于观测到流动不稳定性行为，实验之初，给定回路一个极小流量的扰动。

图5 给出了使用低阶差分方案BDF1+FOU 在不同控制体数下的数值结果。由图可知控制体数低于62 时，数值误差过大，无法预测出不稳定行为；控制体数达到102 时才预测出不稳定性行为，控制体数越多，数值误差越小，预测的不稳定行为振荡周期越短。由于数值误差的阻尼作用，粗糙网格可以稳定物理上不稳定的系统。采用不同计算节点，低阶差分数值结果会出现系统稳定和不稳定2 种相矛盾的情况。因此数值误差对不稳定性行为准确预测有很不利的影响，低阶差分方案BDF1+FOU 精度低，不利于研究流动不稳定行为。

图5 强不稳定模式质量流量：低阶BDF1+FOU方案的计算结果Fig.5 Mass flow of strong instability mode:calculation results of BDF1+FOU scheme

图6 给出了高阶差分方案与低阶差分方案结果的对比。高阶差分方案模拟流动不稳定性行为具有相当大的优势。高阶差分方案TOU_TVD_CFL 只需22 个控制体就可以预测出不稳定行为，因此高阶算法有效减小了数值扩散，提高了计算精度。高阶差分算法采用22 个控制体的数值精度可达到低阶差分算法控制体数为182 的计算精度，振荡周期相近。

图6 强不稳定模式质量流量：高阶和低阶方案的计算结果对比Fig.6 Mass flow of strong instability mode:comparison between low order and high order difference scheme

图7 给出了低阶差分方案BDF1+FOU 计算的流体温度与质量流量关系图，随着控制体数增多，精度越高，混沌行为（Chaotic Behavior）越明显。图8 给出了高阶BDF2+TOU_TVD_CFL 方案计算的流体温度与质量流量关系图。同样可以看出，随着控制体数增多，可以预测到流动反转的混沌演化过程，流体温度与质量流量关系图类似蝴蝶形状。

图7 低阶方案BDF1+FOU 计算的流体温度与质量流量关系图Fig.7 Fluid temperature versus mass flow calculated by low order scheme BDF1+FOU

图8 高阶BDF2+TOU_TVD_CFL 方案计算的流体温度与质量流量关系图Fig.8 Fluid temperature versus mass flow calculated by high order scheme BDF2+TOU_TVD_CFL

针对弱不稳定工况，回路总长为80.2 m。弱不稳定模式时，回路流量不会出现流动反转现象，其维持一个流动方向上的近似等幅振荡。图9 给出了采用低阶方案计算的质量流量随时间变化图。控制体数为722 时质量流量的数值结果依旧出现缓慢的衰减，低阶差分结果误差大。图10 给出了弱不稳定工况下采用高阶差分方案计算的质量流量随时间变化图。高阶算法的结果有效减小人工粘性，抑制振荡衰减，控制体数为102 时就可以模拟出等幅振荡的流动不稳定性行为。

图9 弱不稳定模式质量流量低阶方案的计算结果Fig.9 Calculation results of low order scheme by mass flow of weak instability mode

图10 弱不稳定模式质量流量高阶方案的计算结果Fig.10 Calculation results of high order scheme by mass flow of weak instability mode

针对L=5.2 m 的稳定模式工况，表现为：质量流量初始增加到最大值，然后经过一段时间的衰减型振荡后维持稳定运动。图11 和图12 给出了稳定模式工况下质量流量的计算结果，高阶方案与低阶方案计算结果对比见图13。因为人工粘性的影响，低阶差分预测更小的振荡振幅，更快达到稳定状态；控制体数为102 的低阶结果精度比控制体数为22 的高阶结果精度还低，振荡衰减更快。由于高阶差分对稳态没影响，最终高阶差分预测的稳态流量与低阶结果一致。

图11 稳定模式质量流量低阶方案的计算结果Fig.11 Calculation results of low order scheme by mass flow of stable mode

图12 稳定模式质量流量高阶方案的计算结果Fig.12 Calculation results of high order scheme by mass flow of stable mode

图13 稳定模式质量流量高阶和低阶方案的计算结果对比Fig.13 Comparison between low order and high order scheme by mass flow of stable mode

4 结 语

本文采用两流体双压力两相流模型高精度数值算法对Welander 自然循环的流动不稳定性问题进行计算分析，数值结果表明：高阶差分格式可以以较少的网格精确地预测不稳定性边界和混沌行为，而采用低阶差分在网格不够精细的情况下有可能将自然循环不稳定模式预测为稳定模式，得到错误的结果。因此高精度数值解法能有效减小数值扩散，提高自然循环问题的预测精度。