Science and Technology of Nuclear Installations

Volume 2017, Article ID 4252975, 14 pages

https://doi.org/10.1155/2017/4252975

## Implementation and Comparison of High-Resolution Spatial Discretization Schemes for Solving Two-Fluid Seven-Equation Two-Pressure Model

^{1}State Key Laboratory of Multiphase Flow in Power Engineering, School of Nuclear Science and Technology, Xi’an Jiaotong University, No. 28, Xianning West Road, Xi’an, Shaanxi 710049, China^{2}Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu 610041, China

Correspondence should be addressed to Fei Chao; nc.ude.utjx.uts@iefoahc and Jianqiang Shan; nc.ude.utjx.liam@nahsqj

Received 11 August 2017; Revised 22 October 2017; Accepted 5 November 2017; Published 29 November 2017

Academic Editor: Alejandro Clausse

Copyright © 2017 Pan Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

As compared to the two-fluid single-pressure model, the two-fluid seven-equation two-pressure model has been proved to be unconditionally well-posed in all situations, thus existing with a wide range of industrial applications. The classical 1st-order upwind scheme is widely used in existing nuclear system analysis codes such as RELAP5, CATHARE, and TRACE. However, the 1st-order upwind scheme possesses issues of serious numerical diffusion and high truncation error, thus giving rise to the challenge of accurately modeling many nuclear thermal-hydraulics problems such as long term transients. In this paper, a semi-implicit algorithm based on the finite volume method with staggered grids is developed to solve such advanced well-posed two-pressure model. To overcome the challenge from 1st-order upwind scheme, eight high-resolution total variation diminishing (TVD) schemes are implemented in such algorithm to improve spatial accuracy. Then the semi-implicit algorithm with high-resolution TVD schemes is validated on the water faucet test. The numerical results show that the high-resolution semi-implicit algorithm is robust in solving the two-pressure two-fluid two-phase flow model; Superbee scheme and Koren scheme give two highest levels of accuracy while Minmod scheme is the worst one among the eight TVD schemes.

#### 1. Introduction

In many industrial applications especially in nuclear industries, two-phase flows exist widely and are the most important phenomenon. Accurate calculation of two-fluid flows is a subject of intense interest and of great importance in research topics. There are many important models in literature for describing two-phase flows. Herein, the two-fluid six-equation model seems to be the most complete approximation for two-phase flows. Hence, current main reactor thermal-hydraulics analysis codes such as RALAP5, TRACE, TRAC, and CATHARE are all based on such six-equation model.

However, due to an instantaneous equilibrium pressure assumption, such six-equation model has been proved to be ill-posed, which means that the initial value problem of the six-equation system is nonhyperbolic and could lead to numerical oscillations. From the book of De Bertodano et al. [1], which is of great reference value for investigating the stability of two-fluid model, the oscillations are caused by the Kelvin–Helmholtz instability (KH). When the relative velocity exceeds the Kelvin–Helmholtz instability criterion, the oscillations occur. In this book, De Bertodano et al. [1] show how the KH instability behaves in horizontal stratified flow for a well-posed fixed-flux model with surface tension compared to an ill-posed one without surface tension. For the well-posed model, the short-wavelength ripples die out and the large wave grows with time while for the ill-posed model the short-wavelength has a larger growth rate and dominates the solution after a short time.

For simulating complicated two-phase flow phenomenon of many nuclear power plant accidents, the ill-posedness of the differential equations presents a problem for higher order numerical methods or finer grids. Phasic vapor/liquid velocities are generally different, especially in the (fast) transients of nuclear power plant accidents where two-phase flow phenomenon is very complicated and phasic relative velocity may exceed Kelvin–Helmholtz criterion. To overcome the ill-posed issue of two-fluid single-pressure model and obtain well-posed numerical model, there are three important ideas: implementation of the interfacial pressure term used in CATHARE [2], implementation of the virtual mass force term used in RELAP5 [3], and application of the two-pressure model. The interfacial pressure differential term and the virtual mass force differential term are adding to phasic momentum equations to restore the hyperbolicity. The present authors investigated the ill-posed characteristic and analyze ill-posed regions of the two-fluid single-pressure model and the effect of the virtual mass force and the interfacial pressure on improving the ill-posedness [4]. The results showed that such two-fluid six-equation single-pressure model cannot completely avoid the ill-posedness with the virtual mass force and the interfacial pressure, only the two-fluid two-pressure model in which each phase is assumed to have its own pressure is a well-posed model in all situations.

Many researchers carried out the study of two-pressure model due to the well-posed advantage since 1976. The most important two-pressure model is the two-fluid seven-equation model to which much attention has been devoted [5–14]. Such seven-equation two-pressure model consists of two mass conservation equations, two momentum conservation equations, two energy conservation equations, and a volume fraction transport equation. Using the volume fraction propagation as a basic conservation equation is widely used because it changes the structure of conservation equation system and makes seven-equation model unconditionally hyperbolic (well-posed) in the sense of Hadamard [15]. In recent years, the development of advanced thermal-hydraulics system analysis codes has aroused great interest such as RELAP-7 [16], NEPTUNE project [17], CASL [18], and CATHARE-3 [19]. RELAP7 code developed by the Idaho National Laboratory is ongoing now. It should be noted that RELAP7 code is based on such seven-equation model.

In the last few decades, a volume of work has been conducted on the numerical computation of this seven-equation two-pressure model. Liang et al. [20] and Zein et al. [21] presented the operator splitting approach to decompose the seven-equation model into the hyperbolic operator and the relaxation operators. They used Godunov scheme with the HLLC flux to gain the numerical scheme for solving such seven-equation model. Gallouët et al. [22] proposed two finite volume methods based on Rusanov scheme and Godunov scheme to solve such two-pressure model. The source terms such as relaxation terms, the phase change, and gravity were calculated by a fractional step approach in his scheme. Ambroso et al. [23] constructed a new approximate Riemann solver for the numerical approximation of the seven-equation model. Berry et al. [6] and Abgrall and Saurel [24] developed the discrete equation method (DEM) with a HLLC scheme to solve the seven-equation model. They analyzed two pressure relaxation cases: infinitely fast and bounded with a realistic specific interfacial area. Moreover, Berry et al. took into account the interphase heat and mass transfer in this two-pressure model. The semi-implicit method with the staggered grid mesh is widely used in existing reactor thermal-hydraulics analysis codes (RELAP5, TRAC, TRACE, etc.) due to the advantage of high efficiency and stability. However, there is little information about the semi-implicit scheme to solve the seven-equation two-pressure model in existing public literature. Consequently, more investigations are essential to studying the semi-implicit algorithm for solving such two-pressure model.

In addition to great interest in the development and simulation of advanced two-fluid model, high-order accuracy schemes have also attracted great increasing attention. Many current reactor thermal-hydraulics codes like RELAP5 and CATHARE were developed based on classical first-order upwind scheme. Using such low-order scheme to make the convection terms of conservation equations discrete produces excessive numerical diffusion. This disadvantage has been realized in many nuclear thermal-hydraulics applications such as hydraulic load analysis of loss of coolant accident [25], long term transient natural circulation flow [26, 27], flow instability [28] or stability analysis [29], and boron solute transport [30]. In the past, there are many classical linear schemes like central-difference scheme (CD), QUICK, third-order upwind scheme (TOU), Fromm scheme [31], and second-order upwind scheme (SOU) to reduce high numerical diffusion [32]. These linear schemes are at least of 2nd-order accuracy and they are unbounded. However, Godunov [33] proved that linear unbounded high-order schemes are not mathematically monotonic as compared to 1st-order upwind scheme, allowing unphysical oscillations under some circumstances. Only bounded nonlinear high-order schemes can be monotone and stable. Nonlinear flux limiter schemes which fulfill total variation diminishing (TVD) criteria are one of the most effective methods to achieve high-order accuracy and stability and many researchers have investigated these nonlinear flux limiters. Tiselj and Petelin [34] developed the PDE2 code based on the six-equation model with flux limiter Minmod in achieving second-order accuracy. Wang et al. [29] have implemented high-resolution flux limiters into TRACE code to improve spatial accuracy of convection terms in mass/energy equations and the performance of six nonlinear flux limiters was assessed. Wang et al. [35, 36] also used Lax-Wendroff (L-W) scheme with flux limiters to achieve the 2nd-order accuracy for both spatial discretization and time integration and added the ENO limiter scheme into TRACE for BWR stability analysis. Abu Saleem et al. [37] developed a new flux limiter based on a combination of 1st-order upwind scheme and 3rd-order QUICK scheme to achieve high-resolution of the solver for the two-phase single-pressure model. Zou et al. [30] adapted an existing high-resolution spatial discretization scheme to reduce numerical errors. In his work, the high-resolution scheme was applied for the mass and momentum equations only, and energy equations were neglected. In all the work mentioned above, the high-order schemes are performed for two-fluid six-equation single-pressure model only; few studies have been done on implementing high-resolution scheme in solving two-pressure model.

In the present work, a semi-implicit algorithm using high-resolution TVD schemes is derived to calculate the two-pressure model. The accuracy and robustness of the proposed numerical scheme are validated with the water faucet benchmark test. Eight high-resolution flux limiter schemes, Minmod [38], Superbee [39], Harmonic [40], OSPRE [41], Van Albada [42], SMART [43], Koren [44], and MUSCL [45], are implemented to evaluate the performance of their accuracy. Consequently, this research will potentially lay the foundation for simulating two-phase flows based on the two-pressure model with higher fidelity algorithm. This paper is organized as follows. Section 2 introduces two-fluid two-pressure mathematical model. The numerical algorithm using high-resolution TVD scheme for solving this model is described in Section 3. Next, the validation test for the proposed scheme is presented in Section 4. At last, Section 5 is devoted to the conclusion.

#### 2. Two-Fluid Two-Pressure Mathematical Model

In nuclear industries, the one-dimensional form of two-fluid model is more economical and used widely. The two-fluid seven-equation two-pressure model used in this paper is summarized as follows [16, 22, 23, 46], which allows for a change of the pipeline cross section [34].The gas and liquid volume fractions and are related by the saturation constraint ; the pressure relaxation coefficient stands for the relaxation rate at which the phase pressures equilibrate. As derived in [8], the interfacial pressure and velocity are, respectively, determined byin which The pressure relaxation coefficient can be expressed in the following form (see [22, 23, 46]):where is the pressure relaxation time (s). Here we chose an empirical constant s and the numerical results in this paper show reasonable agreement with the exact/analytical solutions. The net interfacial mass transfer per unit volume can be calculated by [3, 37, 47]

#### 3. Numerical Scheme

In this paper, the finite volume method based on the staggered mesh is used to discretize seven conservation equations of the two-pressure model and the semi-implicit solution algorithm using high-resolution TVD schemes is developed to solve such two-pressure model. Figure 1 shows the staggered mesh schematic. Based on the staggered mesh, the scalar variables such as phasic pressures , specific internal energies , phasic temperatures , and phasic volume fractions are described at the volume centers and the vector quantities (velocities ) are defined at the volume edges.