﻿ Welander自然循环问题高精度数值模拟
 舰船科学技术  2021, Vol. 43 Issue (7): 93-98    DOI: 10.3404/j.issn.1672-7649.2021.07.019 PDF
Welander自然循环问题高精度数值模拟

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

Numerical simulation on Welander natural circulation using high accuracy difference schemes
CHAO Fei1, YANG Wen1, TAI Yun1, QIU Jin-rong1, SHAN Jian-qiang2
1. Wuhan Second Ship Design and Research Institute, Wuhan 430064, China;
2. School of Nuclear Science and Technology, Xi′an Jiaotong University, Xi′an 710049, China
Abstract: The first-order difference scheme is widely employed to discretize basic conservation equations in current main thermal-hydraulics analysis code. However, first-order difference scheme has the characteristics of large numerical truncation error and low calculation accuracy when simulating natural circulation. In this paper, the high-order accuracy numerical algorithm for two-fluid two-pressure model is used to simulate Welander natural circulation problem. Numerical results imply that high-order accuracy numerical algorithm could prevent excessive numerical diffusion effectively and improve the prediction accuracy, which demonstrates the advantage of using high-order schemes.
Key words: thermal-hydraulics analysis code     two-fluid two-pressure model     high-order accuracy numerical algorithm     natural circulation
0 引　言

1 两流体双压力模型

 $A\frac{{\partial {\alpha _g}}}{{\partial t}} + A{v_{\rm{int} }}\frac{{\partial {\alpha _g}}}{{\partial x}} = A\mu \left( {{P_g} - {P_f}} \right) + A\frac{{{\varGamma _g}}}{{{\rho _{\rm{int} }}}}{\text{。}}$ (1)

 ${\alpha _g} + {\alpha _f} = 1,$ (2)

 $\mu = \frac{1}{{{\tau _P}}}\frac{{{\alpha _f}{\alpha _g}}}{{{P_f} + {P_g}}}\text{。}$ (3)

 $\frac{\partial }{{\partial t}}\left( {{\alpha _g}{\rho _g}} \right) + \frac{1}{A}\frac{\partial }{{\partial x}}\left( {{\alpha _g}{\rho _g}{v_g}A} \right) = {\varGamma _g},$ (4)
 $\frac{\partial }{{\partial t}}\left( {{\alpha _f}{\rho _f}} \right) + \frac{1}{A}\frac{\partial }{{\partial x}}\left( {{\alpha _f}{\rho _f}{v_f}A} \right) = - {\varGamma _g}\text{。}$ (5)

 $\begin{split} \frac{\partial }{{\partial t}}\left( {{\alpha _g}{\rho _g}{v_g}A} \right) +\;& \frac{\partial }{{\partial x}}\left( {{\alpha _g}{\rho _g}v_g^2A} \right) = - {\alpha _g}A\frac{{\partial {P_g}}}{{\partial x}} + \\ & ({P_{\rm{int} }}\; \; - {P_g})A\frac{{\partial {\alpha _g}}}{{\partial x}} + {\alpha _g}{\rho _g}{{\rm{g}}_{\rm{x}}}A + {\Gamma _g}{v_{\rm{int} }}A - \\ & \frac{1}{8}{\rho _c}{C_{iD}}{A_{\rm{int} }}\left( {{v_g} - {v_f}} \right)\left| {{v_g} - {v_f}} \right|A - \\ & \frac{1}{2}\frac{{{f_{wall,g}}}}{D}{\rho _g}v_g^{}\left| {v_g^{}} \right|{\alpha _g}A,\\[-15pt] \end{split}$ (6)
 $\begin{split} \frac{\partial }{{\partial t}}\left( {{\alpha _f}{\rho _f}{v_f}A} \right) +\;& \frac{\partial }{{\partial x}}\left( {{\alpha _f}{\rho _f}v_f^2A} \right) = - {\alpha _f}A\frac{{\partial {P_f}}}{{\partial x}} +\\ & ({P_{\rm{int} }} - {P_f})A\frac{{\partial {\alpha _f}}}{{\partial x}} + {\alpha _f}{\rho _f}{{\rm{g}}_{\rm{x}}}A - {\Gamma _g}{v_{\rm{int} }}A + \\ & \frac{1}{8}{\rho _c}{C_{iD}}{A_{\rm{int} }}\left( {{v_g} - {v_f}} \right)\left| {{v_g} - {v_f}} \right|A - \\ & \frac{1}{2}\frac{{{f_{wall,f}}}}{D}{\rho _f}v_f^{}\left| {v_f^{}} \right|{\alpha _f}A\text{。}\\[-15pt] \end{split}$ (7)

 $\begin{split} \frac{\partial }{{\partial t}}\left( {{\alpha _g}{\rho _g}{e_g}} \right) +& \frac{1}{A}\frac{\partial }{{\partial x}}\left( {{\alpha _g}{\rho _g}{e_g}{v_g}A} \right) +\\ &\frac{1}{A}{P_g}\frac{\partial }{{\partial x}}\left( {{\alpha _g}{v_g}A} \right) = - {P_{\rm{int} }}\frac{{\partial {\alpha _g}}}{{\partial t}} +\\ &({P_g} - {P_{\rm{int} }}){v_g}\frac{{\partial {\alpha _g}}}{{\partial x}} + {\varGamma _g}\frac{{{{({v_g} - {v_{\rm{int} }})}^2}}}{2} + {Q_{wg}} +\\ &{H_{ig}}({T^s} - {T_g}) + {\varGamma _g}h_g^* + \frac{1}{2}\frac{{{f_{wall,g}}}}{D}{\rho _g}v_g^2\left| {v_g^{}} \right|{\alpha _g} \end{split},$ (8)
 $\begin{split} \frac{\partial }{{\partial t}}\left( {{\alpha _f}{\rho _f}{e_f}} \right) +& \frac{1}{A}\frac{\partial }{{\partial x}}\left( {{\alpha _f}{\rho _f}{e_f}{v_f}A} \right) +\frac{1}{A}{P_f}\frac{\partial }{{\partial x}}\left( {{\alpha _f}{v_f}A} \right) = \\ & - {P_{\rm{int} }}\frac{{\partial {\alpha _f}}}{{\partial {\rm{t}}}} +({P_f} - {P_{\rm{int} }}){v_f}\frac{{\partial {\alpha _f}}}{{\partial x}} - \\ &{\varGamma _g}\frac{{{{({v_f} - {v_{\rm{int} }})}^2}}}{2} + {Q_{wf}} {H_{if}}({T^s} - {T_f}) - \\ &{\varGamma _g}h_f^* + \frac{1}{2}\frac{{{f_{wall,k}}}}{D}{\rho _f}v_f^2\left| {v_f^{}} \right|{\alpha _f}\text{。}\\[-15pt] \end{split}$ (9)

2 高精度数值算法 2.1 半隐式数值解法

 图 1 交错网格示意图 Fig. 1 Schematic drawing of staggered mesh
2.2 瞬态项二阶差分格式

 $\begin{split} {\left( {\frac{{\partial \varphi }}{{\partial t}}} \right)_{{t_{n{\rm{ + }}1}}}} = &\frac{{\left( {\Delta {t_{old}} + 2\Delta {t_{}}} \right)}}{{\Delta {t_{}}\left( {\Delta {t_{old}} + \Delta {t_{}}} \right)}}\varphi _{}^{n{\rm{ + }}1} -\frac{{\left( {\Delta {t_{old}} + \Delta {t_{}}} \right)}}{{\Delta {t_{old}}\Delta {t_{}}}}\varphi _{}^n +\\ & \frac{{\Delta {t_{}}}}{{\Delta {t_{old}}\left( {\Delta {t_{old}} + \Delta {t_{}}} \right)}}\varphi _{}^{n - 1}\; {\text{。}}\end{split}$ (10)

 \begin{split} &{\left( {\frac{{\partial \varphi }}{{\partial t}}} \right)_{{t_{n + 1}}}} = \frac{{\delta \varphi _{}^{n + 1}}}{{\Delta t}} = \\ &\left\{ \begin{aligned} &{ \dfrac{1}{{\Delta {t_{}}}}\left( {\varphi _{}^{n + 1} - \varphi _{}^n} \right)},\qquad\qquad {\rm{BDF1}}, \\ & \begin{aligned}\dfrac{1}{{\Delta {t_{}}}}\left( \dfrac{{\left( {\Delta {t_{old}} + 2\Delta {t_{}}} \right)}}{{\left( {\Delta {t_{old}} + \Delta {t_{}}} \right)}}\varphi _{}^{n + 1} - \dfrac{{\left( {\Delta {t_{old}} + \Delta {t_{}}} \right)}}{{\Delta {t_{old}}}}\varphi _{}^n +\right.\\ \left.\dfrac{{\left( {\Delta {t_{}}} \right)_{}^2}}{{\Delta {t_{old}}\left( {\Delta {t_{old}} + \Delta {t_{}}} \right)}}\varphi _{}^{n - 1} \right),\;\;\, {\rm{BDF2}}\text{。}\qquad\;\;\end{aligned} \end{aligned} \right. \end{split} (11)
2.3 对流项三阶TVD差分格式

 $\begin{split}&{\left. {\frac{\partial }{{\partial x}}\left( {{\alpha _k}{\rho _k}{e_k}{v_k}A} \right)} \right|_L} =\\ &\quad\frac{{\left( {\dot \alpha _{k,j + 1}^n\dot \rho _{k,j + 1}^n\dot e_{k,j + 1}^n} \right)v_{k,j + 1}^{n + 1}{A_{j + 1}} - \left( {\dot \alpha _{k,j}^n\dot \rho _{k,j}^n\dot e_{k,j}^n} \right)v_{k,j}^{n + 1}{A_j}}}{{\Delta {x_L}}}\end{split},$ (12)

 ${\dot{\phi }}^{n}{}_{j}=\left\{\begin{array}{l}{\phi }_{K}{}^{n}+\dfrac{\Delta {x}_{K}}{2}\psi \left(r\right){\left(\dfrac{\partial \phi }{\partial x}\right)}_{j-1}^{n},{\text{若}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}{v}_{j}>0,\\ {\phi }_{L}{}^{n}-\dfrac{\Delta {x}_{L}}{2}\psi \left(r\right){\left(\dfrac{\partial \phi }{\partial x}\right)}_{j+1}^{n}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}},{\text{若}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}{v}_{j}<0{\text{。}}\end{array}\right.$ (13)

 $r=\left\{\begin{array}{l}{\left(\dfrac{\partial \phi }{\partial x}\right)}_{j}^{n}\Biggr/{\left(\dfrac{\partial \phi }{\partial x}\right)}_{j-1}^{n},{\text{若}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}{v}_{j}>0,\\ {\left(\dfrac{\partial \phi }{\partial x}\right)}_{j}^{n}\Biggr/{\left(\dfrac{\partial \phi }{\partial x}\right)}_{j+1}^{n},{\text{若}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}\rm{\hspace{0.05em}}{v}_{j}<0{\text{。}}\end{array}\right.$ (14)

 $\psi \left( {{r_{}}} \right) = \max \left\{ {0,\min \left[ \begin{gathered} \min \left( {2\frac{{1 - C}}{C}\; \; ,2} \right)\; {r_{}}, \\ \min \left( {2\frac{{1 - C}}{C}\; \; ,2} \right),\frac{{2r + 1}}{3} \end{gathered} \right]} \right\}\text{。}$ (15)

3 计算分析 3.1 Welander自然循环问题

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

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

 图 3 Welander自然循环问题流动不稳定性边界 Fig. 3 Flow instability map of Welander natural circulation
 $\left\{ \begin{gathered} \alpha = \frac{{g\; \beta \; \Delta T\; \left( {L - 2\Delta s} \right)}}{{2(k\; \; \Delta s)_{}^2}}, \\ \varepsilon = \frac{a}{{16}}\left( {\frac{{{\rho _f}\; D\; k\; \Delta s}}{{{\mu _f}}}} \right)_{}^{1 - b}\frac{{R\; L}}{{2k\; \; \Delta s}} \text{。}\\ \end{gathered} \right.$ (16)

 图 4 Welander自然循环节点图 Fig. 4 Nodalization of Welander natural circulation
3.2 数值结果

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

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

 图 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

 图 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

 图 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 结　语

 [1] CHAO F., SHAN J., GOU J., et al.. Study on the algorithm for solving two-fluid seven-equation two-pressure model[J]. Annals of Nuclear Energy, 2018, 111: 379-390. DOI:10.1016/j.anucene.2017.09.021 [2] CHAO F., LIU D., SHAN J., et al.. Development of temporal and spatial high‐order schemes for two‐fluid seven‐equation two‐pressure model and its applications in two‐phase flow benchmark problems[J]. International Journal for Numerical Methods in Fluids, 2018, 88(4): 169-192. DOI:10.1002/fld.4515 [3] WELANDER P. On the oscillatory instability of a differentially heated fluid loop[J]. Journal of Fluid Mechanics, 1967, 29(1): 17-30. DOI:10.1017/S0022112067000606 [4] AMBROSINI W, FERRERI JC. The effect of truncation error on the numerical prediction of linear stability boundaries in a natural circulation single-phase loop[J]. Nuclear Engineering and Design, 1998, 183(1-2): 53-76. DOI:10.1016/S0029-5493(98)00157-5 [5] AMBROSINI W, FERRERI JC. Prediction of stability of one-dimensional natural circulation with a low diffusion numerical scheme[J]. Annals of Nuclear Energy, 2003, 30(15): 1505-1537. DOI:10.1016/S0306-4549(03)00119-1