A new numerical algorithm for two-phase flows drift-flux model with staggered grid in porous media

. FLICA4 is a 3D compressible code dedicated to reactor core analysis. It solves a compressible drift-ﬂux model for two-phase ﬂows in a porous medium [2]. To deﬁne convective ﬂuxes, FLICA4 uses a speciﬁc ﬁnite volume numerical method based on an extension of the Roe’s approximate Riemann colocated solver [3]. Nevertheless, analysis of this method shows that at low Mach number, it is necessary to apply modiﬁcations to the 2D or 3D geometries on a cartesian mesh otherwise this method does not converge to the right solution when the mach number goes to zero [4]. For this reason, we apply a so-called “pressure correction“. Although this correction is necessary to reach the required precision, it may produces some checkerboard oscillations in space in the situations we are interested in, especially in the 1D case. Since these checkerboard oscillations are sometimes critical and may lead to unstable solutions in some cases, we investigate another numerical algorithm to solve this compressible drift-ﬂux model in the low Mach regim. The aim of this work is to propose a new compressible scheme accurate and robust at low Mach number on staggered grid since checkerboard oscillations cannot exist on this type of discretisation [8]. The accuracy and robustness of this new scheme are veriﬁed in low Mach regime with test cases describing a simpliﬁed nuclear core ”Boiling channel”. The behavior of this scheme is also tested in the compressible regime with or without shock waves.


Introduction
FLICA4 is a 3D compressible code dedicated to reactor core analysis. This code solves for a compressible drift-flux model for two-phase flows in a porous medium [2]. To derive convective fluxes, FLICA4 uses a specific finite volume numerical method based on an extension of the Roe's approximate Riemann colocated solver [3]. Nevertheless, an analysis of this type of method shows that at low Mach number, it is necessary to apply modifications to the 2D or 3D geometries on a cartesian mesh otherwise this method does not converge to the right solution when the mach number goes to zero [3]. For this reason, we apply a so-called "pressure correction". Although this correction is necessary to reach the required precision, it may produces some checkerboard oscillations in space, especially in the 1D case.
Since these checkerboard oscillations are sometimes critical and may lead to unstable solutions in some cases, we also investigate another numerical algorithm to solve this compressible drift-flux model in the low Mach regim. The key point is to develope a compressible solver on staggered grid since checkerboard oscillations cannot exist on this type of discretisation [8]. The aim of this work is to present such a compressible scheme and to verify it in the low Mach regime with test cases describing a simplified nuclear core. Then, the behavior of this scheme is tested in the compressible regime with or without shock waves.
The compressible solver on staggered grid that we develop follows the finite volume approach for all 4 balance equations. The time discretization of the equations is based on a semi-implicit scheme. Since the equations are not linear, the solution at each time step is obtained by a Newton-Raphson iterative method. This method gives a linear system of equations for the increments of the principal variables. The chosen solution algorithm [6] consists at first in eliminating the velocity increments as functions of the pressure increments by rewriting the momentum equations. Substituting the velocity increments into the non-linear system gives a system involving only the pressure increments. The successive elimination of the scalar variables other than the pressure variable gives a linear system on the pressure. The resolution of this linear system allows to determine the velocity and the other variables. Preliminary numerical experiments are presented and compared with analytic solutions [7,12]. §4.1 shows that the current numerical scheme computes solutions very close to the analytical solutions in the scope of low Mach number flows (see figures 4, 5 and 6). The numerical method uses nonconservative formulation of the four equations. Therefore, difficulties to get precise shock wave solutions were expected. Figures 7 and 9 indeed show a lack of precision of the new numerical scheme to capture shock waves although it is stable. Nevertheless, the scheme captures compressible regular solutions with accuracy (see figure 8). As a consequence, our scheme is robust and is enough accurate to capture at the same time subsonic compressible solutions and low Mach solutions: this property is the main requirement for our applications.

Porous 4 equations model
FLICA4 is a four equation code. The equations are mixture mass balance, phasic mass balance, mixture momentum balance and mixture energy balance. The 4 equations are averaged in space, time and phase over control volumes. A drift-flux model is used to account for the slip between vapor and liquid phases. The fluid is compressible and the two-phases are assumed to be at the same pressure. One of the phases is assumed to at saturation temperature. The equations written in the non-conservative are given by: Above, the blue terms are given by the user, the red terms are computed with physical modeling and the green terms are defined by the equation of state The nomenclature is given as follows : t is the time, φ is the porosity, ρ is the mixture density, h is the mixture specific enthalpy, − → V is the mixture velocity, − → V r is the relative velocity between the two phases, C is the vapor mass concentration, P is the mixture pressure, e is the mixture internal energy, H v is the steam total enthalpy, H l is the liquid total enthalpy, − → g is the gravity force, Q is the power density,τ is the stress tensor, τ f is the friction forces and q is the heat flux. Contribution of the diffusive terms ∇ · (φK cv ∇C), ∇ ·τ and ∇ · q are neglected.

Full time and space discretization
For space discretization, a staggered grid (see figure 1) is used where scalar variables (pressure, density, enthalpy, etc) are computed at center of the control volumes. Velocity or momentum variables are computedb on volume faces. It differs from a collocated grid, where all the variables are computed at the same position. Using a staggered grid is a simple way to avoid odd-even decoupling between pressure and velocity. Odd-even decoupling is a discretization error that can occur on collocated grids and which leads to checkerboard patterns in the solutions [15].

Discretization of the mixture mass equation
The time discretization of the mixture mass equation is based on a semi-implicit scheme. More precisely, the mass fluxes are approximated on the interface using a donor semi-implicit scheme where the velocity is implicit while the scalar variables are explicit. The space discretization in the three directions of the mass fluxes in the mixture mass equation is done by approximating the density on the cell faces by using a donor formulation.
More precisely, the finite volumes discretization of the mixture mass equation involves its integration in time between t n and t n+1 and in space on an elementary control volume M K c (see figure 2): By using the divergence theorem, equation (3) becomes: or in discrete form: To establish the discrete mixture mass equation, all we need now is to approximate the flux F n+1 σ on the faces of a given cell M K c . To ensure the stability of the numerical scheme, we use a donor formulation of the convection term: In (6) As a result, the discrete form of the mixture mass equation (5) can be written as: Taking into account the state equation (2), equation (7) becomes:

Discretization of the mixture internal energy balance equation
The spatial and temporal discretization of the mixture internal energy balance equation is done by using the approach used for the mixture mass equation.
The finite volumes discretization of the mixture internal energy equation involves its integration in time between t n and t n+1 and in space on an elementary control volume M K c (see figure 2): The discrete form of (9) is: where F n+1 σ , G n+1 σ and H n+1 σ are the approximations of the fluxes on the interface σ at time t n+1 : n σ are defined only at the cells center. On the faces, they are approximated using the following donor formulation: σ is the relative velocity which is given by using the Ishii model [9,10]: C 0 is a parameter that adjusts the mixture velocity and − → V v,lim is the vapor velocity limit. The concentration C n σ and the void fraction α n σ on the interface σ are approximated by also using a donor formulation.
Hence, the discrete form of the mixture internal energy equation (10) can be written as: Taking into account the state equation (2) and the relation (11), equation (12) becomes:

Discretization of the vapor mass equation
The spatial and temporal discretization of the vapor mass balance equation is done by using the approach used for the mixture mass equation.
We integrate the vapor mass equation between the time instants t n and t n+1 on the cell M K c (see figure 2): The discrete form of (14) is: where F n+1 σ and G n+1 σ are the approximations of the fluxes on the interface σ at time t n+1 : n σ are defined only at the cells center. On the faces, they are approximated using the following donor formulation: is calculated using a correlation [2].
As a result, the discrete form of the vapor mass equation (15) can be written as: Taking into account the state equation (2) and the relation (11), the equation (16) becomes:

Discretization of the mixture momentum conservation equation
We solve the mixture momentum equation of (1) in the non-conservative form: The equation (18) can be divided by the porosity which then remains only in the In this article, we assume that the porosity is constant, the above equation (18) becomes: Future work will be conducted to raise this hypothesis and discretize the term In the sequel, we detail the spatial and temporal discretization of the equation (20) projected in the direction e x (the discretization in the directions e y and e z , is obtained by analogy).
The projection of equation (20) in the direction − → e x gives: The finite volumes discretization involves the integration of equation (20) in time between t n and t n+1 and in space on an elementary control volume M K+ u (see figure 3): The discrete form is: We now detail each terms in 2D for simplification reasons.
• Approximation of the scalar variables at the interfaces of M K+ u : We use the same method to determine [ρ r C r (1 − C r )] n σ : At the interfaces σ z + and σ z − , we easily get the same result as at the interfaces σ y + and σ y − if we replace y by z.

• Approximation of friction forces
where τ x f w is the wall friction and τ x f s is the singular friction. We obtain: where K x and f x w are the friction coefficients.
Finally, the discretization of the mixture momentum equation in the direction e x (22) can be written as: By analogy, we obtain the discretization in the directions e y and e z :

Construction of the linear system to solve
Let (S) denote the non-linear system we ought to solve at each physical time step: System (S) is solved by using a Newton-Raphson iterative method which consists in solving a linearisation of (S) at each iteration. More precisely, at each iteration k of the Newton-Raphson algorithm, we solve the linear system: and where k is the number of the iteration, the matrix J is the Jacobian matrix of the system (S) and the vector S is the right-hand-side vector containing the residuals of equations (26) evaluated at the previous iteration. The solution U n+1 is U k when k → +∞.
To solve the non-linear system (26), we compare two different methods: • the "Full Jacobian" method (see §3.2), • the "Pressure-based Solver" method (see §3.3). These methods are used by the CATHARE code to solve the six equations model [6]. The Full Jacobian method is used to deal with the 1D problems since the size of the matrix allows it. However, in the 2D and 3D problems CATHARE code uses the pressure-based method.

The "Full Jacobian" method
This method involves the inversion of the matrix in (26) when its size is reasonable. Actually in this case this method is very useful since it is simple to implement and because it allows the possibility of "impliciting" all the variables in the discretization step which leads to a better resolution.
Despite the efficiency of the full Jacobian method to treat the "small sized" problems, its use is bounded by a limit on the jacobian matrix size. Beyond this limit, the pressure-based method becomes more efficient.

The "pressure-based solver" method
The "semi-implicit" scheme we used to discretize the four equations (1) allows to simplify significantly the terms of the matrix that occur in the momentum equations F 4 , F 5 and F 6 . This will enable the expression of the velocity increments as functions of the pressure increments (see section 3.3.1) and then to eliminate them. In section §3.3.2 we will see how to eliminate all the scalar variables (but the pressure increments) and then to obtain a linear equation where only the pressure increments should occur.

Elimination of the velocity increments
The purpose of this step is to write the velocity increments as a function of the pressure increments. To do so we consider only the momentum equations which corresponds to this "partial" linear system: According to (23),(24) and (25), the functions F 4 , F 5 and F 6 do not depend on mixture enthalpy h and mass vapor concentration C at time step t n+1 . Thus, their derivatives with respect to these variables are equals to zero. In the same way: • the derivatives of F 4 with respect to V y and V z are equal to zero, • the derivatives of F 5 with respect to V x and V z are equal to zero, • the derivatives of F 6 with respect to V x and V y are equal to zero. Thus, system (28) becomes: Moreover, by developing equation (29) in the cells M K+ u , M K+ v and M K+ w , we obtain the following equations system: as functions of the pressure increments.
This step enabled the writing of the velocity increments ∆U 4 , ∆U 5 and ∆U 6 at each face of the mesh as functions of the pressure increments. This will be useful in the next step (see §3.3.2) to establish a linear pressure equation (such as only pressure increments are unknown).

Triangulation
In this step we aim to eliminate the scalar variables increments. To do so, we use the three scalar variables equations (mixture mass, mass vapor concentration and mixture energy), which corresponds to the following system: As we did previously, some simplifications take place (thanks to (8) and (13)) since F 1 and F 2 do not depend on C at t n+1 . Thus, system (31) becomes: At this stage applying the operation (L 1 ) ←− (L 1 ) × ∂F 2 ∂h − (L 2 ) × ∂F 1 ∂h on the system (32) results in: where:

Low Mach regime
In this section, we check the efficiency of the numerical method described in Sections 2and3 by carrying out various 1D tests. For all these tests, we used the stiffiend gas law for the equation of state law [5].

Channel with varying porosity
Test description: Stationnary liquid flow in a 1D channel with varying porosity. The channel is 4.2 m long. The porosity in the middle of the channel (1.4 m¡x¡2.8 m) is set to 0.5 and to 1 elsewhere. The pressure at the outlet of the channel is fixed at 155 bar and the liquid velocity at the inlet of the channel is 1 m s −1 . The friction wall is set to zero. The gravity is neglected. Here, the objective is to test the ability of the spatial discretization to conserve momentum in case of a 1D varying porosity.
Results: Figure 4 presents the velocity and the pressure profiles along the z-axis. Test description: Stationnary liquid flow in a 1D channel with charge loss. The channel is 4.2 m long. The singular charge loss and models friction effects due to the presence of an element (like a mixing grid) at the middle of the core 2.1 m. The pressure at the outlet of the channel is fixed at 155 bar and the liquid velocity at the inlet of the channel is 1 m s −1 . The friction wall is set to zero. The gravity is neglected. We take K = 100 as the charge loss coefficient.
Results: Figure 5 shows that the pressure decreases after the flow passed through the mixing grid. Figure 5. Channel with singular charge loss: Pressure-based solver (red), analytic solution (green)

Boiling channel
Test description: The physical quantities that we use in these tests match the functioning of the PWR (Pressurized Water Reactors). We consider a 4.2 m long channel heated by a uniform thermal flux Q =1.E8 W m −3 on which we impose the following conditions: inlet concentration C = 0, inlet enthalpies h l =1300 kJ kg −1 and h v =2600 kJ kg −1 , inlet velocities u l = u v =1 m s −1 and outlet pressure P s =155 bar. The friction wall is set to zero. The gravity is taken into account.
Results: In the figure 6, we compare the results of homogeneous ( − → V r = 0) model (1) obtained with our pressure solver to an analytical solution obtained with the low Mach mixture model proposed in [5].

Compressible regime
To test our Pressure-based Solver in the compressible regim, we study test cases described in [12] and presented in Table 1. All tests are performed with an ideal gases ρ = γ γ−1 P h , with a constant γ = 1.4. All chosen data consist of two constant states separated by a discontinuity at x = x 0 . The position of the discontinuity is stated in the legend. The spatial domain is 0 ≤ x ≤ 1 and the numerical solution is computed with 100 and 5000 cells.
Test ρ L (kg.m −3 ) V L (m.s −1 ) P L (P a) ρ R (kg.m −3 ) V R (m.s −1 ) P R (P a) For test 1, the Sod's Shock Tube problem [12] is modified slightly. The solution of the problem has a right shock wave, a right travelling contact wave and a left sonic rarefaction wave. The purpose of this test is to assess the entropy satisfaction property of the numerical methods. The results for this test are shown in Figure 7 against the exact results. We can see that the Pressure-based Solver struggles slightly with the internal energy, but otherwise performs very close to the exact Riemann solver with particularly accurate result around the sonic point unlike Roe scheme [12].  Table 1: Pressure-based solver (red and blue lines), Roe solver (sky blue line) and exact (dash) solutions compared at time 0.2 Test 2 consists of two symmetric rarefaction waves and a trivial contact wave, with the star region between the non-linear waves close to vacuum. This problem is a good assessment of the performance of numerical methods for low-density flows. The results for this test are shown in Figure 8 against the exact results and we can see that the accuracy of the numerical results is nearly from those of the exact Riemann solver unlike Roe scheme which fails near low-density flows [13].  Table 1: Numerical (red and blue lines) and exact (dash) solutions compared at time 0.12 The robustness and accuracy of the Pressure-based solver is tested with Test 3. The solution of Test 3 consists of a strong shock wave, a contact surface and a left rarefaction wave. We can see from Figures 9 (100 cells) that the Pressure-based solver struggles slightly with the density, but otherwise performs very close to the exact solution. We can also see that the results obtained with the Pressure-based and Roe solvers are near. However, when the mesh is refined (1000 cells), we see that the Pressure-based solver do not converge to the exact solution. Figure 9. Pressure-based Riemann solver applied to test 3 of Table 1: Pressure-based solver (red and blue lines), Roe solver (sky blue line) and exact (dash) solutions compared at time 0.012