High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow

NUCLEAR ENERGY SCIENCE AND ENGINEERING

High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow

A. Hajizadeh，

H. Kazeminejad，

S. Talebi

Nuclear Science and Techniques

Vol.29, No.8

Article number 114

Published in print 01 Aug 2018

Available online 04 Jul 2018

DOI：10.1007/s41365-018-0456-3

52501

The drift-flux model has a practical importance in two-phase flow analysis. In this study, a finite volume solution is developed for a transient four-equation drift-flux model through the staggered mesh, leading to the development of a fully implicit discretization method. The main advantage of the fully implicit method is its unconditional stability. Newton's scheme is a popular method of choice for the solution of a nonlinear system of equations arising from fully implicit discretization of field equations. However, the lack of convergence robustness and the construction of Jacobian matrix have created several difficulties for the researchers. In this paper, a fully implicit model is developed based on the SIMPLE algorithm for two-phase flow simulations. The drawbacks of Newton's method are avoided in the developed model. Different limiter functions are considered, and the stabilized method is developed under steady and transient conditions. The results obtained by the numerical modeling are in good agreement with the experimental data. As expected, the results prove that the developed model is not restricted by any stability limit.

Fully implicitTwo-phase flowDrift-flux modelPressure-based algorithm

1 Introduction

The two phase modeling of fluid flow in a vertical channel has a great importance in nuclear reactor operational and accident analysis, especially in boiling water reactors. From 1970 onwards, large number of codes, including RELAP [1] and TRACE [2], have been developed for thermo-hydraulic analysis of nuclear reactors. These codes are developed based on the semi-implicit methods. The first order accuracy in time and the stability restriction on time step are two main drawbacks of the semi-implicit approaches [3]. Some nearly implicit methods are developed based on the two-step approach. In the first step, mass and energy equations are evaluated, and then, in the second step, momentum equations are solved using a fully implicit approach [4]. In this method, the implicitness of the solution is evaluated and the time step limitation corresponding to the interphase exchange processes is ignored [4]. However, the time step is still limited to a value of 10 to 20 times greater than the material Courant limit [1].

With the development of computers, taking advantage of more advanced numerical methods has become possible [5]. The recent efforts for implementation of fully implicit methods have demonstrated the success and robustness of this type of procedure. Some advantages of the fully implicit approaches in two phase flow problems, especially during the transients, are as below:

1. An accident or transient can range from very fast to very slow; therefore, development of a numerical model with no restrictions on the time step can be beneficial. The fully implicit approach is a method of choice for this purpose.

2. The temporal accuracy of the fully implicit methods is 2^nd order [6].

3. Because of the existence of interphase transport phenomena (such as mass, momentum, and energy transported in the interface) in two-phase flow conditions, the system of governing equations is highly non-linear and are also coupled with each other. Hence, a fully implicit numerical method is useful because it solves all discretized equations simultaneously.

Using the fully implicit methods for discretization of governing equations, a system of nonlinear equations is obtained. Developing a numerical solution for a nonlinear system of equations has always been a challenge. Research carried out in this regard has relied on the methods such as Newton’s method for linearization of a nonlinear system of equations. Several efforts have been made to apply Newton’s method for simulation of two-phase flows [7-10]. However, the implementation of Newton’s method can present difficulties, some of which are as follows:

1. The problem with convergence is frequently encountered in the Newton’s method [5]. To overcome this problem, a good initial guess and implementation of an appropriate globalization technique are required. This procedure increases the complexity of the problem, especially in the two-phase flow problems where finding a good initial guess requires some primary simulations.

2. The construction of the Jacobian matrix has been identified as a problem in the implementation of Newton’s method. The Jacobian matrix can be obtained by numerical estimation of the derivatives or analytical derivation of the equations. The numerical estimation of the derivatives is extremely time-consuming and requires high memory usage to store the Jacobian matrix. The analytical derivation of the field equations is also cumbersome, especially in the two-phase problems with complex constitutive equations [3].

3. Even with the availability of the Jacobian matrix, solving the linear equations obtained is challenging since the Jacobian matrix arising from the physical model field equations is a sparse matrix [11].

Recently, the Jacobian-free Newton–Krylov (JFNK) method has been developed to avoid the problems related to the construction of the Jacobian matrix. Few attempts have been made to use the JFNK methods for the solution of two-phase flow problems [3,6,12,13]. However, the implementation of this method has several difficulties, including the strong dependence of the success of the JFNK method on the development of an efficient preconditioner. The preconditioning technique depends on the intended application, and is expensive and difficult to build in some cases [5]. Furthermore, the JFNK method applies a matrix-vector product that can introduce some errors to the results.

The aim of the present study is to develop a high order fully implicit model for simulation of two-phase flow problem in a vertical channel. For that purpose, a fully implicit numerical method, which avoids the problems associated with Newton’s method, is developed based on the iterative nature of the SIMPLE method. Therefore, the two-phase flow simulation without the stability limit requirements was carried out. Several TVD schemes, namely the Van Leer, Van Albada, SUPERBEE, UMIST, and Quick schemes, in addition to the central differencing, upwind differencing, linear upwind differentiating, and Quick schemes were implemented on staggered grids. The SIMPLE algorithm is reliable and robust for solving the single-phase flow problems. However, because of the strong coupling between the momentum and the continuity equations and the effect of the interphase mass transfer, convergence problem arises in the two-phase flow [14]. One solution to this difficulty is to write the interfacial terms in the equations, in an explicit manner, to obtain a pseudo-solution that resolves the inter-phase coupling [15]. In this study, some special techniques are used to overcome this difficulty.

There are two choices for two-phase flow modeling including the drift-flux model (DFM) and the two-fluid model. The DFM, originally developed by Zuber and Findlay [16], is an appropriate model for simulation of two-phase systems. This model is widely used in nuclear research; for example, the RETRAN code uses the DFM [17]. Talebi et al [18] studied the transient two-phase flow in a vertical channel using the DFM, while the TAPINS code solved the DFM field equations using the Gauss–Seidel method [19]. Furthermore, the DFM is used for the analysis of subchannels in light water reactors. For example, Khan and Yi [20] used the DFM for the development of the boiling water reactor subchannel analysis code. Since the results of advanced computer code are used in the licensing process [19], the exact solution of the applied model is very important. Therefore, in the present work, high-order finite volume methods developed by the implementation of the high order limiter functions are stabilized using the special techniques, and finally, these developed methods are evaluated.

The DFM field equations and the corresponding constitutive equations are described in Section 2. The numerical model developed in this study is presented in Section 3 and the results are discussed in Section 4. Finally, the conclusion is presented in Section 5.

2 Experimental analysis

2.1 One-Dimensional Drift-Flux Model

2.1.1 Governing Equations

The DFM considers the mixture as a whole instead of two separate phases. One-dimensional transient four-equation DFM, developed by Hibiki and Ishii [21], consists of two continuity equations, namely the mixture momentum equation and the mixture energy equation, presented below [3,18,19].

The mixture continuity equation can be written as:

\frac{\partial ρ_{m}}{\partial t} + \frac{\partial}{\partial z} (ρ_{m} v_{m}) = 0,

(1)

The continuity equation for dispersed phase is

\frac{\partial α_{g} ρ_{g}}{\partial t} + \frac{\partial}{\partial z} (α_{g} ρ_{g} v_{m}) + \frac{\partial}{\partial z} (\frac{α_{g} ρ_{l}}{ρ_{m} ρ_{g}} \bar{V_{dj}}) = Γ,

(2)

where $ρ_{m}$ , $ρ_{g}$ , and $ρ_{l}$ are mixture, vapor, and liquid-phase densities. The drift velocity ( $\bar{V_{dj}}$ ) and the evaporation rate ( $Γ$ ) are obtained by appropriate constitutive equations.

The mixture momentum equation can be written as

\frac{\partial ρ_{m} v_{m}}{\partial t} + \frac{\partial}{\partial z} (ρ_{m} v_{m}^{2}) + \frac{\partial}{\partial z} (\frac{α_{g}}{(1 - α_{g})} \frac{ρ_{l}}{ρ_{m} ρ_{g}} {\bar{V_{dj}}}^{2}) = - \frac{\partial P}{\partial z} - ρ_{m} g_{z} - \frac{f_{m}}{2 D} ρ_{m} v_{m} | v_{m} |,

(3)

where D, $g_{z}$ , and $f_{m}$ denote the channel hydraulic diameter, gravity vector, and friction factor, respectively.

The mixture enthalpy-energy equation can be stated as follows:

\frac{\partial ρ_{m} h_{m}}{\partial t} + \frac{\partial}{\partial z} (ρ_{m} h_{m} v_{m}) = \frac{q_{W}^{""} ξ_{h}}{A} - \frac{\partial}{\partial z} (\frac{α_{g} ρ_{l} ρ_{g}}{ρ_{m}} Δ h_{gl} \bar{V_{dj}}) + \frac{\partial P}{\partial t} + (v_{m} + \frac{α_{g} (ρ_{l} - ρ_{g})}{ρ_{m}} \bar{V_{dj}}) \frac{\partial P}{\partial z},

(4)

where $q_{W}^{""}$ is the wall heat flux, $ξ_{h}$ is the heated perimeter, and $A$ is the channel area. In Eqs. (1) to (4), pressure ( $P$ ), mixture velocity ( $ρ_{m}$ ), mixture enthalpy ( $h_{m}$ ), and void fraction ( $α_{g}$ ) are the primary variables calculated by solving the discretized equations.

2.1.2 Constitutive Equations

By surface averaging, the information about the variables perpendicular to the direction of the flow in a channel is essentially lost. Therefore, in DFM, it is necessary to employ appropriate constitutive equations to predict the state variables, drift velocity, two-phase friction, and evaporation rate.

To calculate the required state variables such as enthalpy and density of the liquid phase and vapor, the relations reported by the International Association for the Properties of Water and Steam (IAPWS) were implemented [22]. There are two options for the calculation of the phase temperatures ( $T_{g}$ and $T_{l}$ ): in the first option, the phase temperatures are considered to be equal, while in the second option, saturation is considered in the gas phase. However, Hirt and colleagues [23] found no noticeable difference between these two. Hence, in the present study, gas-phase was considered to be at saturation temperature. To predict the frictional pressure drop in the two-phase flows, the model developed by Sun and Mishima [24] is considered the best option based on extensive validation conducted by applying the empirical data. In relation to the prediction of the evaporation rate, calculation of the wall evaporation rate and the mass transfer rate are required. The Lahey model is used [25] to determine the wall evaporation rate and the model proposed by Andersen et al [26] is used to calculate the interfacial mass transfer rate. These relationships are the function of the flow regimes; therefore, the flow regime transition criteria provided by Andersen et al [26] are used. The Chexal–Lellouche correlation [27] is applied to obtain drift velocity. Finally, an appropriate model is necessary to calculate the quality. The equilibrium quality is calculated from

x_{eq} = \frac{h_{m} - h_{l}^{s}}{h_{g}^{s} - h_{l}^{s}},

(5)

where $h_{g}^{s}, h_{l}^{s}$ are the vapor and liquid saturation enthalpy. The quality $x_{g}$ is calculated as a function of the equilibrium quality $x_{eq}$ according to the Levy [28] suggestion

x_{g} = x_{eq} + x_{eq} (Z_{OSV}) exp (\frac{x_{eq}}{x_{eq} (Z_{OSV})} - 1),

(6)

The onset of significant void (OSV) takes place where a significant increase in the void fraction occurs ( $Z_{OSV}$ ). The following equation is proposed by Delhaye et al. [29] to calculate $Z_{OSV}$ :

Z_{OSV} = Z_{ONB} + {h_{l}^{} [T_{l} (Z_{OSV})] - h_{l}^{} [T_{l} (Z_{ONB})]} \frac{G D}{4 q_{W}^{""}} .

(7)

The onset of nucleate boiling (ONB) is defined as the point where first bubbles appear on the wall nucleation sites ( $Z_{ONB}$ ). This is calculated by using the following relation proposed by Delhaye et al. [29]:

Z_{ONB} = \frac{G C_{pl} D}{4} [\frac{T^{s} - T_{l,in} + Δ T^{s}}{q_{W}^{""}} - \frac{1}{h_{l0}}] .

(8)

where $h_{l 0}$ is calculated from the Dittus–Boelter correlation and ΔT^s is estimated by the Frost and Dzakowic [30] correlation, which can be stated as follows:

Δ T^{s} = {(\frac{8 Q σ T^{s}}{k_{l} h_{lg} ρ_{g}})}^{0.5} P r_{l} .

(9)

2.2 Numerical Method

2.2.1 Discretization of the Governing Equation

The division of the physical space into several control volumes is the first step in the finite volume method. For this purpose, some nodes are considered in the flow domain and the control volume boundaries are located between the adjacent nodes [31]. Thus, each node is surrounded by a control volume. Then, discretization is performed on the basis of the staggered grid method to analyze the 'checker-board’ pressure field. The staggered grid is shown in Fig. 1. In this nodalization, the pressure, enthalpy, and void fraction are obtained on the nodes inside the control volume and the velocity is calculated in the interfacial nodes. By integrating Eqs. (1)–(4) on the control volumes and by using the fully implicit Euler method for a time, a set of non-linear equations are obtained. For the description of the axial nodes and time step, "i" and "n" are used as the subscript and the superscript in the discretized equations.

Fig. 1

(Color online) Control volume and staggered grid

The mixture continuity equation can be written as:

\frac{ρ_{m}_{i}^{n + 1} - ρ_{m}_{i}^{n}}{Δ t} Δ z - {(ρ_{m} v_{m})}_{i - \frac{1}{2}}^{n + 1} + {(ρ_{m} v_{m})}_{i + \frac{1}{2}}^{n + 1} = 0.

(10)

The continuity equation for the dispersed phase can be stated as

\begin{array}{l} \frac{{(α_{g} ρ_{g})}_{i}^{n + 1} - {(α_{g} ρ_{g})}_{i}^{n}}{Δ t} Δ z - {(ρ_{g} α_{g} v_{m})}_{i - \frac{1}{2}}^{n + 1} + {(ρ_{g} α_{g} v_{m})}_{i + \frac{1}{2}}^{n + 1} \\ + {[\frac{α_{g} ρ_{g} ρ_{l}}{ρ_{m}} \bar{V_{dj}}]}_{i + \frac{1}{2}}^{n + 1} - {[\frac{α_{g} ρ_{g} ρ_{l}}{ρ_{m}} \bar{V_{dj}}]}_{i - \frac{1}{2}}^{n + 1} - Γ_{i}^{n + 1} Δ z = 0. \end{array}

(11)

The mixture momentum equation can be expressed as

\begin{array}{l} \frac{{(ρ_{m} ν_{m})}_{i + \frac{1}{2}}^{n + 1} - {(ρ_{m} ν_{m})}_{i + \frac{1}{2}}^{n}}{Δ t} Δ z - {[ρ_{m} v_{m} v_{m}]}_{i}^{n + 1} + {[ρ_{m} v_{m} v_{m}]}_{i + 1}^{n + 1} + {[P_{i + 1} - P_{i}]}^{n + 1} \\ + {(\frac{f_{m}}{2 D} ρ_{m} {(v_{m})}^{2} Δ z)}_{i + \frac{1}{2}}^{n + 1} + ρ_{m}_{i + \frac{1}{2}}^{n + 1} g_{z} Δ z \\ + {(\frac{α_{g}}{(1 - α_{g})} \frac{ρ_{g} ρ_{l}}{ρ_{m}} {\bar{V}}_{dj}^{2})}_{i + 1}^{n + 1} - {(\frac{α_{g}}{(1 - α_{g})} \frac{ρ_{g} ρ_{l}}{ρ_{m}} {\bar{V}}_{dj}^{2})}_{i}^{n + 1} = 0. \end{array}

(12)

Finally, the mixture enthalpy-energy equation can be written as follows:

(13)

In Eqs. (5)–(8), pressure, enthalpy, and void fraction appear in the staggered grid; however, they should be calculated at the grid points. Moreover, velocity appears in the grid points, which should be obtained at the staggered grid. In Versteeg and Malalasekera's published book [31], several differencing schemes are discussed in detail. As explained in that book, the central differencing method uses the Taylor series with second-order truncation error. It is also mentioned that the results obtained by a numerical method will be physically realistic when the numerical method has three essential properties of conservativeness, boundedness, and transportiveness. The transportiveness property is not present in the central differencing approach. This method also cannot identify the flow direction; therefore, this is an unstable scheme. To ensure the stability and transportiveness, the upwind differencing method is developed by taking into account of the flow direction. Since the upwind differencing scheme is inherently a type of backward differencing method, its corresponding truncation error is of the order of one. The higher-order discretization approaches that consider the flow direction are required to increase the accuracy and to reduce the numerical errors, while remaining stable at the same time. From one point of view, it can be said that the other higher-order methods are formulated as summation of the upwind differencing and an additional function ψ to improve the order of accuracy. Versteeg and Malalasekera [31] explained that this function and the nature of the method were determined by the ratio of the upwind-side gradient to the downwind-side gradient, with the ratio being called "r". To achieve second-order accuracy, it is required that the developed function crosses the point (1,1) in the r–ψ diagram [31]. The central differencing, Quick scheme, linear upwind differencing scheme, and all TVD schemes possess this requirement and therefore, they all are second-order accurate. Because they are based on the upwind differencing, they possess the conservativeness, boundedness, and transportiveness requirement and consider the flow direction, thereby ensuring stability. Using the limiter functions ( $ψ$ ), the relation between the parameters on the staggered grids is obtained. The relations for different parameters are stated below.

Pressure:

\begin{array}{l} r_{i + \frac{1}{2}}^{P} = \frac{P_{i} - P_{i - 1}}{P_{i + 1} - P_{i}} \\ P_{i + \frac{1}{2}} = P_{i} + \frac{1}{2} ψ (r_{i + \frac{1}{2}}^{P}) (P_{i + 1} - P_{i}) . \end{array}

(14)

Enthalpy:

\begin{array}{l} r_{i + \frac{1}{2}}^{h_{m}} = \frac{h_{m}_{i} - h_{m}_{i - 1}}{h_{m}_{i + 1} - h_{m}_{i}} \\ h_{m}_{i + \frac{1}{2}} = h_{m}_{i} + \frac{1}{2} ψ (r_{i + \frac{1}{2}}^{h_{m}}) (h_{m}_{i + 1} - h_{m}_{i}) . \end{array}

(15)

Void fraction:

\begin{array}{l} r_{i + \frac{1}{2}}^{α_{g}} = \frac{α_{g}_{i} - α_{g}_{i - 1}}{α_{g}_{i + 1} - α_{g}_{i}} \\ α_{g}_{i + \frac{1}{2}} = α_{g}_{i} + \frac{1}{2} ψ (r_{i + \frac{1}{2}}^{α_{g}}) (α_{g}_{i + 1} - α_{g}_{i}) . \end{array}

(16)

Mass flux (G):

\begin{array}{l} r_{i}^{G_{m}} = \frac{G_{m}_{i - \frac{1}{2}} - G_{m}_{i - \frac{3}{2}}}{G_{m}_{i + \frac{1}{2}} - G_{m}_{i - \frac{1}{2}}} \\ G_{m}_{i} = G_{m}_{i - \frac{1}{2}} + \frac{1}{2} ψ (r_{i}^{G_{m}}) (G_{m}_{i + \frac{1}{2}} - G_{m}_{i - \frac{1}{2}}) . \end{array}

(17)

Ψ is a flux limiter function. The widely used form of this function is presented in Table 1 [31].

Most frequently used limiter functions

FVM scheme	Limiter function
Central differencing	1
Upwind differencing	0
Linear upwind differentiating	γ
Quick	$\frac{3 + r}{4}$
TVD Van Leer	$\frac{r + \| r \|}{1 + r}$
TVD Van Albada	$\frac{r + r^{2}}{1 + r^{2}}$
TVD Min-Mod	$\max (0, m i n (r, 1))$
TVD SUPERBEE	$m a x (0, m i n (2 r, 1), m i n (r, 2))$
TVD UMIST	$m a x (0, m i n (2 r, \frac{1 + 3 r}{4}, \frac{3 + r}{4}, 2))$
TVD Quick	$m a x (0, m i n (2 r, \frac{3 + r}{4}, 2))$

2.2.2 Present model based on the SIMPLE Algorithm

Four primary variables, namely, pressure, mixture velocity (or mass flux), enthalpy, and void fraction should be determined using the field equations. Other parameters are evaluated using the constitutive equations. The SIMPLE scheme is an iterative method; therefore, each iteration commences with a series of guess values for primary variables. In the first iteration of the first time step, the guess parameters are considered to be equal to the initial and boundary conditions; while in the case of other time steps, the guess parameters in the first iteration are equal to the value of primary variables obtained in the previous time step. For other iterations in each time step, the guess parameters are same as the parameters obtained in the previous iteration. The mixture velocity equation and the mass conservation equation will be used to ensure the coupling between the velocity and the pressure. However, both enthalpy and void fraction should be corrected iteratively by the corresponding equations. By replacing n+1 with k+1 for the new iteration and k for the previous one as described below, a system of linear equations for four primary variables at k+1 is obtained.

Mixture continuity equation

{(ρ_{m} v_{m})}_{i + \frac{1}{2}}^{k + 1} - {(ρ_{m} v_{m})}_{i - \frac{1}{2}}^{k + 1} = \frac{ρ_{m}_{i}^{n} - ρ_{m}_{i}^{k + 1}}{Δ t} Δ z .

(18)

Continuity equation for the dispersed phase

(19)

Mixture momentum equation

(20)

Mixture enthalpy-energy equation

(21)

The velocity and pressure fields are closely interlinked, and thus, must be solved simultaneously. The SIMPLE algorithm is an iterative approach in which pressure and velocity are coupled indirectly. The procedure starts with guessed vales of the pressures and velocity fields. The implicit iterative form of the momentum equation for the velocity at the staggered grid is in the following form:

(22)

where $v_{m}^{*}$ and $P_{}^{*}$ are the guessed velocities and pressure, which should be corrected based on the following correction equations:

v_{m} = v_{m}^{*} + δ v_{m},

(23)

P = P_{}^{*} + δ P,

(24)

where $δ v_{m}$ and $δ P$ are the correction terms for the velocity and the pressure, respectively. In the SIMPLE approach, the correct velocity field is linked to the correct pressure field to give the correcting equation for the velocity. The detailed discussion of the method is provided by Versteeg and Malalasekera [31] Based on that method, the velocity correction is obtained as follows:

v_{m}_{i + \frac{1}{2}} = v_{m}^{*}_{i + \frac{1}{2}} - \frac{1}{(ρ_{m}_{i + \frac{1}{2}}^{k} \frac{Δ z}{Δ t} + {[ρ_{m}_{i + \frac{1}{2}} v_{m}_{i + 1}]}^{k})} (δ P_{i + 1} - δ P_{i}),

(25)

Except for the momentum equation, the velocity field must satisfy the continuity equation as well. The continuity equation is discretized on the scalar control volume shown in Fig. 1. Introducing the equation for the corrected velocities into the discretized continuity equation, the following algebraic equation for the pressure correction is obtained:

A^{P}_{i, i} δ P_{i} + A^{P}_{i, i + 1} δ P_{i + 1} + A^{P}_{i, i - 1} δ P_{i - 1} = S^{P}_{i},

(26)

where

\begin{array}{l} A_{i, i + 1}^{δ P} = - \frac{ρ_{m}_{i + \frac{1}{2}}^{k}}{(ρ_{m}_{i + \frac{1}{2}}^{k} \frac{Δ z}{Δ t} + {[ρ_{m}_{i + \frac{1}{2}} v_{m}_{i + 1}]}^{k})} \\ A_{i, i - 1}^{δ P} = - \frac{ρ_{m}_{i - \frac{1}{2}}^{k}}{(ρ_{m}_{i - \frac{1}{2}}^{k} \frac{Δ z}{Δ t} + {[ρ_{m}_{i - \frac{1}{2}} v_{m}_{i}]}^{k})} \\ A_{i, i}^{δ P} = - A_{i, i + 1}^{δ P} - A_{i, i - 1}^{δ P} \\ f^{δ P} (P, v_{m}, α_{g}, h_{m}) = \frac{ρ_{m}_{i}^{n} - ρ_{m}_{i}^{k}}{Δ t} Δ z . \end{array}

(27)

2.2.3 Special Techniques

The SIMPLE method is an iterative method solving the mixture continuity equation and the mixture momentum equation by using the guess and correct procedure. In this approach, the other equations are solved one by one. By paying attention to the importance of void fraction and its dominant influence on the results, the dispersed continuity equation is solved in the present study first by employing the pressure and velocity values obtained by the guess and correct procedure and estimated enthalpy, which is obtained in the previous iteration. All these parameters have some errors that influence the calculated void fraction in this stage. Since pressure and velocity are corrected during the converged procedure, the error of guess enthalpy creates a visible error in the void fraction. Then, the obtained void fraction with its error will be introduced into the mixture enthalpy equation to determine the enthalpy. The error of guess enthalpy causes a visible error in the void fraction, which, in turn, amplifies the error of obtained enthalpy in consecutive iterations. This phenomenon causes the oscillatory divergence of the model, with oscillations appearing around the desired value. To diminish the oscillations in the present model, two special techniques are required.

1. At first, the dependence of enthalpy and the void fraction is reduced by applying the following equation for liquid enthalpy

h_{l}^{} = \frac{h_{m} - x_{g} h_{g}^{}}{1 - x_{g}},

(28)

where $h_{l}^{}$ and $h_{g}^{}$ are the liquid and vapor phase enthalpy, respectively and x is the quality. The use of Eq. (28) is crucial for stable convergence of the algorithm.

2. Then, based on the oscillatory divergence of the model around an averaged value, void fraction is considered to be the average of the obtained void fraction in the present and the previous iterations.

The stable solution will also allow the use of the high order limiter function for improved accuracy. The procedure of implementation of the SIMPLE algorithm is shown in Fig. 2.

Fig. 2

SIMPLE algorithm with special technique

3 Results and Discussions

There are three general methods for the validation of a numerical model: 1. Compare the results with the experimental data; 2. Compare the results with the results of validated robust software model; 3. Compare the results with the analytical solution. Since the analytical solution of the complex two-phase flow problems is cumbersome and often impossible to derive, the first two methods are widely used. Bartolomei et al. [32,33] have conducted many experimental works, with the derived results used for validation of two-phase numerical models in many research studies. The several cases used in this study are illustrated in Table 2. Their measurements are reported only for the steady state. Since the transient behavior of the results is of great importance for the evaluation of implemented fully implicit method in this study, the RELAP code is used for the transient study. The RELAP code [34] is one of the strong and popular software for thermo-hydraulic analysis of nuclear reactors and is widely used in various design stages of light water reactors. Extensive research in conjunction with the validation of the software has been reported in the literature. The research studies confirm its capabilities in the simulation of boiling water flows [34,35]. For modeling of the two-phase flows, RELAP5/MOD2 employs a full non-equilibrium, six-equation, two-fluid model. For the intended calculations, the software uses the semi-implicit method. The model component for modeling by RELAP5 is presented in Fig. 3. In this study, eight finite volume cells are applied for all the cases.

Experimental conditions for simulation

Case	P (MPa)	Mass flux (kg/m² s)	Wall heat flux (MW/m²)	T_inlet(K)	Reference
1	14.70	2014	1.72	545	[30]
2	10.81	966	1.13	502	[30]
3	6.81	998	0.44	521	[30]
4	4.50	900	0.38	506.14	[29]
5	3.00	900	0.38	481.85	[29]
6	1.50	900	0.38	448.47	[29]
7	7.01	996	1.98	434	[30]

Fig. 3

RELAP nodalization (TDV: Time-Dependent Volume, SJ: Single Junction TDJ: Time-Dependent Junction)

3.1 Effectiveness of Special Techniques

Throughout this section, the following relative iteration tolerance for the pressure velocity correction loop is used:

\frac{{‖ d P ‖}_{2}}{{‖ P ‖}_{2}} \leq 10^{- 7}, \frac{{‖ d G_{m} ‖}_{2}}{{‖ G_{m} ‖}_{2}} \leq 10^{- 4} .

(29)

For the outer loop, the relative iteration tolerance for the void fraction and enthalpy is considered:

\frac{{‖ d h_{m} ‖}_{2}}{{‖ h_{m} ‖}_{2}} \leq 10^{- 7}, \frac{{‖ d α_{g} ‖}_{2}}{{‖ α_{g} ‖}_{2}} \leq 10^{- 2} .

(30)

Furthermore, the SIMPLE scheme is terminated if no convergence is achieved after 200 iterations. Based on the literature, these criteria are appropriate for two-phase flow simulations [36-38].

Fig. 4 represents the convergence behavior of the solution with and without the aforementioned technique (Section 3.2), for Case 3. Based on the figures, the results show an oscillatory convergence around the final results. By implementation of the special techniques, the results are flattened, and final results can be achieved with less number of iterations.

Fig. 4

Convergence behavior of the void fraction for case 3 (ES-T: ESpecial Technique)

As Table 3 demonstrates, the rate of convergence is different for different flux limiter functions. Without the implementation of the special technique, the TVD Min-Mod, linear upwind differentiating, and TVD SUPERBEE schemes did not converge even after 200 iterations. By using the special technique, all the above schemes are stabilized, and the convergence is achieved with less number of iteration and lower CPU time. The upwind differencing scheme is the fastest scheme with the least number of iterations among all methods, while UMIST is the fastest among the TVD methods.

Number of iterations and CPU(s) (case 3)

FVM scheme	Number of iterations (#)		CPU (s)
	SIMPLE(ES-T)	SIMPLE	SIMPLE(ES-T)	SIMPLE
Central differencing	58	61	1.30	1.86
Upwind differencing	14	35	0.29	1.34
Linear upwind differentiating	28	>200	0.61	-
Quick	23	50	0.50	1.79
TVD Van Leer	21	74	0.42	2.53
TVD Van Albada	21	119	0.41	4.43
TVD Min-Mod	25	>200	0.54	-
TVD SUPERBEE	31	>200	0.64	-
TVD UMIST	19	64	0.39	2.42
TVD Quick	23	>200	0.51	-

3.2 Steady State

To verify the capability of the present model, simulation of the two-phase flow within a uniformly heated vertical tube at different conditions with a wide range of pressure and heat flux were performed. The analysis is carried out for the parameters presented in Table 2. Based on the report published by Bartolomei et al. [33], the maximum absolute errors of the measured void fraction are between ±0.04. For the comparison study, the mean absolute errors (MAE) of the calculated results ( $ε_{α}$ ) are obtained by the following equation:

ε_{α} = \frac{1}{N} \sum_{i = 1}^{N} | α_{n u m, i} - α_{\exp, i} |,

(31)

where $α_{n u m}$ and $α_{\exp}$ are the void fractions calculated by the numerical method and experimental measurement, respectively. Bartolomei et al. explained that the scattering of the experimental data could be due to the inaccuracy of maintaining operational parameters and their corresponding measurement errors, and Zou et al added that this could be a result of the stochastic nature of the two-phase flow phenomenon [3,33]. Because of the scattering of the experimental data, disparate values for MAE are obtained. However, Table 4 demonstrates that the maximum averaged error of the results in comparison with the measured data does not exceed the inherent absolute error of the experimental data. In the simulations with eight axial nodes, upwind differencing represents the minimum averaged MAE (i.e. 0.024) between all methods and the UMIST between the TVD schemes.

Mean Absolute error

FVM scheme	MAE
	Case 1	Case 2	Case 3	Case 4	Case 5	Case 6	Average
Central differencing	0.015	0.054	0.016	0.019	0.031	0.064	0.033
Upwind differencing	0.010	0.028	0.015	0.012	0.018	0.061	0.024
Linear upwind differentiating	0.010	0.048	0.011	0.016	0.027	0.064	0.029
Quick	0.015	0.052	0.014	0.018	0.029	0.065	0.032
TVD Van Leer	0.013	0.049	0.012	0.017	0.028	0.065	0.031
TVD Van Albada	0.013	0.049	0.012	0.017	0.028	0.065	0.031
TVD Min-Mod	0.013	0.047	0.025	0.016	0.027	0.064	0.032
TVD SUPERBEE	0.015	0.053	0.013	0.019	0.030	0.065	0.033
TVD UMIST	0.013	0.049	0.012	0.016	0.027	0.065	0.030
TVD Quick	0.015	0.051	0.013	0.018	0.027	0.065	0.032
RELAP	0.025	0.030	0.027	0.024	0.026	0.077	0.035

The order of the numerical schemes together with the constitutive equations (especially the interface heat and mass transfer correlations) determines the accuracy of the results. All the flux limiter functions labeled with TVD, represent a second-order accurate TVD discretization schemes. Therefore, for the rest of this study, TVD UMIST and upwind differencing are selected based on the CPU time and the MAE. To evaluate the accuracy of the numerical schemes, it is required to obtain the analytical solution results, which are cumbersome in the complicated two-phase problems. Because it is expected that the numerical results converge to the analytical results by increasing the number of nodes, an analysis is performed for increasing number of the nodes in Case 3. The void fraction obtained for 160 nodes is considered as the reference value. Then, a relative error corresponding to each number node is calculated as the relative difference between obtained the void fraction and the reference void fraction. As seen in Fig. 5, the TVD UMIST scheme is less dependent on the number of nodes. Furthermore, it is visible that the upwind differencing simulation results approach the TVD UMIST results with the increase in the number of nodes. In all number of nodes, the UMIST as a TVD scheme results in a less relative error, which proves its higher degree of accuracy.

Fig. 5

Void fraction relative error for case 3

Based on Fig. 6, for different values of inlet subcooling, the estimated axial void fraction distribution along the pipe and the results obtained by the RELAP5 code are in reasonable agreement with the measured data. The estimated initial point of vapor generation predicted by the present model is found acceptable for different inlet subcooling, compared to the RELAP5 code and the experimental data. The small deviations are essentially due to the applied models of heat transfer.

Fig. 6

Comparison of the present results with the experimental data and RELAP5 code

Case 7 is provided to evaluate the stability of model in a wide range of flow regimes. The flow regime maps presented by Andersen et al. [26], as shown below, are used.

{\begin{cases} α \leq 0 .3 \begin{matrix} \begin{matrix} \begin{matrix}  \end{matrix} \end{matrix} \end{matrix} bubbly flow \\ 0 .3< α \leq α_{tran} -0 .1 \begin{matrix} churn flow \end{matrix} \\ α_{tran} -0 .1< α \leq α_{tran} \begin{matrix} transition region \end{matrix} \\ α_{tran} < α \leq 1.0 \begin{matrix} \begin{matrix}  \end{matrix} \end{matrix} dispersed annular flow \end{cases} .

(32)

α_{tran} =(1+ \frac{4 ρ_{g}}{ρ_{l}})(\frac{1}{C_{0}})- \frac{4 ρ_{g}}{ρ_{l}} .

(33)

C_{0} {=C}_{B} {(C}_{B} -1) \sqrt{\frac{ρ_{g}}{ρ_{l}}} .

(34)

C_{B} =1 .393-0 .015log(Re)

(35)

Fig. 7 demonstrates the stability of the present model in a wide range of flow regimes. It should be mentioned that based on the present study, all discretization methods are stable except the central differencing scheme, which shows oscillatory instability when a transition between flow regimes occurs.

Fig. 7

Results for different flow regime (Case 7)

3.3 Transient

The accurate prediction of the transient behavior of void fraction through a two-phase system is very important for safety analysis. To analyze the accuracy of the results during a typical transient, it was assumed that the heat flux increases to 440 kW/m² for case 3 and 380 kW/m²for cases 4 and 5 in the duration of 1000 s. In another test case, it was assumed that the mass flux decreases from 2500 kg/m² s to 998 kg/m² s for case 3 and from 2500 kg/m² s to 900 kg/m² s for cases 4 and 5 in the duration of 1000 s. Among the four primary variables in the two-phase flow analysis, the void fraction is the most important parameter with dominant effect on other parameters and consequences of the two-phase transient. The value of the void fraction determines important feedbacks in nuclear reactors and is of great importance in the prediction of critical heat flux in the pressurized water reactors and critical power ratio in boiling water reactors. Fig. 8 shows the void fraction variation at the pipe outlet during the heat flux and the mass flux transient in a vertical boiling channel. As seen in the figure, the results obtained by present model are in a good agreement with the RELAP5 results for the transient simulations in boiling channel, with both following the same trend. By increasing the heat flux, the liquid temperature increases, and first bubbles grow on the nucleation sites of the heated wall at an elevation corresponding to the onset of the nucleate boiling. The present model and the RELAP5 code predict very similar values for the onset of boiling. A consistent trend is observed in the case of decreasing mass flux, where reducing mass flux leads to the increase in temperature and commencement of nucleation boiling. The number of bubbles increases until the bubbly flow is developed in the flow path where randomly distributed small bubbles travel at different velocities. It can be seen from Fig. 8 that when the transient continues, the void fraction increases beyond the value of 0.3, after which the churn flow is expected. The churn flow is developed when the number of bubble increases, the bubbles coalesce, and similar sized stable bubbles are formed.

Fig. 8

Comparison of void fraction versus time with the RELAP5 code

Since the present model solves the drift flux equations by a fully implicit method, it is expected that the time step is not restricted by the Courant–Friedrich–Levy limit. To investigate this, it is assumed that in case 3, the power increases to 440 kW/m² from 400 kW/m² in 10 s. The exit void fraction versus the time step is shown in Fig. 9 to compare the two methods—the present model as a fully implicit method, and the semi-implicit method implemented in the RELAP code. The Courant–Friedrich–Levy limit is also presented in Fig. 9. Based on the two graphs, the following results can be achieved:

Fig. 9

Exit void fraction versus time step, (a) Present model, (b) RELAP

· By changing the time step, it can be seen that the uncertainty of the void fraction obtained by the present model (fully implicit method) is in the order of 1.0e-5, while this amount for the semi-implicit method used in the RELAP code is in the order of 1.0e-3. It means the results are less dependent on the chosen time step in the present model.

· In the present model, the time step is not restricted by the Courant–Friedrich–Levy limit, while the time step in the semi-implicit method used in the RELAP code is limited to maintain the stability. The benefits of a fully implicit method in comparison with the semi-implicit method are that the time step can be selected large enough according to the type of the transient.

· The obtained trend for the void fraction at the last node and the last time step in the present model is a definite trend, making it is possible to easily fit a second-degree function on it, while the trend is not uniform in the case of the semi-implicit method implemented in the RELAP5 code.

To evaluate the stability of the present model when a small time step is required, a typical fast power transient is hypothesized to occur in case 7. The power variation is shown in Fig. 10a. A small time step equal to 0.06 s is considered for this simulation. The obtained void fraction is shown in Fig. 10b. An expected peak on the calculated void fraction as a result of the peak in the power is predicted well by the present model. This peak appears in the transient regime, which is obtained based on a weighted summation of the void fraction calculated by the constitutive equations corresponding to the annular flow and churn flow. It can be concluded that the present model is able to simulate the transients covering a wide range of flow regimes using a small time step.

Fig. 10

(a) Power, (b) void fraction

4 Conclusion

This paper presents a new robust and accurate model for simulation of a two-phase flow in a vertical channel using the SIMPLE algorithm. A number of TVD schemes, namely the Van Leer, Van Albada, SUPERBEE, UMIST, and Quick schemes, in addition to the central differencing, upwind differencing, Linear upwind differentiating, and Quick schemes were implemented on staggered grids. The DFM field equations were discretized by a fully implicit method and a system of nonlinear equations was obtained. Using the SIMPLE algorithm, nonlinear equations were linearized and solved. As the construction of the Jacobian matrix is not required for linearization in the present model, the associated difficulties are not raised here. Furthermore, the convergence of the results, contrary to what is expected in the case of Newton's method, does not depends on the good initial guess. According to the results, the advantage of the present model over the semi-implicit method is that the time step is not restricted by the CFL limit. Therefore, considerably large time step can be selected according to the transient time scale. The obtained results are in good agreement with the experimental data for a wide range of conditions. Among the different finite volume methods, the upwind differencing is the most robust and fast method with minimum MAE. The TVD schemes are more accurate, and among them, the UMIST method can be considered as the optimal scheme based on the convergence rate, CPU time, and MAE.

References

Nuclear Safety Analysis Division. Relap5 / Mod3.3 Code Manual Volume I : Code Structure, System Models, and Solution Methods. I, (2001).

Division of System Analysis. TRACE V5.0 THEORY MANUAL: Field Equations, Solution Methods, and Physical Models. (2008).

L. Zou, H. Zhao, H. Zhang,

Numerical implementation, verification and validation of two-phase flow four-equation drift flux model with Jacobian-free Newton-Krylov method

. Ann. Nucl. Energy 87, 707-719 (2015). doi: 10.1016/j.anucene.2015.07.033

Baidu Scholar

Google Scholar

J.A. Trapp, R.A. Riemke,

A nearly-implicit hydrodynamics numerical scheme for two-phase flows

. J. Comput. Phys. 66 62-82, (1986). doi: 10.1016/0021-9991(86)90054-9

Baidu Scholar

Google Scholar

D.A. Knoll, D.E. Keyes,

Jacobian-free Newton-Krylov methods: a survey of approaches and applications

. J. Comput. Phys. 193, 357-397 (2004). doi: 10.1016/j.jcp.2003.08.010

Baidu Scholar

Google Scholar

M.A. Pope, V.A. Mousseau,

Accuracy and Efficiency of a Coupled Neutronics and Thermal Hydraulics Model

. 41, 885-892 (2008). doi: 10.5516/NET.2009.41.7.885.

Baidu Scholar

Google Scholar

D. Bestion,

The physical closure laws in the CATHARE code

. Nucl. Eng. Des. 124, 229-245 (1990). doi: 10.1016/0029-5493(90)90294-8

Baidu Scholar

Google Scholar

F. Barre, M. Parent, B. Brun,

Advanced numerical methods for thermal hydraulics

. Nucl. Eng. Des. 145, 147-158 (1993). doi: 10.1016/0029-5493(93)90064-G

Baidu Scholar

Google Scholar

C. Frepoli, J.H. Mahaffy, K. Ohkawa,

Notes on the implementation of a fully-implicit numerical scheme for a two-phase three-field flow model

. Nucl. Eng. Des. 225, 191-217 (2003). doi: 10.1016/S0029-5493(03)00159-6

Baidu Scholar

Google Scholar

R.A.A. Saleem, T. Kozlowski,

Development of accurate and stable two-phase two-fluid model solver

. Proc. ICAPP 2014 6-9 (2014). doi: 10.1155/2015/575380

Baidu Scholar

Google Scholar

Y. Saad,

Iterative Methods for Sparse Linear Systems Second Edition Yousef Saad

. 535 (2003). doi: 10.2277/0898715342

Baidu Scholar

Google Scholar

V.A. Mousseau,

Implicitly balanced solution of the two-phase flow equations coupled to nonlinear heat conduction

. J. Comput. Phys. 200, 104-132 (2004). doi: 10.1016/j.jcp.2004.03.009

Baidu Scholar

Google Scholar

A. Ashrafizadeh, Investigation of a Jacobian-free Newton-Krylov solution to multiphase flows. (2014). doi: 10.1002/fld.3999.

T.F. Miller, D.J. Miller,

Fourier analysis of the IPSA / PEA algorithms applied to multiphase flows with mass transfer

. Comput. Fluids 32, 197-221 (2003). doi: 10.1016/S0045-7930(02)00005-1.

Baidu Scholar

Google Scholar

S. Patankar, Numerical heat transfer and fluid flow. (CRC press, 1980).

N. Zuber, J.A. Findlay,

Average Volumetric Concentration in Two-Phase Flow Systems

. J. Heat Transfer 87, 453-468 (1965). doi: 10.1115/1.3689137

Baidu Scholar

Google Scholar

E.D. Hughes, M.P. Paulsen, L. J. Agee,

A drift-flux model of two-phase flow for RETRAN

. Nucl. Technol. 54, 410-421 (1981).

Baidu Scholar

Google Scholar

S. Talebi, H. Kazeminejad, H. Davilu,

A numerical technique for analysis of transient two-phase flow in a vertical tube using the Drift Flux Model

. Nucl. Eng. Des. 242, 316-322 (2012). doi: 10.1016/j.nucengdes.2011.10.038

Baidu Scholar

Google Scholar

Y.G. Lee, G.C. Park,

TAPINS: A thermal-hydraulic system code for transient analysis of a fully-passive integral PWR

. Nucl. Eng. Technol. 45, 439-458 (2013). doi: 10.5516/NET.02.2013.025

Baidu Scholar

Google Scholar

H.J. Khan, H.C. Yi, Subchannel analysis in BWR fuel bundles. 12, 559-572 (1985). doi: 10.1016/0306-4549(85)90029-5.

T. Hibiki, M. Ishii,

One-dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes

. Int. J. Heat Mass Transf. 46, 4935-4948 (2003).

Baidu Scholar

Google Scholar

W. Wagner,

The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use

. J. Phys. Chem. Ref. Data 31, 387 (1999). doi: 10.1063/1.1461829

Baidu Scholar

Google Scholar

C.W. Hirt, T. A. Oliphant, W. C. Rivard et al., SOLA-LOOP: a nonequilibrium, drift-flux code for two-phase flow in networks. (Los Alamos Scientific Lab., 1979).

L. Sun, K. Mishima,

Evaluation analysis of prediction methods for two-phase flow pressure drop in mini-channels

. Int. J. Multiph. Flow 35, 47-54 (2009). doi: 10.1016/j.ijmultiphaseflow.2008.08.003.

Baidu Scholar

Google Scholar

R.T. Lahey,

A mechanistic subcooled boiling model

. In Proceedings of the 6th International Heat Transfer Conference 1, 293-297 (1978).

Baidu Scholar

Google Scholar

J.G.M. Andersen, R. Harrington, B. Hizoum,

COBRAG Subchannel Code: Model Description Report NEDE-32199P

, Revision 1, 80 p (2007).

Baidu Scholar

Google Scholar

B. Chexal,

The Chexal-Lellouche void fraction correlation for generalized applications

. Electr. Power Res. Inst. Rep. EPRI NSAC/139 (1991).

Baidu Scholar

Google Scholar

S. Levy,

Forced convection subcooled boiling—prediction of vapor volumetric fraction

. Int. J. Heat Mass Transf. 10, 951-965 (1967). doi: 10.1016/0017-9310(67)90071-3.

Baidu Scholar

Google Scholar

J.M. Delhaye, F. Maugin, J.M. Ochterbeck,

Void fraction predictions in forced convective subcooled boiling of water between 10 and 18 Mpa

. Int. J. Heat Mass Transf. 47, 4415-4425 (2004). doi: 10.1016/j.ijheatmasstransfer.2004.05.018

Baidu Scholar

Google Scholar

W. Frost, G.S. Dzakowic,

An extension of the method of predictive incipient boiling on commercially finished surfaces

. ASME-AIChE Heat Transfer Conference, Seattle (1967).

Baidu Scholar

Google Scholar

H.K. Versteegand, W. Malalasekera, An Introduction to Computational Fluid Dynamics. Longman Scientific & Technical, New York (1995).

G. Bartolomei, V. Chanturiya,

Experimental study of true void fraction when boiling subcooled water in vertical tubes

. Therm. Eng. 14, 123-128 (1967).

Baidu Scholar

Google Scholar

G.G. Bartolomei, V.G. Brantov, Y.S. Molochnikov et al.,

An experimental investigation of true volumetric vapor content with subcooled boiling in tubes

. Therm. Eng. 29(3), 132-135 (1982).

Baidu Scholar

Google Scholar

A.S. Devkin, A.S. Podosenov,

International Agreement Report RELAP5 / MOD3 Subcooled Boiling Model Assessment. NUREG/IA-0

, (1998).

Baidu Scholar

Google Scholar

S. Hari, Y. A. Hassan,

Refinement of the RELAP5/MOD3. 2 subcooled boiling model for low-pressure conditions

. Proc. 8th Int. Conf. in Nuclear Eng., Baltimore, USA (2000).

Baidu Scholar

Google Scholar

A. Ashrafizadeh, C.B. Devaud, N.U. Aydemir,

A Jacobian-free Newton-Krylov method for thermal hydraulics simulations

. Int. J. Numer. Meth. Fluids 77(10), 590-615 (2015). doi: 10.1002/fld.3999.

Baidu Scholar

Google Scholar

C. Frepoli, J.H. Mahaffy, K. Ohkawa,

Notes on the implementation of a fully implicit numerical scheme for a two-phase three-field flow model

. Nucl. Eng. Des. 225, 191-217 (2003). doi: 10.1016/S0029-5493(03)00159-6

Baidu Scholar

Google Scholar

V.A. Mousseau,

A fully implicit hybrid solution method for a two-phase thermal-hydraulic model

. J. Heat Transfer 127, 531-539 (2005). doi: 10.1115/1.1865223

Baidu Scholar

Google Scholar