# Maximum principle for the mass fraction in a system with two mass balance equations

Gauthier Lazare, Qingqing Feng, Philippe Helluy, Jean-Marc Hérard, Frank

Hülsemann, Stéphane Pujet

► To cite this version:

Gauthier Lazare, Qingqing Feng, Philippe Helluy, Jean-Marc Hérard, Frank Hülsemann, et al. Maximum principle for the mass fraction in a system with two mass balance equations. Comptes Rendus. Mécanique, 2024, 352, pp.81-98. 10.5802/crmeca.244. hal-04479205

HAL Id: hal-04479205

https://hal.science/hal-04479205

Submitted on 4 Apr 2024

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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Maximum principle for the mass fraction in a system with two mass balance equations

Gauthier Lazare, Qingqing Feng, Philippe Helluy, Jean-Marc Hérard, Frank

Hülsemann, Stéphane Pujet

► To cite this version:

Gauthier Lazare, Qingqing Feng, Philippe Helluy, Jean-Marc Hérard, Frank Hülsemann, et al. Maximum principle for the mass fraction in a system with two mass balance equations. Comptes Rendus. Mécanique, 2024, 352, pp.81-98. 10.5802/crmeca.244. hal-04479205

HAL Id: hal-04479205

https://hal.science/hal-04479205v1

Submitted on 4 Apr 2024

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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Maximum principle for the mass fraction in a system with two mass balance equations

*Principe du maximum pour le titre massique d'un système incluant deux bilans de masse.*

Gauthier Lazare<sup>® a, b</sup>, Qingqing Feng<sup>a</sup>, Philippe Helluy<sup>® b</sup>, Jean-Marc Hérard<sup>® a, c</sup>, Frank Hulsemann<sup>® a</sup> and Stéphane Pujet<sup>a</sup>

**Abstract.** Three Finite Volume schemes are proposed in this note to satisfy the maximum principle for the mass fraction y, solution of an unsteady balance equation, including a relative velocity between phases and a source term. The continuous maximum principle is examined first. Then, linear implicit discrete schemes are detailed in a multi-dimensional and unstructured framework.

**Résumé.** Dans cette note, trois schémas Volumes Finis sont proposés pour respecter le principe du maximum du titre massique, solution d'une équation de bilan instationnaire, incluant un déséquilibre en vitesse avec une vitesse relative non nulle et un terme source. Le principe du maximum continu est d'abord étudié puis les schémas discrets linéaires implicites sont détaillés dans un cadre multi-dimensionnel non structuré.

**Keywords.** maximum principle, Finite Volume scheme, two-phase flow, non-equilibrium velocity.

**Mots-clés.** Principe du maximum, schéma Volumes Finis, écoulement diphasique, déséquilibre en vitesse.

**Funding.** The first author is supported by EDF R&D and ANRT through CIFRE contract 2022/1027.

This article is a draft (not yet accepted!)

<footnote><sup>a</sup> EDF R&D Chatou - 6 quai Waltier, 78400, Chatou</footnote>

<footnote><sup>b</sup> IRMA, UMR 7501, 7 rue Descartes, 67000 Strasbourg</footnote>

<footnote><sup>c</sup> I2M - Institut de Mathématiques de Marseille, Aix Marseille Université, France</footnote>

<footnote>*E-mails:* gauthier.lazare@edf.fr (G. Lazare), qingqing.feng@edf.fr (Q. Feng), philippe.helluy@unistra.fr (P. Helluy), jean-marc.herard@edf.fr (J-M. Herard), frank.hulsemann@edf.fr (F. Hulsemann), stephane.pujet@edf.fr (S. Pujet)</footnote>
<figure><img src="image_2.png" /></figure>

# Maximum principle for the mass fraction in a system with two mass balance equations

*Principe du maximum pour le titre massique d'un système incluant deux bilans de masse.*

Gauthier Lazare<sup>® a, b</sup>, Qingqing Feng<sup>a</sup>, Philippe Helluy<sup>® b</sup>, Jean-Marc Hérard<sup>® a, c</sup>, Frank Hulsemann<sup>® a</sup> and Stéphane Pujet<sup>a</sup>

**Abstract.** Three Finite Volume schemes are proposed in this note to satisfy the maximum principle for the mass fraction y, solution of an unsteady balance equation, including a relative velocity between phases and a source term. The continuous maximum principle is examined first. Then, linear implicit discrete schemes are detailed in a multi-dimensional and unstructured framework.

**Résumé.** Dans cette note, trois schémas Volumes Finis sont proposés pour respecter le principe du maximum du titre massique, solution d'une équation de bilan instationnaire, incluant un déséquilibre en vitesse avec une vitesse relative non nulle et un terme source. Le principe du maximum continu est d'abord étudié puis les schémas discrets linéaires implicites sont détaillés dans un cadre multi-dimensionnel non structuré.

**Keywords.** maximum principle, Finite Volume scheme, two-phase flow, non-equilibrium velocity.

**Mots-clés.** Principe du maximum, schéma Volumes Finis, écoulement diphasique, déséquilibre en vitesse.

**Funding.** The first author is supported by EDF R&D and ANRT through CIFRE contract 2022/1027.

This article is a draft (not yet accepted!)

<footnote><sup>a</sup> EDF R&D Chatou - 6 quai Waltier, 78400, Chatou</footnote>

<footnote><sup>b</sup> IRMA, UMR 7501, 7 rue Descartes, 67000 Strasbourg</footnote>

<footnote><sup>c</sup> I2M - Institut de Mathématiques de Marseille, Aix Marseille Université, France</footnote>

<footnote>*E-mails:* gauthier.lazare@edf.fr (G. Lazare), qingqing.feng@edf.fr (Q. Feng), philippe.helluy@unistra.fr (P. Helluy), jean-marc.herard@edf.fr (J-M. Herard), frank.hulsemann@edf.fr (F. Hulsemann), stephane.pujet@edf.fr (S. Pujet)</footnote>
# 1. Introduction

We consider in this note an unsteady mass balance equation for the mass fraction $y$ in a two-phase flow model, including a relative velocity between phases, such as the drift-flux model [1]. The maximum principle for the mass fraction is widely investigated, either for diffusive problems (see among others [2]) or in hybrid convection-diffusion problems (see among others [3-7]). Most of the time, a null relative velocity is considered in the convective flux model. We focus in the sequel on the system involving two mass balance equations with a relative velocity (see also [4, 5]) and a source term.

Once the continuous maximum principle has been examined in section 2, three distinct linear Finite Volume schemes complying with the discrete maximum principle with no (or with a weak) restriction on the time step are proposed in section 3. Eventually, numerical simulations are used to assess the accuracy of the schemes in section 4. The proposed schemes can be used in a broader framework such as two-phase flow models.

# 2. Governing equations for the model

## 2.1. System of equations

A flow characterized by a density $\rho$ and a flow rate $\mathbf{q}$ is studied on a domain $\Omega$, over a time period $[0, T]$. The boundary of $\Omega$ is noted $\Gamma$. The outward unit normal is noted $\mathbf{n}_\Gamma$. Then

$$ \Gamma = \Gamma_+ \cup \Gamma_- \cup \Gamma_w, \quad (1) $$

where:

$$ \begin{aligned} \Gamma_+ &= \{\mathbf{x} \in \Gamma, \mathbf{q}(\mathbf{x}, t) \cdot \mathbf{n}_\Gamma < 0\}, \\ \Gamma_- &= \{\mathbf{x} \in \Gamma, \mathbf{q}(\mathbf{x}, t) \cdot \mathbf{n}_\Gamma > 0\}, \\ \Gamma_w &= \{\mathbf{x} \in \Gamma, \mathbf{q}(\mathbf{x}, t) \cdot \mathbf{n}_\Gamma = 0\}. \end{aligned} \quad (2) $$

This study focuses on the evolution of the mass fraction $y$ of a species in this flow. The species evolves according to a flow rate $\mathbf{q}_g = \mathbf{q} + (1 - y)\mathbf{q}_r$, where the relative velocity between phases $\mathbf{u}_r$ [1] is used in the relative flow rate $\mathbf{q}_r = \rho \mathbf{u}_r$. The mixture velocity is defined by $\mathbf{u} := \mathbf{q}/\rho$. The mass fraction $y$ tends to return to the equilibrium $\bar{y}$ after a characteristic relaxation time $\tau$. With this modeling, the quantity $\rho y$ follows a conservation law, which is associated with the total mass conservation law. These equations are provided below (see also [5] for a similar system).

$$ \begin{aligned} & \frac{\partial \rho}{\partial t} + \underline{\nabla} \cdot \mathbf{q} = 0, \\ & \frac{\partial (\rho y)}{\partial t} + \underline{\nabla} \cdot (y\mathbf{q}) + \underline{\nabla} \cdot (y(1-y)\mathbf{q}_r) = \rho \frac{\bar{y} - y}{\tau}. \end{aligned} \quad (3) $$

Given the definition of the relative flow rate, the boundary $\Gamma$ can also be split according to the sign of the flow rate $\mathbf{q}_\xi = \mathbf{q} + (1 - 2y)\mathbf{q}_r$:

$$ \Gamma = \Gamma_+^\xi \cup \Gamma_-^\xi \cup \Gamma_w^\xi, \quad (4) $$

where:

$$ \begin{aligned} \Gamma_+^\xi &= \{\mathbf{x} \in \Gamma, \mathbf{q}_\xi \cdot \mathbf{n}_\Gamma < 0\}, \\ \Gamma_-^\xi &= \{\mathbf{x} \in \Gamma, \mathbf{q}_\xi \cdot \mathbf{n}_\Gamma > 0\}, \\ \Gamma_w^\xi &= \{\mathbf{x} \in \Gamma, \mathbf{q}_\xi \cdot \mathbf{n}_\Gamma = 0\}. \end{aligned} \quad (5) $$
# 1. Introduction

We consider in this note an unsteady mass balance equation for the mass fraction $y$ in a two-phase flow model, including a relative velocity between phases, such as the drift-flux model [1]. The maximum principle for the mass fraction is widely investigated, either for diffusive problems (see among others [2]) or in hybrid convection-diffusion problems (see among others [3–7]). Most of the time, a null relative velocity is considered in the convective flux model. We focus in the sequel on the system involving two mass balance equations with a relative velocity (see also [4, 5]) and a source term.

Once the continuous maximum principle has been examined in section 2, three distinct linear Finite Volume schemes complying with the discrete maximum principle with no (or with a weak) restriction on the time step are proposed in section 3. Eventually, numerical simulations are used to assess the accuracy of the schemes in section 4. The proposed schemes can be used in a broader framework such as two-phase flow models.

# 2. Governing equations for the model

## 2.1. System of equations

A flow characterized by a density $\rho$ and a flow rate $\mathbf{q}$ is studied on a domain $\Omega$, over a time period $[0, T]$. The boundary of $\Omega$ is noted $\Gamma$. The outward unit normal is noted $\mathbf{n}_{\Gamma}$. Then

$$ \Gamma = \Gamma_+ \cup \Gamma_- \cup \Gamma_w, \quad (1) $$

where:

$$ \begin{aligned} \Gamma_+ &= \{\mathbf{x} \in \Gamma, \mathbf{q}(\mathbf{x}, t) \cdot \mathbf{n}_{\Gamma} < 0\}, \\ \Gamma_- &= \{\mathbf{x} \in \Gamma, \mathbf{q}(\mathbf{x}, t) \cdot \mathbf{n}_{\Gamma} > 0\}, \\ \Gamma_w &= \{\mathbf{x} \in \Gamma, \mathbf{q}(\mathbf{x}, t) \cdot \mathbf{n}_{\Gamma} = 0\}. \end{aligned} \quad (2) $$

This study focuses on the evolution of the mass fraction $y$ of a species in this flow. The species evolves according to a flow rate $\mathbf{q}_g = \mathbf{q} + (1 - y)\mathbf{q}_r$, where the relative velocity between phases $\mathbf{u}_r$ [1] is used in the relative flow rate $\mathbf{q}_r = \rho \mathbf{u}_r$. The mixture velocity is defined by $\mathbf{u} := \mathbf{q}/\rho$. The mass fraction $y$ tends to return to the equilibrium $\bar{y}$ after a characteristic relaxation time $\tau$. With this modeling, the quantity $\rho y$ follows a conservation law, which is associated with the total mass conservation law. These equations are provided below (see also [5] for a similar system).

$$ \begin{aligned} & \frac{\partial \rho}{\partial t} + \underline{\nabla} \cdot \mathbf{q} = 0, \\ & \frac{\partial (\rho y)}{\partial t} + \underline{\nabla} \cdot (y\mathbf{q}) + \underline{\nabla} \cdot (y(1-y)\mathbf{q}_r) = \rho \frac{\bar{y} - y}{\tau}. \end{aligned} \quad (3) $$

Given the definition of the relative flow rate, the boundary $\Gamma$ can also be split according to the sign of the flow rate $\mathbf{q}_{\xi} = \mathbf{q} + (1 - 2y)\mathbf{q}_r$:

$$ \Gamma = \Gamma_+^\xi \cup \Gamma_-^\xi \cup \Gamma_w^\xi, \quad (4) $$

where:

$$ \begin{aligned} \Gamma_+^\xi &= \{\mathbf{x} \in \Gamma, \mathbf{q}_\xi \cdot \mathbf{n}_\Gamma < 0\}, \\ \Gamma_-^\xi &= \{\mathbf{x} \in \Gamma, \mathbf{q}_\xi \cdot \mathbf{n}_\Gamma > 0\}, \\ \Gamma_w^\xi &= \{\mathbf{x} \in \Gamma, \mathbf{q}_\xi \cdot \mathbf{n}_\Gamma = 0\}. \end{aligned} \quad (5) $$
The four parameters $\mathbf{q}, \mathbf{q}_r, \bar{y}$ and $\tau$ are given functions :

* $\mathbf{q}(\mathbf{x}, t)$ : the mixture flow rate, given by a momentum conservation equation or provided as input data,

* $\mathbf{q}_r(\mathbf{x}, t)$ : the relative flow rate, provided as input data. It is often given with a drift-flux closure law related to the other parameters,

* $\bar{y}(\mathbf{x}, t) \in [0, 1]$ : the equilibrium mass fraction, reached after a characteristic relaxation time $\tau(\mathbf{x}, t) > 0$. $\bar{y}$ and $\tau$ are obtained through closure laws.

This system of equations must be completed by initial conditions for the mass fraction and the density: $y(\mathbf{x}, t = 0) = y_0(\mathbf{x})$, $\rho(\mathbf{x}, t = 0) = \rho_0(\mathbf{x})$, and by suitable boundary conditions for the mass fraction (and also for the density on $\Gamma_+$).

## 2.2. Continuous maximum principle

**Property 1 (Continuous maximum principle)**: For the density $\rho(\mathbf{x}, t) > 0$, assume initial condition such that $y(\mathbf{x}, t = 0) \in [0, 1]$ and boundary conditions such that $y(\mathbf{x} \in \Gamma_+, t) \in [0, 1]$, with $\Gamma_+^\xi = \{\mathbf{x} \in \Gamma, (\mathbf{q} + (1-2y)\mathbf{q}_r) \cdot \mathbf{n}_\Gamma < 0\}$. Consider closure laws for source terms such that $\bar{y}(\mathbf{x}, t) \in [0, 1]$ and $\tau(\mathbf{x}, t) > 0$. If the quantity $((1-2y)\bar{\nabla} \cdot \mathbf{q}_r + 2\mathbf{q}_r \cdot \bar{\nabla}y)/\rho$ is bounded on $\Omega \times [0, T]$, then the mass fraction $y(\mathbf{x}, t)$ solution of (3) lies in $[0, 1]$ on $\Omega \times [0, T]$.

<figure><img src="image_3.png" /></figure>

Using the same methodology as [5] and [8], the proof of this property is given in **Appendix A**. We recall that the positivity of the density is ensured as long as the divergence of the mixture velocity $\mathbf{u}$ is bounded and that the incoming flow has a positive density $\rho(\mathbf{x} \in \Gamma_+, t) > 0$. In the following, the focus shifts on the mass fraction governing equation with unknown $y$. A discrete Finite Volume scheme for the density $\rho$ will be assumed (see equation (9)).

**Remark:** In order to ensure the continuous maximum principle for the mass fraction $y$ and the positivity of the density $\rho$, two different input boundaries $\Gamma_+$ and $\Gamma_+^\xi$ are considered in the most general case. For our numerical applications, we will consider co-current flows entering the domain, hence the phase velocities have the same sign, such that

$$ \mathbf{u}_g \cdot \mathbf{u}_l \ge 0 \quad \text{with} \quad \begin{cases} \mathbf{u}_g = \mathbf{u} + (1-y)\mathbf{u}_r \\ \mathbf{u}_l = \mathbf{u} - y\mathbf{u}_r \end{cases}. \qquad (6) $$

Hence, noting that

$$ \begin{aligned} \mathbf{q} &= \rho(y\mathbf{u}_g + (1-y)\mathbf{u}_l), \\ \mathbf{q}_\xi &= \rho((1-y)\mathbf{u}_g + y\mathbf{u}_l), \end{aligned} \qquad (7) $$

we get at once $\mathbf{q}_\xi \cdot \mathbf{q} = \rho^2(y(1-y)(\mathbf{u}_g - \mathbf{u}_l)^2 + \mathbf{u}_g \cdot \mathbf{u}_l) \ge 0$ but also $\Gamma_{lw}^\xi = \Gamma_w$, $\Gamma_+^\xi = \Gamma_+$ and $\Gamma_-^\xi = \Gamma_-$.
The four parameters $\mathbf{q}, \mathbf{q}_r, \bar{y}$ and $\tau$ are given functions :

* $\mathbf{q}(\mathbf{x}, t)$ : the mixture flow rate, given by a momentum conservation equation or provided as input data,

* $\mathbf{q}_r(\mathbf{x}, t)$ : the relative flow rate, provided as input data. It is often given with a drift-flux closure law related to the other parameters,

* $\bar{y}(\mathbf{x}, t) \in [0, 1]$ : the equilibrium mass fraction, reached after a characteristic relaxation time $\tau(\mathbf{x}, t) > 0$. $\bar{y}$ and $\tau$ are obtained through closure laws.

This system of equations must be completed by initial conditions for the mass fraction and the density: $y(\mathbf{x}, t = 0) = y_0(\mathbf{x})$, $\rho(\mathbf{x}, t = 0) = \rho_0(\mathbf{x})$, and by suitable boundary conditions for the mass fraction (and also for the density on $\Gamma_+$).

## 2.2. Continuous maximum principle

**Property 1 (Continuous maximum principle)**: For the density $\rho(\mathbf{x}, t) > 0$, assume initial condition such that $y(\mathbf{x}, t = 0) \in [0, 1]$ and boundary conditions such that $y(\mathbf{x} \in \Gamma_+, t) \in [0, 1]$, with $\Gamma_+^\xi = \{\mathbf{x} \in \Gamma, (\mathbf{q} + (1-2y)\mathbf{q}_r) \cdot \mathbf{n}_\Gamma < 0\}$. Consider closure laws for source terms such that $\bar{y}(\mathbf{x}, t) \in [0, 1]$ and $\tau(\mathbf{x}, t) > 0$. If the quantity $((1-2y)\bar{\nabla} \cdot \mathbf{q}_r + 2\mathbf{q}_r \cdot \bar{\nabla}y)/\rho$ is bounded on $\Omega \times [0, T]$, then the mass fraction $y(\mathbf{x}, t)$ solution of (3) lies in $[0, 1]$ on $\Omega \times [0, T]$.

<figure><img src="image_3.png" /></figure>

Using the same methodology as [5] and [8], the proof of this property is given in **Appendix A**. We recall that the positivity of the density is ensured as long as the divergence of the mixture velocity $\mathbf{u}$ is bounded and that the incoming flow has a positive density $\rho(\mathbf{x} \in \Gamma_+, t) > 0$. In the following, the focus shifts on the mass fraction governing equation with unknown $y$. A discrete Finite Volume scheme for the density $\rho$ will be assumed (see equation (9)).

**Remark:** In order to ensure the continuous maximum principle for the mass fraction $y$ and the positivity of the density $\rho$, two different input boundaries $\Gamma_+$ and $\Gamma_+^\xi$ are considered in the most general case. For our numerical applications, we will consider co-current flows entering the domain, hence the phase velocities have the same sign, such that

$$ \mathbf{u}_g \cdot \mathbf{u}_l \ge 0 \quad \text{with} \quad \begin{cases} \mathbf{u}_g = \mathbf{u} + (1-y)\mathbf{u}_r \\ \mathbf{u}_l = \mathbf{u} - y\mathbf{u}_r \end{cases}. \qquad (6) $$

Hence, noting that

$$ \begin{aligned} \mathbf{q} &= \rho(y\mathbf{u}_g + (1-y)\mathbf{u}_l), \\ \mathbf{q}_\xi &= \rho((1-y)\mathbf{u}_g + y\mathbf{u}_l), \end{aligned} \qquad (7) $$

we get at once $\mathbf{q}_\xi \cdot \mathbf{q} = \rho^2(y(1-y)(\mathbf{u}_g - \mathbf{u}_l)^2 + \mathbf{u}_g \cdot \mathbf{u}_l) \ge 0$ but also $\Gamma_{lw}^\xi = \Gamma_w$, $\Gamma_+^\xi = \Gamma_+$ and $\Gamma_-^\xi = \Gamma_-$.
### 3. Numerical method

The domain $\Omega$ is discretised in $N$ cells: $\Omega = \cup_{i \in [1,N]} \Omega_i$. The time interval $[0, T]$ is also discretised in $N_t$ intervals so that:

$$t^0 = 0; \quad \forall n \in [0, N_t], \quad t^{n+1} = t^n + \Delta t^n \quad \text{and} \quad T = \sum_{i \in [0, N_t]} \Delta t^n. \qquad (8)$$

The continuous equations are integrated on each cell $i$ between time $t^n$ (superscript $n$) and time $t^{n+1}$ (superscript $n+1$). We denote with a superscript $*$ data taken at an intermediate state $t^*$, which is either $t^n$ or $t^{n+1}$. The explicit or implicit choice will be decided according to the wanted properties of the scheme. The following notations are used:

* $\theta_i^n$ the approximate value of $\theta$ on cell $i$ at time step $n$,

* $\omega_i$ the volume of cell $i$,

* $j \in v(i)$ the neighbors of cell $i$,

* $S_{ij}$ the surface of intersection of cell $i$ with its neighbor cell $j$,

* $\mathbf{n}_{ij}$ the normal unit vector of the surface $S_{ij}$ outward cell $i$,

* $q_{ij}^* = \mathbf{q}_{ij}^* \cdot \mathbf{n}_{ij}$ the normal mixture flow rate between cell $i$ and cell $j$,

* $(q_r)_{ij}^* = (q_r)_{ij}^* \cdot \mathbf{n}_{ij}$ the normal relative flow rate between cell $i$ and cell $j$,

* $(q_g)_{ij}^* = q_{ij}^* + (1 - y_{ij}^n)(q_r)_{ij}^*$ the normal gas flow rate between cell $i$ and cell $j$. The explicit choice for $y_{ij}^n$ will result in a linear numerical scheme. The spatial discretisation for the mass fraction $y_{ij}^n$ at interface $ij$ is detailed in the sequel.

Here, $q_{ij}^*$ (respectively $(q_r)_{ij}^*$) is an estimation of the flow rate $\mathbf{q}$ (respectively $\mathbf{q}_r$) at interface $ij$ at time $t^* \in [t^n, t^{n+1}]$.

#### 3.1. Finite Volume Discretisation

Finite Volume methods are well adapted to treat conservation issues and to keep physical quantities in valid bounds. Focusing on the system of mass balances (3), we recall that paper [5] introduces some **non linear implicit** scheme (with respect to $y$), which guarantees the discrete maximum principle with no restriction on the time step, even if some non-zero relative velocity is accounted for. Here, we propose first a **linear implicit** scheme with some restriction on the time step in order to comply with the discrete maximum principle. Then two other **linear implicit** schemes satisfying the maximum principle without any condition on the time step are presented.

The Finite Volume scheme for the total mass conservation is assumed to be:

$$\omega_i(\rho_i^{n+1} - \rho_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* = 0. \qquad (9)$$

Turning to the mass fraction equation, the following choices have been made:

* The unsteady term is decomposed in two parts:

$$\begin{align}
\Delta t^n \int_{\omega_i} \frac{\partial \rho y}{\partial t} d\Omega &= ((\rho y)_i^{n+1} - (\rho y)_i^n)\omega_i \nonumber \\
&= \rho_i^n (y_i^{n+1} - y_i^n)\omega_i + y_i^{n+1} \underbrace{(\rho_i^{n+1} - \rho_i^n)\omega_i}_{\text{using (9)}}, \tag{10} \\
&= -\Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \nonumber
\end{align}$$
### 3. Numerical method

The domain $\Omega$ is discretised in $N$ cells: $\Omega = \cup_{i \in [1, N]} \Omega_i$. The time interval $[0, T]$ is also discretised in $N_t$ intervals so that:

$$t^0 = 0; \quad \forall n \in [0, N_t], \quad t^{n+1} = t^n + \Delta t^n \quad \text{and} \quad T = \sum_{i \in [0, N_t]} \Delta t^n. \qquad (8)$$

The continuous equations are integrated on each cell $i$ between time $t^n$ (superscript $n$) and time $t^{n+1}$ (superscript $n+1$). We denote with a superscript $*$ data taken at an intermediate state $t^*$, which is either $t^n$ or $t^{n+1}$. The explicit or implicit choice will be decided according to the wanted properties of the scheme. The following notations are used:

* $\theta_i^n$ the approximate value of $\theta$ on cell $i$ at time step $n$,

* $\omega_i$ the volume of cell $i$,

* $j \in v(i)$ the neighbors of cell $i$,

* $S_{ij}$ the surface of intersection of cell $i$ with its neighbor cell $j$,

* $\mathbf{n}_{ij}$ the normal unit vector of the surface $S_{ij}$ outward cell $i$,

* $q_{ij}^* = \mathbf{q}_{ij}^* \cdot \mathbf{n}_{ij}$ the normal mixture flow rate between cell $i$ and cell $j$,

* $(q_r)_{ij}^* = (\mathbf{q}_r)_{ij}^* \cdot \mathbf{n}_{ij}$ the normal relative flow rate between cell $i$ and cell $j$,

* $(q_g)_{ij}^* = q_{ij}^* + (1 - y_{ij}^n)(q_r)_{ij}^*$ the normal gas flow rate between cell $i$ and cell $j$. The explicit choice for $y_{ij}^n$ will result in a linear numerical scheme. The spatial discretisation for the mass fraction $y_{ij}^n$ at interface $ij$ is detailed in the sequel.

Here, $q_{ij}^*$ (respectively $(q_r)_{ij}^*$) is an estimation of the flow rate $\mathbf{q}$ (respectively $\mathbf{q}_r$) at interface $ij$ at time $t^* \in [t^n, t^{n+1}]$.

#### 3.1. Finite Volume Discretisation

Finite Volume methods are well adapted to treat conservation issues and to keep physical quantities in valid bounds. Focusing on the system of mass balances (3), we recall that paper [5] introduces some **non linear implicit** scheme (with respect to $y$), which guarantees the discrete maximum principle with no restriction on the time step, even if some non-zero relative velocity is accounted for. Here, we propose first a **linear implicit** scheme with some restriction on the time step in order to comply with the discrete maximum principle. Then two other **linear implicit** schemes satisfying the maximum principle without any condition on the time step are presented.

The Finite Volume scheme for the total mass conservation is assumed to be:

$$\omega_i(\rho_i^{n+1} - \rho_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* = 0. \qquad (9)$$

Turning to the mass fraction equation, the following choices have been made:

* The unsteady term is decomposed in two parts:

$$\begin{align}
\Delta t^n \int_{\omega_i} \frac{\partial \rho y}{\partial t} d\Omega &= (\rho y_i)^{n+1} - (\rho y_i)^n \omega_i \nonumber \\
&= \rho_i^n (y_i^{n+1} - y_i^n) \omega_i + y_i^{n+1} \underbrace{(\rho_i^{n+1} - \rho_i^n)\omega_i}_{\text{using (9)}}, \tag{10} \\
&= -\Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \nonumber
\end{align}$$
• The mass fraction is taken implicit in the source term

$$ \int_{\omega_i} \rho \frac{\bar{y} - y}{\tau} d\Omega \approx \rho_i^{\#} \frac{\bar{y}_i^{\#} - y_i^{n+1}}{\tau_i^{\#}} \omega_i \quad (11) $$

where $t^{\#} \in [t^n, t^{n+1}]$ such that:

$$ \begin{cases} \rho_i^{\#} > 0, \\ \bar{y}_i^{\#} \in [0, 1], \\ \tau_i^{\#} > 0. \end{cases} \quad (12) $$

In the sequel, we will use: $\rho_i^{\#} = \rho_i^n$, $\bar{y}_i^{\#} = \bar{y}_i^n$ and $\tau_i^{\#} = \tau_i^n$.

• Turning to convection contributions, two different methods are used. The first method simply considers the total convection: $\nabla \cdot (yq) + \nabla \cdot (y(1-y)q_r) = \nabla \cdot (yq_g)$. The scheme using this method is nicknamed QG and is detailed in the sequel. It introduces a condition on the time step to preserve the maximum principle. The second method treats the mixture convection and the relative convection separately. Two different schemes are detailed, labeled QRd scheme and QRq scheme. The latter two schemes approximate the non-linear convection so that the problem is well-posed and the discrete mass fraction remains in physical bounds $[0, 1]$ whatever the time step is.

In the sequel, the following notation is used for the sign of $z$: $sg(z) = \begin{cases} 1 & \text{if } z \ge 0 \\ 0 & \text{if } z < 0 \end{cases}$. In order to simplify notations, $sg_{ij} := sg(q_{ij}^*)$, $sg_{ij}^g := sg((q_g)_i^*)$ and $sg_{ij}^r := sg((q_r)_i^*)$ are used.

## 3.2. Global Scheme with QG Scheme

### 3.2.1. Definition of QG scheme

For the gas flow rate convection, a standard implicit upwind scheme is used [9]. With QG scheme, the global scheme is written with the sign convention, for each cell $\Omega_i$:

$$ \rho_i^n \omega_i (y_i^{n+1} - y_i^n) - \Delta t^n \sum_{j \in \nu(i)} S_{ij} q_{ij}^* y_i^{n+1} + \Delta t^n \sum_{j \in \nu(i)} S_{ij} (q_g)_i^* \{ \mathrm{sg}_{ij}^g y_i^{n+1} + (1 - \mathrm{sg}_{ij}^g) y_j^{n+1} \} = \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n} \quad (13) $$

The mass fraction $y_{ij}^n$ at interface $ij$ involved in the gas flow rate $(q_g)_i^* = q_{ij}^* + (1 - y_{ij}^n)(q_r)_i^*$ is approximated here using an upwind scheme based on the mixture flow rate $q_{ij}^*$:

$$ y_{ij}^n \approx \begin{cases} y_i^n & \text{if } q_{ij}^* \ge 0, \\ y_j^n & \text{if } q_{ij}^* < 0. \end{cases} \quad (14) $$

Other consistent choices could be made as long as $y_{ij}^n \in [0, 1]$.
• The mass fraction is taken implicit in the source term

$$ \int_{\omega_i} \rho \frac{\bar{y} - y}{\tau} d\Omega \approx \rho_i^{\#} \frac{\bar{y}_i^{\#} - y_i^{n+1}}{\tau_i^{\#}} \omega_i \quad (11) $$

where $t^{\#} \in [t^n, t^{n+1}]$ such that:

$$ \begin{cases} \rho_i^{\#} > 0, \\ \bar{y}_i^{\#} \in [0, 1], \\ \tau_i^{\#} > 0. \end{cases} \quad (12) $$

In the sequel, we will use: $\rho_i^{\#} = \rho_i^n$, $\bar{y}_i^{\#} = \bar{y}_i^n$ and $\tau_i^{\#} = \tau_i^n$.

• Turning to convection contributions, two different methods are used. The first method simply considers the total convection: $\nabla \cdot (yq) + \nabla \cdot (y(1-y)q_r) = \nabla \cdot (yq_g)$. The scheme using this method is nicknamed QG and is detailed in the sequel. It introduces a condition on the time step to preserve the maximum principle. The second method treats the mixture convection and the relative convection separately. Two different schemes are detailed, labeled QRd scheme and QRq scheme. The latter two schemes approximate the non-linear convection so that the problem is well-posed and the discrete mass fraction remains in physical bounds $[0, 1]$ whatever the time step is.

In the sequel, the following notation is used for the sign of $z$: $\text{sg}(z) = \begin{cases} 1 & \text{if } z \ge 0 \\ 0 & \text{if } z < 0 \end{cases}$. In order to simplify notations, $\text{sg}(q_{ij}^*) := \text{sg}(q_{ij}^*)$, $\text{sg}_{ij}^g := \text{sg}((q_g)_i^*)$ and $\text{sg}_{ij}^r := \text{sg}((q_r)_i^*)$ are used.

## 3.2. Global Scheme with QG Scheme

### 3.2.1. Definition of QG scheme

For the gas flow rate convection, a standard implicit upwind scheme is used [9]. With QG scheme, the global scheme is written with the sign convention, for each cell $\Omega_i$:

$$ \rho_i^n \omega_i (y_i^{n+1} - y_i^n) - \Delta t^n \sum_{j \in \nu(i)} S_{ij} q_{ij}^* y_i^{n+1} + \Delta t^n \sum_{j \in \nu(i)} S_{ij} (q_g)_i^* \{ \text{sg}_{ij}^g y_i^{n+1} + (1 - \text{sg}_{ij}^g) y_j^{n+1} \} = \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n} \quad (13) $$

The mass fraction $y_{ij}^n$ at interface $ij$ involved in the gas flow rate $(q_g)_i^* = q_{ij}^* + (1 - y_{ij}^n)(q_r)_i^*$ is approximated here using an upwind scheme based on the mixture flow rate $q_{ij}^*$:

$$ y_{ij}^n \approx \begin{cases} y_i^n & \text{if } q_{ij}^* \ge 0, \\ y_j^n & \text{if } q_{ij}^* < 0. \end{cases} \quad (14) $$

Other consistent choices could be made as long as $y_{ij}^n \in [0, 1]$.
### 3.2.2. Discrete maximum principle

We examine now if the discrete mass fraction remains within $[0, 1]$ on the $N$ cells with this scheme.

**Property 2 (Maximum principle for the mass fraction with QG scheme):** Assume that the physical parameters are such that: $\tau_k^n > 0$, $\bar{y}_k^n \in [0, 1]$, $k \in [1, N]$. If the initial conditions are such that $\forall k \in [1, N]$, $y_k^n \in [0, 1]$, then the global scheme with QG ensures that $y_i^{n+1}$ remains in $[0, 1]$, $\forall i \in [1, N]$, when the time step $\Delta t^n$ verifies, $\forall i \in [1, N]$:

$$ (1 - y_i^n) + (1 - \bar{y}_i^n) \frac{\Delta t^n}{\tau} + \frac{\Delta t^n}{\rho_i^n \omega_i} \sum_{j \in v(i)} S_{ij} (1 - y_{ij}^n) (q_r)_{ij}^* \geq 0. \quad (15) $$

<figure><img src="image_4.png" /></figure>

**Remark:** Obviously, the condition on the time step is automatically verified when $q_r = 0$.

**Proof.** Using equation (13) obtained for each cell $i \in [1, N]$, the system can be expressed in matrix form.

$$ \mathbf{A} \mathbf{Y}^{n+1} = \mathbf{B} \text{ with } \forall (i, j) \in [1, N]^2, \begin{cases} \mathbf{A} = (a_{ij}), \\ \mathbf{B} = (b_i), \\ \mathbf{Y}^n = (y_i^n). \end{cases} \quad (16) $$

The discrete system for $\hat{\mathbf{Y}} = \mathbf{1} - \mathbf{Y}$ can be expressed too.

$$ \hat{\mathbf{A}} \hat{\mathbf{Y}}^{n+1} = \hat{\mathbf{B}} = \mathbf{A} \times \mathbf{1} - \mathbf{B} \text{ with } \forall (i, j) \in [1, N]^2, \begin{cases} \hat{\mathbf{B}} = (\hat{b}_i) = (\sum_j a_{ij} - b_i), \\ \hat{\mathbf{A}} = (\hat{a}_{ij}) = (a_{ij}), \\ \hat{\mathbf{Y}}^n = (\hat{y}_i^n) = (1 - y_i^n). \end{cases} \quad (17) $$

We also introduce the quantity $\Lambda_i(\mathbf{A})$ on each cell by $\Lambda_i(\mathbf{A}) = |a_{ii}| - \sum_{j \in v(i)} |a_{ij}|$. Note that $\Lambda_i(\mathbf{A}) = \Lambda_i(\hat{\mathbf{A}})$.

The coefficients of the matrix system and the quantity $\Lambda_i$ are

$$ \begin{align} a_{ii} &= \hat{a}_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij}^g - \frac{q_{ij}^*}{\rho_i^n} \right), \notag \\ a_{ij} &= \hat{a}_{ij} = \Delta t^n \frac{S_{ij}}{\omega_i} \frac{(q_g)_{ij}^*}{\rho_i^n} (1 - \mathrm{sg}_{ij}^g), \notag \\ b_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n, \notag \\ \hat{b}_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right). \tag{18} \end{align} $$

Assuming that $a_{ii} > 0$,

$$ \Lambda_i(\mathbf{A}) = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right). $$

**Remark:** The following two formulations are equivalent:

$$ \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right) = \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} (1 - y_{ij}^n) (q_r)_{ij}^* \quad (19) $$
### 3.2.2. Discrete maximum principle

We examine now if the discrete mass fraction remains within $[0, 1]$ on the $N$ cells with this scheme.

**Property 2 (Maximum principle for the mass fraction with QG scheme):** Assume that the physical parameters are such that: $\tau_k^n > 0$, $\bar{y}_k^n \in [0, 1]$, $k \in [1, N]$. If the initial conditions are such that $\forall k \in [1, N]$, $y_k^n \in [0, 1]$, then the global scheme with QG ensures that $y_i^{n+1}$ remains in $[0, 1]$, $\forall i \in [1, N]$, when the time step $\Delta t^n$ verifies, $\forall i \in [1, N]$:

$$ (1 - y_i^n) + (1 - \bar{y}_i^n) \frac{\Delta t^n}{\tau} + \frac{\Delta t^n}{\rho_i^n \omega_i} \sum_{j \in v(i)} S_{ij} (1 - y_{ij}^n) (q_r)_{ij}^* \geq 0. \quad (15) $$

<figure><img src="image_4.png" /></figure>

**Remark:** Obviously, the condition on the time step is automatically verified when $q_r = 0$.

**Proof.** Using equation (13) obtained for each cell $i \in [1, N]$, the system can be expressed in matrix form.

$$ A \mathbf{Y}^{n+1} = B \mathbf{Y}^{n+1} \text{ with } \forall (i, j) \in [1, N]^2, \begin{cases} \mathbf{A} = (a_{ij}), \\ \mathbf{B} = (b_i), \\ \mathbf{Y}^n = (y_i^n). \end{cases} \quad (16) $$

The discrete system for $\hat{\mathbf{Y}} = \mathbf{1} - \mathbf{Y}$ can be expressed too.

$$ \hat{\mathbf{A}} \hat{\mathbf{Y}}^{n+1} = \hat{\mathbf{B}} = \mathbf{A} \times \mathbf{1} - \mathbf{B} \text{ with } \forall (i, j) \in [1, N]^2, \begin{cases} \hat{\mathbf{B}} = (\hat{b}_i) = (\sum_j a_{ij} - b_i), \\ \hat{\mathbf{A}} = (\hat{a}_{ij}) = (a_{ij}), \\ \hat{\mathbf{Y}}^n = (\hat{y}_i^n) = (1 - y_i^n). \end{cases} \quad (17) $$

We also introduce the quantity $\Lambda_i(\mathbf{A})$ on each cell by $\Lambda_i(\mathbf{A}) = |a_{ii}| - \sum_{j \in v(i)} |a_{ij}|$. Note that $\Lambda_i(\mathbf{A}) = \Lambda_i(\hat{\mathbf{A}})$.

The coefficients of the matrix system and the quantity $\Lambda_i$ are

$$ \begin{align} a_{ii} &= \hat{a}_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij}^g - \frac{q_{ij}^*}{\rho_i^n} \right), \notag \\ a_{ij} &= \hat{a}_{ij} = \Delta t^n \frac{S_{ij}}{\omega_i} \frac{(q_g)_{ij}^*}{\rho_i^n} (1 - \mathrm{sg}_{ij}^g), \notag \\ b_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n, \notag \\ \hat{b}_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right). \tag{18} \end{align} $$

Assuming that $a_{ii} > 0$,

$$ \Lambda_i(\mathbf{A}) = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right). $$

**Remark:** The following two formulations are equivalent:

$$ \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right) = \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} (1 - y_{ij}^n) (q_r)_{ij}^* \quad (19) $$
The conditions $\forall (i, j) \in [1, N]^2, a_{ii} > 0, a_{ij} \le 0$ and $\Lambda_i > 0$ ensure that the matrix $\mathbf{A}$ is invertible and its inverse has positive coefficients (see [10, 11]). In order to satisfy the maximum principle, the vectors $\mathbf{B}$ and $\hat{\mathbf{B}}$ also need to fulfill $b_i \ge 0$ and $\hat{b}_i \ge 0$ with $i \in [1, N]$. To summarize, the conditions for the scheme are

$$ \forall i \in [1, N], \left\{ \begin{array}{l} a_{ii} > 0, a_{ij} \le 0, j \in v(i), \\ \Lambda_i(\mathbf{A}) > 0, \\ b_i \ge 0, \hat{b}_i \ge 0. \end{array} \right. \qquad (20) $$

For $\tau_k^n > 0, \bar{y}_k^n \in [0, 1], k \in [1, N]$ and when $y_i^n \in [0, 1], i \in [1, N]$, the maximum principle for the mass fraction with QG scheme is satisfied when the properties (20) are fulfilled. Actually, the condition on $(\hat{b}_i)$ is equivalent to (15). If condition (15) is verified, then the coefficients $(a_{ii})$ and $(\Lambda_i)$ are positive as $a_{ii} \ge \Lambda_i \ge \hat{b}_i$. The other two conditions $a_{ij} \le 0$ and $b_i \ge 0$ are always satisfied. $\square$

### 3.3. Global Scheme with QRd scheme

#### 3.3.1. Definition of QRd scheme

For the mixture flow rate contribution to the convection flux, a standard upwind scheme is once again used as in equation (13).

The global scheme with QRd scheme is written:

$$ \begin{aligned} & \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \{(1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1})\} \\ & + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_i^* \{\mathrm{sg}_{ij}^r y_i^{n+1} (1 - y_j^n) + (1 - \mathrm{sg}_{ij}^r) y_j^n (1 - y_i^{n+1})\} \end{aligned} \qquad (21) \\ & = \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \end{aligned} $$

#### 3.3.2. Discrete maximum principle

We examine now if the discrete mass fraction remains within $[0, 1]$ on the $N$ cells with this scheme.

**Property 3 (Maximum principle for the mass fraction with QRd scheme):** Assume that the physical parameters are such that: $\tau_k^n > 0, \bar{y}_k^n \in [0, 1], \forall k \in [1, N]$. If the initial conditions are such that $\forall k \in [1, N], y_k^n \in [0, 1]$, then the global scheme with QRd ensures that $y_i^{n+1}$ remains in $[0, 1], \forall i \in [1, N]$, whatever the time step $\Delta t^n$ is. $\square$

**Remark:** The scheme is rigorously conservative in space only for the discrete steady states $(y_i^{n+1} = y_i^n, i \in [1, N])$. This is not an issue in practice, when the method is applied to the computations of steady states.

**Proof.** Using the same methodology as for QG scheme, a similar discrete system (noted $d$) is written from (21) for QRd scheme:

$$ \mathbf{A}_d \mathbf{Y}^{n+1} = \mathbf{B}_d \text{ with } \forall (i, j) \in [1, N]^2, \left\{ \begin{array}{l} \mathbf{A}_d = ((a_d)_{ij}), \\ \mathbf{B}_d = ((b_d)_i), \\ \mathbf{Y}^n = (y_i^n). \end{array} \right. \qquad (22) $$
The conditions $\forall (i, j) \in [1, N]^2, a_{ii} > 0, a_{ij} \le 0$ and $\Lambda_i > 0$ ensure that the matrix **A** is invertible and its inverse has positive coefficients (see [10, 11]). In order to satisfy the maximum principle, the vectors **B** and $\hat{\mathbf{B}}$ also need to fulfill $b_i \ge 0$ and $\hat{b}_i \ge 0$ with $i \in [1, N]$. To summarize, the conditions for the scheme are

$$ \forall i \in [1, N], \left\{ \begin{array}{l} a_{ii} > 0, a_{ij} \le 0, j \in v(i), \\ \Lambda_i(\mathbf{A}) > 0, \\ b_i \ge 0, \hat{b}_i \ge 0. \end{array} \right. \qquad (20) $$

For $\tau_k^n > 0, \bar{y}_k^n \in [0, 1], k \in [1, N]$ and when $y_i^n \in [0, 1], i \in [1, N]$, the maximum principle for the mass fraction with QG scheme is satisfied when the properties (20) are fulfilled. Actually, the condition on $(\hat{b}_i)$ is equivalent to (15). If condition (15) is verified, then the coefficients $(a_{ii})$ and $(\Lambda_i)$ are positive as $a_{ii} \ge \Lambda_i \ge \hat{b}_i$. The other two conditions $a_{ij} \le 0$ and $b_i \ge 0$ are always satisfied. $\square$

### 3.3. Global Scheme with QRd scheme

#### 3.3.1. Definition of QRd scheme

For the mixture flow rate contribution to the convection flux, a standard upwind scheme is once again used as in equation (13).

The global scheme with QRd scheme is written:

$$ \begin{aligned} & \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \{(1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1})\} \\ & + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_i^* \{\mathrm{sg}_{ij}^r y_i^{n+1} (1 - y_j^n) + (1 - \mathrm{sg}_{ij}^r) y_j^n (1 - y_i^{n+1})\} \end{aligned} \qquad (21) \\ & = \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \end{aligned} $$

#### 3.3.2. Discrete maximum principle

We examine now if the discrete mass fraction remains within $[0, 1]$ on the $N$ cells with this scheme.

**Property 3 (Maximum principle for the mass fraction with QRd scheme):** Assume that the physical parameters are such that: $\tau_k^n > 0, \bar{y}_k^n \in [0, 1], \forall k \in [1, N]$. If the initial conditions are such that $\forall k \in [1, N], y_k^n \in [0, 1]$, then the global scheme with QRd ensures that $y_i^{n+1}$ remains in $[0, 1]$, $\forall i \in [1, N]$, whatever the time step $\Delta t^n$ is. $\square$

**Remark:** The scheme is rigorously conservative in space only for the discrete steady states $(y_i^{n+1} = y_i^n, i \in [1, N])$. This is not an issue in practice, when the method is applied to the computations of steady states.

**Proof.** Using the same methodology as for QG scheme, a similar discrete system (noted $d$) is written from (21) for QRd scheme:

$$ \mathbf{A}_d \mathbf{Y}^{n+1} = \mathbf{B}_d \text{ with } \forall (i, j) \in [1, N]^2, \left\{ \begin{array}{l} \mathbf{A}_d = ((a_d)_{ij}), \\ \mathbf{B}_d = ((b_d)_i), \\ \mathbf{Y}^n = (y_i^n). \end{array} \right. \qquad (22) $$
Using the blue color to denote the terms linked to QRd scheme, the coefficients of the matrix system and the quantity $\Lambda_i$ are

$$
\begin{align}
(a_d)_{ii} &= (\hat{a}_d)_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{q_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - y_j^n) \right), \\
(a_d)_{ij} &= (\hat{a}_d)_{ij} = \Delta t^n \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}), \\
(b_d)_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - 1) y_j^n \right), \tag{23} \\
(\hat{b}_d)_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij}^r (1 - y_j^n) \right), \\
\Lambda_i(\mathbf{A}_d) &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - y_j^n) \right). \nonumber
\end{align}
$$

For given values $q_{kl}^n$ and $(q_r)_{kl}^n$, $(k, l) \in [1, N]^2$, and for $\tau_k^n > 0$, $\bar{y}_k^n \in [0, 1]$, $k \in [1, N]$, the following conditions are fulfilled when $y_i^n \in [0, 1]$, $i \in [1, N]$:

$$
\forall i \in [1, N], \quad \left\{
\begin{array}{ll}
(a_d)_{ii} > 0, & (a_d)_{ij} \le 0, \quad j \in v(i), \\
\Lambda_i(\mathbf{A}_d) > 0, & \\
(b_d)_i \ge 0, & (b_d)_i \ge 0,
\end{array}
\right.
\qquad (24)
$$

whatever the time step $\Delta t^n$ is.

Once again using the results [10, 11], the matrix $\mathbf{A}_d$ is invertible and its inverse has positive coefficients. Hence the maximum principle for the mass fraction $y$ is satisfied, using QRd scheme, whatever the time step $\Delta t^n$ is. $\square$

3.4. Global Scheme with QRq scheme

3.4.1. Definition of QRq scheme

QRq scheme not only takes into account the sign of the relative flow rate but also the sign of the mixture flow rate inside the relative flow rate contribution to the convection flux. The global QRq scheme is

$$
\begin{align}
& \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \left\{ (1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1}) \right\} \nonumber \\
& \quad + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_{ij}^* \left[ \begin{aligned}[t]
& \mathrm{sg}_{ij} \{ y_i^n (1 - y_i^{n+1}) (1 - \mathrm{sg}_{ij}^r) + y_i^{n+1} (1 - y_i^n) \mathrm{sg}_{ij}^r \} \\
& + (1 - \mathrm{sg}_{ij}) \{ y_j^n (1 - y_j^{n+1}) (1 - \mathrm{sg}_{ij}^r) + y_j^{n+1} (1 - y_j^n) \mathrm{sg}_{ij}^r \}
\end{aligned} \right] \tag{25} \\
&= \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \nonumber
\end{align}
$$
Using the blue color to denote the terms linked to QRd scheme, the coefficients of the matrix system and the quantity $\Lambda_i$ are

$$
\begin{align}
(a_d)_{ii} &= (\hat{a}_d)_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{q_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - y_j^n) \right), \\
(a_d)_{ij} &= (\hat{a}_d)_{ij} = \Delta t^n \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}), \\
(b_d)_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - 1) y_j^n \right), \tag{23} \\
(\hat{b}_d)_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij}^r (1 - y_j^n) \right), \\
\Lambda_i(\mathbf{A}_d) &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - y_j^n) \right). \nonumber
\end{align}
$$

For given values $q_{kl}^n$ and $(q_r)_{kl}^n, (k, l) \in [1, N]^2$, and for $\tau_k^n > 0$, $\bar{y}_k^n \in [0, 1]$, $k \in [1, N]$, the following conditions are fulfilled when $y_i^n \in [0, 1]$, $i \in [1, N]$:

$$
\forall i \in [1, N], \left\{
\begin{array}{ll}
(a_d)_{ii} > 0, & (a_d)_{ij} \le 0, \quad j \in v(i), \\
\Lambda_i(\mathbf{A}_d) > 0, & \\
(b_d)_i \ge 0, & (b_d)_i \ge 0,
\end{array}
\right.
\qquad (24)
$$

whatever the time step $\Delta t^n$ is.

Once again using the results [10, 11], the matrix $\mathbf{A}_d$ is invertible and its inverse has positive coefficients. Hence the maximum principle for the mass fraction $y$ is satisfied, using QRd scheme, whatever the time step $\Delta t^n$ is. $\square$

3.4. Global Scheme with QRq scheme

3.4.1. Definition of QRq scheme

QRq scheme not only takes into account the sign of the relative flow rate but also the sign of the mixture flow rate inside the relative flow rate contribution to the convection flux. The global QRq scheme is

$$
\begin{align}
& \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \left\{ (1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1}) \right\} \nonumber \\
& \quad + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_{ij}^* \left[ \begin{aligned}[t]
& \mathrm{sg}_{ij} \{ y_i^n (1 - y_i^{n+1}) (1 - \mathrm{sg}_{ij}^r) + y_i^{n+1} (1 - y_i^n) \mathrm{sg}_{ij}^r \} \\
& + (1 - \mathrm{sg}_{ij}) \{ y_j^n (1 - y_j^{n+1}) (1 - \mathrm{sg}_{ij}^r) + y_j^{n+1} (1 - y_j^n) \mathrm{sg}_{ij}^r \}
\end{aligned} \right] \tag{25} \\
&= \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \nonumber
\end{align}
$$
3.4.2. Discrete maximum Principle

**Property 4 (Maximum Principle for the mass fraction with QRq scheme):** Assume that the physical parameters are such that: $\tau_k^n > 0, \bar{y}_k^n \in [0,1], k \in [1, N]$. If the initial conditions are such that $k \in [1, N]$, $y_k^n \in [0, 1]$, then QRq scheme ensures that $y_i^{n+1}$ remains in $[0, 1], i \in [1, N]$, whatever the time step $\Delta t^n$ is, provided that the flow rates fulfill the following conditions, when $q_{ij}^* < 0$:

* If $(q_r)_{ij}^* \ge 0$

$$q_{ij}^* + (q_r)_{ij}^* (1 - y_j^n) \le 0. \qquad (26)$$

* Otherwise if $(q_r)_{ij}^* < 0$

$$q_{ij}^* - (q_r)_{ij}^* y_j^n \le 0. \qquad (27)$$

No condition arises when $q_{ij}^* \ge 0$.

<figure><img src="image_5.png" /></figure>

For a co-current flow, these conditions are automatically verified.

**Remark:** Once again, the scheme is rigorously conservative in space only for the steady states $(y_i^{n+1} = y_i^n, i \in [1, N])$.

**Proof.** Using the same methodology as for QRd scheme, the discrete system from (25) is written

$$\mathbf{A}_q \mathbf{Y}^{n+1} = \mathbf{B}_q \text{ with } \forall (i, j) \in [1, N]^2, \begin{cases} \mathbf{A}_q = ((a_q)_{ij}), \\ \mathbf{B}_q = ((b_q)_i), \\ \mathbf{Y}^n = (y_i^n). \end{cases} \qquad (28)$$

Using the blue color to denote the terms linked to QRq scheme, the coefficients of the matrix system and the quantity $\Lambda_i(\mathbf{A}_q)$ are

$$(a_q)_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{q_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij} (\mathrm{sg}_{ij}^r - y_i^n) \right),$$

$$(a_q)_{ij} = \Delta t^n \frac{S_{ij}}{\omega_i} \frac{q_{ij}^*}{\rho_i^n} (1 - \mathrm{sg}_{ij}) + \Delta t^n \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (1 - \mathrm{sg}_{ij}) (\mathrm{sg}_{ij}^r - y_j^n),$$

$$(b_q)_i = y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\mathrm{sg}_{ij}^r - 1) (\mathrm{sg}_{ij}^r y_i^n + (1 - \mathrm{sg}_{ij}^r) y_j^n) \right),$$

$$(\hat{b}_q)_i = (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij}^r (1 - [\mathrm{sg}_{ij}^r y_i^n + (1 - \mathrm{sg}_{ij}^r) y_j^n]) \right).$$

Assuming that $\forall j \in v(i), (a_q)_{ij} < 0$,

$$\Lambda_i(\mathbf{A}_q) = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \left( \mathrm{sg}_{ij}^r - \left[ \mathrm{sg}_{ij}^r y_i^n + (1 - \mathrm{sg}_{ij}^r) y_j^n \right] \right) \right). \qquad (29)$$

Examining the sign of these coefficients raises two conditions on the flow rates.

* When $q_{ij}^* < 0$ and $(q_r)_{ij}^* \ge 0$

$$q_{ij}^* + (q_r)_{ij}^* (1 - y_j^n) \le 0. \qquad (30)$$

* When $q_{ij}^* < 0$ and $(q_r)_{ij}^* < 0$

$$q_{ij}^* - (q_r)_{ij}^* y_j^n \le 0. \qquad (31)$$
### 3.4.2. Discrete maximum Principle

**Property 4 (Maximum Principle for the mass fraction with QRq scheme):** Assume that the physical parameters are such that: $\tau_k^n > 0, \bar{y}_k^n \in [0,1], k \in [1, N]$. If the initial conditions are such that $k \in [1, N]$, $y_k^n \in [0, 1]$, then QRq scheme ensures that $y_i^{n+1}$ remains in $[0, 1], i \in [1, N]$, whatever the time step $\Delta t^n$ is, provided that the flow rates fulfill the following conditions, when $q_{ij}^* < 0$:

* If $(q_r)_{ij}^* \ge 0$

$$q_{ij}^* + (q_r)_{ij}^* (1 - y_j^n) \le 0. \qquad (26)$$

* Otherwise if $(q_r)_{ij}^* < 0$

$$q_{ij}^* - (q_r)_{ij}^* y_j^n \le 0. \qquad (27)$$

No condition arises when $q_{ij}^* \ge 0$.

<figure><img src="image_5.png" /></figure>

For a co-current flow, these conditions are automatically verified.

**Remark:** Once again, the scheme is rigorously conservative in space only for the steady states $(y_i^{n+1} = y_i^n, i \in [1, N])$.

**Proof.** Using the same methodology as for QRd scheme, the discrete system from (25) is written

$$\mathbf{A}_q \mathbf{Y}^{n+1} = \mathbf{B}_q \text{ with } \forall (i, j) \in [1, N]^2, \begin{cases} \mathbf{A}_q = ((a_q)_{ij}), \\ \mathbf{B}_q = ((b_q)_i), \\ \mathbf{Y}^n = (y_i^n). \end{cases} \qquad (28)$$

Using the blue color to denote the terms linked to QRq scheme, the coefficients of the matrix system and the quantity $\Lambda_i(\mathbf{A}_q)$ are

$$(a_q)_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{q_{ij}^*}{\rho_i^n} (\text{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \text{sg}_{ij} (\text{sg}_{ij}^r - y_i^n) \right),$$

$$(a_q)_{ij} = \Delta t^n \frac{S_{ij}}{\omega_i} \frac{q_{ij}^*}{\rho_i^n} (1 - \text{sg}_{ij}) + \Delta t^n \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (1 - \text{sg}_{ij}) (\text{sg}_{ij}^r - y_j^n),$$

$$(b_q)_i = y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} (\text{sg}_{ij}^r - 1) (\text{sg}_{ij}^r y_i^n + (1 - \text{sg}_{ij}^r) y_j^n) \right),$$

$$(\hat{b}_q)_i = (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \text{sg}_{ij}^r (1 - [\text{sg}_{ij}^r y_i^n + (1 - \text{sg}_{ij}^r) y_j^n]) \right).$$

Assuming that $\forall j \in v(i), (a_q)_{ij} < 0$,

$$\Lambda_i(\mathbf{A}_q) = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \frac{(q_r)_{ij}^*}{\rho_i^n} \left( \text{sg}_{ij}^r - \left[ \text{sg}_{ij}^r y_i^n + (1 - \text{sg}_{ij}^r) y_j^n \right] \right) \right). \qquad (29)$$

Examining the sign of these coefficients raises two conditions on the flow rates.

* When $q_{ij}^* < 0$ and $(q_r)_{ij}^* \ge 0$

$$q_{ij}^* + (q_r)_{ij}^* (1 - y_j^n) \le 0. \qquad (30)$$

* When $q_{ij}^* < 0$ and $(q_r)_{ij}^* < 0$

$$q_{ij}^* - (q_r)_{ij}^* y_j^n \le 0. \qquad (31)$$
When $q_{kl}^n$ and $(q_r)_{kl}^n$, $(k, l) \in [1, N]^2$ satisfy the previous conditions (30)-(31), and if $\tau_k^n > 0$, $\bar{y}_k^n \in [0, 1]$, $k \in [1, N]$, the following properties are verified when $y_i^n \in [0, 1]$, $i \in [1, N]$:

$$ \forall i \in [1, N], \left\{ \begin{array}{l} (a_q)_{ii} > 0, (a_q)_{ij} \le 0, j \in \nu(i), \\ \Lambda_i(\mathbf{A}_q) > 0, \\ (b_q)_i \ge 0, (\hat{b}_q)_i \ge 0, \end{array} \right. \qquad (32) $$

**whatever the time step $\Delta t^n$ is.**

Once more, using common results of [10, 11], the matrix $\mathbf{A}_q$ is invertible and its inverse has positive coefficients. Hence the maximum principle for the mass fraction with QRq scheme is verified whatever the time step $\Delta t^n$ is, as long as the conditions (30)-(31) on the flow rates are respected. $\square$

<figure><img src="image_6.png" /></figure>

### 3.5. Boundary Conditions

These Finite Volume cell schemes presented above are only valid for cells which do not contain faces on the boundary. The cells (noted $i$ here) containing at least one face on the boundary are tackled in **Appendix B**. Boundary conditions are considered for three schemes with a given valid mass fraction $y_\infty \in [0, 1]$ outside. The inlet/outlet flow rates for the faces on the boundary are noted $q_{i\infty}$ and $(q_r)_{i\infty}$.

We summarise here the main conclusions of **Appendix B**. If a co-current flow is considered, QRq scheme still satisfies the discrete maximum principle without any condition on the time step. Turning to QG scheme, a slightly different condition on the time step arises on cell $i$ sharing a face with the boundary. This condition must be taken into account as it can be the most constraining one. Eventually, QRd scheme still has no limit on the time step. However, the boundary flux used in the QRd scheme must be handled carefully on an outlet face (with respect to the mixture flow rate $q_{i\infty} > 0$). When the mass fraction is unknown on this face, the expression of the flux should be modified or another boundary flux should be preferred.
When $q_{kl}^n$ and $(q_r)_{kl}^n$, $(k, l) \in [1, N]^2$ satisfy the previous conditions (30)-(31), and if $\tau_k^n > 0$, $\bar{y}_k^n \in [0, 1]$, $k \in [1, N]$, the following properties are verified when $y_i^n \in [0, 1]$, $i \in [1, N]$:

$$ \forall i \in [1, N], \left\{ \begin{array}{l} (a_q)_{ii} > 0, (a_q)_{ij} \le 0, j \in \nu(i), \\ \Lambda_i(A_q) > 0, \\ (b_q)_i \ge 0, (\hat{b}_q)_i \ge 0, \end{array} \right. \qquad (32) $$

whatever the time step $\Delta t^n$ is.

Once more, using common results of [10, 11], the matrix $A_q$ is invertible and its inverse has positive coefficients. Hence the maximum principle for the mass fraction with QRq scheme is verified whatever the time step $\Delta t^n$ is, as long as the conditions (30)-(31) on the flow rates are respected. $\square$

<figure><img src="image_6.png" /></figure>

### 3.5. Boundary Conditions

These Finite Volume cell schemes presented above are only valid for cells which do not contain faces on the boundary. The cells (noted $i$ here) containing at least one face on the boundary are tackled in **Appendix B**. Boundary conditions are considered for three schemes with a given valid mass fraction $y_\infty \in [0, 1]$ outside. The inlet/outlet flow rates for the faces on the boundary are noted $q_{i\infty}$ and $(q_r)_{i\infty}$.

We summarise here the main conclusions of **Appendix B**. If a co-current flow is considered, QRq scheme still satisfies the discrete maximum principle without any condition on the time step. Turning to QG scheme, a slightly different condition on the time step arises on cell $i$ sharing a face with the boundary. This condition must be taken into account as it can be the most constraining one. Eventually, QRd scheme still has no limit on the time step. However, the boundary flux used in the QRd scheme must be handled carefully on an outlet face (with respect to the mixture flow rate $q_{i\infty} > 0$). When the mass fraction is unknown on this face, the expression of the flux should be modified or another boundary flux should be preferred.
## 4. Numerical results

In order to be able to compare the numerical solution of the different schemes, obtained with an analytical solution of (3), a steady one-dimensional test case is considered. The physical parameters ($\rho, q, q_r, \tau, \bar{y}$) are taken uniform and constant, and will be noted with an exponent 0. The simplified system is

$$ \begin{cases} \rho = \rho^0, & q = q^0, \\ \frac{\partial}{\partial x}(q^0 y) + \frac{\partial}{\partial x}(q_r^0 y(1-y)) = \frac{\rho^0 (\bar{y}^0 - y)}{\tau^0}, \end{cases} \quad (33) $$
with $q_r = q_r^0; \tau = \tau^0; \bar{y} = \bar{y}^0$.

Setting $q_{\text{lin}}^0 = 2q_r^0$, $q_{\text{log}}^0 = q^0 + q_r^0 - q_{\text{lin}}^0 \bar{y}^0$, $\gamma^0 = \frac{\rho^0}{\tau^0}$, $y_\infty = \bar{y}^0$, this equation can be rewritten in the following form:

$$ (q_{\text{log}}^0 + q_{\text{lin}}^0 (y_\infty - y)) \frac{\partial y}{\partial x} = \gamma^0 (y_\infty - y). \quad (34) $$

A domain $\Omega = [0, 1]$ $m$ is considered and a Dirichlet boundary condition is used on the left boundary $y(0, t) = y_0$. The analytical solution $y(x)$ reads

$$ q_{\text{lin}}^0 (y - y_0) - q_{\text{log}}^0 \ln \left( \frac{y_\infty - y}{y_\infty - y_0} \right) = \gamma^0 x. \quad (35) $$

For numerical simulations, we use the parameters of Table 1.

<table><caption>**Table 1.** Parameters of the test case</caption><thead><tr><th>Quantity</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>y<sub>∞</sub></td><td>3.04 × 10<sup>-2</sup></td><td>—</td></tr><tr><td>q<sup>0</sup><sub>lin</sub></td><td>-6000</td><td>kg.m<sup>-2</sup>.s<sup>-1</sup></td></tr><tr><td>q<sup>0</sup><sub>log</sub></td><td>682.57</td><td>kg.m<sup>-2</sup>.s<sup>-1</sup></td></tr><tr><td>γ<sup>0</sup></td><td>7 × 10<sup>4</sup></td><td>kg.m<sup>-3</sup>.s<sup>-1</sup></td></tr></tbody></table>

The stationary numerical solution $y_{\text{num}}$ and the exact solution $y_{\text{ex}}$ are computed, using a mesh with $\Delta x = 10^{-3}$ m. Results are given on figure 1 for QRd scheme. Similar solutions are obtained with QG and QRq schemes.

Figure 2 shows the convergence rate with the three schemes for this test case using different mesh sizes: $\Delta x \in [10^{-5}, 10^{-1}]$ m. The $L_2$ norm of the error $\text{err}_2$, computed for each mesh size, is

$$ \text{err}_2 = \frac{\|y_{\text{num}} - y_{\text{ex}}\|_2}{\|y_{\text{ex}}\|_2} \quad \text{with} \quad \|x\|_2 = \sqrt{\sum_{i=1}^{N_x} x_i^2}, \quad N_x = \lfloor \frac{1}{\Delta x} \rfloor. \quad (36) $$

As expected, the three different schemes comply with the convergence rate of first order in space. For a given mesh size, QRq scheme is more accurate than QRd scheme. QG scheme and QRq scheme have similar accuracy for this test case.
## 4. Numerical results

In order to be able to compare the numerical solution of the different schemes, obtained with an analytical solution of (3), a steady one-dimensional test case is considered. The physical parameters ($\rho, q, q_r, \tau, \bar{y}$) are taken uniform and constant, and will be noted with an exponent 0. The simplified system is

$$ \begin{cases} \rho = \rho^0, & q = q^0, \\ \frac{\partial}{\partial x}(q^0 y) + \frac{\partial}{\partial x}(q_r^0 y(1-y)) = \frac{\rho^0 (\bar{y}^0 - y)}{\tau^0}, \end{cases} \quad (33) $$
with $q_r = q_r^0; \tau = \tau^0; \bar{y} = \bar{y}^0$.

Setting $q_{\text{lin}}^0 = 2q_r^0$, $q_{\text{log}}^0 = q^0 + q_r^0 - q_{\text{lin}}^0 \bar{y}^0$, $\gamma^0 = \frac{\rho^0}{\tau^0}$, $y_\infty = \bar{y}^0$, this equation can be rewritten in the following form:

$$ (q_{\text{log}}^0 + q_{\text{lin}}^0 (y_\infty - y)) \frac{\partial y}{\partial x} = \gamma^0 (y_\infty - y). \quad (34) $$
A domain $\Omega = [0, 1]$ $m$ is considered and a Dirichlet boundary condition is used on the left boundary $y(0, t) = y_0$. The analytical solution $y(x)$ reads

$$ q_{\text{lin}}^0 (y - y_0) - q_{\text{log}}^0 \ln \left( \frac{y_\infty - y}{y_\infty - y_0} \right) = \gamma^0 x. \quad (35) $$
For numerical simulations, we use the parameters of Table 1.

<table><caption>**Table 1.** Parameters of the test case</caption><thead><tr><th>Quantity</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>y<sub>∞</sub></td><td>3.04 × 10<sup>-2</sup></td><td>—</td></tr><tr><td>q<sup>0</sup><sub>lin</sub></td><td>-6000</td><td>kg.m<sup>-2</sup>.s<sup>-1</sup></td></tr><tr><td>q<sup>0</sup><sub>log</sub></td><td>682.57</td><td>kg.m<sup>-2</sup>.s<sup>-1</sup></td></tr><tr><td>γ<sup>0</sup></td><td>7 × 10<sup>4</sup></td><td>kg.m<sup>-3</sup>.s<sup>-1</sup></td></tr></tbody></table>

The stationary numerical solution $y_{\text{num}}$ and the exact solution $y_{\text{ex}}$ are computed, using a mesh with $\Delta x = 10^{-3}$ m. Results are given on figure 1 for QRd scheme. Similar solutions are obtained with QG and QRq schemes.

Figure 2 shows the convergence rate with the three schemes for this test case using different mesh sizes: $\Delta x \in [10^{-5}, 10^{-1}]$ m. The $L_2$ norm of the error $\text{err}_2$, computed for each mesh size, is

$$ \text{err}_2 = \frac{\|y_{\text{num}} - y_{\text{ex}}\|_2}{\|y_{\text{ex}}\|_2} \quad \text{with} \quad \|x\|_2 = \sqrt{\sum_{i=1}^{N_x} x_i^2}, \quad N_x = \lfloor \frac{1}{\Delta x} \rfloor. \quad (36) $$

As expected, the three different schemes comply with the convergence rate of first order in space. For a given mesh size, QRq scheme is more accurate than QRd scheme. QG scheme and QRq scheme have similar accuracy for this test case.
<figure><img src="image_7.png" /><figcaption>**Figure 1.** Mass fraction $y$ as a function of $x$ - QRd scheme (blue line) and exact (dotted orange) solutions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 2.** $L_2$ Norm of the error as a function of the mesh size for QG (green), QRd (red) and QRq (orange) schemes.</figcaption></figure>

**5. Conclusion**

Three different linear first-order schemes approximating the non-linear relative velocity term have been presented. QG scheme (13) is always conservative in time and space but a condition appears on the time step to preserve the discrete maximum principle when a non null velocity is considered. On the other side, QRd and QRq schemes (21,25) comply with the maximum principle without any constraint on the time step. However, the latter two schemes are conservative in space only once the steady state solution is reached. This is obviously not an issue for many practical computations which aim at approximating steady states. Eventually, QRq scheme seems to be more accurate than QRd scheme.
<figure><img src="image_7.png" /><figcaption>**Figure 1.** Mass fraction $y$ as a function of $x$ - QRd scheme (blue line) and exact (dotted orange) solutions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 2.** $L_2$ Norm of the error as a function of the mesh size for QG (green), QRd (red) and QRq (orange) schemes.</figcaption></figure>

**5. Conclusion**

Three different linear first-order schemes approximating the non-linear relative velocity term have been presented. QG scheme (13) is always conservative in time and space but a condition appears on the time step to preserve the discrete maximum principle when a non null velocity is considered. On the other side, QRd and QRq schemes (21,25) comply with the maximum principle without any constraint on the time step. However, the latter two schemes are conservative in space only once the steady state solution is reached. This is obviously not an issue for many practical computations which aim at approximating steady states. Eventually, QRq scheme seems to be more accurate than QRd scheme.
**Conflicts of interest**

The authors declare no competing financial interest.

**Dedication**

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

**Acknowledgments**

The first author is supported by EDF R&D and ANRT through CIFRE contract 2022/1027. The authors would like to thank Erwan Le Coupanec for his guidance and help during the writing of this article.
**Conflicts of interest**

The authors declare no competing financial interest.

**Dedication**

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

**Acknowledgments**

The first author is supported by EDF R&D and ANRT through CIFRE contract 2022/1027. The authors would like to thank Erwan Le Coupanec for his guidance and help during the writing of this article.
## Appendix A - Continuous maximum principle

The studied equations are

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot \mathbf{q} = 0, \\ \frac{\partial(\rho y)}{\partial t} + \nabla \cdot (y\mathbf{q}) + \nabla \cdot (y(1-y)\mathbf{q}_r) = \rho \frac{\bar{y}-y}{\tau}. \quad (37) $$

The notation $\xi = y(1 - y)$ is used to study the maximum principle for y. Indeed, if $\xi \ge 0$ then

$$ y(1 - y) \ge 0 \Leftrightarrow y \in [0, 1]. \quad (38) $$

The governing equation for $\xi$ can be obtained from the non-conservative equations of y and $1 - y$ as described below:

$$ \left( \rho \partial_t y + \mathbf{q} \cdot \nabla y + \nabla \cdot (\xi \mathbf{q}_r) \right) = \rho \frac{\bar{y}-y}{\tau} \times (1-y) \\ + \left( \rho \partial_t (1-y) + \mathbf{q} \cdot \nabla (1-y) - \nabla \cdot (\xi \mathbf{q}_r) \right) = \rho \frac{(1-\bar{y})-(1-y)}{\tau} \times y. \quad (39) $$

Then

$$ \rho \partial_t \xi + \mathbf{q} \cdot \nabla \xi + (1-2y)\nabla \cdot (\xi \mathbf{q}_r) = \rho \frac{S-2\xi}{\tau} \quad \text{with } S = y(1-\bar{y}) + \bar{y}(1-y). \quad (40) $$

Using the standard notations $\xi_+ = \max(\xi, 0)$ and $\xi_- = -\min(\xi, 0)$, equation (40) is multiplied by $(-\xi_-)$:

$$ \left( \rho \partial_t \left( \frac{\xi_-^2}{2} \right) + \mathbf{q} \cdot \nabla \left( \frac{\xi_-^2}{2} \right) \right) + (1-2y) \left( \xi_-^2 \nabla \cdot \mathbf{q}_r - \xi_- \mathbf{q}_r \cdot \nabla \xi \right) = \rho \frac{-S\xi_- - 2\xi_-^2}{\tau} \\ \left( \partial_t \left( \rho \frac{\xi_-^2}{2} \right) + \nabla \cdot \left( \mathbf{q} \frac{\xi_-^2}{2} \right) \right) + \nabla \cdot \left( (1-2y)\mathbf{q}_r \frac{\xi_-^2}{2} \right) + \xi_-^2 \left\{ (1-2y)\nabla \cdot \mathbf{q}_r - \frac{1}{2}\nabla \cdot \left( (1-2y)\mathbf{q}_r \right) \right\} = - \frac{\rho(S+2\xi_-)\xi_-}{\tau}. \quad (41) $$

Finally, the equation obtained is

$$ \partial_t \left( \rho \frac{\xi_-^2}{2} \right) + \nabla \cdot \left( \left[ \mathbf{q} + (1-2y)\mathbf{q}_r \right] \frac{\xi_-^2}{2} \right) + \frac{\xi_-^2}{2} \left\{ (1-2y)\nabla \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \nabla (1-2y) \right\} = - \frac{\rho(S+2\xi_-)}{\tau} \xi_- . \quad (42) $$

We define the quantity $E(t) = \int_{\Omega} \rho \frac{\xi_-^2}{2} d\Omega$ and split the outside surface domain according to the sign of $\mathbf{q}_{\xi} = \mathbf{q} + (1-2y)\mathbf{q}_r$:

$$ \begin{aligned} \Gamma_w^{\xi} &= \{\mathbf{x} \in \Gamma, \mathbf{q}_{\xi} \cdot \mathbf{n}_{\Gamma} = 0\}, \\ \Gamma_+^{\xi} &= \{\mathbf{x} \in \Gamma, \mathbf{q}_{\xi} \cdot \mathbf{n}_{\Gamma} < 0\}, \\ \Gamma_-^{\xi} &= \{\mathbf{x} \in \Gamma, \mathbf{q}_{\xi} \cdot \mathbf{n}_{\Gamma} > 0\}. \end{aligned} \quad (43) $$

Then, the quantity $E(t)$ verifies

$$ \frac{dE(t)}{dt} = - \int_{\Gamma} \left[ \left( \mathbf{q} + (1-2y)\mathbf{q}_r \right) \frac{\xi_-^2}{2} \right] \cdot \mathbf{n}_{\Gamma} d\Gamma - \int_{\Omega} \frac{\xi_-^2}{2} \left\{ (1-2y)\nabla \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \nabla (1-2y) \right\} d\Omega \\ - \int_{\Omega} \frac{\rho(S+2\xi_-)}{\tau} \xi_- . \quad (44) $$
## Appendix A - Continuous maximum principle

The studied equations are

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot \mathbf{q} = 0, \\ \frac{\partial(\rho y)}{\partial t} + \nabla \cdot (y\mathbf{q}) + \nabla \cdot (y(1-y)\mathbf{q}_r) = \rho \frac{\bar{y}-y}{\tau}. \quad (37) $$

The notation $\xi = y(1 - y)$ is used to study the maximum principle for y. Indeed, if $\xi \ge 0$ then

$$ y(1 - y) \ge 0 \Leftrightarrow y \in [0, 1]. \quad (38) $$

The governing equation for $\xi$ can be obtained from the non-conservative equations of y and $1 - y$ as described below:

$$ \left( \rho \partial_t y + \mathbf{q} \cdot \nabla y + \nabla \cdot (\xi \mathbf{q}_r) \right) = \rho \frac{\bar{y}-y}{\tau} \times (1-y) \\ + \left( \rho \partial_t (1-y) + \mathbf{q} \cdot \nabla (1-y) - \nabla \cdot (\xi \mathbf{q}_r) \right) = \rho \frac{(1-\bar{y})-(1-y)}{\tau} \times y. \quad (39) $$

Then

$$ \rho \partial_t \xi + \mathbf{q} \cdot \nabla \xi + (1-2y) \nabla \cdot (\xi \mathbf{q}_r) = \rho \frac{S-2\xi}{\tau} \quad \text{with } S = y(1-\bar{y}) + \bar{y}(1-y). \quad (40) $$

Using the standard notations $\xi_+ = \max(\xi, 0)$ and $\xi_- = -\min(\xi, 0)$, equation (40) is multiplied by $(-\xi_-)$:

$$ \left( \rho \partial_t \left( \frac{\xi_-^2}{2} \right) + \mathbf{q} \cdot \nabla \left( \frac{\xi_-^2}{2} \right) \right) + (1-2y) \left( \xi_-^2 \nabla \cdot \mathbf{q}_r - \xi_- \mathbf{q}_r \cdot \nabla \xi \right) = \rho \frac{-S\xi_- - 2\xi_-^2}{\tau} \\ \left( \partial_t \left( \rho \frac{\xi_-^2}{2} \right) + \nabla \cdot \left( \mathbf{q} \frac{\xi_-^2}{2} \right) \right) + \nabla \cdot \left( (1-2y)\mathbf{q}_r \frac{\xi_-^2}{2} \right) + \xi_-^2 \left\{ (1-2y)\nabla \cdot \mathbf{q}_r - \frac{1}{2}\nabla \cdot \left( (1-2y)\mathbf{q}_r \right) \right\} = - \frac{\rho(S+2\xi_-)\xi_-}{\tau}. \quad (41) $$

Finally, the equation obtained is

$$ \partial_t \left( \rho \frac{\xi_-^2}{2} \right) + \nabla \cdot \left( \left[ \mathbf{q} + (1-2y)\mathbf{q}_r \right] \frac{\xi_-^2}{2} \right) + \frac{\xi_-^2}{2} \left\{ (1-2y)\nabla \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \nabla (1-2y) \right\} = - \frac{\rho(S+2\xi_-)}{\tau} \xi_- . \quad (42) $$

We define the quantity $E(t) = \int_{\Omega} \rho \frac{\xi_-^2}{2} d\Omega$ and split the outside surface domain according to the sign of $\mathbf{q}_{\xi} = \mathbf{q} + (1-2y)\mathbf{q}_r$:

$$ \begin{aligned} \Gamma_w^{\xi} &= \{\mathbf{x} \in \Gamma, \mathbf{q}_{\xi} \cdot \mathbf{n}_{\Gamma} = 0\}, \\ \Gamma_+^{\xi} &= \{\mathbf{x} \in \Gamma, \mathbf{q}_{\xi} \cdot \mathbf{n}_{\Gamma} < 0\}, \\ \Gamma_-^{\xi} &= \{\mathbf{x} \in \Gamma, \mathbf{q}_{\xi} \cdot \mathbf{n}_{\Gamma} > 0\}. \end{aligned} \quad (43) $$

Then, the quantity $E(t)$ verifies

$$ \frac{dE(t)}{dt} = - \int_{\Gamma} \left[ \left( \mathbf{q} + (1-2y)\mathbf{q}_r \right) \frac{\xi_-^2}{2} \right] \cdot \mathbf{n}_{\Gamma} d\Gamma - \int_{\Omega} \frac{\xi_-^2}{2} \left\{ (1-2y)\nabla \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \nabla (1-2y) \right\} d\Omega \\ - \int_{\Omega} \frac{\rho(S+2\xi_-)}{\tau} \xi_- . \quad (44) $$
**Assumptions:**

(1) Relevant mass fraction on the boundary conditions: $\xi_{-}(\mathbf{x} \in \Gamma_{+}^{\xi}, t) = 0$.

(2) Mass fraction for the initial condition such that: $\xi_{-}(\mathbf{x} \in \Omega, t = 0) = 0$ and thus $E(t = 0) = 0$.

(3) $\frac{1}{\rho}[(1-2y)\underline{\nabla} \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \underline{\nabla}(1-2y)] \in \mathcal{L}^{\infty}(\Omega, [0, T])$.

(4) Equilibrium mass fraction such that: $\bar{y} \in [0, 1]$ implying $(S + 2\xi_{-}) \ge 0$ (proof of lemma below).

(5) Positive relaxation time scale $\tau > 0$ and density $\rho > 0$.

Thus we have

$$ \frac{dE(t)}{dt} \le - \int_{\Omega} \rho \frac{\xi_{-}^2}{2} \left[ \frac{1}{\rho} \left( (1-2y)\underline{\nabla} \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \underline{\nabla}(1-2y) \right) \right] d\Omega. \quad (45) $$

Assuming the assumptions previously presented,

$$ \frac{dE(t)}{dt} \le \left\| \frac{1}{\rho} \left[ (1-2y)\underline{\nabla} \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \underline{\nabla}(1-2y) \right] \right\|_{\infty} E(t). \quad (46) $$

Using the Grönwall's inequality [12] and $E(0) = 0$, it enables to conclude that $E(t) = 0$ on $[0,T]$.

**Lemma 1.** If $\bar{y} \in [0, 1]$ then $S + 2\xi_{-} \ge 0$ with $S = \bar{y}(1 - y) + y(1 - \bar{y})$

**Proof.** Let's assume that $\bar{y} \in [0, 1]$.

**Case $\xi_{-} = 0$**

If $\xi_{-} = 0$, then $y \in [0, 1]$ implying that $S \ge 0$.

As a consequence, we have $S + 2\xi_{-} \ge 0$

**Case $\xi_{-} \ne 0$**

$$ \xi = -\xi_{-} = \text{thus} \begin{cases} y < 0 \\ \text{or } y > 1. \end{cases} \quad (47) $$

If $y < 0$, using $S = \bar{y} + y(1 - 2\bar{y})$, the following inequality is verified:

$$ \begin{array}{l} y \le S \le 1 - y \\ \Rightarrow \quad \begin{array}{l} Sy \ge y(1 - y) = \xi \\ \text{and} \\ S(1 - y) \ge y(1 - y) = \xi. \end{array} \\ \text{So} \quad S(y + 1 - y) = S \ge 2\xi = -2\xi_{-}. \end{array} \quad (48) $$

If $y > 1$, the same proof can be used in a symmetric way. Finally, we have $S + 2\xi_{-} \ge 0$ for $y \in \mathbb{R}$. $\square$
**Assumptions:**

(1) Relevant mass fraction on the boundary conditions: $\xi_{-}(\mathbf{x} \in \Gamma_{+}^{\xi}, t) = 0$.

(2) Mass fraction for the initial condition such that: $\xi_{-}(\mathbf{x} \in \Omega, t = 0) = 0$ and thus $E(t = 0) = 0$.

(3) $\frac{1}{\rho}[(1-2y)\underline{\nabla} \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \underline{\nabla}(1-2y)] \in \mathcal{L}^{\infty}(\Omega, [0, T])$.

(4) Equilibrium mass fraction such that: $\bar{y} \in [0, 1]$ implying $(S + 2\xi_{-}) \ge 0$ (proof of lemma below).

(5) Positive relaxation time scale $\tau > 0$ and density $\rho > 0$.

Thus we have

$$ \frac{dE(t)}{dt} \le - \int_{\Omega} \rho \frac{\xi_{-}^2}{2} \left[ \frac{1}{\rho} \left( (1-2y)\underline{\nabla} \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \underline{\nabla}(1-2y) \right) \right] d\Omega. \quad (45) $$

Assuming the assumptions previously presented,

$$ \frac{dE(t)}{dt} \le \left\| \frac{1}{\rho} \left[ (1-2y)\underline{\nabla} \cdot \mathbf{q}_r - \mathbf{q}_r \cdot \underline{\nabla}(1-2y) \right] \right\|_{\infty} E(t). \quad (46) $$

Using the Grönwall's inequality [12] and $E(0) = 0$, it enables to conclude that $E(t) = 0$ on $[0,T]$.

**Lemma 1.** If $\bar{y} \in [0, 1]$ then $S + 2\xi_{-} \ge 0$ with $S = \bar{y}(1 - y) + y(1 - \bar{y})$

**Proof.** Let's assume that $\bar{y} \in [0, 1]$.

**Case $\xi_{-} = 0$**

If $\xi_{-} = 0$, then $y \in [0, 1]$ implying that $S \ge 0$.

As a consequence, we have $S + 2\xi_{-} \ge 0$

**Case $\xi_{-} \ne 0$**

$$ \xi = -\xi_{-} = \text{thus} \begin{cases} y < 0 \\ \text{or } y > 1. \end{cases} \quad (47) $$

If $y < 0$, using $S = \bar{y} + y(1 - 2\bar{y})$, the following inequality is verified:

$$ \begin{array}{l} y \le S \le 1 - y \\ \Rightarrow \quad \begin{array}{l} S y \ge y(1 - y) = \xi \\ \text{and} \\ \text{So} \quad S(y + 1 - y) = S \ge 2\xi = -2\xi_{-}. \end{array} \end{array} \quad (48) $$

If $y > 1$, the same proof can be used in a symmetric way. Finally, we have $S + 2\xi_{-} \ge 0$ for $y \in \mathbb{R}$. $\square$
## Appendix B - Boundary conditions

As shown on Figure 3, we consider a cell i containing several faces shared with $j \in v(i)$ neighboring cells, a wall face (noted w) and an inlet/outlet face (noted with an index $i\infty$). The inlet/outlet face has a surface $S_{i\infty}$, a given valid mass fraction $y_{\infty} \in [0, 1]$ outside and an outward unit normal $n_{\infty}$. The flow rates on this face are noted

$$ 
\begin{aligned}
q_{i\infty} &= \mathbf{q}_{i\infty} \cdot \mathbf{n}_{\infty}, \\
(q_r)_{i\infty} &= (\mathbf{q}_r)_{i\infty} \cdot \mathbf{n}_{\infty}, \\
(q_g)_{i\infty} &= q_{i\infty} + (1 - y_{i\infty})(q_r)_{i\infty},
\end{aligned} 
\quad (49) $$

where

$$ y_{i\infty} = \begin{cases} y_i^n & \text{si } q_{i\infty} \ge 0, \\ y_{\infty} & \text{si } q_{i\infty} < 0. \end{cases} \quad (50) $$

Only one inlet/outlet face is considered here but the discussion is valid for several ones. For quantities on the inlet/outlet boundary face, the instants considered are not specified here for the sake of readability. These values are given as data, hence any consistent formula can be used, with $y_{\infty} \in [y_{\infty}^n, y_{\infty}^{n+1}]$. Conditions (20),(24) and (32) are computed in order to verify the discrete maximum principle for Cell i. Blue terms are terms added due to the inlet/outlet boundary face.

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Cell i with $(k, l, m)$ neighboring inner cells, a wall boundary face and an inlet/outlet boundary face.</figcaption></figure>
## Appendix B - Boundary conditions

As shown on Figure 3, we consider a cell i containing several faces shared with $j \in v(i)$ neighboring cells, a wall face (noted w) and an inlet/outlet face (noted with an index $i\infty$). The inlet/outlet face has a surface $S_{i\infty}$, a given valid mass fraction $y_{\infty} \in [0, 1]$ outside and an outward unit normal $n_{\infty}$. The flow rates on this face are noted

$$ 
\begin{align} 
q_{i\infty} &= \mathbf{q}_{i\infty} \cdot \mathbf{n}_{\infty}, \nonumber \\ 
(q_r)_{i\infty} &= (\mathbf{q}_r)_{i\infty} \cdot \mathbf{n}_{\infty}, \tag{49} \\ 
(q_g)_{i\infty} &= q_{i\infty} + (1 - y_{i\infty})(q_r)_{i\infty}, \nonumber 
\end{align} 
$$

where

$$ y_{i\infty} = \begin{cases} y_i^n & \text{si } q_{i\infty} \ge 0, \\ y_{\infty} & \text{si } q_{i\infty} < 0. \end{cases} \tag{50} $$

Only one inlet/outlet face is considered here but the discussion is valid for several ones. For quantities on the inlet/outlet boundary face, the instants considered are not specified here for the sake of readability. These values are given as data, hence any consistent formula can be used, with $y_{\infty} \in [y_{\infty}^n, y_{\infty}^{n+1}]$. Conditions (20),(24) and (32) are computed in order to verify the discrete maximum principle for Cell i. Blue terms are terms added due to the inlet/outlet boundary face.

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Cell i with $(k, l, m)$ neighboring inner cells, a wall boundary face and an inlet/outlet boundary face.</figcaption></figure>
5.1. Boundary condition for QG scheme

For Cell i from Figure 3, QG scheme is written

$$
\begin{align}
\rho_i^n \omega_i (y_i^{n+1} - y_i^n) - \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* y_i^{n+1} - \Delta t^n S_{i\infty} q_{i\infty} y_i^{n+1} & \nonumber \\
& + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_g)_{ij}^* \{ \mathrm{sg}_{ij}^g y_i^{n+1} + (1 - \mathrm{sg}_{ij}^g) y_j^{n+1} \} & \tag{51} \\
& + \Delta t^n S_{i\infty} (q_g)_{i\infty} \{ \mathrm{sg}_{i\infty}^g y_i^{n+1} + (1 - \mathrm{sg}_{i\infty}^g) y_\infty \} & \nonumber \\
& = \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. & \nonumber
\end{align}
$$

Using this scheme, coefficients arising from conditions (20) are

$$
\begin{align}
a_{ii} &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} \mathrm{sg}_{ij}^g - \frac{q_{ij}^*}{\rho_i^n} \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty}}{\omega_i} \left( \frac{(q_g)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty}^g - \frac{q_{i\infty}}{\rho_i^n} \right), \nonumber \\
a_{ij} &= \Delta t^n \frac{S_{ij}}{\omega_i} \frac{(q_g)_{ij}^*}{\rho_i^n} (1 - \mathrm{sg}_{ij}^g), \nonumber \\
b_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n - \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_g)_{i\infty}}{\rho_i^n} (1 - \mathrm{sg}_{i\infty}^g) y_\infty, \nonumber \\
\hat{b}_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty}}{\omega_i} \left( \frac{(q_g)_{i\infty}}{\rho_i^n} (y_\infty + (1 - y_\infty)\mathrm{sg}_{i\infty}^g) - \frac{q_{i\infty}}{\rho_i^n} \right). \tag{52}
\end{align}
$$

Assuming that $a_{ii} > 0$,

$$
\Lambda_i(\mathbf{A}) = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right) \\
+ \Delta t^n \frac{S_{i\infty}}{\omega_i} \left( \frac{(q_g)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty}^g - \frac{q_{i\infty}}{\rho_i^n} \right).
$$

The three conditions $a_{ii} > 0$, $\Lambda_i > 0$ and $\hat{b}_i \ge 0$ are not always fulfilled. Using the formula $(q_g)_{i\infty}(y_\infty + (1-y_\infty)\mathrm{sg}_{i\infty}^g) \le \mathrm{sg}_{i\infty}^g (q_g)_{i\infty}$, it can be stated that $a_{ii} \ge \Lambda_i \ge \hat{b}_i$. The only remaining condition on the time step is once again on $\hat{b}_i$. This condition, which is slightly different from (15), should be monitored as it can become the most constraining one.

5.2. Boundary condition for QRd scheme
5.1. Boundary condition for QG scheme

For Cell i from Figure 3, QG scheme is written

$$
\begin{align}
\rho_i^n \omega_i (y_i^{n+1} - y_i^n) - \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* y_i^{n+1} - \Delta t^n S_{i\infty} q_{i\infty} y_i^{n+1} & \nonumber \\
& + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_g)_{ij}^* \{ \operatorname{sg}_{ij}^g y_i^{n+1} + (1 - \operatorname{sg}_{ij}^g) y_j^{n+1} \} & \tag{51} \\
& + \Delta t^n S_{i\infty} (q_g)_{i\infty} \{ \operatorname{sg}_{i\infty}^g y_i^{n+1} + (1 - \operatorname{sg}_{i\infty}^g) y_\infty \} & \nonumber \\
& = \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. & \nonumber
\end{align}
$$

Using this scheme, coefficients arising from conditions (20) are

$$
\begin{align}
a_{ii} &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} \operatorname{sg}_{ij}^g - \frac{q_{ij}^*}{\rho_i^n} \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty}}{\omega_i} \left( \frac{(q_g)_{i\infty}}{\rho_i^n} \operatorname{sg}_{i\infty}^g - \frac{q_{i\infty}}{\rho_i^n} \right), \nonumber \\
a_{ij} &= \Delta t^n \frac{S_{ij}}{\omega_i} \frac{(q_g)_{ij}^*}{\rho_i^n} (1 - \operatorname{sg}_{ij}^g), \nonumber \\
b_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n - \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_g)_{i\infty}}{\rho_i^n} (1 - \operatorname{sg}_{i\infty}^g) y_\infty, \nonumber \\
\hat{b}_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty}}{\omega_i} \left( \frac{(q_g)_{i\infty}}{\rho_i^n} (y_\infty + (1 - y_\infty) \operatorname{sg}_{i\infty}^g) - \frac{q_{i\infty}}{\rho_i^n} \right). \tag{52}
\end{align}
$$

Assuming that $a_{ii} > 0$,

$$
\Lambda_i(\mathbf{A}) = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \sum_{j \in v(i)} \frac{S_{ij}}{\omega_i} \left( \frac{(q_g)_{ij}^*}{\rho_i^n} - \frac{q_{ij}^*}{\rho_i^n} \right) \\
+ \Delta t^n \frac{S_{i\infty}}{\omega_i} \left( \frac{(q_g)_{i\infty}}{\rho_i^n} \operatorname{sg}_{i\infty}^g - \frac{q_{i\infty}}{\rho_i^n} \right).
$$

The three conditions $a_{ii} > 0$, $\Lambda_i > 0$ and $\hat{b}_i \ge 0$ are not always fulfilled. Using the formula $(q_g)_{i\infty}(y_\infty + (1-y_\infty)\operatorname{sg}_{i\infty}^g) \le \operatorname{sg}_{i\infty}^g (q_g)_{i\infty}$, it can be stated that $a_{ii} \ge \Lambda_i \ge \hat{b}_i$. The only remaining condition on the time step is once again on $\hat{b}_i$. This condition, which is slightly different from (15), should be monitored as it can become the most constraining one.

5.2. Boundary condition for QRd scheme
For Cell *i* from Figure 3, QRd scheme is written

$$
\begin{align}
& \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \{ (1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1}) \} \nonumber \\
& \quad + \Delta t^n S_{i\infty} q_{i\infty} (1 - \mathrm{sg}_{i\infty}) (y_\infty - y_i^{n+1}) \nonumber \\
& \quad + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_{ij}^* \{ \mathrm{sg}_{ij}^r y_i^{n+1} (1 - y_j^n) + (1 - \mathrm{sg}_{ij}^r) y_j^n (1 - y_i^{n+1}) \} \nonumber \\
& \quad + \Delta t^n S_{i\infty} (q_r)_{i\infty} \{ \mathrm{sg}_{i\infty}^r y_i^{n+1} (1 - y_\infty) + (1 - \mathrm{sg}_{i\infty}^r) y_\infty (1 - y_i^{n+1}) \} \nonumber \\
&= \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \tag{53}
\end{align}
$$

Using this scheme, coefficients arising from conditions (24) are

$$
\begin{align}
(a_d)_{ii} &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij}^r - y_j^n) \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty}^r - y_\infty), \tag{54} \\
(a_d)_{ij} &= \Delta t^n \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}), \nonumber \\
(b_d)_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij}^r - 1) y_j^n \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty} - 1) y_\infty + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty}^r - 1) y_\infty, \tag{55} \\
(\hat{b}_d)_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} \mathrm{sg}_{ij}^r (1 - y_j^n) \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty} - 1) (1 - y_\infty) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} \mathrm{sg}_{i\infty}^r (1 - y_\infty), \nonumber \\
\Lambda_i(\mathbf{A}_d) &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij}^r - y_j^n) \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty}^r - y_\infty). \nonumber
\end{align}
$$

Once more, no condition on the time step arises for the QRd scheme. However, it should be
handled carefully for an outlet boundary (with respect to mixture flow rate $q_{i\infty} > 0$). Indeed,
whatever the sign of the relative flow rate ($q_r)_{i\infty$}, the flux requires a given value of $y_\infty$. If the
value of the outlet mass fraction $y_\infty$ is unknown, a different consistent formula should be used to
approximate $y_\infty \in [0, 1]$ (for instance $y_\infty = y_i^n$). Otherwise the flux should be modified.

5.3. Boundary condition for QRq scheme
For Cell i from Figure 3, QRd scheme is written

$$
\begin{align}
& \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \{ (1 - \operatorname{sg}_{ij}) (y_j^{n+1} - y_i^{n+1}) \} \nonumber \\
& \quad + \Delta t^n S_{i\infty} q_{i\infty} (1 - \operatorname{sg}_{i\infty}) (y_\infty - y_i^{n+1}) \nonumber \\
& \quad + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_{ij}^* \{ \operatorname{sg}_{ij}^r y_i^{n+1} (1 - y_j^n) + (1 - \operatorname{sg}_{ij}^r) y_j^n (1 - y_i^{n+1}) \} \nonumber \\
& \quad + \Delta t^n S_{i\infty} (q_r)_{i\infty} \{ \operatorname{sg}_{i\infty}^r y_i^{n+1} (1 - y_\infty) + (1 - \operatorname{sg}_{i\infty}^r) y_\infty (1 - y_i^{n+1}) \} \nonumber \\
&= \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \tag{53}
\end{align}
$$

Using this scheme, coefficients arising from conditions (24) are

$$
\begin{align}
(a_d)_{ii} &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (\operatorname{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\operatorname{sg}_{ij}^r - y_j^n) \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty}^r - y_\infty), \tag{54} \\
(a_d)_{ij} &= \Delta t^n \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (1 - \operatorname{sg}_{ij}), \nonumber \\
(b_d)_i &= y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\operatorname{sg}_{ij}^r - 1) y_j^n \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty} - 1) y_\infty + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty}^r - 1) y_\infty, \tag{55} \\
(\hat{b}_d)_i &= (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} \operatorname{sg}_{ij}^r (1 - y_j^n) \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty} - 1) (1 - y_\infty) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} \operatorname{sg}_{i\infty}^r (1 - y_\infty), \nonumber \\
\Lambda_i(\mathbf{A}_d) &= 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\operatorname{sg}_{ij}^r - y_j^n) \right) \nonumber \\
&\quad + \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} (\operatorname{sg}_{i\infty}^r - y_\infty). \nonumber
\end{align}
$$

Once more, no condition on the time step arises for the QRd scheme. However, it should be
handled carefully for an outlet boundary (with respect to mixture flow rate $q_{i\infty} > 0$). Indeed,
whatever the sign of the relative flow rate ($q_r$)$_{i\infty}$, the flux requires a given value of $y_\infty$. If the
value of the outlet mass fraction $y_\infty$ is unknown, a different consistent formula should be used to
approximate $y_\infty \in [0, 1]$ (for instance $y_\infty = y_i^n$). Otherwise the flux should be modified.

5.3. Boundary condition for QRq scheme
For Cell i from Figure 3, QRq scheme is written

$$
\begin{align}
& \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \{(1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1})\} \nonumber \\
& \quad + \Delta t^n S_{i\infty} q_{i\infty} (1 - \mathrm{sg}_{i\infty}) (y_\infty - y_i^{n+1}) \nonumber \\
& \quad + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_{ij}^* \left[ \mathrm{sg}_{ij} \{y_i^n (1 - y_i^{n+1}) (1 - \mathrm{sg}_{ij}^r) + y_i^{n+1} (1 - y_i^n) \mathrm{sg}_{ij}^r\} \right. \nonumber \\
& \qquad \left. + (1 - \mathrm{sg}_{ij}) \{y_j^n (1 - y_j^{n+1}) (1 - \mathrm{sg}_{ij}^r) + y_j^{n+1} (1 - y_j^n) \mathrm{sg}_{ij}^r\} \right] \nonumber \\
& \quad + \Delta t^n S_{i\infty} (q_r)_{i\infty} \left[ \mathrm{sg}_{i\infty} \{y_i^n (1 - y_i^{n+1}) (1 - \mathrm{sg}_{i\infty}^r) + y_i^{n+1} (1 - y_i^n) \mathrm{sg}_{i\infty}^r\} \right. \nonumber \\
& \qquad \left. + (1 - \mathrm{sg}_{i\infty}) \{y_\infty (1 - y_\infty)\} \right] \nonumber \\
&= \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \tag{55}
\end{align}
$$

Using this scheme, coefficients arising from conditions (32) are

$$
(a_q)_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} \mathrm{sg}_{ij} (\mathrm{sg}_{ij}^r - y_i^n) \right) \\
+ \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} \mathrm{sg}_{i\infty} (\mathrm{sg}_{i\infty}^r - y_i^n),
$$

$$
(a_q)_{ij} = \Delta t^n \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}) + \Delta t^n \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}) (\mathrm{sg}_{ij}^r - y_j^n),
$$

$$
(b_q)_i = y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij}^r - 1) (\mathrm{sg}_{ij} y_i^n + (1 - \mathrm{sg}_{ij}) y_j^n) \right) \\
+ \Delta t^n \frac{S_{i\infty}}{\omega_i} (\mathrm{sg}_{i\infty} - 1) \left( \frac{q_{i\infty} + (1 - y_\infty) (q_r)_{i\infty}}{\rho_i^n} \right) y_\infty \\
+ \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_r)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty} (\mathrm{sg}_{i\infty}^r - 1) y_i^n,
$$

$$
\begin{align*}
(\hat{b}_q)_i ={}& (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} \mathrm{sg}_{ij}^r (1 - [\mathrm{sg}_{ij} y_i^n + (1 - \mathrm{sg}_{ij}) y_j^n]) \right) \\
& + \Delta t^n \frac{S_{i\infty}}{\omega_i} (\mathrm{sg}_{i\infty} - 1) \frac{q_{i\infty} - y_\infty (q_r)_{i\infty}}{\rho_i^n} (1 - y_\infty) \\
& + \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_r)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty} \mathrm{sg}_{i\infty}^r (1 - y_i^n).
\end{align*}
$$

Assuming that $\forall j \in v(i), (a_q)_{ij} \le 0$,

$$
\begin{align}
\Lambda_i(\mathbf{A}_q) = 1 &+ \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} \left( \mathrm{sg}_{ij}^r - \left[ \mathrm{sg}_{ij} y_i^n + (1 - \mathrm{sg}_{ij}) y_j^n \right] \right) \right) \nonumber \\
&+ \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{q_{i\infty}}{\rho_i^n} (\mathrm{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_r)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty} (\mathrm{sg}_{i\infty}^r - y_i^n). \tag{56}
\end{align}
$$
For Cell *i* from Figure 3, QRq scheme is written

$$
\begin{align}
& \rho_i^n \omega_i (y_i^{n+1} - y_i^n) + \Delta t^n \sum_{j \in v(i)} S_{ij} q_{ij}^* \{(1 - \mathrm{sg}_{ij}) (y_j^{n+1} - y_i^{n+1})\} \nonumber \\
& \quad + \Delta t^n S_{i\infty} q_{i\infty} (1 - \mathrm{sg}_{i\infty}) (y_\infty - y_i^{n+1}) \nonumber \\
& \quad + \Delta t^n \sum_{j \in v(i)} S_{ij} (q_r)_i^* \left[ \mathrm{sg}_{ij} \{y_i^n (1 - y_i^{n+1})(1 - \mathrm{sg}_{ij}^r) + y_i^{n+1} (1 - y_i^n) \mathrm{sg}_{ij}^r\} \right. \nonumber \\
& \qquad \left. + (1 - \mathrm{sg}_{ij}) \{y_j^n (1 - y_j^{n+1})(1 - \mathrm{sg}_{ij}^r) + y_j^{n+1} (1 - y_j^n) \mathrm{sg}_{ij}^r\} \right] \nonumber \\
& \quad + \Delta t^n S_{i\infty} (q_r)_{i\infty} \left[ \mathrm{sg}_{i\infty} \{y_i^n (1 - y_i^{n+1})(1 - \mathrm{sg}_{i\infty}^r) + y_i^{n+1} (1 - y_i^n) \mathrm{sg}_{i\infty}^r\} \right. \nonumber \\
& \qquad \left. + (1 - \mathrm{sg}_{i\infty}) \{y_\infty (1 - y_\infty)\} \right] \nonumber \\
&= \Delta t^n \rho_i^n \omega_i \frac{\bar{y}_i^n - y_i^{n+1}}{\tau_i^n}. \tag{55}
\end{align}
$$

Using this scheme, coefficients arising from conditions (32) are

$$
(a_q)_{ii} = 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij} - 1) + \sum_{j \in v(i)} \frac{S_{ij} (q_r)_i^*}{\omega_i \rho_i^n} \mathrm{sg}_{ij} (\mathrm{sg}_{ij}^r - y_i^n) \right) \\
+ \Delta t^n \frac{S_{i\infty} q_{i\infty}}{\omega_i \rho_i^n} (\mathrm{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty} (q_r)_{i\infty}}{\omega_i \rho_i^n} \mathrm{sg}_{i\infty} (\mathrm{sg}_{i\infty}^r - y_i^n),
$$

$$
(a_q)_{ij} = \Delta t^n \frac{S_{ij} q_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}) + \Delta t^n \frac{S_{ij} (q_r)_{ij}^*}{\omega_i \rho_i^n} (1 - \mathrm{sg}_{ij}) (\mathrm{sg}_{ij}^r - y_j^n),
$$

$$
(b_q)_i = y_i^n + \frac{\Delta t^n}{\tau_i^n} \bar{y}_i^n + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_i^*}{\omega_i \rho_i^n} (\mathrm{sg}_{ij}^r - 1) (\mathrm{sg}_{ij} y_i^n + (1 - \mathrm{sg}_{ij}) y_j^n) \right) \\
+ \Delta t^n \frac{S_{i\infty}}{\omega_i} (\mathrm{sg}_{i\infty} - 1) \left( \frac{q_{i\infty} + (1 - y_\infty)(q_r)_{i\infty}}{\rho_i^n} \right) y_\infty \\
+ \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_r)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty} (\mathrm{sg}_{i\infty}^r - 1) y_i^n,
$$

$$
\begin{align*}
(\hat{b}_q)_i ={}& (1 - y_i^n) + \frac{\Delta t^n}{\tau_i^n} (1 - \bar{y}_i^n) + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_i^*}{\omega_i \rho_i^n} \mathrm{sg}_{ij}^r (1 - [\mathrm{sg}_{ij} y_i^n + (1 - \mathrm{sg}_{ij}) y_j^n]) \right) \\
& + \Delta t^n \frac{S_{i\infty}}{\omega_i} (\mathrm{sg}_{i\infty} - 1) \frac{q_{i\infty} - y_\infty (q_r)_{i\infty}}{\rho_i^n} (1 - y_\infty) \\
& + \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_r)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty} \mathrm{sg}_{i\infty}^r (1 - y_i^n).
\end{align*}
$$

Assuming that $\forall j \in v(i), (a_q)_{ij} \le 0$,

$$
\begin{align}
\Lambda_i(\mathbf{A}_q) ={}& 1 + \frac{\Delta t^n}{\tau_i^n} + \Delta t^n \left( \sum_{j \in v(i)} \frac{S_{ij} (q_r)_i^*}{\omega_i \rho_i^n} \left( \mathrm{sg}_{ij}^r - \left[ \mathrm{sg}_{ij} y_i^n + (1 - \mathrm{sg}_{ij}) y_j^n \right] \right) \right) \nonumber \\
& + \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{q_{i\infty}}{\rho_i^n} (\mathrm{sg}_{i\infty} - 1) + \Delta t^n \frac{S_{i\infty}}{\omega_i} \frac{(q_r)_{i\infty}}{\rho_i^n} \mathrm{sg}_{i\infty} (\mathrm{sg}_{i\infty}^r - y_i^n). \tag{56}
\end{align}
$$
The same conditions as (30) and (31) on the flow rates for neighboring inner cells again arise from condition $a_{ij}^q \le 0$. Similar conditions appear for the boundary face from condition $b_i^q \ge 0$ and $\hat{b}_i^q \ge 0$. When $q_{i\infty} < 0$, the inlet boundary flow rates must verify:

$$ 
\begin{aligned}
q_{i\infty} + (1 - y_{\infty})(q_r)_{i\infty} &\le 0, \\
q_{i\infty} - (q_r)_{i\infty}y_{\infty} &\le 0.
\end{aligned} 
\qquad (57) 
$$

When (57) is verified, no condition on the time step appears. Conditions (57) are automatically verified when a co-current flow is considered, because in this case, the mixture flow rate $q_{i\infty}$, the gas flow rate $(q_g)_{i\infty} = q_{i\infty} + (1 - y_{\infty})(q_r)_{i\infty}$ and the liquid flow rate $(q_l)_{i\infty} = q_{i\infty} - y_{\infty}(q_r)_{i\infty}$ have the same sign.

When (57) is not verified, a condition on the time step appears from condition $b_i^q \ge 0$ or from condition $\hat{b}_i^q \ge 0$.
The same conditions as (30) and (31) on the flow rates for neighboring inner cells again arise from condition $a_{ij}^q \le 0$. Similar conditions appear for the boundary face from condition $b_i^q \ge 0$ and $\hat{b}_i^q \ge 0$. When $q_{i\infty} < 0$, the inlet boundary flow rates must verify:

$$ 
\begin{aligned}
q_{i\infty} + (1 - y_{\infty})(q_r)_{i\infty} &\le 0, \\
q_{i\infty} - (q_r)_{i\infty}y_{\infty} &\le 0.
\end{aligned} 
\qquad (57) 
$$

When (57) is verified, no condition on the time step appears. Conditions (57) are automatically verified when a co-current flow is considered, because in this case, the mixture flow rate $q_{i\infty}$, the gas flow rate $(q_g)_{i\infty} = q_{i\infty} + (1 - y_{\infty})(q_r)_{i\infty}$ and the liquid flow rate $(q_l)_{i\infty} = q_{i\infty} - y_{\infty}(q_r)_{i\infty}$ have the same sign.

When (57) is not verified, a condition on the time step appears from condition $b_i^q \ge 0$ or from condition $\hat{b}_i^q \ge 0$.
References

[1] N. Zuber, J. Findlay, "Average Volumetric Concentration in Two-Phase Flow Systems", *J. Heat Transfer* **87** (1965), no. 4, p. 453-468, https://doi.org/10.1115/1.3689137.

[2] R. Eymard, G. Henry, R. Herbin, F. Hubert, R. Klöfkorn, G. Manzini, "3D Benchmark on Discretization Schemes for Anisotropic Diffusion Problems on General Grids", in *Finite Volumes for Complex Applications VI Problems & Perspectives* (Berlin, Heidelberg), Springer Berlin Heidelberg, 2011, p. 895-930.

[3] P. Frolkovič, "Maximum principle and local mass balance for numerical solutions of transport equation coupled with variable density flow", *Acta Mathematica Universitatis Comenianae. New Series* **67** (1998).

[4] L. Gastaldo, R. Herbin, J.-C. Latché, "An unconditionally stable finite element-finite volume pressure correction scheme for the drift-flux model", *ESAIM: Mathematical Modelling and Numerical Analysis* **44** (2010), no. 2 (en), p. 251-287, http://www.numdam.org/articles/10.1051/m2an/2010002/.

[5] L. Gastaldo, R. Herbin, J.-C. Latché, "A discretization of phase mass balance in fractional step algorithms for the drift-flux model", *IMA Journal of Applied Mathematics* **31** (2011), no. 1, p. 116-146, https://hal.science/hal-00308838.

[6] B. Larrouturou, "How to preserve the mass fractions positivity when computing compressible multi-component flows", Research Report RR-1080, INRIA, 1989, https://inria.hal.science/inria-00075479.

[7] K. Lipnikov, D. Svyatskiy, Y. Vassilevski, "Minimal stencil finite volume scheme with the discrete maximum principle", *Russian Journal of Numerical Analysis and Mathematical Modelling* **27** (2012), no. 4, p. 369-386, https://doi:10.1515/rnam-2012-0020.

[8] R. Lewandowski, B. Mohammadi, "Existence and positivity results for the $\varphi-\theta$ and a modified $k-\epsilon$ two-equation turbulence models", *Mathematical Models and Methods in Applied Sciences* **03** (1993), no. 02, p. 195-215, https://doi.org/10.1142/S0218202593000114.

[9] R. Eymard, T. Gallouët, R. Herbin, "Finite volume methods", in *Solution of Equation in Rn (Part 3), Techniques of Scientific Computing (Part 3)*, Handbook of Numerical Analysis, vol. 7, Elsevier, 2000, p. 713-1018, https://www.sciencedirect.com/science/article/pii/S1570865900070058.

[10] P. G. Ciarlet, "Discrete Maximum Principle for Finite-Difference Operators.", *Aequationes mathematicae* **4** (1970), p. 338-352, http://eudml.org/doc/136095.

[11] R. S. Varga, *Matrix Iterative Analysis*, Springer Berlin, Heidelberg, 2009.

[12] T. H. Gronwall, "Note on the Derivatives with Respect to a Parameter of the Solutions of a System of Differential Equations", *Annals of Mathematics* **20** (1919), no. 4, p. 292-296, http://www.jstor.org/stable/1967124.
References

[1] N. Zuber, J. Findlay, "Average Volumetric Concentration in Two-Phase Flow Systems", *J. Heat Transfer* **87** (1965), no. 4, p. 453-468, https://doi.org/10.1115/1.3689137.

[2] R. Eymard, G. Henry, R. Herbin, F. Hubert, R. Klöfkorn, G. Manzini, "3D Benchmark on Discretization Schemes for Anisotropic Diffusion Problems on General Grids", in *Finite Volumes for Complex Applications VI Problems & Perspectives* (Berlin, Heidelberg), Springer Berlin Heidelberg, 2011, p. 895-930.

[3] P. Frolkovič, "Maximum principle and local mass balance for numerical solutions of transport equation coupled with variable density flow", *Acta Mathematica Universitatis Comenianae. New Series* **67** (1998).

[4] L. Gastaldo, R. Herbin, J.-C. Latché, "An unconditionally stable finite element-finite volume pressure correction scheme for the drift-flux model", *ESAIM: Mathematical Modelling and Numerical Analysis* **44** (2010), no. 2 (en), p. 251-287, http://www.numdam.org/articles/10.1051/m2an/2010002/.

[5] L. Gastaldo, R. Herbin, J.-C. Latché, "A discretization of phase mass balance in fractional step algorithms for the drift-flux model", *IMA Journal of Applied Mathematics* **31** (2011), no. 1, p. 116-146, https://hal.science/hal-00308838.

[6] B. Larrouturou, "How to preserve the mass fractions positivity when computing compressible multi-component flows", Research Report RR-1080, INRIA, 1989, https://inria.hal.science/inria-00075479.

[7] K. Lipnikov, D. Svyatskiy, Y. Vassilevski, "Minimal stencil finite volume scheme with the discrete maximum principle", *Russian Journal of Numerical Analysis and Mathematical Modelling* **27** (2012), no. 4, p. 369-386, https://doi:10.1515/rnam-2012-0020.

[8] R. Lewandowski, B. Mohammadi, "Existence and positivity results for the $\varphi-\theta$ and a modified $k-\epsilon$ two-equation turbulence models", *Mathematical Models and Methods in Applied Sciences* **03** (1993), no. 02, p. 195-215, https://doi.org/10.1142/S0218202593000114.

[9] R. Eymard, T. Gallouët, R. Herbin, "Finite volume methods", in *Solution of Equation in Rn (Part 3), Techniques of Scientific Computing (Part 3)*, Handbook of Numerical Analysis, vol. 7, Elsevier, 2000, p. 713-1018, https://www.sciencedirect.com/science/article/pii/S1570865900070058.

[10] P. G. Ciarlet, "Discrete Maximum Principle for Finite-Difference Operators.", *Aequationes mathematicae* **4** (1970), p. 338-352, http://eudml.org/doc/136095.

[11] R. S. Varga, *Matrix Iterative Analysis*, Springer Berlin, Heidelberg, 2009.

[12] T. H. Gronwall, "Note on the Derivatives with Respect to a Parameter of the Solutions of a System of Differential Equations", *Annals of Mathematics* **20** (1919), no. 4, p. 292-296, http://www.jstor.org/stable/1967124.