2. 西安交通大学 核科学与技术学院 西安 710049
2. School of Nuclear Science and Technology, Xi'an Jiaotong University, Xi'an 710049, China
RELAP5程序长期以来在反应堆系统的热工水力及安全分析领域都有着广泛的应用，包括乏燃料池(Spent fuel pool)的安全设计[1]。但RELAP5的基本偏微分方程组包含虚的特征根，因此不具备双曲性，这种两流体模型不具备适定性[2-3]。在实际应用中可能导致由不适定性引起的数值解发散，出现不真实的计算结果。在数值离散方面，RELAP5作为系统程序采用了稀疏的交错网格，动量守恒方程中的对流项采用一阶迎风方法进行离散。这有利于增强程序的数值稳定性和提高计算效率，然而也使得数值扩散过大，影响计算精度。
本文基于RELAP5/MOD3程序，从改善计算稳定性和提高数值精度方面对其进行性能改进并且通过实验验证。在计算稳定性改进方面，分析了动量守恒方程中不同的虚拟质量力形式以及界面压力项对RELAP5六方程特征根的影响，修正了RELAP5原本的虚拟质量力形式和增加界面压力项使其特征根全部为实数来增强模型的双曲性。设计一个空泡份额微扰随时间变化的数值实验来验证改进后的RELAP5程序在计算稳定性方面的性能；通过一个相沉降数值实验来验证改进后的RELAP5程序预测结果在不同网格划分方案下的行为。在数值精度改进方面，Wang等[4]研究了二阶总变差减小(Total variation diminish,TVD)格式在系统程序TRACE中的实现。TRACE是系统程序TRAC和RELAP5组合后的版本。发现二阶TVD格式比一阶迎风格式的计算更准确，且与一阶迎风格式一样高效和稳定。张小英等[5]开发了二阶精度的一维两流体模型程序TFIT。该程序未考虑能量方程，用于绝热不可压缩问题。通过对Ransom水龙头问题进行对比计算，二阶TVD格式比一阶精度格式的计算准确度明显提高，而且同样保持良好的稳定性。本文采用结合了Minmod flux limiter的二阶TVD格式代替RELAP5程序中对流项的一阶迎风离散格式，用于实际工程计算。基于Ransom水龙头数值实验和实际工程试验Marviken CFT 15对改进后的程序进行验证，评估对数值扩散问题的改善。
1 计算稳定性改进 1.1 RELAP5两流体模型的非双曲性为分析RELAP5两流体模型的非双曲性，将其6个守恒方程改写为如下的矩阵方程形式：
| $A(\phi )\frac{\partial \phi }{\partial t}+B(\phi )\frac{\partial \phi }{\partial x}+C(\phi )=0$ | (1) |
式中：$\phi ={{(P\text{,}\ {{\alpha }_{g}}\text{,}\ {{\upsilon }_{g}}\text{,}\ {{\upsilon }_{f}}\text{,}\ {{U}_{g}}\text{,}\ {{U}_{f}})}^{T}}$表示RELAP5的直接求解变量，分别为压力、空泡份额、气相速度、液相速度、气相内能和液相内能；A、B和C为包含未知量的矩阵[6×6]。设${{\lambda }_{i}}\text{(}i=1\text{,}\,6\text{)}$为式(1)的特征根，它们满足以下特征方程：
| $det(B-{{\lambda }_{i}}A)=0$ | (2) |
式中：det表示求矩阵的行列式。当${{\lambda }_{i}}\text{(}i=1\text{,}\,6\text{)}$全部为实数时，式(1)方程组具有双曲性；当${{\lambda }_{i}}\text{(}i=1\text{,}6\text{)}$至少有一个为包含虚部的复数时，方程组不具备双曲性。式(2)为6次多项式方程，通过演算直接获取${{\lambda }_{i}}$的解析解将十分困难。通过简化，以下情况表明RELAP5六方程两流体模型确实存在虚的特征根，从而不是双曲性的。
假设液相不可压缩，气相可压缩，则除了重根${{\lambda }_{\text{1,2}}}={{\upsilon }_{g}}\text{,}{{\upsilon }_{f}}$，还包含以下一元四次方程的4个根：
| $\begin{align} & \frac{{{\alpha }_{g}}}{c_{g}^{2}}(\frac{P}{{{\rho }_{f}}}\frac{\partial {{\rho }_{f}}}{\partial {{U}_{f}}}-{{\rho }_{f}}){{({{\lambda }_{i}}-{{\upsilon }_{g}})}^{2}}{{({{\lambda }_{i}}-{{\upsilon }_{f}})}^{2}}+ \\ & \quad \begin{matrix} {} & {} \\ \end{matrix}{{\rho }_{g}}{{\alpha }_{f}}{{({{\lambda }_{i}}-{{\upsilon }_{g}})}^{2}}+{{\rho }_{f}}{{\alpha }_{g}}{{({{\lambda }_{i}}-{{\upsilon }_{f}})}^{2}}=0 \\ \end{align}$ | (3) |
只有以下条件成立时[5]，特征根才全部为实数：
| ${{\upsilon }_{g}}={{\upsilon }_{f}}$ | (4) |
或：
| ${{({{\upsilon }_{g}}-{{\upsilon }_{f}})}^{2}}\ge c_{g}^{2}\frac{{{[{{({{\rho }_{f}}{{\alpha }_{g}})}^{1/3}}+{{({{\rho }_{g}}{{\alpha }_{f}})}^{1/3}}]}^{3}}}{{{\alpha }_{g}}({{\rho }_{f}}-P\text{/}{{\rho }_{f}}\cdot \partial {{\rho }_{f}}\text{/}\partial {{U}_{f}})}$ | (5) |
式中：$c_{g}^{2}={1}/{\text{(}\partial {{\rho }_{g}}\text{/}\partial P\text{)}}\;$，$c_{f}^{2}=1/\text{(}\partial {{\rho }_{f}}\text{/}\partial P\text{)}$，二者具有当地音速的意义。当气、液两相速度大小相差较小时，式(5)将不再成立，因此方程组将包含虚数根，不具有双曲性。为改善RELAP5方程组的双曲性，可在动量守恒方程中根据物理意义合理地添加微分项，例如虚拟质量力或者界面压力项，以改变特征方程的结构形式，从而避免虚数特征根的出现。
1.2 虚拟质量力的影响虚拟质量力[6]在多相流问题中刻画了离散相流体在连续相流体中运动时，带动其周围的连续相流体运动所需要的额外力。RELAP5采用如下形式的虚拟质量力：
| $\begin{align} & F_{vm}^{g}=-{{C}_{vm}}{{\rho }_{m}}{{\alpha }_{g}}{{\alpha }_{f}}[\frac{\partial ({{\upsilon }_{g}}-{{\upsilon }_{f}})}{\partial t}+{{\upsilon }_{g}}\frac{\partial ({{\upsilon }_{g}}-{{\upsilon }_{f}})}{\partial x}+ \\ & \ \ \ \ \ \ \ \ \ (\mu -2)({{\upsilon }_{g}}-{{\upsilon }_{f}})\frac{\partial {{\upsilon }_{g}}}{\partial x}+(1-\mu )({{\upsilon }_{g}}-{{\upsilon }_{f}})\frac{\partial {{\upsilon }_{f}}}{\partial x}] \\ \end{align}$ | (6) |
式中：$F_{vm}^{g}$为气相动量守恒方程中的虚拟质量力项；Cvm为虚拟质量力系数。在RELAP5程序中，实际采用的虚拟质量力对应于μ=1且忽略空间导数项，此时须满足以下的条件才能保证所有的特征根均为实数：
| $\begin{align} & C_{vm}^{2}\rho _{m}^{2}{{\text{(}{{\alpha }_{f}}{{\upsilon }_{g}}+{{\alpha }_{g}}{{\upsilon }_{f}}\text{)}}^{2}}\text{/}4+ \\ & \quad {{C}_{vm}}{{\rho }_{m}}\text{(}{{\upsilon }_{g}}-{{\upsilon }_{f}}\text{)(}{{\rho }_{f}}{{\upsilon }_{f}}-{{\rho }_{g}}{{\upsilon }_{g}}\text{)}{{\alpha }_{g}}{{\alpha }_{f}}- \\ & \quad {{\rho }_{g}}{{\rho }_{f}}{{\alpha }_{g}}{{\alpha }_{f}}{{\text{(}{{\upsilon }_{g}}-{{\upsilon }_{f}}\text{)}}^{2}}\ge 0 \\ \end{align}$ | (7) |
式中：${{\rho }_{m}}={{\rho }_{g}}{{\alpha }_{g}}+{{\rho }_{f}}{{\alpha }_{f}}$为两相混合密度。可见如果不存在虚拟质量力或者虚拟质量力系数过小，则式(7)不成立，在某些条件下改进前的RELAP5模型确实存在不适定的问题。当保留式(6)中的空间偏微分项时，分析特征根的形式变得更加复杂。分别定义如下参数：
| $\begin{align} & {{a}_{m}}={{\left( \frac{{{\rho }_{g}}{{\alpha }_{f}}+{{\rho }_{f}}{{\alpha }_{g}}+{{C}_{vm}}{{\rho }_{m}}}{{{\rho }_{g}}{{\rho }_{f}}+{{C}_{vm}}\rho _{m}^{2}} \right)}^{1/2}}\times \\ & \ {{\left[ \frac{(\rho _{g}^{2}-P\frac{\partial {{\rho }_{g}}}{\partial {{U}_{g}}})(\rho _{f}^{2}-P\frac{\partial {{\rho }_{f}}}{\partial {{U}_{f}}})}{\frac{{{\rho }_{f}}{{\alpha }_{f}}}{c_{f}^{2}}(\rho _{g}^{2}-P\frac{\partial {{\rho }_{g}}}{\partial {{U}_{g}}})+\frac{{{\rho }_{g}}{{\alpha }_{g}}}{c_{g}^{2}}(\rho _{f}^{2}-P\frac{\partial {{\rho }_{f}}}{\partial {{U}_{f}}})} \right]}^{\frac{1}{2}}} \end{align}$ | (8) |
| $\varepsilon =\frac{{{\upsilon }_{r}}}{{{a}_{m}}}=\frac{{{\upsilon }_{g}}-{{\upsilon }_{f}}}{{{a}_{m}}}$ | (9) |
式中：am具有气液混合物等效音速的意义。如果ε在数量级上足够小，即气、液相的速度大小相对值远小于混合物音速时，经分析须满足如下的条件才能保证所有的特征根均为实数：
| ${{C}_{vm}}\ge {2{{({{\alpha }_{g}}{{\alpha }_{f}}{{\rho }_{g}}{{\rho }_{f}})}^{1/2}}}/{{{\rho }_{m}}}\;$ | (10) |
在大部分情况下，ε都可以被认为是数量级很小的量。当气、液相速度差与音速相当(${{\upsilon }_{r}}\ge 0.3{{a}_{m}}$)时，Tiselj等[7]认为无论虚拟质量力系数如何取值也无法使得方程组具备双曲性。本文在采用式(10)的虚拟质量力系数的基础上，进一步研究界面压力项对两流体模型适定性的影响。
1.3 界面压力项的影响在RELAP5的模型假设中，气、液和相间界面的压力相等，因此也称为单压力模型，这造成了守恒方程中某些重要信息的缺失，是导致模型不具备适定性的重要因素。在另一个著名的反应堆热工水力系统分析程序CATHARE中，通过在动量守恒方程中引入相间界面压力项[8]来考虑各相压力的差异，即：
| $F_{p}^{g}=-\Delta {{P}_{i}}\frac{\partial {{\alpha }_{g}}}{\partial x}$ | (11) |
| $F_{p}^{f}=-\Delta {{P}_{i}}\frac{\partial {{\alpha }_{f}}}{\partial x}$ | (12) |
式中：$\Delta {{P}_{i}}$为界面压力项系数。界面压力项在动量方程中引入了空泡份额的空间偏微分项。同时存在虚拟质量力和界面压力项时，方程组的特征方程化解为：
| $\begin{align} & \alpha _{g}^{2}{{(1-{{\alpha }_{g}})}^{2}}({{\upsilon }_{g}}-{{\lambda }_{i}})({{\upsilon }_{f}}-{{\lambda }_{i}})[{{a}_{1}}{{({{\upsilon }_{g}}-{{\lambda }_{i}})}^{2}}{{({{\upsilon }_{f}}-{{\lambda }_{i}})}^{2}}+ \\ & \ {{b}_{1}}({{\upsilon }_{f}}-{{\lambda }_{i}}){{({{\upsilon }_{g}}-{{\lambda }_{i}})}^{3}}+{{c}_{1}}({{\upsilon }_{g}}-{{\lambda }_{i}}){{({{\upsilon }_{f}}-{{\lambda }_{i}})}^{3}}+ \\ & \ {{d}_{1}}{{({{\upsilon }_{g}}-{{\lambda }_{i}})}^{2}}+{{e}_{1}}{{({{\upsilon }_{f}}-{{\lambda }_{i}})}^{2}}+\ {{f}_{1}}({{\upsilon }_{g}}-{{\lambda }_{i}})({{\upsilon }_{f}}-{{\lambda }_{i}})+ \\ & \ {{g}_{4}}{{({{\upsilon }_{g}}-{{\lambda }_{i}})}^{2}}+{{h}_{4}}{{({{\upsilon }_{f}}-{{\lambda }_{i}})}^{2}}+{{w}_{4}}]=0 \\ \end{align}$ | (13) |
式中：a1、b1、c1、d1、f1、g4、h4、w4为系数，其中仅g4、h4、w4包含界面压力项系数$\Delta {{P}_{i}}$。与虚拟质量力的情况不同的是，界面压力项使一元四次方程中新增了负的独立项w4，从几何的角度来看这相当于将一元四次方程所代表的曲线沿着负y轴方向移动了|w4|单位长度。因此要求在保持方程双曲性的同时，界面压力项的值不能过大，否则将掩盖真实的物理过程，令计算出现不真实的解。无论通过直接求解析解或者作图法，判断式(13)的特征根是否全部为实数都十分困难。这里假设两相的可压缩性均可忽略，则方程大为简化，判断以下一元二次方程是否有实根：
| $\begin{align} & {{\rho }_{g}}{{\alpha }_{f}}{{({{\lambda }_{i}}-{{\upsilon }_{g}})}^{2}}+{{\rho }_{f}}{{\alpha }_{g}}{{({{\lambda }_{i}}-{{\upsilon }_{f}})}^{2}}+ \\ & \ \ \quad \ {{C}_{vm}}{{\rho }_{m}}({{\lambda }_{i}}-{{\upsilon }_{g}})({{\lambda }_{i}}-{{\upsilon }_{f}})-\Delta {{P}_{i}}=0 \\ \end{align}$ | (14) |
式(14)判别式非负，求得界面压力项系数须满足如下关系式：
| $\Delta {{P}_{i}}\ge \frac{{{\alpha }_{g}}{{\alpha }_{f}}{{\rho }_{g}}{{\rho }_{f}}-0.25C_{vm}^{2}\rho _{m}^{2}}{{{\alpha }_{f}}{{\rho }_{g}}+{{\alpha }_{g}}{{\rho }_{f}}+{{C}_{vm}}{{\rho }_{m}}}{{({{\upsilon }_{g}}-{{\upsilon }_{f}})}^{2}}$ | (15) |
从物理量纲来看$\Delta {{P}_{i}}$具有动压的意义。当虚拟质量力的系数为零时，式(15)与CATHARE程序中$\Delta {{P}_{i}}$的形式一致。当虚拟质量力的系数满足式(10)时，即${{\alpha }_{g}}{{\alpha }_{f}}{{\rho }_{g}}{{\rho }_{f}}-0.25C_{vm}^{2}\rho _{m}^{2}\le 0$时，无需界面压力项方程组也能满足双曲性要求，否则界面压力项能够对虚拟质量力的效果起到补充作用（例如在水平分层流动中虚拟质量力系数为零的情况下）。虽然虚拟质量力和界面压力项的目的都与改善RELAP5两流体模型的适定性有关，但是二者之间存在一定的区别：虚拟质量力在两相流中是具有对应物理意义的，而界面压力修正项则缺乏坚实的物理根据。它虽然考虑到单压力模型中界面压力与各相平均压力之间并不相等，但是真正的界面压力非常复杂，因此采用经验关系式近似。不可把式(11)和(12)看作是真实物理过程中的界面压力，可以理解为使两流体模型的方程组具有双曲性而人工添加的微分项，所以为防止过分强调其影响，一般在合适的取值范围内取最小的值。在程序的实际实施过程中，还通常对$\Delta {{P}_{i}}$乘以一个稍大于1的系数ξ。此外，RELAP5在水平分层流动时考虑到了界面压力与各相平均压力的不同而使用了这种压力修正项，CATHARE则在所有流型中予以使用，以改善基本守恒方程的双曲性。
2 计算精度改进RELAP5采用了一阶精度的数值方案进行对流项的离散，大量实验和工程应用证明这对大部分系统瞬态行为的预测能够达到相应的精度要求，同时也简化了程序的编程难度，提高了计算效率（特别是在进行不确定性分析时）。RELAP5采取的是交错网格方案，质量和能量守恒方程中的对流项采用一阶迎风格式[9]进行离散：
| $\overset{\bullet }{\mathop{{{\phi }_{j}}}}\,=\frac{1}{2}({{\phi }_{K}}+{{\phi }_{L}})+\frac{1}{2}\frac{\left| {{\upsilon }_{j}} \right|}{{{\upsilon }_{j}}}({{\phi }_{K}}-{{\phi }_{L}}),\ \ {{\upsilon }_{j}}\ne 0$ | (16) |
式中：${{\dot{\phi }}_{j}}$也被称为受赠项(Donored quantities)，定义在标量控制体单元的面上；$\phi _{K}^{{}}$和${{\phi }_{L}}$为标量控制体单元中心的平均量；${{\upsilon }_{j}}$为流体速度。一阶迎风格式无条件稳定，这有助于保持RELAP5的计算稳定性。但同时该格式也存在明显的数值扩散问题，尤其在稀疏的网格方案中更加无法准确捕捉流场中的不连续性特征。如果完全采用高阶离散格式，虽然能够提高流动平滑区域的计算精度，但是在流场中的不连续点附近计算容易变得不稳定，数值解呈现出振荡或“毛刺”的行为。除了常见的对流项离散格式，如一阶迎风、二阶迎风和QUICK格式等，混合型的离散格式（TVD flux limiter格式）在空气动力学的数值模拟中已经取得成功的应用，近年来也越来越多地出现在能源化工等计算流体动力学(Computational Fluid Dynamics,CFD)模拟领域[10]。线性的TVD最多表现为一阶精度，与非线性Flux limiter结合的TVD，则既具有高阶精度又尽可能避免出现计算发散。TVD flux limiter格式通常用于无粘单相可压缩流问题，本文尝试将该方法应用于反应堆热工水力系统中常见的气液两相流问题，在间断附近利用一阶精度格式的特点提高计算稳定性，在光滑区域利用高精度格式的特点提高程序的计算精度。
采用Waterson[11]关于Flux limiter的通用概念，对流项的离散格式可表达为低阶（一阶迎风）和高阶格式的组合：
| $\overset{\bullet }{\mathop{\phi _{j}^{n}}}\,=\left\{ \begin{matrix} \phi _{K}^{n}+\frac{\Delta {{x}_{K}}}{2}\psi (r_{j,j-1}^{n})\left( \frac{\partial \phi }{\partial x} \right)_{j-1}^{n}\ \ \ \ \ \ \ \ \upsilon _{j}^{n}>0 \\ \phi _{L}^{n}+\frac{\Delta {{x}_{L}}}{2}\psi (r_{j,j+1}^{n})\left( \frac{\partial \phi }{\partial x} \right)_{j+1}^{n}\ \ \ \ \ \ \ \ \upsilon _{j}^{n} ＜0 \\ \end{matrix} \right.$ | (17) |
| $\begin{matrix} r_{j,j-1}^{n}={{{\left( \frac{\partial \phi }{\partial x} \right)}_{j}}}/{{{\left( \frac{\partial \phi }{\partial x} \right)}_{j-1}}}\; & \upsilon _{j}^{n}>0 \\ \end{matrix}$ | (18) |
| $\begin{matrix} r_{j,j+1}^{n}={{{\left( \frac{\partial \phi }{\partial x} \right)}_{j}}}/{{{\left( \frac{\partial \phi }{\partial x} \right)}_{j+1}}\ }\; & \upsilon _{j}^{n} ＜0 \\ \end{matrix}$ | (19) |
式中：$\Delta {{x}_{K}}$($\Delta {{x}_{L}}$)为标量控制体单元的长度；$\partial \phi /\partial x$ 为标量控制体单元界面上的参数梯度；$\psi (r)$为Flux limiter函数；$r_{j,j-1}^{n}$($r_{j,j+1}^{n}$)衡量控制体单元界面j上参数的平滑程度。Flux limiter函数的具体形式将会对程序的计算精度和稳定性产生重要影响。合适的Flux limiter函数形式应该符合Flux limiter和TVD格式的初衷，即提高程序计算精度的同时，保证程序具有合理的稳定性。Flux limiter的形式也不宜过于复杂，以致降低整个程序的运行效率。出于实际计算中的稳定性和精度考虑，本文采用Minmod flux limiter函数：$\psi (r)=\text{max }\!\![\!\!\text{ 0,}\,\text{min(}r\text{,1) }\!\!]\!\!\text{ }$。
3 改进措施的验证 3.1 计算稳定性改进措施验证数值实验在RELAP5程序的开发与验证过程中起到重要作用。采用类似文献[8]中的空泡份额微扰实验对稳定性改进措施进行验证。在长度5 m的水平管道中，流动着0.1 MPa的饱和蒸汽和饱和水，速度分别保持0.1 m·s-1和1.0 m·s-1。除了距离入口2.0 m处的初始空泡份额为0.5008，其他位置的初始空泡份额均为0.5。进行数值模拟时，整个管道划分为99个控制体，不计壁面摩擦。在这种情况下，由于数值扩散的作用，空泡份额微扰将随时间而衰减。但是对于不具备双曲性的模型，扰动可能迅速增长，最终计算结果发散。在本例中虚拟质量力的效应相当微弱，因此对于改进前的RELAP5，式(7)可能无法得到满足，如图 1(a)所示，改进前的RELAP5计算得到的空泡份额值迅速发散。如图 1(b)所示，对于改进后的RELAP5在虚拟质量力和界面压力项的共同作用下模型保持了双曲性，微扰的增长受到抑制，保持了计算的稳定性。图 2为改进后的RELAP5在不考虑虚拟质量力时计算得到的空泡份额分布随时间的变化，可见微扰迅速衰减。本例说明虚拟质量力以及界面压力项的作用在程序中得到了实现，改进后的RELAP5相比改进前的RELAP5计算稳定性明显提高。
|
图 1 空泡份额分布随时间变化 (a) 改进前的RELAP5，(b) 改进后的RELAP5 Figure 1 Time-dependent distribution of the void fraction (a) Default RELAP5,(b) Improved RELAP5 |
|
图 2 改进后的RELAP5不考虑虚拟质量力时空泡份额分布随时间的变化 Figure 2 Time-dependent distribution of the void fraction without virtual mass term of improved RELAP5. |
考虑更复杂的情况，采用如图 3所示的液相沉降数值实验进行验证。利用改进后的RELAP5分别计算了控制体划分数目(Control Volume,CV)为25、49和99的情况，并且与改进前的RELAP5进行比较。图 4给出了改进前/后的RELAP5计算得到控制体1中的空泡份额计算值随时间的变化。随着控制体数目增多，控制体单元尺寸减小，改进前的RELAP5模型非双曲性、不适定的缺点表现得更明显。在精细的网格下更容易出现不稳定性，计算过程中也出现了不真实的结果。而改进后的RELAP5模型适定性得以改善，即使计算网格加密计算结果也不容易出现不稳定的行为，达到了改进的目的。
|
图 3 沉降数值实验示意图 Figure 3 Sketch of the phase sedimentation numerical test. |
|
图 4 控制体1内的空泡份额随时间变化 (a) 改进前的RELAP5，(b) 改进后的RELAP5 Figure 4 Variation of the void fraction in CV #1 with time. (a) Default RELAP5,(b) Improved RELAP5 |
采用Ransom[12]提出的水龙头(Water faucet)数值实验来对本文所作的精度改进进行验证。Water faucet问题常被用于验证两相流求解器的数值性能。如图 5所示，在长度12 m、直径1 m的竖直圆管中初始时刻充满常温常压静止的空气，初始时刻圆管中还存在一股射流，其速度处处为10 m·s-1，且轴向的空泡份额均为0.2。管道入口空泡份额和液体速度分别保持为0.2 m·s-1和10 m·s-1，出口保持常压。计算开始后，水柱各部分在重力的作用下加速运动。不计管道壁面摩擦，整个过程绝热。Ransom在假设液相不可压缩的基础上推导了该问题的精确解：
| ${{\alpha }_{g}}=\left\{ \begin{matrix} {{\alpha }_{g0}}{{\text{(}1+2gz\upsilon _{f0}^{-2}\text{)}}^{-1/2}} & z ＜{{\upsilon }_{f0}}t+g{{t}^{2}}\text{/}2 \\ {{\alpha }_{g0}} & {{\upsilon }_{f0}}t+g{{t}^{2}}\text{/}2\le z\le H \\ \end{matrix} \right.$ | (20) |
式中：${{\alpha }_{g0}}$为初始空泡份额；${{\alpha }_{g}}$为空泡份额；z为距管道入口的距离，m；${{\upsilon }_{f0}}$为液相初始速度，m·s-1；H为管道总长度。
|
图 5 Ransom水龙头数值实验示意图 Figure 5 Sketch of the Ransom water faucet numerical test. |
分别划分96和120个控制体时，改进前/后计算得到的空泡份额分布与分析解的对比如图 6所示。与改进前的RELAP5一阶迎风格式(Upwind Difference,UD)相比，随着控制体数目的增多，二阶TVD格式(Minmod)更接近精确解(Analytic)。一阶迎风格式在控制体单元尺寸过大的情况下，数值扩散更加明显。从本次数值实验可以看出，改进后的RELAP5二阶精度格式体现了程序捕捉流场间断的优势。
|
图 6 空泡份额沿流动方向上的分布 (a) 96个控制体，(b) 120个控制体 Figure 6 Variation of the void fraction along the flow (a) 96 CVs are employed,(b) 120 CVs are employed |
采用Marviken CFT 15全尺寸大型临界流实验来进一步验证改进措施[13]。实验中存放气液混合物的压力容器高24.55 m，平均内径5.2 m，排出管长6.3 m，内径0.72 m，喷嘴直径0.5 m，长径比为3.6。采用的控制体划分方案如图 7所示，采用Henry-Fauske临界流喷放模型。
|
图 7 Marviken CFT 15控制体划分方案 Figure 7 Nodalization scheme for the Marviken CFT 15 experiment. |
图 8给出了改进前/后的RELAP5计算得到的CFT 15破口流量与实验数据的比较。从图 8中可见，改进后的RELAP5计算得到的破口流量计算值与实验数据符合得很好，只在两相喷放阶段出现了一定的偏差。如图 8(b)所见，改进后的RELAP5由于精度提高，计算稳定性增强，能够更好地预测从破口喷放的流量随时间的变化。
|
图 8 改进前的RELAP5 (a)和改进后的RELAP5 (b)计算得到的破口流量 Figure 8 Break flow rate calculated by the default RELAP5 (a) and improved RELAP5 (b). |
为改善RELAP5两流体模型非双曲性的缺点，通过分析特征根得到虚拟质量力和界面压力项，用于保证模型应满足的条件；为改善RELAP5两流体模型仅具有一阶精度的不足，将二阶TVD格式应用到程序中，同时为兼顾计算稳定性选择了能够获得二阶精度、形式简单的Minmod 函数作为Flux limiter替代RELAP5的一阶迎风格式。实验验证结果表明，对RELAP5的改进措施是可行、合理的，达到了改进程序性能的目的，为更好地发挥系统程序在安全分析领域的作用具有积极意义，对改善其他最佳估算程序的性能也具有借鉴意义。
| [1] |
Ge L, Wang H T, Zhang G L, et al. Thermal-hydraulic design and transient analysis of passive cooling system for CPR1000 spent fuel storage pool[J].
Nuclear Science and Techniques, 2016, 27 (1) : 8 .
DOI: 10.1007/s41365-016-0017-6 ( 0)
|
| [2] |
Lyczkowski R W, Gidaspow D, Solbrig C W, et al. Characteristics and stability analyses of transient one-dimensional two-phase flow equations and their finite difference approximations[J].
Nuclear Science and Engineering, 1978, 66 (3) : 378 –396.
DOI: 10.13182/NSE78-4 (0)
|
| [3] |
Lax P D. On the nation of hyperbolicity[J].
Communications on Pure and Applied Mathematics, 1980, 33 (3) : 395 –397.
DOI: 10.1002/cpa.3160330309 (0)
|
| [4] |
Wang D, Mahaffy J H, Staudenmeier J, et al. Implementation and assessment of high-resolution numerical methods in TRACE[J].
Nuclear Engineering and Design, 2013, 263 (0) : 327 –341.
DOI: 10.1016/j.nucengdes.2013.05.015 (0)
|
| [5] |
张小英, 丁斐, 陈佳跃. 反应堆一维两流体模型二阶精度数值解法研究[J].
核动力工程, 2013, 34 (4) : 27 –32.
DOI: 10.3969/j.issn.0258-0926.2013.04.007DOI:10.3969/j.issn.0258-0926.2013.04.007 ZHANG Xiaoying, DING Fei, CHEN Jiayue. Study on numerical solution of the two dimensional two fluid model of the reactor[J]. Nuclear Power Engineering, 2013, 34 (4) : 27 –32. DOI: 10.3969/j.issn.0258-0926.2013.04.007DOI:10.3969/j.issn.0258-0926.2013.04.007 ( 0)
|
| [6] |
Drew D, Cheng L, Lahey R. The analysis of virtual mass effects in two-phase flow[J].
International Journal of Multiphase Flow, 1979, 5 (4) : 233 –242.
DOI: 10.1016/0301-9322(79)90023-5 (0)
|
| [7] |
Tiselj I, Petelin S. Modelling of two-phase flow with second-order accurate scheme[J].
Journal of Computational Physics, 1997, 136 (2) : 503 –521.
DOI: 10.1006/jcph.1997.5778 (0)
|
| [8] |
Hogon L, Lee U, Kim K, et al. Application of hyperbolic two-fluids equations to reactor safety code[J].
Nuclear Engineering and Technology, 2003, 35 (1) : 45 –54.
(0)
|
| [9] |
RELAP5 Code Development Team. RELAP5/MOD3 code manual Vol I:code structure, system models, and solution methods[R]. Idaho Falls:Idaho National Laboratory, 2001
(0)
|
| [10] |
Versteeg H K, Malalasekera W.
An introduction to computational fluid dynamics:the finite volume method[M]. Essex, UK: Pearson Education, 2007 .
(0)
|
| [11] |
Waterson N P, Deconinck H. Design principles for bounded higher-order convection schemes-a unified approach[J].
Journal of Computational Physics, 2007, 224 (1) : 182 –207.
DOI: 10.1016/j.jcp.2007.01.021 (0)
|
| [12] |
Ransom V. Numerical benchmark tests[A]. Hewitt G, Delhay J, Zuber N. Multiphase science and technology[M]. Washington DC:Hemisphere Publishing Corporation, 1987
(0)
|
| [13] |
Kim K, Kim H J. Assessment of RELAP5/MOD2 critical flow model using Marviken test data 15 and 24[R]. Washington, DC:Office of Nuclear Regulatory Research, U.S. Nuclear Regulatory Commission, 1992
(0)
|