一种稳定性增强及高精度数值方法在RELAP5中的实现与评价

李 伟 苏光辉 秋穗正 田文喜 陈 平 李垣明

1（中国核动力研究设计院 核反应堆系统设计技术重点实验室 成都 610213）

2（西安交通大学 核科学与技术学院 西安 710049）

一种稳定性增强及高精度数值方法在RELAP5中的实现与评价

李 伟1苏光辉2秋穗正2田文喜2陈 平1李垣明1

1（中国核动力研究设计院 核反应堆系统设计技术重点实验室 成都 610213）

2（西安交通大学 核科学与技术学院 西安 710049）

在计算稳定性的改进方面，修正RELAP5程序中的虚拟质量力(Virtual mass force)形式，同时添加了新的界面压力项(Interface pressure)；在计算精度的改进方面，采用具有总变差减小(Total variation diminish, TVD)特点的高精度通量限制器(Flux limiter)方法取代RELAP5程序原本的一阶迎风方法来离散质量和能量守恒方程中的对流项。在模拟水平管道内空泡份额微扰随时间发展的数值实验中，相比改进前的RELAP5，改进后的RELAP5计算得到的微扰幅度并未增长；在模拟液相沉降的数值实验中，改进前的RELAP5程序计算得到了不真实的空泡份额分布，而改进后的RELAP5在不同的网格数量下能够得到收敛的稳定解。对Ransom 水龙头数值实验和Marviken CFT 15大破口喷放实验的计算表明，改进后的具有二阶TVD格式的RELAP5程序能够得到更接近实验数据的计算结果。

RELAP5，稳定性，精确性，双曲性，虚拟质量力，界面压力项

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.1RELAP5两流体模型的非双曲性

为分析RELAP5两流体模型的非双曲性，将其6个守恒方程改写为如下的矩阵方程形式：

表示RELAP5的直接求解变量，分别为压力、空泡份额、气相速度、液相速度、气相内能和液相内能；A、B和C为包含未知量的矩阵[6×6]。设为式(1)的特征根，它们满足以下特征方程：

只有以下条件成立时[5]，特征根才全部为实数：

或：

1.2虚拟质量力的影响

虚拟质量力[6]在多相流问题中刻画了离散相流体在连续相流体中运动时，带动其周围的连续相流体运动所需要的额外力。RELAP5采用如下形式的虚拟质量力：

式中：ma具有气液混合物等效音速的意义。如果ε在数量级上足够小，即气、液相的速度大小相对值远小于混合物音速时，经分析须满足如下的条件才能保证所有的特征根均为实数：

在大部分情况下，ε都可以被认为是数量级很小的量。当气、液相速度差与音速相当时，Tiselj等[7]认为无论虚拟质量力系数如何取值也无法使得方程组具备双曲性。本文在采用式(10)的虚拟质量力系数的基础上，进一步研究界面压力项对两流体模型适定性的影响。

1.3界面压力项的影响

在RELAP5的模型假设中，气、液和相间界面的压力相等，因此也称为单压力模型，这造成了守恒方程中某些重要信息的缺失，是导致模型不具备适定性的重要因素。在另一个著名的反应堆热工水力系统分析程序CATHARE中，通过在动量守恒方程中引入相间界面压力项[8]来考虑各相压力的差异，即：

式中：iPΔ为界面压力项系数。界面压力项在动量方程中引入了空泡份额的空间偏微分项。同时存在虚拟质量力和界面压力项时，方程组的特征方程化解为：

式(14)判别式非负，求得界面压力项系数须满足如下关系式：

从物理量纲来看iPΔ具有动压的意义。当虚拟质量力的系数为零时，式(15)与CATHARE程序中的形式一致。当虚拟质量力的系数满足式(10)时，即时，无需界面压力项方程组也能满足双曲性要求，否则界面压力项能够对虚拟质量力的效果起到补充作用（例如在水平分层流动中虚拟质量力系数为零的情况下）。虽然虚拟质量力和界面压力项的目的都与改善RELAP5两流体模型的适定性有关，但是二者之间存在一定的区别：虚拟质量力在两相流中是具有对应物理意义的，而界面压力修正项则缺乏坚实的物理根据。它虽然考虑到单压力模型中界面压力与各相平均压力之间并不相等，但是真正的界面压力非常复杂，因此采用经验关系式近似。不可把式(11)和(12)看作是真实物理过程中的界面压力，可以理解为使两流体模型的方程组具有双曲性而人工添加的微分项，所以为防止过分强调其影响，一般在合适的取值范围内取最小的值。在程序的实际实施过程中，还通常对iPΔ乘以一个稍大于1的系数ξ。此外，RELAP5在水平分层流动时考虑到了界面压力与各相平均压力的不同而使用了这种压力修正项，CATHARE则在所有流型中予以使用，以改善基本守恒方程的双曲性。

2 计算精度改进

RELAP5采用了一阶精度的数值方案进行对流项的离散，大量实验和工程应用证明这对大部分系统瞬态行为的预测能够达到相应的精度要求，同时也简化了程序的编程难度，提高了计算效率（特别是在进行不确定性分析时）。RELAP5采取的是交错网格方案，质量和能量守恒方程中的对流项采用一阶迎风格式[9]进行离散：

采用Waterson[11]关于Flux limiter的通用概念，对流项的离散格式可表达为低阶（一阶迎风）和高阶格式的组合：

式中：KxΔ(LxΔ)为标量控制体单元的长度；/xφ∂∂为标量控制体单元界面上的参数梯度；()rψ为Flux limiter函数衡量控制体单元界面j上参数的平滑程度。Flux limiter函数的具体形式将会对程序的计算精度和稳定性产生重要影响。合适的Flux limiter函数形式应该符合Flux limiter和TVD格式的初衷，即提高程序计算精度的同时，保证程序具有合理的稳定性。Flux limiter的形式也不宜过于复杂，以致降低整个程序的运行效率。出于实际计算中的稳定性和精度考虑，本文采用Minmod flux limiter函数

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) 改进后的RELAP5Fig.1 Time-dependent distribution of the void fraction.(a) Default RELAP5, (b) Improved RELAP5

图2 改进后的RELAP5不考虑虚拟质量力时空泡份额分布随时间的变化Fig.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 沉降数值实验示意图Fig.3 Sketch of the phase sedimentation numerical test.

图4 控制体1内的空泡份额随时间变化 (a) 改进前的RELAP5，(b) 改进后的RELAP5Fig.4 Variation of the void fraction in CV #1 with time. (a) Default RELAP5, (b) Improved RELAP5

3.2计算精度改进措施验证

采用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在假设液相不可压缩的基础上推导了该问题的精确解：

式中：0gα为初始空泡份额；gα为空泡份额；z为距管道入口的距离，m；0fν为液相初始速度，m·s-1；H为管道总长度。

图5 Ransom水龙头数值实验示意图Fig.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个控制体Fig.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控制体划分方案Fig.7 Nodalization scheme for the Marviken CFT 15 experiment.

图8给出了改进前/后的RELAP5计算得到的CFT 15破口流量与实验数据的比较。从图8中可见，改进后的RELAP5计算得到的破口流量计算值与实验数据符合得很好，只在两相喷放阶段出现了一定的偏差。如图8(b)所见，改进后的RELAP5由于精度提高，计算稳定性增强，能够更好地预测从破口喷放的流量随时间的变化。

图8 改进前的RELAP5 (a)和改进后的RELAP5 (b)计算得到的破口流量Fig.8 Break flow rate calculated by the default RELAP5 (a) and improved RELAP5 (b).

4 结语

为改善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

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

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

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

5 张小英, 丁斐, 陈佳跃. 反应堆一维两流体模型二阶精度数值解法研究[J]. 核动力工程, 2013,34(4): 27-32. DOI: 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.007

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

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

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

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

10 Versteeg H K, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method[M]. Essex, UK: Pearson Education, 2007

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

12 Ransom V. Numerical benchmark tests[A]. Hewitt G, Delhay J, Zuber N. Multiphase science and technology[M]. Washington DC: Hemisphere Publishing Corporation, 1987

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

Implementation and assessment of a stability-enhancing and high-resolution numerical scheme in RELAP5

LI Wei1SU Guanghui2QIU Suizheng2TIAN Wenxi2CHEN Ping1LI Yuanming1
1(Science Technology on Reactor System Design Technology Laboratory,Nuclear Power Institute of China,Chengdu 610213,China)
2(School of Nuclear Science and Technology,Xi’an Jiaotong University,Xi’an 710049,China)

Background：RELAP5 has been widely used in the field of nuclear thermal-hydraulics and safety analysis. However, improvement issues in terms of numerical stability and accuracy exist due to its basic single-pressure two-fluid model.Purpose：This study aims to enhance the numerical stability and numerical accuracy of RELAP5 by modifying the mathematical structure of RELAP5 basic model.Methods：Firstly, the default mathematical expression of the virtual mass force was modified in the RELAP5 code. Secondly, an interface pressure term was incorporated in the code, and a total variation diminish (TVD)-type flux limiter method was implemented in the RELAP5 code in place of the default 1st-order upwind scheme for the convective terms in mass and energy conversation equations.Results：In a numerical test where the evolvements of void fraction perturbation were concerned, the modified RELAP5 predicted no growth of the perturbation amplitude in contrast to a rapidly developing divergence for the default RELAP5. In addition, for the phase segmentation problem, convergence was accomplished on finer mesh by the modified RELAP5, while non-physical distribution of the void fraction was rendered by the default RELAP5 especially for finer mesh. When applied to the Ransom water faucet numerical test and Marviken CFT 15 experiments, the improved RELAP5 with 2nd-order of accuracy was more reliable in that thepredictions were closer to the experimental results.Conclusion：Compared with the original RELAP5, the modified or improved RELAP5 showed evident numerical stability and accuracy during the assessment in this work.

RELAP5, Stability, Accuracy, Hyperbolicity, Virtual mass, Interface pressure

TL33

10.11889/j.0253-3219.2016.hjs.39.110601

李伟，男，1986年出生，2015年于西安交通大学核科学与技术专业获博士学位，主要从事核燃料设计、反应堆热工水力及安全分析First author: LI Wei, male, born in 1986, graduated from Xi’an Jiaotong University with a doctoral degree in 2015, major in nuclear science and technology, mainly focusing on nuclear fuel research and design, reactor thermal-hydraulic and safety analysis

2016-05-30，

2016-09-01