Numerical simulation of a barotropic two-phase flow model with miscible phases

Jean Bussac, Khaled Saleh

► To cite this version:

Jean Bussac, Khaled Saleh. Numerical simulation of a barotropic two-phase flow model with miscible phases. International Conference on Finite Volumes for Complex Applications, Oct 2023, Strasbourg, France. 10.1007/978-3-031-40860-1_7. hal-04038701

HAL Id: hal-04038701

https://hal.science/hal-04038701

Submitted on 21 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Numerical simulation of a barotropic two-phase flow model with miscible phases

Jean Bussac¹ and Khaled Saleh²

¹ Nantes Université, CNRS
Laboratoire de Mathématiques Jean Leray, LMJL

44000 Nantes, France jean.bussac@univ-nantes.fr,

² Université de Lyon, CNRS UMR 5209, Université Lyon 1
Institut Camille Jordan
69622 Villeurbanne, France

**Abstract.** This paper addresses the numerical approximation of a compressible barotropic two-phase flow model with miscible conditions. The first phase is a liquid while the second phase corresponds to a gaseous mixture, which contains two components that share the same volume (for example vapor and air). The fluid dynamic is depicted by a Baer-Nunziato-type model, involving non conservative coupling terms. The core of the paper is the simulation of this model by a Suliciu relaxation scheme initially designed for immiscible mixtures. The numerical results illustrate the convergence of the scheme, its robustness for low volume fraction regimes and its computational cost efficiency.

**Keywords:** multiphase flow, relaxation scheme, miscible conditions

## 1 Introduction

This work concerns the simulation of the barotropic three-field two-phase Baer-Nunziato model, studied in [6]. This model corresponds to a mixture of a liquid phase (l) plus a gaseous phase containing two miscible gases: vapor (v) and an inert gas (g). In the nuclear industry, such flows may appear in accidental configurations such as the Loss of Coolant Accident (LOCA) [10]. In the one-dimensional framework, it is a seven equation system of partial differential equations involving non-conservative coupling terms, which reads

$$ 
\begin{align} 
\partial_t \alpha_l + u_l \partial_x \alpha_l &= 0, \nonumber \\
\partial_l (\alpha_l \rho_l) + \partial_x (\alpha_l \rho_l u_l) &= 0, \nonumber \\
\partial_l (\alpha_l \rho_l u_l) + \partial_x (\alpha_l \rho_l u_l^2 + \alpha_l p_l(\rho_l)) + p_v(\rho_v) \partial_x \alpha_v + p_g(\rho_g) \partial_x \alpha_g &= 0, \nonumber \\
\text{and for } k = v, g: \tag{1}
\end{align} 
$$

$$ 
\begin{align*}
\partial_t (\alpha_k \rho_k) + \partial_x (\alpha_k \rho_k u_k) &= 0, \\
\partial_t (\alpha_k \rho_k u_k) + \partial_x (\alpha_k \rho_k u_k^2 + \alpha_k p_k(\rho_k)) - p_k(\rho_k) \partial_x \alpha_k &= 0.
\end{align*} 
$$

For $k \in \{l, v, g\}$, $\rho_k > 0$ and $u_k \in \mathbb{R}$ are respectively the density and velocity of field k. The pressure of field k only depends on the corresponding field density and is given by a barotropic equation of state (see [2, 6]) $p_k: \rho_k \mapsto p_k(\rho_k)$ and we assume that $p'_k(\rho_k) > 0$.
Finally $\alpha_k \in (0, 1)$ is the statistical fraction of field $k$. Since the vapor field and the inert gas share the same gaseous phase we have the saturation constraint

$$ \alpha_v = \alpha_g = 1 - \alpha_l. \qquad (2) $$

**Remark 1.** A mixture of three immiscible phases can be modeled by a similar set of PDEs. Assuming that the fields $v$ and $g$ do not share the same phase, the phase fraction of vapor is computed by the additional transport equation:

$$ \partial_t \alpha_v + u_l \partial_x \alpha_v = 0 $$

and $\alpha_g$ is determined by the following saturation constraint:

$$ \alpha_l + \alpha_v + \alpha_g = 1. \qquad (3) $$

**Remark 2.** Source terms can be added to the model (1) to account for mass and thermodynamical transfers.

**Proposition 1.** Denote $\mathcal{U} := (\alpha_l, \alpha_l \rho_l, \alpha_l \rho_l u_l, \alpha_v \rho_v, \alpha_v \rho_v u_v, \alpha_g \rho_g, \alpha_g \rho_g u_g)$. System (1) is weakly hyperbolic and admits the following real eigenvalues

$$ \sigma_1(\mathcal{U}) = u_l, \quad \sigma_{2,3}(\mathcal{U}) = u_l \pm c_l, \quad \sigma_{4,5}(\mathcal{U}) = u_v \pm c_v, \quad \sigma_{6,7}(\mathcal{U}) = u_g \pm c_g, $$

where $c_k(\rho_k) = \sqrt{p'_k(\rho_k)}$ is the speed of sound of field $k$. All the characteristic fields are genuinely non-linear excepting $\sigma_1$ which is linearly degenerate. The set of corresponding right eigenvectors spans $\mathbb{R}^7$ if and only if

$$ \alpha_l \neq 0, \quad \alpha_v = \alpha_g \neq 0, \quad |u_k - u_l| \neq c_k(\rho_k) \text{ for } k = g, v. \qquad (4) $$

**Remark 3.** $\alpha_k = 0$ is to be understood in the sense $\alpha_k \to 0$. These low fraction regimes are quite important in practical applications. It is already known in immiscible situations that the Suliciu relaxation scheme provides very good results in these regimes compared to standard schemes [4, 8]. The aim of this work is to extend the Suliciu relaxation scheme to the model (1).

**Remark 4.** The loss of hyperbolicity (4) can be due to the presence of a nearly vanishing phase $\alpha_k \to 0$, or to an interaction between the coupling wave $u_l$ and an acoustic wave $u_k \pm c_k$, $k = g, v$. This second resonant configuration is unlikely to happen in the concerned industrial applications (multiphase flows in nuclear reactors), where the flows have strongly subsonic relative velocities, i.e. a relative Mach number much smaller than one:

$$ M_k = \frac{|u_l - u_k|}{c_k(\rho_k)} \ll 1 \text{ for } k = g, v. $$

Let the specific total energy of field $k$ be defined by $E_k = e_k + u_k^2/2$, where the specific internal energy of the field $k$ is such that $e'_k(\rho_k) = p_k(\rho_k)/\rho_k^2$. One can prove that the function $\mathcal{U} \mapsto (\alpha_k \rho_k E_k)(\mathcal{U})$ is (non strictly) convex (see [9]). The entropy weak solutions of system (1) are the weak solutions that satisfy the following inequality in the weak sense:

$$ \partial_t \left( \sum_{k=l,g,v} \alpha_k \rho_k E_k \right) + \partial_x \left( \sum_{k=l,g,v} (\alpha_k \rho_k E_k u_k + \alpha_k p_k(\rho_k) u_k) \right) \le 0. \qquad (5) $$

For smooth solutions of (1), (5) is an equality.
## 2 Suliciu relaxation scheme

The solutions of (1) are formally recovered as the limit of the solutions of the following enlarged system, involving the additional unknowns $\mathcal{T}_k$, in the regime $\epsilon \to 0$:

$$
\begin{align}
& \partial_t \alpha_l + u_l \partial_x \alpha_l = 0, \nonumber \\
& \partial_t (\alpha_l \rho_l) + \partial_x (\alpha_l \rho_l u_l) = 0, \nonumber \\
& \partial_t (\alpha_l \rho_l u_l) + \partial_x (\alpha_l \rho_l u_l^2 + \alpha_l \pi_v (\tau_l, \mathcal{T}_l)) \nonumber \\
& \qquad + \pi_v (\tau_v, \mathcal{T}_v) \partial_x \alpha_v + \pi_g (\tau_g, \mathcal{T}_g) \partial_x \alpha_g = 0, \nonumber \\
& \partial_t (\alpha_l \rho_l \mathcal{T}_l) + \partial_x (\alpha_l \rho_l \mathcal{T}_l u_l) = \frac{1}{\epsilon} \alpha \rho_l (\tau_l - \mathcal{T}_l), \tag{6} \\
& \text{and for } k=v,g: \nonumber \\
& \partial_t (\alpha_k \rho_k) + \partial_x (\alpha_k \rho_k u_k) = 0, \nonumber \\
& \partial_t (\alpha_k \rho_k u_k) + \partial_x (\alpha_k \rho_k u_k^2 + \alpha_k \pi_k (\tau_k, \mathcal{T}_k)) \nonumber \\
& \qquad - \pi_k (\tau_k, \mathcal{T}_k) \partial_x \alpha_k = 0, \nonumber \\
& \partial_t (\alpha_k \rho_k \mathcal{T}_k) + \partial_x (\alpha_k \rho_k \mathcal{T}_k u_k) = \frac{1}{\epsilon} \alpha \rho_k (\tau_k - \mathcal{T}_k). \nonumber
\end{align}
$$

where, denoting $\tau_k = 1/\rho_k$ the specific volume of field k, we have introduced the linearized pressure laws

$$ \pi_k(\tau_k, \mathcal{T}_k) = P_k(\mathcal{T}_k) + a_k^2(\mathcal{T}_k - \tau_k), \quad k = l, g, v, $$

where $\tau \mapsto P_k(\tau) = p_k(\tau^{-1})$ and $a_k \in \mathbb{R}^+$ suitably chosen.

The Suliciu relaxation scheme consists in building an exact Riemann solver for the convective part of the enlarged system (6), which is (relatively) easy to compute because the system has only linearly degenerate fields. Such an exact Riemann solution was first computed in [4, 3] for two-phase flows and has been extended in [8] for immiscible multiphase barotropic flows. In [8], the transport equations for $\alpha_l$ and $\alpha_v$ are solved and by the relation (3) $\alpha_g$ is computed by enforcing the following initialization:

$$ \alpha_g(t=0) := 1 - \alpha_l(t=0) - \alpha_v(t=0). $$

The extension to the barotropic three-field two-phase Baer-Nunziato model (1) is rather straightforward. It consists in solving the transport equation for $\alpha_l$ while $\alpha_g = \alpha_v$ are computed by enforcing the initialization:

$$ \alpha_g(t=0) = \alpha_v(t=0) := 1 - \alpha_l(t=0). $$

Since in both the miscible and immiscible cases, the phase fractions $\alpha_k$ satisfy the same transport equation:

$$ \partial_t \alpha_k + u_l \partial_x \alpha_k = 0, $$

the expression of the overall Riemann solution is unchanged.

At the numerical level, a fractional step method is used. The first step is a time-advancing step using the exact solution of the Riemann problem for the convective part of (6). The second step consists in an instantaneous relaxation towards the equilibrium system (1) by imposing $\mathcal{T}_k = \tau_k$ in the solution obtained by the first step. The whole scheme can be written in the form of a non conservation finite volume scheme (see [8] for details).
**Proposition 2 (Properties of the Suliciu relaxation scheme).**

- Under a CFL condition [8, (51)], the scheme preserves positive values of the phase fractions and densities.

- The discretizations of the partial masses and total mixture momentum are conservative.

- Under Whitham's condition on the parameters $a_k$ (see [8, (16)], [1]), the numerical solution satisfies a discrete counterpart to the energy inequality (5).

## 3 Numerical results

The relaxation scheme is compared to the non conservative Rusanov scheme [5], which is (to our knowledge) the only available scheme to this multiphase flow model. In the following we denote $U = (\alpha_l, \rho_l, u_l, \rho_v, u_v, \rho_g, u_g)$ the vector of non conservative unknowns. We present two test-cases and plot the approximated solutions for a 100 cell mesh, which corresponds to a realistic $10^6$ cell mesh in a three-dimensional industrial context. The CFL number is set to 0.45.

### 3.1 Test-case 1: pure contact discontinuity

The initial data is the following

$$ 
\begin{aligned}
U_L &= (0.1, 3.0, 0.2, 0.3, 1.0, 0.4, 1.0), \\
U_R &= (0.4, 2.0724212862957074, 0.2, 0.28339834471128433, \\
&\qquad 1.4702967632600488, 0.34584269153214842, 1.5879142504747157).
\end{aligned} 
\quad (7) 
$$

We present the approximated solution in Figure 1. The $L^1$-errors with respect to the space step $\Delta x$ and with respect to the computational cost are plotted on Figure 2. The classical sonic wave oscillations are smaller for the relaxation scheme. The contact discontinuity is much less smeared with the relaxation scheme than with Rusanov's scheme which suggests better accuracy. Indeed, Figure 2 shows that it requires a hundred times more CPU time to obtain the same error with Rusanov's scheme (using more refined meshes) than with the relaxation scheme.

### 3.2 Test-case 2: vanishing-phase test-case with $u_v \pm c_v$ shocks

In this test-case, the left initial data has no liquid phase, which corresponds to $\alpha_l = 0$ (in practice $\alpha_l = 10^{-10}$). We consider a shock-contact-shock wave configuration for the vapor field, while the other fields are only subject to the $u_l$ contact. The wave structure is given in Figure 3. The initial data is given by:

$$ 
\begin{aligned}
U_L &= (1.0 \cdot 10^{-10}, 1.0, 0.3, 0.8, 0.5, 0.6, 0.5), \\
&\quad -0.66086132955857257, 0.54048985831522922, 1.4101040855979630), \\
U_R &= (0.8, 2.1958716914805883, 0.3, 0.9, \\
&\quad -0.49798559025814237, 0.54048985831522922, 1.4101040855979630).
\end{aligned} 
\quad (8) 
$$
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Test-case 1: solution for 100 cells at final time $t = 0.08$. Contact wave of velocity $u_l = 0.2$. Initial data given in (7).</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Test-case 1: $L^1$-error with respect to $\Delta x$ and $L^1$-error with respect to computational cost (in seconds). The $u_g$ plot is not represented since it is very similar to the $u_v$ plot.</figcaption></figure>
The approximated solutions on Figure 3 give an interesting comparison. On the one hand, the relaxation scheme follows the contact discontinuity much better than the Rusanov scheme and does not oscillate on the left constant state, which is a low fraction regime. On the other hand, close intermediate states of the vapor field are rather well captured by the relaxation scheme, even for a 100 cell mesh. On the contrary, observe that the Rusanov scheme misses the $U_2$ state.

## 4 Conclusion

This work confirms the convincing behavior of this relaxation scheme already known for immiscible mixtures. For a given error on the solution, it performs much better than the Rusanov scheme (a factor of order $10^2$), which is a discriminating feature in an industrial approach. Moreover, the low fraction regimes are approximated with robustness. In future works, it would be interesting to consider the thermodynamical relaxation by adding source terms to the convective system, in order to use the relaxation scheme in a global context. The three-phase four-field case [7] can be investigated.

## References

1. Bouchut, F.: Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources. Springer Science & Business Media (2004)

2. Coquel F., Gallouët T., Hérard J.-M., Seguin N.: Closure laws for a two-fluid two-pressure model. Comptes Rendus. Mathématique, Académie des sciences (Paris), 334, 10, 927-932 (2002). https://hal.archives-ouvertes.fr/hal-01484345/file/CRAS2PHI.pdf

3. Coquel F., Hérard J.-M., Saleh K.: A positive and entropy-satisfying finite volume scheme for the Baer-Nunziato model. Journal of Computational Physics 330, 401-435 (2017) Elsevier. https://hal-polytechnique.archives-ouvertes.fr/hal-01261458v4

4. Coquel F., Hérard J.-M., Saleh K., Seguin N.: A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model. ESAIM: Mathematical Modelling and Numerical Analysis 48, 165-206 (2014) EDP Sciences. https://hal.archives-ouvertes.fr/hal-00795568v1

5. Gallouët T., Hérard, J.-M., Seguin N.: Numerical modeling of two-phase flows using the two-fluid two-pressure approach. Mathematical Models and Methods in Applied Sciences 14, 05, 663-700 (2004) World Scientific. https://hal.archives-ouvertes.fr/hal-00003327/document

6. Hérard J.-M., Mathis H.: A three-phase flow model with two miscible phases. ESAIM: Mathematical Modelling and Numerical Analysis 53, 4, 1373-1389 (2019) EDP Sciences. https://hal.archives-ouvertes.fr/hal-01976938

7. Hérard J.-M., Hurisse O. Quibel L.: A four-field three-phase flow model with both miscible and immiscible components. ESAIM: Mathematical Modelling and Numerical Analysis 55, S251-S278 (2021). https://hal.archives-ouvertes.fr/hal-02432793

8. Saleh K.: A relaxation scheme for a hyperbolic multiphase flow model. Part I: barotropic eos. ESAIM: Mathematical Modelling and Numerical Analysis, Apr (2019) EDP Sciences. https://hal.archives-ouvertes.fr/hal-01737681v2/file/Triph_Saleh_V3.pdf

9. Saleh K., Seguin N.: Some mathematical properties of a barotropic multiphase flow model. ESAIM Proc. Surveys 69, 70-78 (2020) EDP Sciences. https://doi.org/10.1051/proc/202069070

10. U.S. NRC: Glossary: Loss of Coolant Accident (LOCA). http://www.nrc.gov/reading-rm/basic-ref/glossary/loss-of-coolant-accident-loca.html
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Test-case 2: solutions for 100 cells at final time $t = 0.02$. Contact wave of velocity $u_l = 0.3$, $u_v - c_v$-shock of velocity -3.3155813734376762 and $u_v + c_v$-shock of velocity 3.0689467170729379. Initial data given in (8).</figcaption></figure>