- Email: [email protected]

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

A solver for the two-phase two-ﬂuid model based on high-resolution total variation diminishing scheme Rabie A. Abu Saleem a,∗ , Tomasz Kozlowski b , Rijan Shrestha c a

Nuclear Engineering Department, Jordan University of Science and Technology, P.O. Box 3030, Irbid 22110, Jordan Department of Nuclear, Plasma and Radiological Engineering, University of Illinois at Urbana-Champaign, 216 Talbot Laboratory, 104 S. Wright St., Urbana, IL 61801, USA c Portland Technology Development, Intel Corporation, 2501 NW 229th Ave Hillsboro OR 97124, USA b

h i g h l i g h t s • • • •

The two-ﬂuid model and the challenges associated with its numerical modeling are investigated. A high-order solver based on ﬂux limiter schemes and the theta method was developed. The solver was compared to existing thermal hydraulics codes used in nuclear industry. The solver was shown to handle fast transients with discontinuities and phase change.

a r t i c l e

i n f o

Article history: Received 29 September 2015 Received in revised form 6 March 2016 Accepted 11 March 2016

a b s t r a c t Finite volume techniques with staggered mesh are used to develop a new numerical solver for the one-dimensional two-phase two-ﬂuid model using a high-resolution, Total Variation Diminishing (TVD) scheme. The solver is implemented to analyze numerical benchmark problems for veriﬁcation and testing its abilities to handle discontinuities and fast transients with phase change. Convergence rates are investigated by comparing numerical results to analytical solutions available in literature for the case of the faucet ﬂow problem. The solver based on a new TVD scheme is shown to exhibit higher-order of accuracy compared to other numerical schemes. Mass errors are also examined when phase change occurs for the shock tube problem, and compared to those of the 1st-order upwind scheme implemented in the nuclear thermal-hydraulics code TRACE. The solver is shown to exhibit numerical stability when applied to problems with discontinuous solutions and results of the new solver are free of spurious oscillations. © 2016 Elsevier B.V. All rights reserved.

1. Introduction As a subject of intense interest in many engineering systems, the study of two-phase ﬂow is of great importance in applied research because it is involved in many industrial applications. In the nuclear industry the phenomenon of two-phase ﬂow plays a crucial role, because water in its liquid and gaseous phases is used as both a coolant and a moderator in many types of reactor cores. It also appears in other mechanical parts of the nuclear reactor, such as heat exchangers, condensers and turbines. Consequently, achieving the optimal design for both operation and safety requires a solid understanding of the fundamental aspects

∗ Corresponding author. E-mail addresses: [email protected] (R.A. Abu Saleem), [email protected] (T. Kozlowski), [email protected] (R. Shrestha). http://dx.doi.org/10.1016/j.nucengdes.2016.03.015 0029-5493/© 2016 Elsevier B.V. All rights reserved.

of two-phase ﬂow physics and its mathematical models. Modeling of two-phase ﬂows entails several difﬁculties; difﬁculties in the mathematical model used to govern the evolution of different properties in space and time, difﬁculties in physical models describing interfacial interactions, difﬁculties in the numerical methods used to solve the model and difﬁculties in building robust solvers for the algebraic equations arising from strongly non-linear discrete equations. In the last few decades, many researchers dedicated signiﬁcant effort to overcome the problems associated with modeling of twophase ﬂow. Because most of the legacy nuclear thermal-hydraulics codes were based on one-dimensional models, the majority of the work targeting the nuclear industry was in that area (Mousseau, 2004; Fullmer et al., 2013; Sokolowski and Koszela, 2012). Other work was performed in terms of multidimensional analysis (2D, 3D) (Balcazar et al., 2014; Bonometti and Magnaudet, 2007). Multidimensional analysis of the two-phase ﬂow allows solving for

256

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263

different ﬂow regimes. This entails interface capturing between the two phases using different techniques, such as the Level-Set (LS) method, the Front-Tracking (FT) method and the Volume-Of-Fluid (VOF) method. Due to its simplicity and ease of application, the VOF is the most popular technique used in CFD commercial codes like ANSYS FLUENT, STAR-CD and OpenFOAM (Ansys Inc., 2006; Lee et al., 2012). The difﬁculty of physical modeling of the two-phase ﬂow arises from the existence of moving and deforming interfaces between the two phases. Fluid properties are discontinuous at these interfaces and ﬂow ﬁelds are complicated, there have been several mathematical models proposed for the two phase ﬂow, one of these models is the two-ﬂuid model used in most nuclear simulators. Numerical difﬁculties depend on many factors, one of which is the desired order of accuracy. In the nuclear industry, most thermal-hydraulics codes like RELAP5 (DSR, 2010) and TRACE (DSA, 2010) employ 1st-order numerical schemes for spatial and temporal discretization. 1st-order schemes are known to possess an inherent numerical viscosity which bestows the dissipation property to the numerical solution. Such dissipation is essential to handle the solution near discontinuities. However, it causes smearing of the solution that can mask physical discontinuities and interfaces (Sweby, 1984). For higher order linear schemes, which are non-monotone, the numerical solution is always expected to pose spurious oscillations near discontinuities. This fact comes as a result of the Godunov’s order barrier theorem (Godunov, 1957) which states that a linear numerical scheme for solving the advection equation, having the property of not generating new extrema (monotone scheme), can be at most 1st-order accurate. This paper is aimed at investigating the one-dimensional twoﬂuid model and the challenges associated with its numerical modeling. For the numerical simulation of the problem, non-linear high-order schemes are investigated as a potential replacement for the 1st-order schemes used in nuclear thermal-hydraulic codes. Analysis for numerical properties is carried out, and a numerical solver based on this analysis is developed for the six-equation two-phase two-ﬂuid model. We present the discrete equations of the model along with its closure relations. Temporal discretization is based on the theta method with different schemes implemented for spatial discretization including the new scheme developed in (Abu Saleem and Kozlowski, 2015). The solver is then tested and veriﬁed using numerical benchmarks and analytical solutions.

2. Mathematical model The mathematical model for the two-phase ﬂow is obtained by applying conservation relations to mass, momentum and energy of different phases separately. The model employed by RELAP5 and TRACE is based on the standard one-dimensional two-ﬂuid model (DSR, 2010; DSA, 2010). This model uses a single pressure for both phases. The following two equations represent the conservation of mass for gas and liquid phases, respectively:

• g is the interfacial mass transfer per unit volume. It is equal in magnitude and opposite in sign between the two phases. g is given by (Mousseau, 2004): g = −

Hig aint (Ts − Tg ) + Hil aint (Ts − Tl ) h∗g − h∗l

In this equation: • • • •

aint is the interfacial area between the two phases. Tg and Tl are gas and liquid temperatures, respectively. Ts is saturation temperature. Hig and Hil are the heat transfer coefﬁcients of gas and liquid phases with the interface.

The phasic mass transfer enthalpies of the two phases (h∗g and h∗l ) are calculated as follows:

h∗g

=

h∗l

=

hgs

if g > 0

hg

Otherwise

hl

if g > 0

hls

Otherwise

(1)

∂(˛l l ) ∂(˛l l ul ) + = −g ∂t ∂x

(2)

(4)

hg and hl are the phasic enthalpies of gas and liquid, respectively. hgs and hls are the phasic saturation enthalpies. Conservation of momentum is also applied to each phase separately to obtain the following two equations for the gas and liquid phases, respectively: ˛g g ∂(ug )

∂t

+ ˛g g ug

∂(ug ) ∂p + ˛g − ˛g g G ∂x ∂x

= −aint FI ug − ul (ug − ul ) − g (uint − ug ) ˛l l ∂(ul )

∂t

+ ˛l l ul

(5)

∂(ul ) ∂p + ˛l − ˛l l G ∂x ∂x

= aint FI ug − ul (ug − ul ) + g (uint − ul )

(6)

In the equations above G is the gravity acceleration and p is the pressure of the system. The term with the interfacial friction coefﬁcient (FI) accounts for momentum losses due to interfacial friction. These losses are equal in magnitude and opposite in sign for the two phases. Terms including in the momentum equations account for the momentum lost or gained by the new mass appearing at the interfacial velocity (uint ). It is assumed there is no momentum storage at the interface, therefore liquid and vapor interface velocities are equal to each other. Conservation of energy is also applied to each phase separately. This yields the following two equations for gas and liquid phases, respectively:

∂(˛g g eg ) ∂(˛g g eg ug ) ∂˛g ∂(˛g ug ) + +p +p ∂t ∂x ∂t ∂x = Hig aint (Ts − Tg ) + g h∗g

∂(˛g g ) ∂(˛g g ug ) + = g ∂t ∂x

(3)

(7)

∂(˛l l el ) ∂(˛l l el ul ) ∂˛ ∂(˛l ul ) + +p l +p = Hil aint (Ts − Tl ) − g h∗l ∂t ∂x ∂t ∂x (8)

In these equations: • ˛g and ˛l are void fractions for gas and liquid phases, respectively. • g and l , are gas and liquid densities. • ug and ul are gas and liquid velocities.

In these equations, eg and el are the speciﬁc energies for gas and liquid phases, respectively. Finally, conservation of volume is applied: ˛g + ˛l = 1

(9)

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263

The constitutive relations for this system of equations are presented in a discrete form in the next section. This leads to a closed mathematical system of equations.

257

n+1

n+1 Fmg = (˛)g,i+1/2 un+1 (u˙ n+1 − u˙ n+1 )+˛ ¯ n+1 (pn+1 − pn+1 ) i g,i g,i+1/2 g,i+1 g,i+1/2 i+1 n+1 n+1 (un+1 − un+1 ) − x(˛)g,i+1/2 G + x¯ g,i+1/2 int ,i+1/2 g,i+1/2

3. Discretized equations In this section we introduce the discrete equations for the twoﬂuid model presented in previous section. Numerical methods used to obtain the discrete equations are ﬁnite volume based on a staggered mesh. In this conﬁguration scalar properties (pressure, energies, densities and void fraction) are deﬁned at cell centers and vector properties (velocities) are deﬁned at cell faces. Fig. 1 shows a schematic of the cell conﬁguration for this type of discretization. Based on the above, we write the discrete equations in a residual form, similar to Mousseau’s work (Mousseau, 2004). In the following set of equations superscripts “n” and “n + 1” represent the old and new time steps, respectively, and subscripts containing “i” determine the spatial position. The temporal discretization is based on the theta method formulation. Different schemes can be obtained for different values of . Explicit scheme is obtained for = 0, the implicit scheme for = 1 and the Crank–Nicolson scheme for = 1/2. The theta method is unconditionally stable for any choice of ∈ [1/2, 1] (Higham, 2000; Farago, 2013).

(16) n+1

n+1 Fml = (˛)l,i+1/2 un+1 (u˙ n+1 − u˙ n+1 )+˛ ¯ n+1 (pn+1 − pn+1 ) l,i i l,i+1/2 l,i+1 l,i+1/2 i+1 n+1 n+1 (un+1 − un+1 ) − x(˛)l,i+1/2 G − x¯ g,i+1/2 int ,i+1/2 l,i+1/2

rescl =

(10)

x n+1 n+1 n (˛ − ˛nl,i l,i ) + (1 − )Fcln + Fcln+1 t l,i l,i

(11)

(17) 3 Residuals from the conservation of energy are calculated at cell centers and given as: n+1/2

reseg =

xpi x n+1 n+1 n+1 n n (˛ e −˛ng,i g,i eg,i )+ t g,i g,i g,i t

n+1 Fcg =

n+1

·

˛

g,i+1/2

·

un+1 − ˛ g,i+1/2

(18)

g,i−1/2

un+1 g,i−1/2

n+1/2

resel =

Fcln+1 =

·

n+1

˛

l,i+1/2

·

n+1

un+1 − ˛ l,i+1/2

l,i−1/2

n+1 − xg,i

n+1 + xg,i

2 Residuals from the conservation of momentum are calculated at cell faces, and given as: resmg =

t

n+1

·

n+1

+ pn+1 (˛˙ n+1 un+1 − ˛˙ n+1 un+1 ) i g,i+1/2 g,i+1/2 g,i−1/2 g,i−1/2

Feln+1

=

·

n+1

(˛e)l,i+1/2 un+1 l,i+1/2

·

n+1

(20)

− (˛e)l,i−1/2 un+1 l,i−1/2

+ pn+1 (˛˙ n+1 un+1 − ˛˙ n+1 un+1 ) i l,i+1/2 l,i+1/2 l,i−1/2 l,i−1/2 n+1 n+1 n+1 ∗ − xHiln+1 an+1 (Ts,i − Tl,i ) + xg,i h int

(21)

n n+1 − ung,i+1/2 ) + (1 − )Fmg + Fmg

(14)

resml =

·

(˛e)g,i+1/2 un+1 − (˛e)g,i−1/2 un+1 g,i+1/2 g,i−1/2

n+1 n+1 n+1 n+1 n+1 ∗ − xHig aint (Ts,i − Tg,i ) − xg,i h

(13)

n+1/2 x (˛)g,i+1/2 (un+1 g,i+1/2

(19)

where Feg and Fel are the energy numerical ﬂuxes given by:

un+1 l,i−1/2

− ˛nl,i ) (˛n+1 l,i

+ (1 − )Feln + Feln+1

(12)

xpi x n+1 n+1 n+1 n n (˛ e − ˛nl,i l,i el,i ) + t l,i l,i l,i t

n+1 Feg =

n+1

(˛n+1 − ˛ng,i ) g,i

n n+1 + Feg + (1 − )Feg

where Fcg and Fcl are the mass numerical ﬂuxes given by:

− x(a¯ n+1 FI n+1 )(un+1 − un+1 ) un+1 − un+1 int g,i+1/2 l,i+1/2 g,i+1/2 l,i+1/2

1 Residuals from the conservation of mass are calculated at cell centers and given as: x n+1 n+1 n n n+1 (˛ rescg = − ˛ng,i g,i ) + (1 − )Fcg + Fcg t g,i g,i

+ x(a¯ n+1 FI n+1 )(un+1 − un+1 ) un+1 − un+1 int g,i+1/2 l,i+1/2 g,i+1/2 l,i+1/2

n+1/2 x n+1 n (˛)l,i+1/2 (un+1 − unl,i+1/2 ) + (1 − )Fml + Fml l,i+1/2 t (15)

where Fmg and Fml are the momentum numerical ﬂuxes given by:

Fig. 1. Schematic of the staggered cell conﬁguration.

A quantity with an over-bar is an arithmetically averaged quantity: ¯ i+1/2 = ˚i + ˚i+1 ˚ 2

(22)

Because scalar properties are deﬁned at cell centers, a value for these properties is needed at cell interfaces in order to calculate numerical ﬂuxes, i.e. the variables with a “dot” in Eqs. (12), (13), (20) and (21). Dotted quantities are called donored quantities in RELAP5 nomenclature. The same holds for momentum numerical ﬂuxes in Eqs. (16) and (17), where velocities are needed at cell centers. Order of accuracy for our discretization depends on the method used to calculate the donored quantities in numerical ﬂuxes. There are many numerical schemes one can apply to calculate these donored quantities. For this paper, results from 1st-order upwind scheme and ﬂux limiter schemes (including the new scheme developed in (Abu Saleem and Kozlowski, 2015)) are shown. This will be discussed in Section 5.

258

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263

4. Closure equations

5.1. First-order upwind

In order to close the discrete system of equations discussed in previous section, equations of state are needed. These equations relate densities and temperatures of the two phases to other state variables (pressure and internal energy). The linearized equations of state for densities are:

For a scalar quantity ˚ (this can be , e, p or ˛) calculated at the center of the ith cell, donored quantities at cell faces are deﬁned as follows:

n+1 n g,i = g,i +

∂g ∂p

∂l ∂p

n+1 n l,i = l,i +

n

(pn+1 − pni ) + i

n

(pn+1 − pni ) + i

∂g ∂eg

∂l ∂el

n

n+1 n (eg,i − eg,i )

n+1 n (el,i − el,i )

(24)

n+1 Tg,i

=

n Tg,i

+

n+1 n Tl,i = Tl,i +

∂T ∂p

n+1 n Ts,i = Ts,i +

n

∂T ∂p

∂T ∂p

(pn+1 − pni ) i

n

(pn+1 i

− pni ) +

n

(pn+1 − pni ) + i

∂Tg ∂eg

∂Tl ∂el

n

n+1 (eg,i

n − eg,i )

n+1 n (el,i − el,i )

hn+1 = hnls,i + ls,i

n+1 hn+1 = eg,i + g,i

n+1 hn+1 = el,i + l,i

∂hgs ∂p

∂hls ∂p

(27)

(pn+1 − pni ) i

(28)

n

(ui−1/2 + ui+1/2 ) +

(pn+1 − pni ) i

pn+1 i n+1 g,i

pn+1 i n+1 l,i

ui+1/2 ui+1/2

1 = 2

(ui−1/2 − ui+1/2 )

(gi + gi+1 ) +

ui+1/2

n

(32)

(33)

ui+1/2

(gi − gi+1 )

(26)

n

(˚i − ˚i+1 )

ui+1/2

As an example, donored gas density (˙ gi+1/2 ), and donored gas velocity u˙ gi are given as: ˙ gi+1/2

for saturation enthalpies and speciﬁc enthalpies of the two phases: hn+1 = hngs,i + gs,i

1 u˙ i = 2

(25)

(˚i + ˚i+1 ) +

Similarly, for a velocity u calculated at the face of the ith cell, the donored quantity at the cell center is deﬁned as:

for saturation temperature and phasic temperatures of gas and liquid:

˙ i+1/2 ˚

ui+1/2

(23)

n

1 = 2

(29)

(30)

(31)

All derivatives in Eqs. (23)–(29) are calculated numerically using water properties. Water properties used in this paper are based on the International Association for Properties of Water and Steam Industrial Formulation 1997 (IAPWS IF-97) standard. XSteam, an implementation MATLAB code of the IAPWS IF-97 standard formulation by Magnus Holmgren (Holmgren, 2006), is used to obtain all necessary water properties. The closure relations discussed here can be substituted directly into Eqs. (10)–(21) to reduce the number of unknown variables to six, namely: pressure, void fraction, gas and liquid velocities and internal energies of the two phases.

5. Calculation of donored quantities As mentioned earlier, the order of accuracy of our discretized equations depends on the numerical scheme used to calculate donored quantities in the numerical ﬂuxes. There are many numerical schemes one can apply to calculate these donored quantities. In this paper, results for the 1st-order upwind scheme and ﬂux limiter schemes (including the new scheme developed in (Abu Saleem and Kozlowski, 2015)) are shown.

1 u˙ gi = 2

(ugi−1/2 + ugi+1/2 ) +

ui+1/2 ui+1/2

(ugi−1/2 + ugi+1/2 )

Since the 1st-order upwind scheme is a monotone scheme, we expect the solution due to this discretization to be free of nonphysical numerical oscillations. However, the solution is expected to be smeared due to excessive dissipation of the scheme. 5.2. Flux limiter schemes The ﬂux limiter approach is based on a combination of a monotone 1st-order scheme and a linear higher-order accurate scheme. The two schemes are combined in a way such that the 1st-order scheme is used near discontinuities to prevent spurious oscillations, and the high-order scheme is used elsewhere to obtain the highest order of accuracy possible. The weighted interpolation between the two schemes depends on the smoothness of the solution. Because the smoothness parameter is calculated for a given variable, and depends on the characteristics of the solution, it is hard to implement the ﬂux limiter approach to a system of nonlinear equations like the two-ﬂuid model. One way to implement the ﬂux limiter method to a system of non-linear equations is to approach the problem in an empirical manner. Here we consider the empirical approach suggested by Torro (Toro, 2009). His approach for a system of m variables is based on calculating the smoothness parameters for each variable (q = qk , k = 1, 2, . . ., m) in the system: L ri+1/2 =

R ri+1/2 =

qi−1/2 qi+1/2 qi+3/2 qi+1/2

=

qi − qi−1 qi+1 − qi

(34)

=

qi+2 − qi+1 qi+1 − qi

(35)

Then, the value for the ﬂux limiter function is calculated as follows: L R ϕ(r) = mink [min(ϕ(ri+1/2 ), ϕ(ri+1/2 ))]

(36)

The ﬂux limiter ϕ in Eq. (36) is a function of the smoothness parameter. There are many different forms of ﬂux limiters. In this

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263

259

paper, the new scheme developed in (Abu Saleem and Kozlowski, 2015) is investigated. For this scheme, the ﬂux limiter is given by:

ϕ(r)NEW

⎧ 8r 9−ı ⎪ 0

(37)

ı+3

In this equation ı is weighting parameter such that ı ∈ [0, 1]. The development of this scheme was based on a combination of the 1st-order upwind scheme and the 3rd-order QUICK scheme (Abu Saleem and Kozlowski, 2015). Other classical limiters are implemented for comparison (Waterson and Deconinck, 2007; Kadalbajoo and Kumar, 2006). These limiters are the Minbee limiter function given by:

ϕ(r)MINBEE =

r

0≤r≤1

1

if r ≥ 1

(38)

And the Superbee limiter function given by:

⎧ 1 ⎪ 2r if 0 ≤ r ≤ ⎪ ⎪ 2 ⎪ ⎪ ⎨ 1 1 if ≤ r ≤ 1 ϕ(r)SUPERBEE = 2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ r if 1 ≤ r ≤ 2 2

(39)

ϕ(r)

=

0 min

if r ≤ 0 2 1+ 2, , , r r+

∈ [0, ∞]

(40)

This scheme was based on a combination of the 1st-order upwind scheme and the 2nd-order upwind scheme (Kadalbajoo and Kumar, 2006). Results from Eq. (36) are then applied to all m ﬂux components. For the two-phase problem we have the following variables:

⎛˛ ⎞ g

⎜p⎟ ⎜ ⎟ ⎜u ⎟ ⎜ g⎟ ⎟ q=⎜ ⎜ ul ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ eg ⎠

n+1 el,i for i ∈ [1, 2, 3, . . ., N]. These equations need to be solved simultaneously by a non-linear algebraic equations solver. In this paper, a solver by Kelley (Kelley, 2003) is used. The solver uses Newton’s method with a direct factorization of the Jacobian. Because the goal of this research is improvement of the numerical method rather than the non-linear solver, we use the non-linear solver without any modiﬁcations.

7. Numerical results

if r ≥ 2

Both limiters were based on a combination of the 1storder upwind scheme and the 2nd-order Lax–Wendroff scheme (Waterson and Deconinck, 2007). The last limiter considered in this paper is Kumar’s limiter function given by: KUMAR

Fig. 2. Illustration of the faucet ﬂow problem.

(41)

el Finally, the residual form of the equations can be written using the following form for the ﬂux terms: F = F 1st + ϕ(r)(F High − F 1st ) where F1st is the numerical ﬂux using the monotone 1st-order upwind scheme, FHigh is the numerical ﬂux using the higher order accurate scheme. 6. Non-linear solver For a discrete domain of N cells, Eqs. (10), (11), (14), (15), (18) and (19) form a system of 6N non-linear algebraic equations with 6N , pn+1 , un+1 , un+1 , en+1 and unknowns. The unknowns are ˛n+1 g,i i g,i+1/2 l,i+1/2 g,i

7.1. Faucet ﬂow problem In this problem, the ﬂuid accelerates through the spatial domain under the effect of gravity, and the initial proﬁle of the void fraction is moved out of the domain under this acceleration. With the convection of the void fraction proﬁle, there is also a moving discontinuity in the proﬁle. This feature of the void fraction proﬁle makes it a very good test for the dissipative properties of the numerical scheme and its stability near discontinuities. This problem has an analytical solution for the case of negligible pressure variation and no interfacial and wall interactions. The analytical solution for the void fraction is given by (Qin, 1992):

˛g (x, t) =

⎧ 0 ⎨ 1 − ˛l ul ⎩

2Gx + ul

1 − ˛0l

2

x ≤ ul t +

Gt 2 2

(42)

Otherwise

The physical domain of this problem consists of a vertical tube containing a uniform column of liquid moving in a vapor annulus. Fluids are initially at saturation and inlet velocities are kept constant. Outlet pressure is also maintained constant at atmospheric and a gravity ﬁeld (G) is applied at the start of the simulation. The problem is illustrated in Fig. 2, and initial and boundary conditions are shown in Table 1. The numerical schemes discussed earlier are applied to discretize the two-ﬂuid model as presented in previous sections of this paper. The analytical solution in Eq. (42) assumes no phase transfer, hence g in the model was set to zero. Table 1 Initial and boundary conditions of the faucet ﬂow problem. Tube length L Initial and boundary liquid velocities ul , uin l Initial and boundary gas velocities ug , uin g Initial and boundary volume fraction of liquid ˛0l , ˛in l Initial and boundary pressure p0 , pout 0 Gravity acceleration G

1m 10 m/s 10 m/s 0.8 1.5 bar 9.8 m/s2

260

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263 0.24

0.24 N=40 N=80 N=160 N=320 Exact Solution

0.235

0.23

0.225

Gas void fraction

Gas void fraction

0.23

θ = 1 (Implicit) θ = 5/6 θ = 2/3 θ = 1/2 (Crank−Nicolson) Exact solution

0.235

0.22 0.215 0.21 0.205

0.225 0.22 0.215 0.21 0.205

0.2 0.2

0.195

0

0.2

0.4

0.6

0.8

1

0.195 0

Distance (m)

0.2

0.3

0.4

0.5

0.6

Numerical results from this simulation are shown for the 1storder upwind scheme and ﬂux limiter schemes, including the new scheme derived in (Abu Saleem and Kozlowski, 2015). 7.1.1. 1st-order upwind scheme For this set of results we consider two cases: the implicit discretization ( = 1) and the Crank–Nicolson discretization ( = 1/2). Simulation time was set to 0.045 s with a time step of size 10−4 s. All results for the 1st-order upwind scheme are characterized by a lack of any spurious oscillations. However, the solution is smeared as compared to the results of ﬂux limiter schemes. Fig. 3 shows results for the 1st-order upwind scheme with implicit method for different mesh sizes: N = 40, 80,160 and 320. Fig. 4 shows results for the Crank–Nicolson method with the same mesh sizes. It can be observed that Crank–Nicolson method is less dissipative than the implicit method. This result can also be deduced from the modiﬁed equation analysis of the numerical scheme. 7.1.2. Flux-limiter schemes As discussed earlier, application of this method for a system of non-linear equations is based on an empirical approach. The smoothness parameters deﬁned in Eqs. (34) and (35) were calculated for all the variables in Eq. (41). It was found that using smoothness parameters based on void fraction (˛) leads to the best

0.8

0.9

1

Fig. 5. Numerical results for new scheme with ı = 1 using different values of , N = 320.

results. Hence all results from ﬂux limiter schemes shown in this section are based on void fraction. Different combinations of limiters and values were investigated to reach the best numerical solution. In practice, different limiters led to results with little qualitative difference. A scheme based on the Minbee limiter was found to be the most dissipative among ﬂux limiter schemes. Changing the value of resulted in slightly different results, with implicit method ( = 1) more dissipative than the Crank–Nicolson method ( = 1/2). The numerical solutions were free of spurious oscillations for all cases, as expected. Fig. 5 shows numerical results due to the new scheme derived in (Abu Saleem and Kozlowski, 2015) with ı = 1 for different choices of . It shows better results for the Crank–Nicolson method ( = 1/2), where dissipation is minimal. Fig. 6 shows numerical results for four different schemes: 1storder upwind scheme, Minbee scheme, Superbee scheme and the new scheme with ı = 1. All results are shown with = 1/2. Results from the 1st-order upwind scheme are the most smeared as compared to other schemes. Results from the Minbee limiter seem to be more smeared than the other two limiters, and the new scheme yields numerical solution that is most accurate.

0.24

0.24 N=40 N=80 N=160 N=320 Exact Solution

0.235

0.7

Distance (m)

Fig. 3. Results for the 1st-order-upwind with implicit temporal discretization.

1st order upwind scheme Minbee flux limiter scheme Superbee flux limiter scheme New scheme with δ=1 Exact solution

0.235

0.23

0.23

0.225

0.225

Gas void fraction

Gas void fraction

0.1

0.22 0.215 0.21 0.205

0.22 0.215 0.21 0.205

0.2 0.2

0.195 0

0.2

0.4 0.6 Distance (m)

0.8

1

0.195

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Distance (m) Fig. 4. Results for the 1st-order-upwind with Crank–Nicolson temporal discretization.

Fig. 6. Numerical results for different ﬂux limiter functions = 1/2 and N = 320.

1

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263 Table 3 Initial conditions of the shock tube problem.

−5

−5.5

−6

log(L)

261

−6.5

Left side of the pipe

Right side of the pipe

˛Lg = 0.25 pL = 20 MPa egL = 2824 kJ/kg elL = 1311 kJ/kg vLg = 0.0 m/s vLl = 0.0 m/s

˛Rg = 0.10 pR = 10 MPa egR = 2836 kJ/kg elR = 1330 kJ/kg vRg = 0.0 m/s vRl = 0.0 m/s

−7 7

2.2

−7.5

−5.5

−5

−4.5

−4

−3.5

−3

−2.5

−2

log(dx) Fig. 7. Convergence of numerical solution for schemes with different limiters. Table 2 Convergence of numerical solution for schemes with different limiters. Numerical scheme

Order of the scheme

1st-order upwind Minbee limiter Superbee limiter Kumar limiter New limiter with ı = 1 New limiter with ı = 0

0.52 0.67 0.74 0.73 0.78 0.80

1 qN (xi ) − qE (xi ) N

1.8

1.6

1.4

1.2

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Distance (m)

The exact solution in Eq. (42) was used to determine the accuracy of numerical schemes. Accuracy was assessed on the basis of the L1 norm, deﬁned as: N

L1 =

2

Pressure (Pa)

−8.5 −6

TRACE New solver

1st order upwind scheme Minbee flux limiter scheme Superbee flux limiter scheme Kumars flux limiter scheme New scheme with δ=0

−8

x 10

(43)

i=1

where N is the number of cells and qN (xi ) is the numerical solution at the ith cell. Convergence rate is represented by a log-log plot of the norm of the error versus mesh size. Fig. 7 shows the results of the space convergence study. The 1st-order upwind scheme converges with the order of 0.52, the Minbee limiter converges with the order of 0.67 and the new scheme with ı = 1 converges with the order of 0.78. The complete results are shown in Table 2.

Fig. 9. Pressure proﬁle for shock tube problem at time = 0.5 ms.

step of 10−5 s was used. Results of the solver based on the new numerical scheme derived in (Abu Saleem and Kozlowski, 2015) are compared to those from the nuclear thermal-hydraulics code TRACE. Fig. 9 shows the pressure proﬁle in the pipe after 0.5 ms of the transient. A formation of a new constant state in the middle of the pipe can be observed. This constant state corresponds to an increasing pressure in the right half of the pipe, and decreasing pressure in the left half. Consequently, phase change happens in both halves, resulting in increasing gas void fraction in the left half, and decreasing gas void fraction in the right half, this is shown in Fig. 10. Velocities and internal energies of the two phases are

7.2. The shock tube problem

TRACE New solver 0.25

Gas void fraction

A schematic of the shock tube problem is shown in Fig. 8. Results from this problem are aimed to show the solver’s ability to handle large differences between temperatures of the two phases. The setup of this problem is a pipe containing water in its two different phases (gas and liquid). Initially the left and right halves of the tube have two different states at saturation and separated by a diaphragm. The left and right states for this problem are deﬁned in Table 3. Both ends of the pipe are kept closed, and the diaphragm is removed at the beginning of the simulation to start the transient. The domain of the problem was divided into 40 cells, and a time

0.3

0.2

0.15

0.1

0.05

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Distance (m) Fig. 8. Schematic of the shock tube problem.

Fig. 10. Gas void fraction proﬁle for shock tube problem at time = 0.5 ms.

1

262

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263 140

0.3

Gas velocity − TRACE Liquid velocity − TRACE Gas velocity − New solver Liquid velocity − New solver

120

N=20 N=40 N=80

0.25

Void fraction

Velocity (m/s)

100

80

60

0.2

0.15

40

0.1 20

0.05 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

0.2

0.4

1

0.6

0.8

1

Distance (m)

Distance (m)

Fig. 13. Void fraction proﬁle for shock tube problem at time = 0.5 ms for different mesh sizes.

shown in Figs. 11 and 12, respectively. It can be seen that the results in these ﬁgures show a qualitatively good agreement between the new solver and TRACE. The solution proﬁle for void fraction was investigated for different mesh sizes (N = 20, 40 and 80). This is shown in Fig. 13. Time convergence was also investigated by changing the time step of the simulation. Three time steps were considered (dt = 10−4 , 2.5 × 10−5 and 6.25 × 10−6 s). Results of this study for void fraction are shown in Fig. 14. The two-phase shock tube problem is also used to compare the 1st-order upwind scheme to the new scheme developed in (Abu Saleem and Kozlowski, 2015). This is done by comparing the mass error of the solver due to each scheme. Because both ends of the pipe are kept closed, the total mass of the system is constant. Any change in the total mass of the system is considered an error. This error can be related to many sources, like the numerical scheme, the non-linear solver, or any other source in modeling the problem. To account for the error due to the numerical scheme, the mass error was observed by changing only the mesh size for the simulation. Three mesh sizes were examined: 40, 80 and 160 cells, corresponding to spatial step sizes of: 0.025, 0.0125 and 0.00625 m,

0.3

dt=1E−4 sec dt=2.5E−5 sec dt=6.25E−6 sec

0.25

Void fraction

Fig. 11. Velocity proﬁles for shock tube problem at time = 0.5 ms.

0.2

0.15

0.1

0.05

0

0.2

0.4

0.6

0.8

1

Distance (m) Fig. 14. Void fraction proﬁle for shock tube problem at time = 0.5 ms for different time steps.

6

2.8

x 10

−12.5

Liquid − TRACE Liquid − New solver Gas − TRACE Gas − New solver

−14

2.2 2 1.8

−14.5 −15 −15.5 −16

1.6

−16.5 −17

1.4 1.2

suggested scheme 1st−order upwind

−13.5

2.4

Log(MassError)

Internal energy (J/Kg)

2.6

−13

−17.5 −5.2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

−5

−4.8

−4.6

−4.4

−4.2

−4

Log(dx)

Distance (m) Fig. 15. Convergence rate for the shock tube problem. Fig. 12. Internal energy proﬁles for shock tube problem at time = 0.5 ms.

−3.8

−3.6

R.A. Abu Saleem et al. / Nuclear Engineering and Design 301 (2016) 255–263

respectively. Accuracy was assessed by calculating the percentage mass error as follows: Mass Error = 100 ∗

|M f − M i | Mi

263

to 20% (Brooks et al., 2012). These aspects of the problem should be tackled intensively to achieve better modeling of two-phase ﬂows.

(44)

where Mf is the total mass of the system at the end of the transient and Mi is the total mass of the system at the beginning of the transient. The convergence rate is represented by a log-log plot of the mass error versus mesh size. Fig. 15 shows the results of the spatial convergence of both schemes. We can observe that the solver with the new scheme exhibits better rate of convergence as compared to the 1st-order upwind scheme. 8. Summary and conclusions A high-order solver based on ﬂux limiter schemes and the theta method was developed for solving the two-ﬂuid model. This approach was compared to the 1st-order upwind scheme used in most existing codes in nuclear industry. The new two-phase twoﬂuid model solver was tested for two benchmark problems: the faucet ﬂow problem and the shock tube problem. The analytical solution for void fraction in the faucet ﬂow problem poses a propagation of discontinuity. This is useful to test the different schemes for their abilities to handle discontinuities, and compare their TVD and dissipation properties. For the 1st-order upwind scheme, results were smeared due to numerical dissipation; this dissipation was minimum for the case of = 1/2 (Crank–Nicolson method). Numerical results for the solver with ﬂux limiter schemes showed better agreement with the analytical solution and an improved rate of convergence. Numerical solutions were less smeared and posed no oscillations near the discontinuity. Most satisfactory results were attained by using the new scheme developed in (Abu Saleem and Kozlowski, 2015) and the theta method with = 1/2. The shock tube problem showed the solver’s capability to handle fast transients with discontinuities and phase change. Results from this test were compared to those from the U.S. NRC TRACE code, and good qualitative agreement between the two results was observed. Errors in the system’s total mass were calculated for two numerical schemes: 1st-order upwind scheme and the new scheme developed in (Abu Saleem and Kozlowski, 2015). The new scheme showed better reduction in mass error as compared to the 1st-order upwind scheme. It is important to mention that there are other sources of errors and uncertainties in the modeling of two-phase ﬂows. For the twoﬂuid model, for example, errors in the closure relations can be up

References Abu Saleem, R., Kozlowski, T., 2015. Development of high-resolution total variation diminishing scheme for linear hyperbolic problems. J. Comput. Eng. 2015, 10, AN:575380. Ansys Inc., 2006. Fluent 6. 3 User’s Guide. Balcazar, N., Jofre, L., Lehmkuhl, O., Castro, J., Rigola, J., 2014. A ﬁnite-volume/levelset method for simulating two-phase ﬂows on unstructured grids. Int. J. Multiph. Flow 64, 55–72. Bonometti, T., Magnaudet, J., 2007. An interface-capturing method for incompressible two-phase ﬂows. Validation and application to bubble dynamics. Int. J. Multiph. Flow 33, 109–133. Brooks, C.S., Hibiki, T., Ishii, M., 2012. Interfacial drag force in one-dimensional twoﬂuid model. Prog. Nucl. Energy 61, 57–68. Division of Safety Analysis, Ofﬁce of Nuclear Regulatory Research, US Nuclear Regulatory Commission, Washington, DC 20555-0001, June 2010. TRACE V5.0 Theory Manual. Field Equations, Solution Methods, and Physical Models. Division of Systems Research, Ofﬁce of Nuclear Regulatory Research, US Nuclear Regulatory Commission, Washington, DC 20555-0001, October 2010. RELAP5/MOD3.3 Code Manual. Farago, I., 2013. Convergence and stability constant of the theta method. In: Proc. Of Conference of Applications of Mathematics, Prague. Fullmer, W.D., Lopez De Bertodano, M.A., Zhang, X., 2013. Veriﬁcation of a high-order ﬁnite difference scheme for the one-dimensional two-ﬂuid model. J. Comput. Multiph. Flows 5, 139–156. Godunov, S.K., 1957. A Difference method for the numerical calculation of discontinuous solution of hydrodynamic equations. Mat. Sb. 47, 271–306, Translated as JPRS 7225 by U.S. Dept. of Commerce, pages 36, 41 and 106 (1960). Higham, D.J., 2000. Mean-square and asymptotic stability of the stochastic theta method. SIAM J. Numer. Anal. 38, 753–769. Holmgren, M., 2006. XSteam – Properties of Water and Steam for MatLab, Retrieved 03.11.12 from http://www.x-eng.com/. Kadalbajoo, M.K., Kumar, R., 2006. A high resolution total variation diminishing scheme for hyperbolic conservation law and related problems. Appl. Math. Comput. 175, 1556–1573. Kelley, C.T., 2003. Solving Nonlinear Equations with Newton’s Method, 1st edition. Society for Industrial and Applied Mathematics, pp. 27–55. Lee, M.S., Aute, V., Riaz, A., Radermacher, R., 2012. A review on direct two-phase, phase change ﬂow simulation methods and their applications. In: Int Refrig. Air Cond. Conf. Mousseau, V.A., 2004. Implicitly balanced solution of the two-phase ﬂow equations coupled to nonlinear heat conduction. J. Comput. Phys. 200, 104–132. Qin, H.Q., 1992. Numerical benchmark test 2.1: faucet ﬂow. Multiph. Sci. Technol. 6, 577–590. Sokolowski, L., Koszela, Z., 2012. RELAP5 capability to predict pressure wave propagation phenomena in single- and two-phase ﬂow conditions. J. Power Technol. 92, 150–165. Sweby, P.K., 1984. High resolution schemes using ﬂux limiters for hyperbolic conservation laws. SIAM J. Numer. Anal. October (21), 995–1011. Toro, E.F., 2009. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, vol. 480., 3rd edition Springer, Dordrecht/Heidelberg, pp. 512–514. Waterson, N.P., Deconinck, H., 2007. Design principles for bounded higher-order convection schemes – a uniﬁed approach. J. Comput. Phys. 224, 182–207.