%~Mouliné par MaN_auto v.0.32.0 2024-03-22 09:16:12 \documentclass[CRMECA, Unicode, XML, screen,biblatex, published]{cedram} %% %\TopicEN{Spontaneous articles} %\TopicFR{Articles spontanés} % \addbibresource{crmeca20230998.bib} \newcommand{\bbR}{\mathbb{R}} %%MATHBF \newcommand{\Q}{\mathbf{q}} \newcommand{\Qr}[0]{\bfq_r} \newcommand{\Qm}[0]{\bfq} \newcommand{\Qg}[0]{\bfq_g} \newcommand{\bfn}{\mathbf{n}} \newcommand{\bfu}{\mathbf{u}} \newcommand{\bfq}{\mathbf{q}} \newcommand{\bfx}{\mathbf{x}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfB}{\mathbf{B}} \newcommand{\clL}{\mathcal{L}} \newcommand{\ybin}[0]{\bar{y}_i^n} \newcommand{\qm}[0]{q_{ij}^*} \newcommand{\qr}[0]{(q_r)_{ij}^*} \newcommand{\qg}[0]{(q_g)_{ij}^*} \newcommand{\qminf}[0]{q_{i\infty}} \newcommand{\qrinf}[0]{(q_r)_{i\infty}} \newcommand{\qginf}[0]{(q_g)_{i\infty}} \newcommand{\yinf}[0]{y_\infty} \newcommand{\Dt}[0]{\Delta t^n} \newcommand{\Dx}[0]{\Delta x} \newcommand{\yin}[0]{y_i^n} \newcommand{\yjn}[0]{y_{j}^n} \newcommand{\yinp}[0]{y_i^{n+1}} \newcommand{\yjnp}[0]{y_{j}^{n+1}} \newcommand{\Min}[0]{\rho_{i}^n} %\newcommand{\Minp}[0]{\rho_{i}^{n+1}} \newcommand{\sg}[0]{\text{sg}_{ij}} \newcommand{\sgt}[0]{\text{sg}^r_{ij}} \newcommand{\sgg}[0]{\text{sg}^g_{ij}} \newcommand{\sgginf}[0]{\text{sg}^g_{i\infty}} \newcommand{\sginf}[0]{\text{sg}_{i\infty}} \newcommand{\sgtinf}[0]{\text{sg}^r_{i\infty}} \newcommand{\Oi}[0]{\omega_i} \newcommand{\Sij}[0]{S_{ij}} \newcommand{\Sinf}[0]{S_{i\infty}} \newcommand{\sumj}[0]{\sum_{j\,\in\,v(i)}} \newcommand{\xv}[0]{\mathbf{x}} \newcommand{\nv}[0]{\bfn_{ij}} \newcommand{\dunderline}[1]{\mathbf{#1}} \newcommand{\DIV}[1]{\underline{\nabla}\cdot #1} \newcommand{\DIVb}[1]{\underline{\nabla}\cdot\big(#1\big)} \newcommand{\DIVB}[1]{\underline{\nabla}\cdot\Big(#1\Big)} \newcommand{\inv}[1]{\frac{1}{#1}} \newcommand{\pdt}[1]{\frac{\partial #1}{\partial t}} \newcommand{\pdx}[1]{\frac{\partial}{\partial x}#1} \newcommand{\pder}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\der}[2]{\frac{d #1}{d #2}} \newcommand{\red}[2][red]{{\color{#1} #2}} \newcommand{\orange}[1]{{\color{orange} #1}} %\newcommand{\green}[1]{{\color{teal} #1}} \newcommand{\blue}[1]{{\color{blue} #1}} \newcommand{\eq}[1]{\begin{equation}#1\end{equation}} \newcommand{\eqsplit}[1]{ \begin{equation} \begin{split}#1 \end{split} \end{equation} } \newcommand{\dif}[1]{\underline{\nabla}{#1}} \newcommand{\difb}[1]{\underline{\nabla}\,\cdot{\big(#1\big)}} \newcommand{\difB}[1]{\underline{\nabla}{\,\cdot\Big(#1\Big)}} \newcommand{\dift}[1]{\partial_t #1} %\newcommand{\diftb}[1]{\partial_t \left(#1\right)} \newcommand{\qxi}[0]{\bfq_\xi} \newcommand{\xim}[0]{\xi_-} \newcommand{\ximc}[0]{\frac{\xim^2}{2}} \newcommand{\xip}[0]{\xi_+} %\newcommand{\ym}[0]{y_\xi} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \graphicspath{{./figures/}} \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} %\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red} \newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}} \newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel} \newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel} \newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel} \newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel} \let\oldhat\hat \renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}} \let\oldforall\forall \renewcommand*{\forall}{\mathrel{\oldforall}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \title{Maximum principle for the mass fraction in a system with two mass balance equations} \alttitle{Principe du maximum pour le titre massique d'un système incluant deux bilans de masse} \altkeywords{Principe du maximum, schéma Volumes~Finis, écoulement diphasique, déséquilibre en vitesse} \keywords{maximum principle, Finite Volume scheme, two-phase flow, non-equilibrium velocity} \author{\firstname{Gauthier} \lastname{Lazare\IsCorresp}\CDRorcid{0009-0000-8991-8851}} \email[G. Lazare]{gauthier.lazare@edf.fr} \address{EDF R\&D Chatou - 6 quai Waltier, 78400, Chatou, France.} \address{IRMA, UMR 7501, 7 rue Descartes, 67000 Strasbourg, France.} \thanks{The first author is supported by EDF R\&D and ANRT through CIFRE contract 2022/1027.} \author{\firstname{Qingqing} \lastname{Feng}} \email[Q. Feng]{qingqing.feng@edf.fr} \address[1]{EDF R\&D Chatou - 6 quai Waltier, 78400, Chatou, France.} \author{\firstname{Philippe} \lastname{Helluy}\CDRorcid{0000-0002-1077-1867}} \email[P. Helluy]{philippe.helluy@unistra.fr} \address[2]{IRMA, UMR 7501, 7 rue Descartes, 67000 Strasbourg, France.} \author{\firstname{Jean-Marc} \lastname{Hérard}\CDRorcid{0000-0001-9099-6072}} \email[J-M. Herard]{jean-marc.herard@edf.fr} \address[1]{EDF R\&D Chatou - 6 quai Waltier, 78400, Chatou, France.} \address{I2M - Institut de Mathématiques de Marseille, Aix Marseille Université, France.} \author{\firstname{Frank} \lastname{Hulsemann}\CDRorcid{0000-0001-5736-4532}} \email[F. Hulsemann]{frank.hulsemann@edf.fr} \address[1]{EDF R\&D Chatou - 6 quai Waltier, 78400, Chatou, France.} \author{\firstname{Stéphane} \lastname{Pujet}} \email[S. Pujet]{stephane.pujet@edf.fr} \address[1]{EDF R\&D Chatou - 6 quai Waltier, 78400, Chatou, France.} \begin{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. \end{abstract} \begin{altabstract} 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é. \end{altabstract} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} \maketitle %\newpage \section{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~\cite{ZuberFindlay}. The maximum principle for the mass fraction has been widely investigated, either for diffusive problems (see among others~\cite{fvca}) or in hybrid convection-diffusion problems (see among others~\cite{Frolkovic, gastaldo2, gastaldo:hal-00308838, larrouturou:inria-00075479, LipnikovSvyatskiyVassilevsk}). Most of the time, a null relative velocity is considered in the convective flux model. In the sequel, we focus on a system involving two mass balance equations with a relative velocity (see also~\cite{gastaldo2, gastaldo:hal-00308838}) and a source term. Once the continuous maximum principle has been examined in Section~\ref{sec:1}, 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~\ref{sec:2}. Eventually, numerical simulations were used to assess the accuracy of the schemes described in Section~\ref{sec:3}. The proposed schemes can be used in a broader framework, such as two-phase flow models. \section{Governing equations for the model}\label{sec:1} \subsection{System of equations} A flow characterized by a density $\rho$ and a flow rate $\Qm$ is studied on a domain $\Omega$, over a period $[0,T]$. The boundary of $\Omega$ is noted $\Gamma$. The outward unit normal is noted $\bfn_{\Gamma}$. Then \begin{equation} \Gamma= \Gamma_+ \cup \Gamma_- \cup\Gamma_w, \end{equation} where : \begin{equation} \begin{split} \Gamma_+ &= \left\{\xv\in\Gamma, \Qm(\xv,t)\cdot\bfn_{\Gamma} < 0\right\},\\ \Gamma_- &= \left\{\xv\in\Gamma, \Qm(\xv,t)\cdot\bfn_{\Gamma} > 0\right\}, \\ \Gamma_w &= \left\{\xv\in\Gamma, \Qm(\xv,t)\cdot\bfn_{\Gamma} = 0\right\}. \end{split} \end{equation} 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 $\Qg = \Qm + (1-y)\Qr$, where the relative velocity between phases $\bfu_r$~\cite{ZuberFindlay} is used in the relative flow rate $\Qr=\rho \bfu_r$. The mixture velocity is defined by $\bfu := \Qm/\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 \hspace{-0.05cm}y$ follows a conservation law that is associated with the total mass conservation law. These equations are provided below (see also~\cite{gastaldo:hal-00308838} for a similar system). \begin{equation}\label{eq:equations_continues} \begin{split} \pdt{\rho}+\DIV{\Qm} &= 0, \\ \pdt{(\rho y)} + \DIVb{y\Qm}+\DIVB{y\left(1-y\right)\Qr}&= \displaystyle\rho\frac{\bar{y}-y}{\tau}. \end{split} \end{equation} Given the definition of the relative flow rate, the boundary $\Gamma$ can also be split according to the sign of the flow rate \mbox{$\qxi = \Qm+(1-2y)\Qr$}: \begin{equation} \Gamma = \Gamma_+^\xi\cup \Gamma_-^\xi\cup\Gamma_w^\xi, \end{equation} where : \begin{equation} \begin{split} \Gamma_+^\xi & = \left\{\xv\in\Gamma,~\bfq_\xi\cdot\bfn_\Gamma < 0\right\}, \\ \Gamma_-^\xi & = \left\{\xv\in\Gamma,~\bfq_\xi\cdot\bfn_\Gamma > 0\right\},\\ \Gamma_w^\xi & = \left\{\xv\in\Gamma,~\bfq_\xi\cdot\bfn_\Gamma = 0\right\}. \end{split} \end{equation} The four parameters $\Qm,\Qr,\bar{y}$ and $\tau$ are given functions: \begin{itemize} \item $\Qm(\xv,t)$: the mixture flow rate, given by a momentum conservation equation or provided as input data, \item $\Qr(\xv,t)$: the relative flow rate, provided as input data. It is often given with a drift-flux closure law related to the other parameters, \item $\bar{y}(\xv,t)\in[0,1]$: the equilibrium mass fraction, reached after a characteristic relaxation time $\tau(\xv,t)>0$. $\bar{y}$ and $\tau$ are obtained through closure laws. \end{itemize} This system of equations must be completed by the initial conditions for the mass fraction and density:~\mbox{$y(\xv,t=0) = y_0(\xv),$} \mbox{$ \rho(\xv,t=0) = \rho_0(\xv)$}, and by suitable boundary conditions for the mass fraction (and also for the density on $\Gamma_+$). \subsection{Continuous maximum principle} \paragraph{Property 1 (Continuous maximum principle):} For the density $\rho(\xv,t)>0$, assume initial condition such that $y(\xv,t=0)\in [0,1]$ and boundary conditions such that $y(\xv\in \Gamma_+^{\xi},t)\in[0,1]$, with $\Gamma_+^{\xi}=\{\xv\in\Gamma,(\Qm+(1-2y)\Qr)\cdot\bfn_\Gamma < 0\}$. Consider closure laws for source terms such that $\bar{y}(\xv,t)\in[0,1]$ and $\tau(\xv,t)>0$. If the quantity $((1-2y)\DIV{\Qr}+2\Qr\cdot\underline{\nabla}y)/\rho$ is bounded on $\Omega\times[0,T]$, then the mass fraction $y(\xv,t)$ solution of~\eqref{eq:equations_continues} lies in $[0,1]$ on $\Omega\times[0,T]$. %{\flushright$\square$\\ } Using the same methodology as~\cite{gastaldo:hal-00308838} and~\cite{lewandowski}, the proof of this property is given in Appendix~\ref{an:pmc}. We recall that the positivity of the density is ensured as long as the divergence of the mixture velocity \emph{u} is bounded and the incoming flow has a positive density $\rho(\bfx\in\Gamma_+,t)>0$. In the following, the focus shifts to the mass fraction governing equation with unknown $y$. A discrete Finite Volume scheme for the density $\rho$ will be assumed (see equation~\eqref{eq:discr_masse}). \begin{rema} To ensure the continuous maximum principle for the mass fraction $y$ and 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 consider co-current flows entering the domain; hence, the phase velocities have the same sign, such that \begin{equation} \mathbf{u_g}\cdot\mathbf{u_l}\geq 0 \text{ with } \left\{ \begin{array}{rcl} \bfu_g &=& \bfu + (1-y)\mathbf{u_r}\\ \mathbf{u_l} &=& \bfu -y\mathbf{u_r} \end{array}\right. . \end{equation} Hence, noting that \begin{equation} \begin{split} \bfq &= \rho\big(y\bfu_g+(1-y)\bfu_l\big),\\ \mathbf{q_\xi} &= \rho\big((1-y)\bfu_g+y\bfu_l\big), \end{split} \end{equation} we obtain at once $\Qm_\xi\cdot\Qm = \rho^2\big(y(1-y)(\bfu_g-\bfu_l)^2+ \bfu_g\cdot \bfu_l\big)\geq 0$ but also $\Gamma_w^\xi=\Gamma_w$, $\Gamma_+^\xi=\Gamma_+$ and $\Gamma_-^\xi=\Gamma_-$. \end{rema} \section{Numerical method}\label{sec:2} The domain $\Omega$ is discretised in $N$ cells: $\Omega=\displaystyle\cup_{i\,\in\,\llbracket 1,N \rrbracket} \Omega_i$. The time interval $[0,T]$ is also discretised in $N_t$ intervals so that: \begin{equation} t^0 = 0~;~\forall n \in [0,N_t], t^{n+1} = t^n + \Dt~\text{and}~T = \sum_{i\,\in\,[0,N_t]} \Dt. \end{equation} 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 is decided according to the desired properties of the scheme. The following notations are used: \begin{itemize} \item $\theta_i^n$ the approximate value of $\theta$ on cell i at time step $n$, \item $\Oi$ the volume of cell $i$, \item $j\in v(i)$ the neighbors of cell $i$, \item $S_{ij}$ the surface of the intersection of cell $i$ with its neighbor cell $j$, \item $\nv$ the normal unit vector of the surface $\Sij$ outward cell $i$, \item $\qm =\Qm_{ij}^*\cdot\nv$ the normal mixture flow rate between cells $i$ and $j$, \item $\qr=(\Qr)_{ij}^*\cdot\nv$ the normal relative flow rate between cells $i$ and $j$, \item $\qg = \qm+(1-y_{ij}^n)\qr$ the normal gas flow rate between cells $i$ and $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 $i\hspace{-0.08cm}j$ is detailed in the sequel. \end{itemize} Here, $\Qm^*_{ij}$ (respectively $(\Qr)^*_{ij}$) is an estimation of the flow rate $\Qm$ (respectively $\Qr$) at interface $i\hspace{-0.08cm}j$ at time $t^* \in [t^n,t^{n+1}]$. \subsection{Finite Volume~Discretisation} Finite Volume methods are well adapted to treat conservation issues and maintain physical properties within valid bounds. Focusing on the system of mass balances~\eqref{eq:equations_continues}, we recall that~\cite{gastaldo:hal-00308838} introduced a \emph{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 first propose a \emph{linear implicit} scheme with some restrictions on the time step to comply with the discrete maximum principle. Then, two other \emph{linear implicit} schemes satisfying the maximum principle without any condition on the time step are presented. The Finite Volume scheme for total mass conservation is assumed to be \begin{equation}\label{eq:discr_masse} \Oi\left(\rho_i^{n+1}-\rho_i^n\right) + \Dt \sumj \Sij\qm = 0 . \end{equation} Turning to the mass fraction equation, the following choices were made : \begin{itemize} \item The unsteady term is decomposed into two parts: \begin{equation} \begin{split} \Dt\int_{\Oi} \pdt{\rho y} d\Omega &= \Big(\left(\rho y\right)^{n+1}_i-\left(\rho y\right)^n_i\Big)\Oi \\ &= \rho_i^n \left(\yinp-\yin\right)\Oi + \yinp \underbrace{\left(\rho_i^{n+1}-\rho_i^n\right)\Oi}_{\displaystyle=-\Dt\sumj \Sij\qm}\hspace{-0.3cm}\text{, using~\eqref{eq:discr_masse}.} \end{split} \end{equation} \item The mass fraction is implicit in the source term \begin{equation} \int_{\Oi} \rho \frac{\bar{y}-y}{\tau} d\Omega\approx \rho_i^\# \frac{\bar{y}_i^\#-y_i^{n+1}}{\tau_i^\#}\Oi \end{equation} where $t^\# \in [t^n,t^{n+1}]$ such that \begin{equation} \left\{ \begin{split}\rho_i^\# &> 0, \\ \bar{y}^\#_i &\in [0,1], \\ \tau_i^\# &> 0. \end{split}\right. \end{equation} In the sequel, we will use: $\rho_i^\# = \rho_i^n, \bar{y}^\#_i = \bar{y}^n_i$ and $ \tau_i^\#=\tau_i^n$. \item Turning to convection contributions, two different methods are used. The first method simply considers the total convection: $\DIVb{y\Qm}+\DIVb{y(1-y)\Qr} = \DIVb{y\Qg}$. 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 mixture convection and relative convection separately. Two different schemes are detailed: the QRd 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 within physical bounds $[0,1]$ regardless of the time step. \end{itemize} In the sequel, the following notation is used for the sign of \[ z: \displaystyle\text{sg}(z) = \left\{ \begin{array}{rl} 1 & \text{ if } z\geq 0 \\ 0 & \text{ if } z<0 \end{array}\right.\,. \] To simplify notations, \[ \sg := \text{sg}\left(\qm\right), \sg^g := \text{sg}\left(\qg\right)\; \text{and}\; \sgt := \text{sg}\left(\qr\right)\; \text{are used}. \] \subsection{Global Scheme with QG Scheme} \subsubsection{Definition of QG scheme} For gas flow rate convection, a standard implicit upwind scheme is used~\cite{EYMARD2000713}. With QG scheme, the global scheme is written with the sign convention for each cell $\Omega_i$: \begin{equation}\label{eq:discretisation_generale_qg} \begin{split} &\Min\Oi\left(\yinp-\yin\right) -\Dt\sumj\Sij\qm~\yinp\\ &+\Dt\sumj\Sij\qg\left\{\sgg\yinp+\left(1-\sgg\right)\yjnp\right\} = \Dt\Min\Oi\frac{\ybin-\yinp}{\tau^n_i}. \end{split} \end{equation} The mass fraction $y_{ij}^n$ at interface $i\hspace{-0.08cm}j$ involved in the gas flow rate $\qg=\qm+(1-y_{ij}^n)\qr$ is approximated using an upwind scheme based on the mixture flow rate $\qm$: \begin{equation} y_{ij}^n\approx\left\{ \begin{array}{rcl} y_i^n &\text{ if }& \qm\geq 0,\\ y_j^n &\text{ if }& \qm<0. \end{array}\right. \end{equation} Other consistent choices can be made as long as $y_{ij}^n\in [0,1]$. %\newpage \subsubsection{Discrete maximum principle}\label{sec:pm_qg} We now examine whether if the discrete mass fraction remains within $[0,1]$ on the $N$ cells using this scheme. \paragraph{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\llbracket 1,N \rrbracket$. If the initial conditions are that $\forall k\in \llbracket 1,N \rrbracket, y_k^n\in[0,1]$, then the global scheme with QG ensures that $\yinp$ remains in $[0,1]$, $\forall i \in\llbracket 1,N \rrbracket$, when the time step $\Dt$ verifies, $\forall i \in \llbracket 1,N \rrbracket$: \begin{equation}\label{eq:condition_qg_dt} \left(1-\yin\right)+ \left(1-\bar{y}^n_i\right)\frac{\Dt}{\tau} + \frac{\Dt}{\rho_i^n \Oi} \sumj \Sij \left(1-y_{ij}^n\right)\qr \geq 0 . \end{equation} %{\flushright$\square$\\ } \begin{rema} The condition on the time step is automatically verified when $\Qr = \mathbf{0}$. \end{rema} \begin{proof} Using equation~\eqref{eq:discretisation_generale_qg} obtained for each cell $i\in \llbracket 1,N \rrbracket$, the system can be expressed in matrix form. \begin{equation} \dunderline{A}~\bfY^{n+1} = \bfB\text{ with } \forall (i,j)\in \llbracket 1,N\rrbracket^2,~\left\{ \begin{array}{rlc} \dunderline{A}&=&\left(a_{ij}\right), \\ \bfB &=&\left(b_i\right), \\ \bfY^n & = & \left(\yin\right). \end{array}\right. \end{equation} The discrete system for $\widehat{\bfY} =\mathbf{1}-\bfY$ can also be expressed as follows: \begin{equation} \dunderline{\widehat{A}}~\widehat{\bfY}^{n+1}=\widehat{\bfB} =\dunderline{A}\times\mathbf{1}-\bfB \text{ with } \forall (i,j)\in\llbracket 1,N\rrbracket^2, \left\{ \begin{array}{rlcll} \mathbf{\widehat{B}} &=& \left(\widehat{b}_i\right) &=& \left(\sum_{j} a_{ij}-b_i\right), \\ ~\dunderline{\widehat{A}} &=& \left(\widehat{a}_{ij}\right) &=& \left(a_{ij}\right), \\ ~\mathbf{\widehat{Y}}^{n} & =& \left(\widehat{y}_i^n\right) &=& \left(1-\yin\right). \end{array}\right. \end{equation} We also introduce the quantity $\Lambda_i(\dunderline{A})$ on each cell by $ \Lambda_i(\dunderline{A}) = |a_{ii}|-\sumj|a_{ij}|$. Note that $\Lambda_i(\dunderline{A}) = \Lambda_i(\dunderline{\widehat{A}})$. The coefficients of the matrix system and the quantity $\Lambda_i$ are \begin{equation}\label{sys:qg} \begin{split} a_{ii} &= \widehat{a}_{ii} = 1+\frac{\Dt}{\tau_i^n} +\Dt\sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}\sgg-\frac{\qm}{\rho_i^n}\right),\\ a_{ij} &= \widehat{a}_{ij} = \Dt~\frac{\Sij}{\Oi}\frac{\qg}{\rho_i^n} \left(1-\sgg\right), \\ b_{i} &= \yin + \frac{\Dt}{\tau_i^n}\ybin, \\ \widehat{b}_i &= \left(1-\yin\right)+\frac{\Dt}{\tau_i^n}\left(1-\ybin\right) +\Dt\sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}-\frac{\qm}{\rho_i^n}\right). \\ \text{Assuming that } &a_{ii}>0, \\ \Lambda_i(\dunderline{A}) & = 1 + \frac{\Dt}{\tau_i^n} + \Dt\sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}-\frac{\qm}{\rho_i^n}\right). \end{split} \end{equation} \begin{rema} These two formulations are equivalent: \begin{equation} \sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}-\frac{\qm}{\rho_i^n}\right) = \sumj\frac{\Sij}{\Oi}\left(1-y_{ij}^n\right)\qr. \end{equation} \end{rema} The conditions $\forall (i,j) \in \llbracket 1,N \rrbracket^2, a_{ii}>0, a_{ij}\leq 0$ and $\Lambda_i> 0$ ensure that matrix $\dunderline{A}$ is invertible and its inverse has positive coefficients (see~\cite{Ciarlet1970,Varga}). To satisfy the maximum principle, vectors $\bfB$ and $\mathbf{\widehat{B}}$ also need to fulfill $b_i \geq 0$ and $\widehat{b}_i\geq 0$ with $i\in \llbracket 1,N\rrbracket$. To summarize, the conditions for the scheme are \begin{equation}\label{properties_qg} \forall i\in \llbracket 1,N \rrbracket,\left\{ \begin{array}{ccccccl} a_{ii} &>& 0 &,& a_{ij} &\leq& 0, j \in v(i),\\ \Lambda_i(\dunderline{A}) &>& 0, && && \\ b_i &\geq & 0 &,& \widehat{b}_i &\geq& 0. \end{array}\right. \end{equation} For $\tau_k^n>0$, $\bar{y}_k^n\in[0,1], k \in\llbracket 1,N \rrbracket$ and when $\yin\in[0,1], i\in\llbracket 1,N \rrbracket$, the maximum principle for the mass fraction with QG scheme is satisfied when properties~\eqref{properties_qg} are fulfilled. Actually, the condition on $(\hat{b}_i)$ is equivalent to~\eqref{eq:condition_qg_dt}. If condition~\eqref{eq:condition_qg_dt} is verified, the coefficients $(a_{ii})$ and $(\Lambda_i)$ are positive as $a_{ii}\geq\Lambda_i\geq\widehat{b}_i$. The other two conditions $a_{ij}\leq 0$ and $b_i\geq 0$ are always satisfied. \end{proof} \subsection{Global Scheme with QRd scheme} \subsubsection{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~\eqref{eq:discretisation_generale_qg}. The global scheme with QRd scheme is written: \begin{equation}\label{eq:discretisation_generale_qrd} \begin{split} \Min\Oi\left(\yinp-\yin\right) &+\Dt\sumj\Sij\qm\left\{\left(1-\sg\right)\left(\yjnp-\yinp\right)\right\}\\ &+\Dt\sumj\Sij\qr\left\{\sgt\yinp\left(1-\yjn\right)+\left(1-\sgt\right)\yjn\left(1-\yinp\right)\right\} \\ &= \Dt\Min\Oi\frac{\ybin-\yinp}{\tau^n_i}. \end{split} \end{equation} \subsubsection{Discrete maximum principle}\label{sec:pm_qrd} Now we examine whether the discrete mass fraction remains within $[0,1]$ on the $N$ cells using this scheme. \paragraph{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\llbracket 1,N \rrbracket$. If the initial conditions are such that $\forall k\in \llbracket 1,N \rrbracket, y_k^n\in[0,1]$, then the global scheme with QRd ensures that $\yinp$ remains in $[0,1]$, $\forall i \in\llbracket 1,N \rrbracket$, whatever the time step $\Dt$ is. %{\flushright$\square$\\ } %\emph{Remark :} \begin{rema} The scheme is rigorously conservative in space only for the discrete steady states ($\yinp = \yin, i\in \llbracket 1,N \rrbracket$). This is not an issue in practice, when the method is applied to the computations of steady states. \end{rema} \begin{proof} Using the same methodology as for QG scheme, a similar discrete system (noted $d$) is written from~\eqref{eq:discretisation_generale_qrd} for QRd scheme: \begin{equation} \dunderline{A}_d~\bfY^{n+1} = \bfB_d\text{ with } \forall (i,j)\in \llbracket 1,N\rrbracket^2,~\left\{ \begin{array}{rlc}\dunderline{A}_d&=&\big((a_d)_{ij}\big), \\ \bfB_d &=&\big((b_d)_i\big), \\ \bfY^n & = & \left(\yin\right). \end{array}\right. \end{equation} For QRd scheme, the coefficients of the matrix system and the quantity $\Lambda_i$ are \begin{equation}\label{sys:qrd} \begin{split} (a_d)_{ii} &= (\widehat{a}_d)_{ii} = 1+\frac{\Dt}{\tau_i^n} +\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n}\left(\sg-1\right) +\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-\yjn\right)\right),\\ (a_d)_{ij} &= \left(\widehat{a}_d\right)_{ij} = \Dt~\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n} \left(1-\sg\right), \\ (b_d)_{i} &= \yin + \frac{\Dt}{\tau_i^n}\ybin + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-1\right)\yjn\right), \\ \left(\widehat{b}_d\right)_i &= \left(1-\yin\right)+\frac{\Dt}{\tau_i^n}\left(1-\ybin\right) +\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\sgt\left(1-\yjn\right) \right), \\ \Lambda_i\left(\dunderline{A}_d\right) & = 1 + \frac{\Dt}{\tau_i^n} + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-\yjn\right)\right). \end{split} \end{equation} For given values $q_{kl}^n$ and $(q_r)_{kl}^n, (k,l)\in \llbracket 1,N \rrbracket^2$, and for $\tau_k^n>0$, $\bar{y}_k^n\in[0,1], k \in\llbracket 1,N \rrbracket$, the following conditions are fulfilled when $\yin\in[0,1], i\in\llbracket 1,N \rrbracket$: \begin{equation}\label{properties_qrd} \forall i\in \llbracket 1,N \rrbracket,\left\{ \begin{array}{ccccccl} (a_d)_{ii} &>& 0 &,& (a_d)_{ij} &\leq& 0, j \in v(i),\\ \Lambda_i(\dunderline{A}_d) &>& 0, && && \\ (b_d)_i &\geq & 0 &,& \left(\widehat{b}_d\right)_i &\geq& 0, \end{array}\right. \end{equation} \emph{whatever the time step $\Dt$ is.} Once again, using the results~\cite{Ciarlet1970,Varga}, the matrix $\dunderline{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 $\Dt$ is. \end{proof} \subsection{Global Scheme with QRq scheme} \subsubsection{Definition of QRq scheme} QRq scheme not only considers 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{multline}\label{eq:discretisation_generale_qrq} \Min\Oi\left(\yinp-\yin\right) +\Dt\sumj\Sij\qm\left\{(1-\sg)\left(\yjnp-\yinp\right)\right\}\\ \begin{aligned} &+\Dt\sumj\Sij\qr{\footnotesize \left[ \begin{split} &\sg \left\{\yin \left(1-\yinp\right)\left(1-\sgt\right) + \yinp \left(1-\yin\right)\sgt\right\} \\ &+\left(1-\sg\right)\left\{\yjn \left(1-\yjnp\right)\left(1-\sgt\right) + \yjnp \left(1-\yjn\right)\sgt\right\} \end{split}\right]} \\ &= \Dt\Min\Oi\frac{\ybin-\yinp}{\tau^n_i}. \end{aligned} \end{multline} %\newpage \subsubsection{Discrete maximum Principle} \paragraph{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\llbracket 1,N \rrbracket$. If the initial conditions are such that $ k\in \llbracket 1,N \rrbracket, y_k^n\in[0,1]$, then QRq scheme ensures that $\yinp$ remains in $[0,1]$, $ i \in\llbracket 1,N \rrbracket$, whatever the time step $\Dt$ is, provided that the flow rates fulfill the following conditions, when $\qm<0$: \begin{itemize} \item If $\qr\geq 0$ \begin{equation} \qm + \qr\left(1-\yjn\right)\leq 0. \end{equation} \item Otherwise, if $\qr< 0$ \begin{equation} \qm-\qr\yjn\leq 0. \end{equation} \end{itemize} No condition arises when $\qm\geq 0$. %{\flushright$\square$\\ } For a co-current flow, these conditions are automatically verified. \begin{rema} Once again, the scheme is rigorously conservative in space only for steady states $(\yinp = \yin, i \in \llbracket 1,N \rrbracket)$. \end{rema} \begin{proof} Using the same methodology as for QRd scheme, the discrete system from~\eqref{eq:discretisation_generale_qrq} is written as \begin{equation} \dunderline{A}_q~\bfY^{n+1} = \bfB_q\text{ with } \forall (i,j)\in \llbracket 1,N \rrbracket^2,~\left\{ \begin{array}{rlc}\dunderline{A}_q &=\big((a_q)_{ij}\big)~, \\ \bfB_q &=\big((b_q)_i\big)~, \\ \bfY^n & = \big(\yin\big). \end{array}\right. \end{equation} Using blue to denote the terms linked to QRq scheme, the coefficients of the matrix system and the quantity $\Lambda_i(\dunderline{A}_q)$ are \begin{equation} \begin{split} (a_q)_{ii} &= 1 +\frac{\Dt}{\tau_i^n}+\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n}\left(\sg-1\right) +\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\sg\left(\sgt-\yin\right)\right),\\ (a_q)_{ij} &= \Dt~\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n} \left(1-\sg\right) + \Dt\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(1-\sg\right)\left(\sgt-\yjn\right),\\ (b_q)_{i} &= \yin +\frac{\Dt}{\tau_i^n}\ybin + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-1\right)\left(\sg\yin+\left(1-\sg\right)\yjn\right)\right), \\ \left(\widehat{b}_q\right)_i &= \left(1-\yin\right)+\frac{\Dt}{\tau_i^n}\left(1-\ybin\right)+\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\sgt\left(1-\left[\sg\yin+\left(1-\sg\right)\yjn)\right]\right) \right). \\ \text{Assuming that } &\forall j\in v(i), (a_q)_{ij}<0, \\ \Lambda_i\left(\dunderline{A}_q\right) &= 1 + \frac{\Dt}{\tau_i^n} + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-\left[\sg\yin+\left(1-\sg\right)\yjn\right]\right)\right). \end{split} \end{equation} Examining the sign of these coefficients raises two conditions for the flow rates. \begin{itemize} \item When $\qm < 0$ and $\qr\geq 0$ \begin{equation}\label{eq:condition_1} \qm + \qr\left(1-\yjn\right)\leq 0. \end{equation} \item When $\qm < 0$ and $\qr< 0$ \begin{equation}\label{eq:condition_2} \qm-\qr\yjn\leq 0. \end{equation} \end{itemize} When $q_{kl}^n$ and $(q_r)_{kl}^n, (k,l)\in \llbracket 1,N \rrbracket^2$ satisfy the previous conditions~\eqref{eq:condition_1}-\eqref{eq:condition_2}, and if $\tau_k^n>0$, $\bar{y}_k^n\in[0,1], k \in\llbracket 1,N \rrbracket$, the following properties are verified when $\yin\in[0,1], i\in\llbracket 1,N \rrbracket$: \begin{equation}\label{properties_qrq} \forall i\in \llbracket 1,N \rrbracket,\left\{ \begin{array}{ccccccl} \left(a_q\right)_{ii} &>& 0 &,& \left(a_q\right)_{ij} &\leq& 0, j \in v(i),\\ \Lambda_i\left(\dunderline{A}_q\right) &>& 0, && && \\ \left(b_q\right)_i &\geq & 0 &,& \left(\widehat{b}_q\right)_i &\geq& 0, \end{array}\right. \end{equation} \emph{whatever the time step $\Dt$ is.} Once more, using common results of~\cite{Ciarlet1970, Varga}, the matrix $\dunderline{A}_q$ is invertible and its inverse has positive coefficients. Hence, the maximum principle for the mass fraction with QRq scheme is verified regardless of the time step $\Dt$, as long as the conditions~\eqref{eq:condition_1}-\eqref{eq:condition_2} on the flow rates are respected. \end{proof} \subsection{Boundary Conditions} These Finite Volume cell schemes presented above are only valid for cells that do not contain faces on the boundary. The cells (noted $i$ here) containing at least one face on the boundary are presented in Appendix~\ref{appb}. Boundary conditions are considered for three schemes with a given valid mass fraction $\yinf\in[0,1]$ outside. The inlet/outlet flow rates for the faces on the boundary are noted $\qminf$ and $\qrinf$. We summarise here the main conclusions of Appendix~\ref{appb}. 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 considered because 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 the outlet face (with respect to the mixture flow rate $\qminf> 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. \section{Numerical results}\label{sec:3} To be able to compare the numerical solution of the different schemes obtained using an analytical solution of~\eqref{eq:equations_continues}, a steady one-dimensional test case is considered. The physical parameters $(\rho,q,q_r,\tau,\bar{y})$ are uniform and constant, and are noted with an exponent $0$. The simplified system~is \begin{equation}\label{eq:4eme_equation_1D_simplifiee} \left\{ \begin{split} &\rho = \rho^0, q = q^0, \\ &\pdx{\left(q^0 y\right)}+\pdx{\left(q_r^0y\left(1-y\right)\right)} =\frac{\rho^0\left(\bar{y}^{0}-y\right)}{\tau^0}, \end{split}\right. \end{equation} with $q_r=q_r^0;~\tau = \tau^0;~\bar{y} = \bar{y}^0$. Setting $q_{lin}^0 = 2 q_r^0,~q_{log}^0 = q^0+q_r^0-q_{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: \begin{equation}\label{eq:4eme_equation_analytique} \left(q_{log}^0+q_{lin}^0\left(y_\infty- y\right)\right)\pder{y}{x}=\gamma^0\left(y_\infty-y\right). \end{equation} 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 $\blue{y}(\orange{x})$ reads \begin{equation}\label{eq:4eme_equation_analytique_solution} q_{lin}^0(\blue{y}-y_0)-q_{log}^0\ln{\left(\frac{y_\infty-\blue{y}}{y_\infty-y_0}\right)}=\gamma^0 \orange{x}. \end{equation} For numerical simulations, we use the parameters in Table~\ref{tab:param_cas_tests_param_constants}. \renewcommand{\arraystretch}{1.2} \begin{table}[h] \caption{Parameters of the test case} %{C{3.5cm}|C{3.5cm}|C{3cm}} \begin{tabular}{c|c|c} \hline Quantity & Value & Unit \\ \hline $y_\infty$ & $3.04\times 10^{-2}$ & $-$ \\ \hline $q_{lin}^0$ & $-6000$ & $kg.m^{-2}.s^{-1}$\\ \hline $q_{log}^0$ & $682.57$ & $kg.m^{-2}.s^{-1}$\\ \hline $\gamma^0$ & $7\times 10^4$ & $kg.m^{-3}.s^{-1}$\\ \hline \end{tabular}\label{tab:param_cas_tests_param_constants} \end{table} The stationary numerical solution $y_{num}$ and the exact solution $y_{ex}$ are computed, using a mesh with $\Dx = 10^{-3}$ m. Results are given in Figure~\ref{fig:sol2} for QRd scheme. Similar solutions are obtained using QG and QRq schemes. \begin{figure}[h] \centering \includegraphics[width=0.9\textwidth]{solution_logct.pdf} \caption{Mass fraction $y$ as a function of $x$ - QRd scheme (blue line) and exact (dotted orange) solutions.}\label{fig:sol2} \end{figure} Figure~\ref{fig:convergence} shows the convergence rate with the three schemes for this test case using different mesh sizes: $\Dx \in[10^{-5},10^{-1}]$ m. The $L_2$ norm of the error $\text{err}_2$, computed for each mesh size, is \begin{equation} \text{err}_{2} = \frac{\left\|y_{num}-y_{ex}\right\|_{2}}{\left\|y_{ex}\right\|_2} \text{ with } \|x\|_2 = \sqrt{\sum_{i=1}^{N_x} x^2},\ N_x = \left\lfloor\frac{1}{\Dx}\right\rfloor. \end{equation} 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. \begin{figure}[h] \centering \includegraphics[width=0.9\textwidth]{convergence_cas2_3s.pdf} \caption{$L_2$ Norm of the error as a function of mesh size for QG (green), QRd (red), and QRq (orange) schemes.}\label{fig:convergence} \end{figure} \section{Conclusion} Three different linear first-order schemes that approximate the nonlinear relative velocity term have been presented. QG scheme~\eqref{eq:discretisation_generale_qg} 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. Conversely, QRd and QRq schemes (\eqref{eq:discretisation_generale_qrd},\eqref{eq:discretisation_generale_qrq}) 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 that aim at approximating steady states. Eventually, QRq scheme seems to be more accurate than QRd scheme. \subsection{Boundary condition for QG scheme} For Cell $i$ from Figure~\ref{scheme}, QG scheme is written \begin{equation} \begin{split} \Min\Oi\left(\yinp-\yin\right) &-\Dt\sumj\Sij\qm~\yinp - \Dt\Sinf\qminf\yinp\\ &+\Dt\sumj\Sij\qg\left\{\sgg\yinp+\left(1-\sgg\right)\yjnp\right\}\\ & + \Dt\Sinf\qginf\Big\{\sgginf\yinp + \left(1-\sgginf\right)\yinf\Big\} \\ &= \Dt\Min\Oi\frac{\ybin-\yinp}{\tau^n_i}. \end{split} \end{equation} Using this scheme, the coefficients arising from conditions~\eqref{properties_qg} are \begingroup \allowdisplaybreaks \begin{equation} \begin{split} a_{ii} =& 1+\frac{\Dt}{\tau_i^n} +\Dt\sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}\sgg-\frac{\qm}{\rho_i^n}\right) \\ &+ \Dt \frac{\Sinf}{\Oi}\left(\frac{\qginf}{\rho_i^n}\sgginf-\frac{\qminf}{\rho_i^n}\right),\\ a_{ij} =& \Dt~\frac{\Sij}{\Oi}\frac{\qg}{\rho_i^n} \left(1-\sgg\right), \\ b_{i} =& \yin + \frac{\Dt}{\tau_i^n}\ybin - \Dt\frac{\Sinf}{\Oi}\frac{\qginf}{\rho_i^n}\left(1-\sgginf\right)\yinf, \\ \widehat{b}_i =& \left(1-\yin\right)+\frac{\Dt}{\tau_i^n}\left(1-\ybin\right) +\Dt\sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}-\frac{\qm}{\rho_i^n}\right) \\ &+\Dt\frac{\Sinf}{\Oi}\left(\frac{\qginf}{\rho_i^n}\left(\yinf+\left(1-\yinf\right)\sgginf\right)-\frac{\qminf}{\rho_i^n}\right). \\ \text{Assuming that } a_{ii}>0, \\ \Lambda_i\left(\dunderline{A}\right) & = 1 + \frac{\Dt}{\tau_i^n} + \Dt\sumj\frac{\Sij}{\Oi}\left(\frac{\qg}{\rho_i^n}-\frac{\qm}{\rho_i^n}\right)\\ &+\Dt\frac{\Sinf}{\Oi}\left(\frac{\qginf}{\rho_i^n}\sgginf-\frac{\qminf}{\rho_i^n}\right). \end{split} \end{equation} \endgroup The three conditions $a_{ii}>0$, $\Lambda_i>0$ and $\widehat{b}_i\geq 0$ are not always fulfilled. Using the formula \mbox{$\qginf (\yinf + (1-\yinf)\sgginf)\leq \sgginf\qginf$}, it can be stated that $a_{ii}\geq \Lambda_i\geq \widehat{b}_i$. The only remaining condition on the time step is again on $\widehat{b}_i$. This condition, which is slightly different from~\eqref{eq:condition_qg_dt}, should be monitored because it can become the most constraining one. \subsection{Boundary condition for QRd scheme} For Cell $i$ from Figure~\ref{scheme}, QRd scheme is written \begin{equation} \begin{split} \Min\Oi\left(\yinp-\yin\right) &+\Dt\sumj\Sij\qm\left\{\left(1-\sg\right)\left(\yjnp-\yinp\right)\right\} \\ &+ \Dt\Sinf\qminf\left(1-\sginf\right)\left(\yinf-\yinp\right)\\ &+\Dt\sumj\Sij\qr\left\{\sgt\yinp\left(1-\yjn\right)+\left(1-\sgt\right)\yjn\left(1-\yinp\right)\right\} \\ &+\Dt\Sinf\qrinf\Big\{\sgtinf\yinp\left(1-\yinf\right)+\left(1-\sgtinf\right)\yinf\left(1-\yinp\right)\Big\} \\ &= \Dt\Min\Oi\frac{\ybin-\yinp}{\tau^n_i}. \end{split} \end{equation} Using this scheme, the coefficients arising from conditions~\eqref{properties_qrd} are \begin{equation} \begin{split} (a_d)_{ii} &= 1+\frac{\Dt}{\tau_i^n} +\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n}\left(\sg-1\right) +\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-\yjn\right)\right)\\ & \quad + \Dt\frac{\Sinf}{\Oi}\frac{\qminf}{\rho_i^n}\left(\sginf - 1\right) + \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\Big(\sgtinf-\yinf\Big),\\ (a_d)_{ij} =& \Dt~\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n} \left(1-\sg\right), \\ (b_d)_{i} =& \yin + \frac{\Dt}{\tau_i^n}\ybin + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-1\right)\yjn\right)\\ &+ \Dt~\frac{\Sinf}{\Oi}\frac{\qminf}{\rho_i^n} \left(\sginf-1\right)\yinf + \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\left(\sgtinf-1\right)\yinf, \\ (\widehat{b}_d)_i =& \left(1-\yin\right)+\frac{\Dt}{\tau_i^n}\left(1-\ybin\right) +\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\sgt\left(1-\yjn\right) \right)\\ &+ \Dt\frac{\Sinf}{\Oi}\frac{\qminf}{\rho_i^n}\left(\sginf - 1\right)\left(1-\yinf\right) + \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\sgtinf\left(1-\yinf\right), \\ \Lambda_i(\dunderline{A}_d) =& 1 + \frac{\Dt}{\tau_i^n} + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-\yjn\right)\right)\\ &+ \Dt\frac{\Sinf}{\Oi}\frac{\qminf}{\rho_i^n}\left(\sginf - 1\right) + \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\Big(\sgtinf-\yinf\Big). \end{split} \end{equation} Once again, 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 $\qminf>0$). Indeed, regardless of the sign of the relative flow rate $\qrinf$, the flux requires a given value of $\yinf$. If the value of the outlet mass fraction $\yinf$ is unknown, a different consistent formula should be used to approximate $\yinf\in[0,1]$ (for instance $\yinf = \yin$). Otherwise, the flux should be modified. \subsection{Boundary condition for QRq scheme} For Cell $i$ from Figure~\ref{scheme}, QRq scheme is written \begin{multline} \Min\Oi\left(\yinp-\yin\right) +\Dt\sumj\Sij\qm\left\{\left(1-\sg\right)\Big(\yjnp-\yinp\Big)\right\}\\ \begin{aligned} &+\Dt\Sinf\qminf\left(1-\sginf\right)\Big(\yinf-\yinp\Big)\\ &+\Dt\sumj\Sij\qr{\footnotesize \left[ \begin{split} &\sg \Big\{\yin \left(1-\yinp\right)\left(1-\sgt\right) + \yinp \left(1-\yin\right)\sgt\Big\} \\ +& \left(1-\sg\right)\Big\{\yjn \left(1-\yjnp\right)\left(1-\sgt\right) + \yjnp \left(1-\yjn\right)\sgt\Big\} \end{split}\right]} \\ &+\Dt\Sinf\qrinf{\footnotesize \left[ \begin{array}{c} \sginf \Big\{\yin \left(1-\yinp\right)\left(1-\sgtinf\right) + \yinp \left(1-\yin\right)\sgtinf\Big\}\\ +\left(1-\sginf\right)\Big\{\yinf \left(1-\yinf\right)\Big\} \end{array}\right]}\\ &= \Dt\Min\Oi\frac{\ybin-\yinp}{\tau^n_i}. \end{aligned} \end{multline} Using this scheme, the coefficients arising from conditions~\eqref{properties_qrq} are \begin{equation} \begin{split} (a_q)_{ii} =& 1 +\frac{\Dt}{\tau_i^n}+\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n}\left(\sg-1\right) +\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\sg\left(\sgt-\yin\right)\right)\\ &+\Dt\frac{\Sinf}{\Oi}\frac{\qminf}{\rho_i^n}\left(\sginf-1\right) + \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\sginf\left(\sgtinf-\yin\right),\\ (a_q)_{ij} =& \Dt~\frac{\Sij}{\Oi}\frac{\qm}{\rho_i^n} \left(1-\sg\right) + \Dt\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(1-\sg\right)\left(\sgt-\yjn\right),\\ (b_q)_{i} =& \yin +\frac{\Dt}{\tau_i^n}\ybin + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\left(\sgt-1\right)\Big(\sg\yin+\left(1-\sg\right)\yjn\Big)\right) \\ &+ \Dt\frac{\Sinf}{\Oi}(\sginf-1)\left(\frac{\qminf+\left(1-\yinf\right)\qrinf}{\rho_i^n}\right)\yinf\\ &+\Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\sginf\left(\sgtinf-1\right)\yin,\\ (\widehat{b}_q)_i =& \left(1-\yin\right)+\frac{\Dt}{\tau_i^n}\left(1-\ybin\right)+\Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\sgt\Big(1-\left[\sg\yin+\left(1-\sg\right)\yjn)\right]\Big) \right)\\ +&\Dt\frac{\Sinf}{\Oi}\left(\sginf-1\right)\frac{\qminf-\yinf\qrinf}{\rho_i^n}\left(1-\yinf\right)\\ +& \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\sginf\sgtinf\left(1-\yin\right).\\ \text{Assuming that } &\forall j\in v(i), (a_q)_{ij}\leq 0, \\ \Lambda_i(\dunderline{A}_q) =& 1 + \frac{\Dt}{\tau_i^n} + \Dt\left(\sumj\frac{\Sij}{\Oi}\frac{\qr}{\rho_i^n}\Big(\sgt-\left[\sg\yin+\left(1-\sg\right)\yjn\right]\Big)\right)\\ &+\Dt\frac{\Sinf}{\Oi}\frac{\qminf}{\rho_i^n}\left(\sginf-1\right) + \Dt\frac{\Sinf}{\Oi}\frac{\qrinf}{\rho_i^n}\sginf\left(\sgtinf-\yin\right). \end{split} \end{equation} The same conditions as~\eqref{eq:condition_1} and~\eqref{eq:condition_2} on the flow rates for neighboring inner cells again arise from condition $a_{ij}^q\leq 0$. Similar conditions appear for the boundary face from condition $b_i^q\geq 0$ and $\widehat{b}_i^q\geq 0$. When $\qminf < 0$, the inlet boundary flow rates must be coherent with \begin{equation}\label{eq:condition_qrq_boundary} \begin{split} \qminf + \left(1-\yinf\right)\qrinf &\leq 0,\\ \qminf - \qrinf\yinf &\leq 0. \end{split} \end{equation} When~\eqref{eq:condition_qrq_boundary} is verified, no condition on the time step appears. Conditions~\eqref{eq:condition_qrq_boundary} are automatically verified when a co-current flow is considered, because in this case, the mixture flow rate $\qminf$, the gas flow rate $\qginf=\qminf + (1-\yinf)\qrinf$ and the liquid flow rate $(q_l)_{i\infty}=\qminf-\yinf\qrinf$ have the same sign. When~\eqref{eq:condition_qrq_boundary} is not verified, a condition on the time step appears from condition $b_i^q\geq 0$ or from condition $\widehat{b}_i^q\geq 0$. \section*{Conflicts of interest} The authors declare no competing financial interest. \section*{Dedication} This manuscript was written with the contributions of all authors. All authors have approved the final version of the manuscript. \section*{Acknowledgments} The authors would like to thank Erwan Le Coupanec for his guidance and help during the writing of this article. \appendix \def\theequation{\Alph{section}\arabic{equation}} \setcounter{equation}{0} \def\thetable{\Alph{section}\arabic{table}} %\setcounter{figure}{0} \section{Continuous maximum principle}\label{an:pmc} The study equations are \begin{equation} \begin{split} \pdt{\rho}+\DIV{\Qm} &= 0, \\ \pdt{\left(\rho y\right)} + \DIVb{y\Qm}+\DIVB{y\left(1-y\right)\Qr}&= \displaystyle\rho\frac{\bar{y}-y}{\tau}. \end{split} \end{equation} The notation $\xi = y(1-y)$ is used to study the maximum principle for y. Indeed, if $\xi\geq 0$ then \begin{equation} y(1-y)\geq 0 \Leftrightarrow y \in [0,1]. \end{equation} The governing equation for $\xi$ can be obtained from the non-conservative equations of $y$ and $1-y$ as follows: \begin{equation}\label{eq:xi_init} \begin{gathered} \left(\rho\dift{y}+\Qm\cdot\dif{y}+\difb{\xi\Qr} = \rho\frac{\bar{y}-y}{\tau} \right) \times \left(1-y\right)\\ + \left(\rho\dift{(1-y)}+\Qm\cdot\dif{\big(1-y\big)}-\difb{\xi\Qr} = \rho\frac{(1-\bar{y})-(1-y)}{\tau} \right) \times y. \end{gathered} \end{equation} Then \begin{equation} \rho\dift{\xi}+\Qm\cdot\dif{\xi}+(1-2y)\difb{\xi\Qr}=\rho\frac{S-2\xi}{\tau} \text{ with } S = y(1-\bar{y})+\bar{y}(1-y). \end{equation} Using the standard notations $\xip = \max(\xi,0)$ and $\xim = -\min(\xi,0)$, equation~\eqref{eq:xi_init} is multiplied by $(-\xim)$: \begin{equation} \begin{split} &\left(\rho\dift{\left(\ximc\right)}+\Qm\cdot\dif{\left(\ximc\right)}\right)+(1-2y)\Big(\xim^2\dif{\,\cdot\,\Qr}-\xim\Qr\cdot\dif{\xi}\Big) =\rho\frac{-S\xim-2\xim^2}{\tau}\\ &\left(\dift{\left(\rho\ximc\right)}+\left(\Qm\ximc\right)\right)+\left(\left(1-2y\right)\Qr\ximc\right)+\xim^2\left\{\left(1-2y\right)\dif{\,\cdot\,\Qr}-\frac{1}{2}\difB{\left(1-2y\right)\Qr}\right\}\\ &\mkern440mu=-\frac{\rho (S+2\xim)\xim}{\tau}. \end{split} \end{equation} Finally, the equation obtained is \begin{equation} \dift{\left(\rho\ximc\right)}+\dif{\cdot\left(\left[\Qm+\left(1-2y\right)\Qr\right]\ximc\right)}+\ximc\Big\{\left(1-2y\right)\dif{\,\cdot\,\Qr}-\Qr\cdot\dif{\left(1-2y\right)}\Big\}=-\frac{\rho (S+2\xim)}{\tau}\xim. \end{equation} We define the quantity $\displaystyle E(t)= \int_\Omega \rho \ximc d\Omega$ and split the outside surface domain according to the sign of $\qxi = \Qm+(1-2y)\Qr$: \begin{equation} \begin{split} \Gamma_w^\xi & = \left\{\xv\in\Gamma,~\bfq_\xi\cdot\bfn_\Gamma = 0\right\},\\ \Gamma_+^\xi & = \left\{\xv\in\Gamma,~\bfq_\xi\cdot\bfn_\Gamma < 0\right\}, \\ \Gamma_-^\xi & = \left\{\xv\in\Gamma,~\bfq_\xi\cdot\bfn_\Gamma > 0\right\}. \end{split} \end{equation} Then, the quantity $E(t)$ verifies \begin{equation} \begin{split} \der{E(t)}{t} = &- \int_\Gamma\left[(\Qm+\left(1-2y\right)\Qr)\ximc\right]\cdot\bfn_\Gamma d\Gamma -\int_\Omega \ximc \Big\{\left(1-2y\right)\dif{\,\cdot\,\Qr}-\Qr\cdot\dif{\left(1-2y\right)}\Big\}d\Omega\\ &- \int_\Omega \frac{\rho (S+2\xim)}{\tau}\xim. \end{split} \end{equation} %\newpage \begin{enonce}[remark]{Assumptions}\leavevmode \begin{enumerate} \item Relevant mass fraction on the boundary conditions: $\xim(\xv\in\Gamma_+^\xi,t)=0$. \item Mass fraction for the initial condition such that: $\xim(\xv\in\Omega,t=0)=0$ and thus $E(t=0)=0$. \item $\inv{\rho}[(1-2y)\dif{\,\cdot\,\Qr}-\Qr\cdot\dif{(1-2y)}]\in \clL^\infty(\Omega,[0,T])$. \item Equilibrium mass fraction such that $\bar{y}\in[0,1]$ implying, $(S+2\xim) \geq 0$ (proof of lemma below). \item Positive relaxation time scale $\tau>0$ and density $\rho>0$. \end{enumerate} \end{enonce} Thus, we have \begin{equation} \der{E(t)}{t} \leq -\int_\Omega \rho\ximc \left[\inv{\rho}\Big(\left(1-2y\right)\dif{\,\cdot\,\Qr}-\Qr\cdot\dif{\big(1-2y\big)}\Big)\right]d\Omega. \end{equation} Assuming the assumptions previously presented, \begin{equation} \der{E(t)}{t} \leq \left\|\inv{\rho}\Big[\left(1-2y\right)\dif{\,\cdot\,\Qr}-\Qr\cdot\dif{\big(1-2y\big)}\Big]\right\|_\infty E(t). \end{equation} Using the Gr\"onwall's inequality~\cite{c680dbbe-119c-31e0-a97a-d790b679674f} and $E(0)=0$, it enables to conclude that $E(t)=0$ on [0,T]. \vspace*{4pt} %\begin{center} %\rule{15cm}{1pt} %\end{center} \begin{lemm} If $\bar{y}\in[0,1]$ then $S+2\xim\geq 0$ with $S=\bar{y}(1-y)+y(1-\bar{y})$ \end{lemm} \begin{proof} Let us assume that $\bar{y}\in[0,1]$. \emph{Case $\xim=0$} \\ If $\xim = 0$, then $y\in [0,1]$, implying that $S\geq 0$. \\ Consequently, we have $S+2\xim\geq 0$ \emph{Case $\xim \neq 0$} \begin{equation} \xi= - \xim = \text{ thus } \left\{ \begin{split} y &< 0 \\ \text{or } y &> 1. \end{split}\right. \end{equation} If $y<0$, using $S=\bar{y}+y(1-2\bar{y})$, the following inequality is verified: %\begin{array}{lrcccl} %&y&\leq& S &\leq& 1-y \\ %\Rightarrow & Sy &\geq& y(1-y)=\xi \\ %\text{ and } &S(1-y)&\geq& y(1-y)=\xi. \\ %\text{ So } & S(y+1-y) = S &\geq& 2\xi = -2\xim. %\end{array} \begin{equation} \begin{split} y &\leq S \leq 1-y \\ \Rightarrow Sy &\geq y(1-y) =\xi \\ \text{ and } S(1-y) &\geq y(1-y) =\xi. \\ \text{ So } S(y+1-y) &= S \geq 2\xi = -2\xim. \end{split} \end{equation} If $y>1$, the same proof can be used in a symmetric manner. Finally, we have $S + 2\xim \geq 0$ for $y\in\mathbb{R}$. \end{proof} \section{Boundary conditions}\label{appb} As shown in Figure~\ref{scheme}, 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 $\Sinf$, a given valid mass fraction $\yinf\in[0,1]$ outside and an outward unit normal $\bfn_\infty$. The flow rates on this face are noted \begin{equation} \begin{split} \qminf &= \Qm_{i\infty}\cdot\bfn_\infty,\\ \qrinf &= (\Qr)_{i\infty}\cdot\bfn_\infty,\\ \qginf &= \qminf + (1-y_{i\infty})\qrinf, \end{split} \end{equation} where \begin{equation} y_{i\infty} =\left\{ \begin{array}{lcl} y_i^n &\text{ si }& \qminf \geq 0, \\ y_\infty &\text{ si }& \qminf<0. \end{array}\right. \end{equation} Only one inlet/outlet face is considered here, but the discussion is valid for several other faces. 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 \mbox{$\yinf\in[\yinf^n,\yinf^{n+1}]$}. Conditions~\eqref{properties_qg},\eqref{properties_qrd} and~\eqref{properties_qrq} are computed to verify the discrete maximum principle for Cell $i$. Blue terms are terms added due to the inlet/outlet boundary face. \tikzstyle{line} = [line width = 2,-] \tikzstyle{separation} = [line width = 2,densely dashed,-] \tikzstyle{arrow} = [line width = 2,<->,>=stealth] \tikzstyle{arrow_simp} = [line width = 2,->,>=stealth] \begin{figure}[h] \centering \begin{tikzpicture}[scale=0.8,font=\LARGE] \draw[arrow_simp,black] (2.0,0.0) -- (2.0,-0.7) node[below] {$\bfn_\infty$}; \draw[fill=black!15] (0.0,0.0) -- (-1.5,2.5) -- (1.0,6.0) -- (6.0,3.5) -- (4.0,0.0); \draw[line, color=blue] (0.0,0.0) -- (4.0,0.0); \draw[line, color=brown] (-1.5,2.5) -- (1.0,6.0) -- (6.0,3.5) -- (4.0,0.0); \draw[line width = 4, color=black] (0.0,0.0) -- (-1.5,2.5); \draw (2.0,3.0) node[black] {\emph{i}}; \draw (-1.0,5.0) node[black] {\emph{k}}; \draw (3.95,5.7) node[black] {\emph{l}}; \draw (5.8,1.2) node[black] {\emph{m}}; \draw (-1.5,1.0) node[black] {\emph{w}}; \draw (2.0,-1.5) node[blue] {$\mathbf{ y_\infty}$}; \draw[arrow, blue] (1.0,-0.7) -- (1.0,0.7); \draw[arrow, blue] (3.0,-0.7) -- (3.0,0.7); \draw (0.8,-0.55) node[blue, below] {$\bfq_{i\infty}$}; \draw (3.3,-0.55) node[blue, below] {$(\bfq_r)_{i\infty}$} ; \draw[arrow,black] (0.0,0.8) -- (4.0,0.8) node[midway, above] {{$S_{i\infty}$}}; \draw[arrow_simp,brown] (0.38,4.10)-- (-0.59,4.80); \draw[arrow_simp,brown] (-0.88,4.40)-- (0.09,3.70); \draw[arrow_simp,brown] (3.46,4.10)-- (3.99,5.17); \draw[arrow_simp,brown] (3.54,5.40)-- (3.01,4.33); \draw[arrow_simp,brown] (4.36,1.83)-- (5.40,1.24); \draw[arrow_simp,brown] (5.64,1.67)-- (4.60,2.26); \draw[separation, color=black] (-1.5,2.5) -- (-2.2,2.5); \draw[separation, color=black] (1.0,6.0) -- (1.0,6.7); \draw[separation, color=black] (6.0,3.5) -- (6.7,3.5); \draw[separation, color=black] (0.0,0.0) -- (-0.5,-0.5); \draw[separation, color=black] (4.0,0.0) -- (4.5,-0.5); \end{tikzpicture} \caption{Cell $i$ with $(k,l,m)$ neighboring inner cells, a wall boundary face, and an inlet/outlet boundary face.}\label{scheme} \end{figure} %\nocite{*} %\bibliographystyle{crunsrt} %\bibliography{Lazare} \printbibliography \end{document}