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

High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow 1 first-author 2 3 corresp Corresponding author: sa.talebi@aut.ac.ir Corresponding author: sa.talebi@aut.ac.ir sa.talebi@aut.ac.ir Research School of Reactor, Nuclear Science and Technology Research Institute, P.O. Box 11365-3486, Tehran, Islamic Republic of Iran. Research School of Reactor, Nuclear Science and Technology Research Institute, P.O. Box 11365-3486, Tehran, Islamic Republic of Iran. Research School of Radiation Applications, Nuclear Science and Technology Research Institute, P.O. Box 11365-3486, Tehran, Islamic Republic of Iran. Research School of Radiation Applications, Nuclear Science and Technology Research Institute, P.O. Box 11365-3486, Tehran, Islamic Republic of Iran. Department of Energy Engineering and Physics, Amirkabir University of Technology (Tehran Polytechnic), 424 Hafez Avenue, P.O. Box 15875-4413, Tehran, Islamic Republic of Iran. Department of Energy Engineering and Physics, Amirkabir University of Technology (Tehran Polytechnic), 424 Hafez Avenue, P.O. Box 15875-4413, Tehran, Islamic Republic of Iran. 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 implicit Two-phase flow Drift-flux model Pressure-based algorithm 1 Introduction 1 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 1 ] and TRACE [ 2 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 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 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 4 ]. However, the time step is still limited to a value of 10 to 20 times greater than the material Courant limit [ 1 1 ]. With the development of computers, taking advantage of more advanced numerical methods has become possible [ 5 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 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 7 - 10 10 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 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 ]. 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 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 3 , 6 6 , 12 12 , 13 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 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 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 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 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 17 ]. Talebi et al [ 18 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 19 ]. Furthermore, the DFM is used for the analysis of subchannels in light water reactors. For example, Khan and Yi [ 20 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 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 1 2.1 One-Dimensional Drift-Flux Model 2 2.1.1 Governing Equations 3 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 21 ], consists of two continuity equations, namely the mixture momentum equation and the mixture energy equation, presented below [ 3 3 , 18 18 , 19 19 ]. The mixture continuity equation can be written as: <math display="block" id="M1"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> The continuity equation for dispersed phase is <math display="block" id="M2"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow></mfrac><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mi>Γ</mi><mo>,</mo></mrow></math> where <math id="M3"><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></math> , <math id="M4"><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow></math> , and <math id="M5"><mrow><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow></math> are mixture, vapor, and liquid-phase densities. The drift velocity ( <math id="M6"><mrow><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow></math> ) and the evaporation rate ( <math id="M7"><mi>Γ</mi></math> ) are obtained by appropriate constitutive equations. The mixture momentum equation can be written as <math display="block" id="M8"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub><msup><mrow/><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mfrac><mrow><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow></mfrac><msup><mrow><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mfrac><mrow><mo>∂</mo><mtext>P</mtext></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>−</mo><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>g</mi><mtext>z</mtext></msub><mo>−</mo><mfrac><mrow><msub><mi>f</mi><mtext>m</mtext></msub></mrow><mrow><mn>2</mn><mi>D</mi></mrow></mfrac><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub><mrow><mo>|</mo><mrow><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>|</mo></mrow><mo>,</mo></mrow></math> where D , <math id="M9"><mrow><msub><mi>g</mi><mtext>z</mtext></msub></mrow></math> , and <math id="M10"><mrow><msub><mi>f</mi><mtext>m</mtext></msub></mrow></math> denote the channel hydraulic diameter, gravity vector, and friction factor, respectively. The mixture enthalpy-energy equation can be stated as follows: <math display="block" id="M11"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>h</mi><mtext>m</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>h</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msubsup><mi>q</mi><mtext>W</mtext><mrow><mtext>""</mtext></mrow></msubsup><msub><mi>ξ</mi><mtext>h</mtext></msub></mrow><mi>A</mi></mfrac><mo>−</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>l</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></mfrac><mi>Δ</mi><msub><mi>h</mi><mrow><mtext>gl</mtext></mrow></msub><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><mo>∂</mo><mtext>P</mtext></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mrow><mo>(</mo><mrow><msub><mi>v</mi><mtext>m</mtext></msub><mo>+</mo><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>l</mtext></msub><mo>−</mo><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></mfrac><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><mtext>P</mtext></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M12"><mrow><msubsup><mi>q</mi><mtext>W</mtext><mrow><mtext>""</mtext></mrow></msubsup></mrow></math> is the wall heat flux, <math id="M13"><mrow><msub><mi>ξ</mi><mtext>h</mtext></msub></mrow></math> is the heated perimeter, and <math id="M14"><mi>A</mi></math> is the channel area. In Eqs. (1) to (4), pressure ( <math id="M15"><mi>P</mi></math> ), mixture velocity ( <math id="M16"><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></math> ), mixture enthalpy ( <math id="M17"><mrow><msub><mi>h</mi><mtext>m</mtext></msub></mrow></math> ), and void fraction ( <math id="M18"><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow></math> ) are the primary variables calculated by solving the discretized equations. 2.1.2 Constitutive Equations 3 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 22 ]. There are two options for the calculation of the phase temperatures ( <math id="M19"><mrow><msub><mi>T</mi><mi>g</mi></msub></mrow></math> and <math id="M20"><mrow><msub><mi>T</mi><mi>l</mi></msub></mrow></math> ): 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 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 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 25 ] to determine the wall evaporation rate and the model proposed by Andersen et al [ 26 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 26 ] are used. The Chexal–Lellouche correlation [ 27 27 ] is applied to obtain drift velocity. Finally, an appropriate model is necessary to calculate the quality. The equilibrium quality is calculated from <math display="block" id="M21"><mrow><msub><mi>x</mi><mrow><mtext>eq</mtext></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>h</mi><mtext>m</mtext></msub><mo>−</mo><msubsup><mi>h</mi><mtext>l</mtext><mtext>s</mtext></msubsup></mrow><mrow><msubsup><mi>h</mi><mtext>g</mtext><mtext>s</mtext></msubsup><mo>−</mo><msubsup><mi>h</mi><mtext>l</mtext><mtext>s</mtext></msubsup></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M22"><mrow><msubsup><mi>h</mi><mtext>g</mtext><mtext>s</mtext></msubsup><mo>,</mo><mtext> </mtext><msubsup><mi>h</mi><mtext>l</mtext><mtext>s</mtext></msubsup></mrow></math> are the vapor and liquid saturation enthalpy. The quality <math id="M23"><mrow><msub><mi>x</mi><mtext>g</mtext></msub></mrow></math> is calculated as a function of the equilibrium quality <math id="M24"><mrow><msub><mi>x</mi><mrow><mtext>eq</mtext></mrow></msub></mrow></math> according to the Levy [ 28 28 ] suggestion <math display="block" id="M25"><mrow><msub><mi>x</mi><mtext>g</mtext></msub><mo>=</mo><msub><mi>x</mi><mrow><mtext>eq</mtext></mrow></msub><mo>+</mo><msub><mi>x</mi><mrow><mtext>eq</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>Z</mi><mrow><mtext>OSV</mtext></mrow></msub></mrow><mo>)</mo></mrow><mtext>exp</mtext><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>x</mi><mrow><mtext>eq</mtext></mrow></msub></mrow><mrow><msub><mi>x</mi><mrow><mtext>eq</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>Z</mi><mrow><mtext>OSV</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> The onset of significant void (OSV) takes place where a significant increase in the void fraction occurs ( <math id="M26"><mrow><msub><mi>Z</mi><mrow><mtext>OSV</mtext></mrow></msub></mrow></math> ). The following equation is proposed by Delhaye et al. [ 29 29 ] to calculate <math id="M27"><mrow><msub><mi>Z</mi><mrow><mtext>OSV</mtext></mrow></msub></mrow></math> : <math display="block" id="M28"><mrow><msub><mi>Z</mi><mrow><mtext>OSV</mtext></mrow></msub><mo>=</mo><msub><mi>Z</mi><mrow><mtext>ONB</mtext></mrow></msub><mo>+</mo><mrow><mo>{</mo><mrow><msubsup><mi>h</mi><mtext>l</mtext><mtext> </mtext></msubsup><mrow><mo>[</mo><mrow><msub><mi>T</mi><mtext>l</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>Z</mi><mrow><mtext>OSV</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>−</mo><msubsup><mi>h</mi><mtext>l</mtext><mtext> </mtext></msubsup><mrow><mo>[</mo><mrow><msub><mi>T</mi><mtext>l</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>Z</mi><mrow><mtext>ONB</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow><mo>}</mo></mrow><mfrac><mrow><mi>G</mi><mi>D</mi></mrow><mrow><mn>4</mn><msubsup><mi>q</mi><mtext>W</mtext><mrow><mtext>""</mtext></mrow></msubsup></mrow></mfrac><mo>.</mo></mrow></math> The onset of nucleate boiling (ONB) is defined as the point where first bubbles appear on the wall nucleation sites ( <math id="M29"><mrow><msub><mi>Z</mi><mrow><mtext>ONB</mtext></mrow></msub></mrow></math> ). This is calculated by using the following relation proposed by Delhaye et al. [ 29 29 ]: <math display="block" id="M30"><mrow><msub><mi>Z</mi><mrow><mtext>ONB</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mi>G</mi><msub><mi>C</mi><mrow><mtext>pl</mtext></mrow></msub><mi>D</mi></mrow><mn>4</mn></mfrac><mrow><mo>[</mo><mrow><mfrac><mrow><msup><mi>T</mi><mtext>s</mtext></msup><mo>−</mo><msub><mi>T</mi><mrow><mtext>l,in</mtext></mrow></msub><mo>+</mo><mi>Δ</mi><msup><mi>T</mi><mtext>s</mtext></msup></mrow><mrow><msubsup><mi>q</mi><mtext>W</mtext><mrow><mtext>""</mtext></mrow></msubsup></mrow></mfrac><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>h</mi><mrow><mtext>l0</mtext></mrow></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> where <math id="M31"><mrow><msub><mi>h</mi><mrow><mtext>l</mtext><mn>0</mn></mrow></msub></mrow></math> is calculated from the Dittus–Boelter correlation and ΔT s is estimated by the Frost and Dzakowic [ 30 30 ] correlation, which can be stated as follows: <math display="block" id="M32"><mrow><mi>Δ</mi><msup><mi>T</mi><mtext>s</mtext></msup><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mn>8</mn><mi>Q</mi><mi>σ</mi><msup><mi>T</mi><mtext>s</mtext></msup></mrow><mrow><msub><mi>k</mi><mtext>l</mtext></msub><msub><mi>h</mi><mrow><mtext>lg</mtext></mrow></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>0.5</mn></mrow></msup><mi>P</mi><msub><mi>r</mi><mtext>l</mtext></msub><mo>.</mo></mrow></math> 2.2 Numerical Method 2 2.2.1 Discretization of the Governing Equation 3 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 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 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. The mixture continuity equation can be written as: <math display="block" id="M33"><mrow><mfrac><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mi>i</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mi>i</mi><mi>n</mi></msubsup></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mtext>Δ</mtext><mi>z</mi><mo>−</mo><msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mi>v</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mi>v</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mn>0.</mn></mrow></math> The continuity equation for the dispersed phase can be stated as <math display="block" id="M34"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mi>i</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mi>i</mi><mi>n</mi></msubsup></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mtext>Δ</mtext><mi>z</mi><mo>−</mo><msubsup><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>)</mo></mrow><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>)</mo></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></mfrac><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>]</mo></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></mfrac><mover accent="true"><mrow><msub><mi>V</mi><mrow><mtext>dj</mtext></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>]</mo></mrow><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>Γ</mi><mi>i</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mtext>Δ</mtext><mi>z</mi><mo>=</mo><mn>0.</mn></mtd></mtr></mtable></math> The mixture momentum equation can be expressed as <math display="block" id="M35"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><msubsup><mrow><mo stretchy="false">(</mo><msub><mi>ρ</mi><mi>m</mi></msub><msub><mi>ν</mi><mi>m</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mrow><mo stretchy="false">(</mo><msub><mi>ρ</mi><mi>m</mi></msub><msub><mi>ν</mi><mi>m</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mtext>Δ</mtext><mi>z</mi><mo>−</mo><msubsup><mrow><mo>[</mo><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>]</mo></mrow><mi>i</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mrow><mo>[</mo><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>]</mo></mrow><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msup><mrow><mo>[</mo><mrow><msub><mi>P</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>P</mi><mi>i</mi></msub></mrow><mo>]</mo></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msup></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>f</mi><mtext>m</mtext></msub></mrow><mrow><mn>2</mn><mi>D</mi></mrow></mfrac><msub><mi>ρ</mi><mtext>m</mtext></msub><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>v</mi><mtext>m</mtext></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mtext>Δ</mtext><mi>z</mi></mrow><mo>)</mo></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>+</mo><msub><mi>ρ</mi><mtext>m</mtext></msub><msubsup><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><msub><mi>g</mi><mtext>z</mtext></msub><mtext>Δ</mtext><mi>z</mi></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mfrac><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></mfrac><msubsup><mrow><mover accent="true"><mi>V</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mtext>dj</mtext></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>α</mi><mi mathvariant="normal">g</mi></msub></mrow><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mfrac><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>m</mtext></msub></mrow></mfrac><msubsup><mrow><mover accent="true"><mi>V</mi><mo stretchy="true">¯</mo></mover></mrow><mrow><mtext>dj</mtext></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mi>i</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mn>0.</mn></mtd></mtr></mtable></math> Finally, the mixture enthalpy-energy equation can be written as follows: http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330319&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330314&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330319&type= 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 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 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 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 ( <math id="M36"><mi>ψ</mi></math> ), the relation between the parameters on the staggered grids is obtained. The relations for different parameters are stated below. Pressure: <math display="block" id="M37"><mtable columnalign="left"><mtr><mtd><msubsup><mi>r</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>P</mi></msubsup><mo>=</mo><mfrac><mrow><msub><mi>P</mi><mi>i</mi></msub><mo>−</mo><msub><mi>P</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mrow><msub><mi>P</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>P</mi><mi>i</mi></msub></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>P</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><msub><mi>P</mi><mi>i</mi></msub><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>ψ</mi><mrow><mo>(</mo><mrow><msubsup><mi>r</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>P</mi></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>P</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>P</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> Enthalpy: <math display="block" id="M38"><mtable columnalign="left"><mtr><mtd><msubsup><mi>r</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><msub><mi>h</mi><mtext>m</mtext></msub></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mi>i</mi></msub><mo>−</mo><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mrow><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mi>i</mi></msub></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mi>i</mi></msub><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>ψ</mi><mrow><mo>(</mo><mrow><msubsup><mi>r</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><msub><mi>h</mi><mtext>m</mtext></msub></mrow></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>h</mi><mtext>m</mtext></msub><msub><mrow/><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> Void fraction: <math display="block" id="M39"><mtable columnalign="left"><mtr><mtd><msubsup><mi>r</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mi>i</mi></msub><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mi>i</mi></msub></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mi>i</mi></msub><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>ψ</mi><mrow><mo>(</mo><mrow><msubsup><mi>r</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub><msub><mrow/><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> Mass flux (G): <math display="block" id="M40"><mtable columnalign="left"><mtr><mtd><msubsup><mi>r</mi><mi>i</mi><mrow><msub><mi>G</mi><mtext>m</mtext></msub></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>3</mn><mn>2</mn></mfrac></mrow></msub></mrow><mrow><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mi>i</mi></msub><mo>=</mo><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>ψ</mi><mrow><mo>(</mo><mrow><msubsup><mi>r</mi><mi>i</mi><mrow><msub><mi>G</mi><mtext>m</mtext></msub></mrow></msubsup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>G</mi><mi mathvariant="normal">m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>G</mi><mtext>m</mtext></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> Ψ is a flux limiter function. The widely used form of this function is presented in Table 1 Table 1 [ 31 31 ]. 1 1 FVM scheme 1 1 Limiter function 1 1 Central differencing 1 1 1 1 1 Upwind differencing 1 1 0 1 1 Linear upwind differentiating 1 1 γ 1 1 Quick 1 1 <math id="M41"><mrow><mfrac><mrow><mn>3</mn><mo>+</mo><mi>r</mi></mrow><mn>4</mn></mfrac></mrow></math> 1 1 TVD Van Leer 1 1 <math id="M42"><mrow><mfrac><mrow><mi>r</mi><mo>+</mo><mrow><mo>|</mo><mi>r</mi><mo>|</mo></mrow></mrow><mrow><mn>1</mn><mo>+</mo><mi>r</mi></mrow></mfrac></mrow></math> 1 1 TVD Van Albada 1 1 <math id="M43"><mrow><mfrac><mrow><mi>r</mi><mo>+</mo><msup><mi>r</mi><mn>2</mn></msup></mrow><mrow><mn>1</mn><mo>+</mo><msup><mi>r</mi><mn>2</mn></msup></mrow></mfrac></mrow></math> 1 1 TVD Min-Mod 1 1 <math id="M44"><mrow><mi>max</mi><mrow><mo>(</mo><mrow><mn>0</mn><mo>,</mo><mi>m</mi><mi>i</mi><mi>n</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> 1 1 TVD SUPERBEE 1 1 <math id="M45"><mrow><mi>m</mi><mi>a</mi><mi>x</mi><mrow><mo>(</mo><mrow><mn>0</mn><mo>,</mo><mi>m</mi><mi>i</mi><mi>n</mi><mrow><mo>(</mo><mrow><mn>2</mn><mi>r</mi><mo>,</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>,</mo><mi>m</mi><mi>i</mi><mi>n</mi><mrow><mo>(</mo><mrow><mi>r</mi><mo>,</mo><mn>2</mn></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> 1 1 TVD UMIST 1 1 <math id="M46"><mrow><mi>m</mi><mi>a</mi><mi>x</mi><mrow><mo>(</mo><mrow><mn>0</mn><mo>,</mo><mi>m</mi><mi>i</mi><mi>n</mi><mrow><mo>(</mo><mrow><mn>2</mn><mi>r</mi><mo>,</mo><mfrac><mrow><mn>1</mn><mo>+</mo><mn>3</mn><mi>r</mi></mrow><mn>4</mn></mfrac><mo>,</mo><mfrac><mrow><mn>3</mn><mo>+</mo><mi>r</mi></mrow><mn>4</mn></mfrac><mo>,</mo><mn>2</mn></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> 1 1 TVD Quick 1 1 <math id="M47"><mrow><mi>m</mi><mi>a</mi><mi>x</mi><mrow><mo>(</mo><mrow><mn>0</mn><mo>,</mo><mi>m</mi><mi>i</mi><mi>n</mi><mrow><mo>(</mo><mrow><mn>2</mn><mi>r</mi><mo>,</mo><mfrac><mrow><mn>3</mn><mo>+</mo><mi>r</mi></mrow><mn>4</mn></mfrac><mo>,</mo><mn>2</mn></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> 2.2.2 Present model based on the SIMPLE Algorithm 3 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 <math display="block" id="M48"><mrow><msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mi>v</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mi>v</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mi>i</mi><mi>n</mi></msubsup><mo>−</mo><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mi>i</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mtext>Δ</mtext><mi>z</mi><mtext> </mtext><mo> </mo><mo>.</mo></mrow></math> Continuity equation for the dispersed phase http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330335&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330332&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330335&type= Mixture momentum equation http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330345&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330341&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330345&type= Mixture enthalpy-energy equation http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330355&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330351&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330355&type= 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: http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330369&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330365&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10330369&type= where <math id="M49"><mrow><msubsup><mi>v</mi><mi>m</mi><mtext>*</mtext></msubsup></mrow></math> and <math id="M50"><mrow><msubsup><mi>P</mi><mrow/><mtext>*</mtext></msubsup></mrow></math> are the guessed velocities and pressure, which should be corrected based on the following correction equations: <math display="block" id="M51"><mrow><msub><mi>v</mi><mi>m</mi></msub><mo>=</mo><msubsup><mi>v</mi><mi>m</mi><mo>*</mo></msubsup><mo>+</mo><mi>δ</mi><msub><mi>v</mi><mi>m</mi></msub><mo>,</mo></mrow></math> <math display="block" id="M52"><mrow><mi>P</mi><mo>=</mo><msubsup><mi>P</mi><mo> </mo><mo>*</mo></msubsup><mo>+</mo><mi>δ</mi><mi>P</mi><mo>,</mo></mrow></math> where <math id="M53"><mrow><mi>δ</mi><msub><mi>v</mi><mi>m</mi></msub></mrow></math> and <math id="M54"><mrow><mi>δ</mi><mi>P</mi></mrow></math> 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 31 ] Based on that method, the velocity correction is obtained as follows: <math display="block" id="M55"><mrow><msub><mi>v</mi><mi>m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><msubsup><mi>v</mi><mi>m</mi><mo>*</mo></msubsup><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><mfrac><mn>1</mn><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>k</mi></msubsup><mfrac><mrow><mtext>Δ</mtext><mi>z</mi></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mo>+</mo><msup><mrow><mrow><mo>[</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><msub><mi>v</mi><mi>m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mo>]</mo></mrow></mrow><mi>k</mi></msup></mrow><mo>)</mo></mrow></mrow></mfrac><mrow><mo>(</mo><mrow><mi>δ</mi><msub><mtext>P</mtext><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><mi>δ</mi><msub><mtext>P</mtext><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> 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 Fig. 1 . Introducing the equation for the corrected velocities into the discretized continuity equation, the following algebraic equation for the pressure correction is obtained: <math display="block" id="M56"><mrow><msup><mi>A</mi><mi>P</mi></msup><msub><mrow/><mrow><mi>i</mi><mo>,</mo><mi>i</mi></mrow></msub><mi>δ</mi><msub><mtext>P</mtext><mi>i</mi></msub><mo>+</mo><msup><mi>A</mi><mtext>P</mtext></msup><msub><mrow/><mrow><mi>i</mi><mo>,</mo><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mi>δ</mi><msub><mtext>P</mtext><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>+</mo><msup><mi>A</mi><mtext>P</mtext></msup><msub><mrow/><mrow><mi>i</mi><mo>,</mo><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mi>δ</mi><msub><mtext>P</mtext><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>=</mo><msup><mi>S</mi><mtext>P</mtext></msup><msub><mrow/><mi>i</mi></msub><mo>,</mo></mrow></math> where <math display="block" id="M57"><mtable columnalign="left"><mtr><mtd><msubsup><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>i</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>δ</mi><mi>P</mi></mrow></msubsup><mo>=</mo><mo>−</mo><mfrac><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>k</mi></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>k</mi></msubsup><mfrac><mrow><mtext>Δ</mtext><mi>z</mi></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mo>+</mo><msup><mrow><mrow><mo>[</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><msub><mi>v</mi><mi>m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mo>]</mo></mrow></mrow><mi>k</mi></msup></mrow><mo>)</mo></mrow></mrow></mfrac></mtd></mtr><mtr><mtd><msubsup><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mrow><mi>δ</mi><mi>P</mi></mrow></msubsup><mo>=</mo><mo>−</mo><mfrac><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>k</mi></msubsup></mrow><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>k</mi></msubsup><mfrac><mrow><mtext>Δ</mtext><mi>z</mi></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mo>+</mo><msup><mrow><mrow><mo>[</mo><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msub><mrow/><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><msub><mi>v</mi><mi>m</mi></msub><msub><mrow/><mi>i</mi></msub></mrow><mo>]</mo></mrow></mrow><mi>k</mi></msup></mrow><mo>)</mo></mrow></mrow></mfrac></mtd></mtr><mtr><mtd><msubsup><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>i</mi></mrow><mrow><mi>δ</mi><mi>P</mi></mrow></msubsup><mo>=</mo><mo>−</mo><msubsup><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>i</mi><mo>+</mo><mn>1</mn></mrow><mrow><mi>δ</mi><mi>P</mi></mrow></msubsup><mo>−</mo><msubsup><mi>A</mi><mrow><mi>i</mi><mo>,</mo><mi>i</mi><mo>−</mo><mn>1</mn></mrow><mrow><mi>δ</mi><mi>P</mi></mrow></msubsup></mtd></mtr><mtr><mtd><msup><mi>f</mi><mrow><mi>δ</mi><mi>P</mi></mrow></msup><mrow><mo>(</mo><mrow><mi>P</mi><mo>,</mo><msub><mi>v</mi><mi>m</mi></msub><mo>,</mo><msub><mi>α</mi><mi>g</mi></msub><mo>,</mo><msub><mi>h</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mi>i</mi><mi>n</mi></msubsup><mo>−</mo><msub><mi>ρ</mi><mi>m</mi></msub><msubsup><mrow/><mi>i</mi><mi>k</mi></msubsup></mrow><mrow><mtext>Δ</mtext><mi>t</mi></mrow></mfrac><mtext>Δ</mtext><mi>z</mi><mo>.</mo></mtd></mtr></mtable></math> 2.2.3 Special Techniques 3 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 <math display="block" id="M58"><mrow><msubsup><mi>h</mi><mi>l</mi><mtext> </mtext></msubsup><mo>=</mo><mfrac><mrow><msub><mi>h</mi><mi>m</mi></msub><mo>−</mo><msub><mi>x</mi><mi>g</mi></msub><msubsup><mi>h</mi><mi>g</mi><mtext> </mtext></msubsup></mrow><mrow><mn>1</mn><mo>−</mo><msub><mi>x</mi><mi>g</mi></msub></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M59"><mrow><msubsup><mi>h</mi><mi>l</mi><mtext> </mtext></msubsup></mrow></math> and <math id="M60"><mrow><msubsup><mi>h</mi><mi>g</mi><mtext> </mtext></msubsup></mrow></math> 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 . 3 Results and Discussions 1 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 32 , 33 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 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 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 34 , 35 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 Fig. 3 . In this study, eight finite volume cells are applied for all the cases. 1 1 Case 1 1 P (MPa) 1 1 Mass flux (kg/m 2 s) 1 1 Wall heat flux (MW/m 2 ) 1 1 T inlet (K) 1 1 Reference 1 1 1 1 1 14.70 1 1 2014 1 1 1.72 1 1 545 1 1 [ 30 30 ] 1 1 2 1 1 10.81 1 1 966 1 1 1.13 1 1 502 1 1 [ 30 30 ] 1 1 3 1 1 6.81 1 1 998 1 1 0.44 1 1 521 1 1 [ 30 30 ] 1 1 4 1 1 4.50 1 1 900 1 1 0.38 1 1 506.14 1 1 [ 29 29 ] 1 1 5 1 1 3.00 1 1 900 1 1 0.38 1 1 481.85 1 1 [ 29 29 ] 1 1 6 1 1 1.50 1 1 900 1 1 0.38 1 1 448.47 1 1 [ 29 29 ] 1 1 7 1 1 7.01 1 1 996 1 1 1.98 1 1 434 1 1 [ 30 30 ] 3.1 Effectiveness of Special Techniques 2 Throughout this section, the following relative iteration tolerance for the pressure velocity correction loop is used: <math display="block" id="M61"><mrow><mfrac><mrow><msub><mrow><mrow><mo>‖</mo><mrow><mi>d</mi><mi>P</mi></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow><mrow><msub><mrow><mrow><mo>‖</mo><mi>P</mi><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow></mfrac><mo>≤</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>7</mn></mrow></msup><mo>,</mo><mfrac><mrow><msub><mrow><mrow><mo>‖</mo><mrow><mi>d</mi><msub><mi>G</mi><mi>m</mi></msub></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow><mrow><msub><mrow><mrow><mo>‖</mo><mrow><msub><mi>G</mi><mi>m</mi></msub></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow></mfrac><mo>≤</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>4</mn></mrow></msup><mo>.</mo></mrow></math> For the outer loop, the relative iteration tolerance for the void fraction and enthalpy is considered: <math display="block" id="M62"><mrow><mfrac><mrow><msub><mrow><mrow><mo>‖</mo><mrow><mi>d</mi><msub><mi>h</mi><mi>m</mi></msub></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow><mrow><msub><mrow><mrow><mo>‖</mo><mrow><msub><mi>h</mi><mi>m</mi></msub></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow></mfrac><mo>≤</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>7</mn></mrow></msup><mo>,</mo><mfrac><mrow><msub><mrow><mrow><mo>‖</mo><mrow><mi>d</mi><msub><mi>α</mi><mi>g</mi></msub></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow><mrow><msub><mrow><mrow><mo>‖</mo><mrow><msub><mi>α</mi><mi>g</mi></msub></mrow><mo>‖</mo></mrow></mrow><mn>2</mn></msub></mrow></mfrac><mo>≤</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>2</mn></mrow></msup><mo>.</mo></mrow></math> 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 36 - 38 38 36-38 ]. Fig. 4 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. As Table 3 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. 1 1 FVM scheme 2 1 Number of iterations (#) 2 1 CPU (s) 1 1   1 1 SIMPLE(ES-T) 1 1 SIMPLE 1 1 SIMPLE(ES-T) 1 1 SIMPLE 1 1 Central differencing 1 1 58 1 1 61 1 1 1.30 1 1 1.86 1 1 Upwind differencing 1 1 14 1 1 35 1 1 0.29 1 1 1.34 1 1 Linear upwind differentiating 1 1 28 1 1 >200 1 1 0.61 1 1 - 1 1 Quick 1 1 23 1 1 50 1 1 0.50 1 1 1.79 1 1 TVD Van Leer 1 1 21 1 1 74 1 1 0.42 1 1 2.53 1 1 TVD Van Albada 1 1 21 1 1 119 1 1 0.41 1 1 4.43 1 1 TVD Min-Mod 1 1 25 1 1 >200 1 1 0.54 1 1 - 1 1 TVD SUPERBEE 1 1 31 1 1 >200 1 1 0.64 1 1 - 1 1 TVD UMIST 1 1 19 1 1 64 1 1 0.39 1 1 2.42 1 1 TVD Quick 1 1 23 1 1 >200 1 1 0.51 1 1 - 3.2 Steady State 2 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 Table 2 . Based on the report published by Bartolomei et al. [ 33 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 ( <math id="M63"><mrow><msub><mi>ε</mi><mi>α</mi></msub></mrow></math> ) are obtained by the following equation: <math display="block" id="M64"><mrow><msub><mi>ε</mi><mi>α</mi></msub><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mrow><mo>|</mo><mrow><msub><mi>α</mi><mrow><mi>n</mi><mi>u</mi><mi>m</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>−</mo><msub><mi>α</mi><mrow><mi>exp</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow><mo>|</mo></mrow></mrow></mstyle><mo>,</mo></mrow></math> where <math id="M65"><mrow><msub><mi>α</mi><mrow><mi>n</mi><mi>u</mi><mi>m</mi></mrow></msub></mrow></math> and <math id="M66"><mrow><msub><mi>α</mi><mrow><mi>exp</mi></mrow></msub></mrow></math> 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 3 , 33 33 ]. Because of the scattering of the experimental data, disparate values for MAE are obtained. However, Table 4 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. 1 1 FVM scheme 7 1 MAE 1 1   1 1 Case 1 1 1 Case 2 1 1 Case 3 1 1 Case 4 1 1 Case 5 1 1 Case 6 1 1 Average 1 1 Central differencing 1 1 0.015 1 1 0.054 1 1 0.016 1 1 0.019 1 1 0.031 1 1 0.064 1 1 0.033 1 1 Upwind differencing 1 1 0.010 1 1 0.028 1 1 0.015 1 1 0.012 1 1 0.018 1 1 0.061 1 1 0.024 1 1 Linear upwind differentiating 1 1 0.010 1 1 0.048 1 1 0.011 1 1 0.016 1 1 0.027 1 1 0.064 1 1 0.029 1 1 Quick 1 1 0.015 1 1 0.052 1 1 0.014 1 1 0.018 1 1 0.029 1 1 0.065 1 1 0.032 1 1 TVD Van Leer 1 1 0.013 1 1 0.049 1 1 0.012 1 1 0.017 1 1 0.028 1 1 0.065 1 1 0.031 1 1 TVD Van Albada 1 1 0.013 1 1 0.049 1 1 0.012 1 1 0.017 1 1 0.028 1 1 0.065 1 1 0.031 1 1 TVD Min-Mod 1 1 0.013 1 1 0.047 1 1 0.025 1 1 0.016 1 1 0.027 1 1 0.064 1 1 0.032 1 1 TVD SUPERBEE 1 1 0.015 1 1 0.053 1 1 0.013 1 1 0.019 1 1 0.030 1 1 0.065 1 1 0.033 1 1 TVD UMIST 1 1 0.013 1 1 0.049 1 1 0.012 1 1 0.016 1 1 0.027 1 1 0.065 1 1 0.030 1 1 TVD Quick 1 1 0.015 1 1 0.051 1 1 0.013 1 1 0.018 1 1 0.027 1 1 0.065 1 1 0.032 1 1 RELAP 1 1 0.025 1 1 0.030 1 1 0.027 1 1 0.024 1 1 0.026 1 1 0.077 1 1 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 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. Based on Fig. 6 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. 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 26 ], as shown below, are used. <math display="block" id="M67"><mrow><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mtext>  </mtext><mi>α</mi><mo>≤</mo><mtext>0</mtext><mtext>.3 </mtext><mtable><mtr><mtd><mrow><mtable><mtr><mtd><mrow><mtable><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr></mtable></mrow></mtd><mtd><mrow/></mtd></mtr></mtable></mrow></mtd><mtd><mrow/></mtd></mtr></mtable><mtext>bubbly flow</mtext></mtd></mtr><mtr><mtd><mtext>  0</mtext><mtext>.3<</mtext><mi>α</mi><mo>≤</mo><msub><mi>α</mi><mrow><mtext>tran</mtext></mrow></msub><mtext>-0</mtext><mtext>.1 </mtext><mtable><mtr><mtd><mrow/></mtd><mtd><mrow><mtext>churn flow</mtext></mrow></mtd></mtr></mtable></mtd></mtr><mtr><mtd><mtext> </mtext><msub><mi>α</mi><mrow><mtext>tran</mtext></mrow></msub><mtext>-0</mtext><mtext>.1<</mtext><mi>α</mi><mo>≤</mo><msub><mi>α</mi><mrow><mtext>tran</mtext></mrow></msub><mtext> </mtext><mtable><mtr><mtd><mrow/></mtd><mtd><mrow><mtext>transition region</mtext></mrow></mtd></mtr></mtable></mtd></mtr><mtr><mtd><mtext> </mtext><msub><mi>α</mi><mrow><mtext>tran</mtext></mrow></msub><mo><</mo><mi>α</mi><mo>≤</mo><mn>1.0</mn><mtext> </mtext><mtable><mtr><mtd><mrow><mtable><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr></mtable></mrow></mtd><mtd><mrow/></mtd></mtr></mtable><mtext>dispersed annular flow</mtext></mtd></mtr></mtable></mrow><mo>.</mo></mrow></math> <math display="block" id="M68"><mrow><msub><mi>α</mi><mrow><mtext>tran</mtext></mrow></msub><mtext>=(1+</mtext><mfrac><mrow><mtext>4</mtext><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow></mfrac><mtext>)(</mtext><mfrac><mtext>1</mtext><mrow><msub><mtext>C</mtext><mtext>0</mtext></msub></mrow></mfrac><mtext>)-</mtext><mfrac><mrow><mtext>4</mtext><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow></mfrac><mo>.</mo></mrow></math> <math display="block" id="M69"><mrow><msub><mtext>C</mtext><mtext>0</mtext></msub><msub><mrow><mtext>=C</mtext></mrow><mtext>B</mtext></msub><msub><mrow><mtext>(C</mtext></mrow><mtext>B</mtext></msub><mtext>-1)</mtext><msqrt><mrow><mfrac><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mrow><msub><mi>ρ</mi><mtext>l</mtext></msub></mrow></mfrac></mrow></msqrt><mtext> </mtext><mo>.</mo></mrow></math> <math display="block" id="M70"><mrow><msub><mtext>C</mtext><mtext>B</mtext></msub><mtext>=1</mtext><mtext>.393-0</mtext><mtext>.015log(Re)</mtext></mrow></math> Fig. 7 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. 3.3 Transient 2 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 2 for case 3 and 380 kW/m 2 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 2 s to 998 kg/m 2 s for case 3 and from 2500 kg/m 2 s to 900 kg/m 2 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 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 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. 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 2 from 400 kW/m 2 in 10 s. The exit void fraction versus the time step is shown in Fig. 9 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 Fig. 9 . Based on the two graphs, the following results can be achieved: · 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 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 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. 4 Conclusion 1 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 ). Relap5 / Mod3.3 Code Manual Volume I : Code Structure, System Models, and Solution Methods. I Division of System Analysis . TRACE V5.0 THEORY MANUAL: Field Equations, Solution Methods, and Physical Models . ( 2008 ). TRACE V5.0 THEORY MANUAL: Field Equations, Solution Methods, and Physical Models 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 http://doi.org/10.1016/j.anucene.2015.07.033 Numerical implementation, verification and validation of two-phase flow four-equation drift flux model with Jacobian-free Newton-Krylov method 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 http://doi.org/10.1016/0021-9991(86)90054-9 A nearly-implicit hydrodynamics numerical scheme for two-phase flows 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 http://doi.org/10.1016/j.jcp.2003.08.010 Jacobian-free Newton-Krylov methods: a survey of approaches and applications 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 http://doi.org/10.5516/NET.2009.41.7.885 . Accuracy and Efficiency of a Coupled Neutronics and Thermal Hydraulics Model 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 http://doi.org/10.1016/0029-5493(90)90294-8 The physical closure laws in the CATHARE code 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 http://doi.org/10.1016/0029-5493(93)90064-G Advanced numerical methods for thermal hydraulics 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 http://doi.org/10.1016/S0029-5493(03)00159-6 Notes on the implementation of a fully-implicit numerical scheme for a two-phase three-field flow model 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 http://doi.org/10.1155/2015/575380 Development of accurate and stable two-phase two-fluid model solver Y. Saad , Iterative Methods for Sparse Linear Systems Second Edition Yousef Saad . 535 ( 2003 ). doi: 10.2277/0898715342 http://doi.org/10.2277/0898715342 Iterative Methods for Sparse Linear Systems Second Edition Yousef Saad 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 http://doi.org/10.1016/j.jcp.2004.03.009 Implicitly balanced solution of the two-phase flow equations coupled to nonlinear heat conduction A. Ashrafizadeh , Investigation of a Jacobian-free Newton-Krylov solution to multiphase flows . ( 2014 ). doi: 10.1002/fld.3999 http://doi.org/10.1002/fld.3999 . Investigation of a Jacobian-free Newton-Krylov solution to multiphase flows 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 http://doi.org/10.1016/S0045-7930(02)00005-1 . Fourier analysis of the IPSA / PEA algorithms applied to multiphase flows with mass transfer S. Patankar , Numerical heat transfer and fluid flow . ( CRC press , 1980 ). Numerical heat transfer and fluid flow 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 http://doi.org/10.1115/1.3689137 Average Volumetric Concentration in Two-Phase Flow Systems 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 ). A drift-flux model of two-phase flow for RETRAN 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 http://doi.org/10.1016/j.nucengdes.2011.10.038 A numerical technique for analysis of transient two-phase flow in a vertical tube using the Drift Flux Model 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 http://doi.org/10.5516/NET.02.2013.025 TAPINS: A thermal-hydraulic system code for transient analysis of a fully-passive integral PWR H.J. Khan , H.C. Yi , Subchannel analysis in BWR fuel bundles . 12 , 559 - 572 ( 1985 ). doi: 10.1016/0306-4549(85)90029-5 http://doi.org/10.1016/0306-4549(85)90029-5 . Subchannel analysis in BWR fuel bundles 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 ). One-dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes 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 http://doi.org/10.1063/1.1461829 The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use 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 ). SOLA-LOOP: a nonequilibrium, drift-flux code for two-phase flow in networks 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 http://doi.org/10.1016/j.ijmultiphaseflow.2008.08.003 . Evaluation analysis of prediction methods for two-phase flow pressure drop in mini-channels R.T. Lahey , A mechanistic subcooled boiling model . In Proceedings of the 6th International Heat Transfer Conference 1 , 293 - 297 ( 1978 ). A mechanistic subcooled boiling model J.G.M. Andersen , R. Harrington , B. Hizoum , COBRAG Subchannel Code: Model Description Report NEDE-32199P , Revision 1 , 80 p ( 2007 ). COBRAG Subchannel Code: Model Description Report NEDE-32199P B. Chexal , The Chexal-Lellouche void fraction correlation for generalized applications . Electr. Power Res. Inst. Rep. EPRI NSAC/ 139 ( 1991 ). The Chexal-Lellouche void fraction correlation for generalized applications 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 http://doi.org/10.1016/0017-9310(67)90071-3 . Forced convection subcooled boiling—prediction of vapor volumetric fraction 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 http://doi.org/10.1016/j.ijheatmasstransfer.2004.05.018 Void fraction predictions in forced convective subcooled boiling of water between 10 and 18 Mpa 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 ). An extension of the method of predictive incipient boiling on commercially finished surfaces H.K. Versteegand , W. Malalasekera , An Introduction to Computational Fluid Dynamics . Longman Scientific & Technical , New York ( 1995 ). An Introduction to Computational Fluid Dynamics G. Bartolomei , V. Chanturiya , Experimental study of true void fraction when boiling subcooled water in vertical tubes . Therm. Eng . 14 , 123 - 128 ( 1967 ). Experimental study of true void fraction when boiling subcooled water in vertical tubes 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 ). An experimental investigation of true volumetric vapor content with subcooled boiling in tubes A.S. Devkin , A.S. Podosenov , International Agreement Report RELAP5 / MOD3 Subcooled Boiling Model Assessment. NUREG/IA-0 , ( 1998 ). International Agreement Report RELAP5 / MOD3 Subcooled Boiling Model Assessment. NUREG/IA-0 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 ). Refinement of the RELAP5/MOD3. 2 subcooled boiling model for low-pressure conditions 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 http://doi.org/10.1002/fld.3999 . A Jacobian-free Newton-Krylov method for thermal hydraulics simulations 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 http://doi.org/10.1016/S0029-5493(03)00159-6 Notes on the implementation of a fully implicit numerical scheme for a two-phase three-field flow model 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 http://doi.org/10.1115/1.1865223 A fully implicit hybrid solution method for a two-phase thermal-hydraulic model No acknowledgments. 10.1007/s41365-018-0456-3 2017-05-05 2017-08-13 2017-08-16 2018-07-04 2018-08-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2018 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-018-0456-3 http://dx.doi.org/10.1007/s41365-018-0456-3 . 2018 Nuclear Science and Techniques 08 29