\documentclass[twoside, 12pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{amssymb}
% 
%\usepackage{graphicx}
\usepackage{psfrag}
%\usepackage{amsmath}
%\usepackage{amsfonts}
\usepackage{eufrak}
\usepackage{eucal}
\usepackage{color}
%\usepackage{subcaption}
\usepackage{graphicx}
\usepackage{grffile}
%\usepackage{float}
\usepackage{float}
\setcounter{MaxMatrixCols}{10}
%TCIDATA{OutputFilter=Latex.dll}
%TCIDATA{Version=5.00.0.2552}
%TCIDATA{<META NAME="SaveForMode" CONTENT="1">}
%TCIDATA{LastRevised=Monday, January 14, 2013 07:41:47}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}

\thispagestyle{empty} \pagestyle{myheadings} \markboth{Souheyla Ouffa  and Djamila Seba  }{TVD WAF scheme with PVRS Riemann Solver for the Drift-Flux Equations }
\newtheorem{theorem}{Theorem}[section]
\newtheorem{acknowledgement}[theorem]{Acknowledgement}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{axiom}[theorem]{Axiom}
\newtheorem{case}[theorem]{Case}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conclusion}[theorem]{Conclusion}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{criterion}[theorem]{Criterion}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{solution}[theorem]{Solution}
\newtheorem{summary}[theorem]{Summary}
\newenvironment{proof}[1][Proof]{\noindent\textbf{#1.} }{\ \rule{0.5em}{0.5em}}
\textheight=8in \textwidth=6.1in \oddsidemargin=0.15in
\evensidemargin=0.15in \numberwithin{equation}{section} \textwidth
16cm \textheight 22.1cm
\topmargin -1cm
\oddsidemargin 0.5cm
\baselineskip=15.5cm

\begin{document}

\begin{center}
	{\Large \textbf{TVD WAF scheme with PVRS Riemann Solver for the Drift-Flux Equations of Two-Phase Flows under Isothermal conditions}}\\[0pt]
		 	\vskip1cm \textbf{Souheyla Ouffa $^{a}$ and Djamila Seba $^{a}$  }
 \\ [0pt]$^{a}$ Dynamic of Engines and Vibroacoustic Laboratory, FSI
UMBB, Boumerdes, Algeria \\[0pt]
e-mail: s.ouffa@univ-boumerdes.dz\\[0pt]
e-mail:  seba@univ-boumerdes.dz \\[0pt]
	
\end{center}
\begin{abstract}
This paper  concentrated on the extension of the TVD version of  Weighted Average Flux (WAF) scheme for the numerical simulation of two-phase gas-liquid flow  for drift-flux model under isotermal condition. The model becomes a hyperbolic system of conservation laws with realistic closure relations where both phases are strongly coupled during their motion. Based on a PVRS  Riemann solver for the drift-flux model, the model is numerically solved by the WAF scheme. The accuracy of the proposed numerical algorithm is verified by comparing it with other numerical methods compared to the exact solutions through a set of different test cases mentioned in the literature ~\cite{dia15}. Simulation results appear  appropriate agreement of tvd WAF scheme together with the exact solutions. we note that the behaviour of WAF scheme is encouraging, more accurate and fast than other numerical methods.  
\end{abstract}
\noindent\textbf{Key words:
	WAF scheme; approximation ; drift-flux model; Riemann problem; numerical simulation\\
}
\section{Introduction}
The drift model that governs the equations presented in this paper is simples of the two-fluid model type. Moreover, the model has simpler expression form with low computational cost for predicting basic flow properties. the drift-flux model is frequently considered for the numerical approximation of practical cases related to two-phase flows. As regard as to the numerical approximations of such models, the model requires advanced and specialised methods due to their mathematical and physical. The model equations are solved exactely by [S. Kuila, T. R. Sekhar and D. Zeidan] in the literature \cite{dia16, dia15}, such the analytical solution suggested here is based on the Riemann problem and  is constructed under the assumption of isothermal conditions.  In \cite{dia16, dia15}, Godunov methods of centred-type are extended to both isothermal and isentropic gas-liquid two-phase flow employing the drift-flux model and validated with the developed exact Riemann solvers. The objective of this paper is to extend another high-order Godunov method of single-phase flow to the isotermal drift-flux model. In this context, the Weighted Average Flux (WAF) approximation proposed initially in \cite{torob}  is extended for simulating compressible gas and liquid two-phase flows. It is worth to note
that the WAF scheme is applied to approximate shock wave propagations appearing in two-phase flows governed by the two-fluid model in \cite{uk1, uk2, kr1, kr2}. This scheme is completely different from  the numerical methods which have been used before to resolve the drift flow model, and almost similar to what was used in \cite{sdj}. The WAF scheme for the isentropic drift-flux model of compressible
two-phase flows. The WAF scheme is a generalisation of the widely known Godunov first-order upwind scheme for the non-linear conservation laws.  On this basis, we propose to approximate this model with the  TVD version of WAF scheme based on the PVRS approximation Riemann solver originally used in literature \cite{torob}. It is a clear scheme that has second-order accuracy in space and time. Moreover, the feature of the scheme is that it is based on the fact that all the distinct fields are included. In the following sections, the used WAF scheme is presented in detail with version TVD and make a description of the user model. In the last section, the results obtained are presented.

\section{Numerical resolution approach}
The drift-flux model can be written in the following form without taking into account the presence source terms ~\cite{dia16, dia15}:
\begin{align}
	\frac{\partial \mathbb{U}}{\partial t} + \frac{\partial \mathbb{F(U)}}{\partial x} = 0, \quad x\in\mathbb{R} \quad {\text{and}} \quad  t\in\mathbb{R}^{+},
	\label{eq1}
\end{align}
in which the state vector of conservative variables $\mathbb{U}$ and the fluxes vector $\mathbb{F}$ are:
\begin{align}
	{\mathbb{U}} =
	\begin{pmatrix}
		\begin{array}{c}
			\alpha_{g}\rho_{g}\\
			\alpha_{l}\rho_{l} \\
			(\alpha_{g}\rho_{g} + \alpha_{l}\rho_{l}) u
		\end{array}
	\end{pmatrix}\quad {\mbox{and}} \quad
	{\mathbb{F}(\mathbb{U})} =
	\begin{pmatrix}
		\begin{array}{c}
			\alpha_{g}\rho_{g}u\\
			\alpha_{l}\rho_{l}u\\
			(\alpha_{g} \rho_{g} + \rho_{l}\alpha_{l})u + P
		\end{array}
	\end{pmatrix},
	\nonumber
\end{align}
which constitute the gas and liquid mass balance equations,  as each component is isothermal,  the mixture pressure   [Kuila et al. (2015)] given by

\begin{align}
	P = a_{g}^{2} P_{g} + a_{l}^{2} P_{l},
	\label{p}
\end{align}
The system \eqref{eq1} in a nonconservative form writen as
\begin{align}
	\frac{\partial \mathbb{V}}{\partial t} + \mathbb{A}\frac{\partial \mathbb{V}}{\partial x} = 0,
	\label{eeq1}
\end{align}
with\begin{equation}
	\mathbb{V}=\left( \begin{array}{c}
		\alpha_{g}\rho_{g}\\ \alpha_{l}\rho_{l}\\u
	\end{array}\right)\quad  and \quad \mathbb{A(V)}=\begin{pmatrix}
		u & 0&\alpha_{g}\rho_{g}\\0& u&\alpha_{l}\rho_{l} \\ \dfrac{a_{g}^{2}}{\alpha_{g}\rho_{g}+\alpha_{l}\rho_{l}} &\dfrac{a_{l}^{2}}{\alpha_{g}\rho_{g}+\alpha_{l}\rho_{l}} &u
	\end{pmatrix}
	\label{PV2}
\end{equation}
It has been shown in the literature \cite{dia15} that system must be strictly deterministic with a complete eigenstructure for an isothermal state. Where eigenvalues are given by:
\begin{align}
	\lambda_{1} = u - 
	\sqrt{\frac{a_{g}^{2} P_{g} +  a_{l}^{2} P_{l} }
		{\alpha_{g} \rho_{g} + \alpha_{l} \rho_{l}}}, 
	\quad \lambda_{2} = u \quad {\mbox{and}} \quad
	\lambda_{3} = u +
	\sqrt{\frac{a_{g}^{2} P_{g} + a_{l}^{2} P_{l} }
		{\alpha_{g} \rho_{g} + \alpha_{l} \rho_{l}}}\cdot
	\label{egn}
\end{align}
the numerical precision of system \eqref{eq1} can be approximated by the conservative and explicit finite volume , and taking into consideration the initial and boundary conditions within a control volume $[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}] \times [t^{n}, t^{n+1}]$ as :
\begin{align}
	{{\mathbb{U}}}_{i}^{n+1}  = {\mathbb{U}}_{i}^{n} +
	\frac{\Delta t}{\Delta x} \Big[{\mathbb{F}}_{{i} - \frac{1}{2}} -
	{\mathbb{F}}_{{i} + \frac{1}{2}} \Big],
	\label{fv1}
\end{align} 
in which $\mathbb{U}_{i}^{n}$  is the space average approximation at time n in the ith computational cell  and $\mathbb{F}_{i+\frac{1}{2}}^{n}$   is the numerical  inter cells flux at time n in the ith computational cell. The distance between two cross sections is $\Delta x$ and $\Delta t$ is the time step calculated using: 
\begin{align}
	\Delta t  = \mathcal{C} \ \frac{\Delta x}{{\lambda}^{^{(n)}}_{_{\text{max}}}},\quad  0\leq \mathcal{C} \leq1,
	\label{fv3}
\end{align}
where $\mathcal{C}$ is the Courant-Friedrichs-Lewy (CFL) number. ${{\lambda}^{^{(n)}}_{_{\text{max}}}}$ given by 
\begin{align}
	{\lambda}^{^{(n)}}_{_{\text{max}}} =
	\max_{i} \big\{|\lambda_{i}|\big\},
	\label{fv4}
\end{align}
is the maximum wave speed at the current time level and $\lambda_{i}$ are eigenvalues corresponding to sound waves.

For the method using , The numerical fluxes at the cell face are computed by solving the local Riemann problem for  system \eqref{eq1} with the  initial data :
\begin{align}
	\mathbb{U}\big(x, t = 0\big) = \left\{ 
	\begin{array}{lcr}
		\mathbb{U}_{i}^{n} \quad \quad & \text{if} &\quad \quad x\le x_{0},\\
		\\
		\mathbb{U}_{i+1}^{n} & \text{if} & x>x_{0},
	\end{array} 
	\right.
	\label{fv2}
\end{align}
in which $x_{0}$ is the initial data discontinuity. 
The WAF scheme~\cite{torob} is a second-order approach of accuracy in time and space without data reconstruction, which is developed in the framework of conservative systems.
In this scheme, the intercell numerical fluxes are calculated through the following integral average form of $\mathbb{F}(\mathbb{U})$ at the half-time step~\cite{torob}:
\begin{align}
	\mathbb{F}_{i + \frac{1}{2}}^{\texttt{{WAF}}} = \frac{1}{\Delta x} 
	\int^{\frac{\Delta x}{2}}_{\frac{-\Delta x}{2}} \mathbb{F}\Bigg[\mathbb{U}_{i + \frac{1}{2}}\bigg(x, \frac{\Delta t}{2}\bigg)\Bigg] dx,
	\label{waf1}
\end{align}
in which $\mathbb{U}_{i + \frac{1}{2}}$ is the solution of the Riemann problem with the initial states~\eqref{fv2}. This solution is based on the information about the complete wave structure of system~\eqref{eq1} making use of the wave jumps. The solution then is based on primitive variable PVRS approximate Riemann solvers for the drift-flux model. For solving linear hyperbolic systems with constant coefficients and and apply  the pvrs rieman solver  was originally mentioned in the literary \cite{torob}, we assume that $\mathbb{A}(\mathbb{V})$  to be constant and the initial data $\mathbb{V}_{L}$,  $\mathbb{V}_{R}$ and the solution are close to a constant state $\bar{\mathbb{V}}$, so 
\begin{align}
	\bar{\mathbb{A}}=\mathbb{A}(\bar{\mathbb{V}})
	\label{Fffv2}
\end{align}
Then  the approximation of the Riemann problem for system  treat  as the system in conservative form by 
\begin{align}
	\frac{\partial \mathbb{V}}{\partial t} + \frac{\partial \mathbb{F(V)}}{\partial x} = 0, 	\label{ffv2} \\ \nonumber \\
	\quad \mathbb{F(V)}\equiv \bar{\mathbb{A}}\mathbb{V}
\nonumber
\end{align}
and within  $\bar{\lambda}_{i}$ 
\begin{align}
	\Delta \mathbb{F}= \bar{\lambda}_{i}\Delta \mathbb{V}
	\label{fkfv2}
\end{align}
Direct application of waves 1 and 3 and after some algebraic manipulation we get the complete solution.
The WAF integral average~\eqref{waf1} can be written in the form of wave structure as:
\begin{align}
	\mathbb{F}_{i + \frac{1}{2}}^{\texttt{{WAF}}} =
	\frac{1}{2}\Big(\mathbb{F}_{i} + \mathbb{F}_{i + 1} \Big)  - \frac{1}{2} \sum_{k=1}^{N}\mbox{sign}(\mathcal{C}_{k})\, \phi_{i + \frac{1}{2}}^{(k)} \, \Delta \mathbb{F}_{i + \frac{1}{2}}^{(k)}.
	\label{wfa2}
\end{align}
where 
\begin{align}
	\phi_{i + \frac{1}{2}}^{(k)}= \phi_{i + \frac{1}{2}}(\phi(r^{k}, |\mathcal{C}|)\quad  and \quad \mathcal{C}_{k} = \lambda_{k} \frac{\Delta t}{\Delta x},
	\label{wfa3}
\end{align}
$\mathcal{C}_{k}$ is the wave CFL number for wave $k$ of speed $\lambda_{k}$, $\phi_{i + \frac{1}{2}}^{(k)}$  is the WAF limiter  employed to apply the Total Variation Diminishing (TVD) condition to the flux during the computation and to suppress any spurious oscillations near discontinuities and high gradients,   Although there are many choices for limiters, it is found that the SUPERBEE limiter of Roe \cite{roe} is the best suited to the current investigation. The WAF limiter, where the SUPERBE function is applied in this work, is given as follow:
\begin{eqnarray}
	\phi(r^{k}, |\mathcal{C}|) = 
	\left\{
	\begin{array}{c c r}
		1                            &   \mbox{if} &  r^{k}\le 0,\\
		1 - 2(1 - |\mathcal{C}_{k}|)\, r^{k} &   \mbox{if} &    0\leq r^{k}\le \frac{1}{2},\\
		|\mathcal{C}_{k}| &   \mbox{if} & \frac{1}{2}\leq r^{k}\le 1\\
		1 - (1 - |\mathcal{C}_{k}|)\, r^{k} &   \mbox{if} &    1\leq r^{k}\le 2,\\
		2|\mathcal{C}_{k}| - 1 &   \mbox{if} &    r^{k}\ge 2,
	\end{array}
	\right.
	\label{wfa5}
\end{eqnarray}
in which 
\begin{align}
	r^{(k)} = \left\lbrace \begin{array}{lll}
		\dfrac{\Delta \chi_{i- \frac{1}{2}}^{k}}{\Delta \chi_{i+ \frac{1}{2}}^{k}} & &\mathcal{C}_{k}> 0 ,\\
		\dfrac{\Delta \chi_{i+ \frac{3}{2}}^{k}}{\Delta \chi_{i+ \frac{1}{2}}^{k}}& &\mathcal{C}_{k}<0
	\end{array}
	\right.
	\label{wfa6}
\end{align}
where the choice of $\chi$ is liquid or gas density as follow:
\begin{align}
	\chi_{g} = \alpha_{g} \rho_{g} \quad {\mbox{or}} \quad \chi_{l} = \alpha_{l} \rho_{l}.
	\label{wfa9}
\end{align}
\section{Computational Results}

The objective of this paper is to extend a weighted average flow (WAF)scheme of the drift flow model of two phase  with using  PVRS approximation of  Riemann solver. To demonstrate the numerical accuracy and computational efficiency of the current WAF scheme, some  published  test cases  have been carried out are performed. \cite{dia16}, where the initial data of this cases are shown  in table 1. \begin{table}[ht!]
	\caption{Initial  data for Test 1 as in ~\cite{dia15}.}
	\begin{center}
		\begin{tabular}{cl}
			\hline
			Test &
			\begin{tabular}{cllllllllllllllll}
				& &&&&  $\rho_{1}$ && $\rho_{2}$ &&& $u$&&& $a_{g}$&& $a_{l}$
			\end{tabular} \\
			\hline
			Test 1 &
			\begin{tabular}{clllll}
				Left data & 50.0&  1000.0 & -100.00&500.0& 100.0 \\
				Right data & 50.0&  1000.0 & 100.00&500.0 &100.0
			\end{tabular}
			\\ 
			\hline 
			Test 2 &
			\begin{tabular}{clllllllll}
				Left data & 6.0& & 1.0 && 0.0&&10.0 &&7.0 \\
				Right data & 1.0&&  1.0 && 0.0&&10.0&& 7.0
			\end{tabular}\\ \hline 
			Test 3 &
			\begin{tabular}{clllllll}
				Left data & 0.1&&  0.1 && 200.0& 700.0 &300\\
				Right data & 0.1&&  0.1 && -200.0& 700.0 &300
			\end{tabular}\\ \hline
			Test 4 &
			\begin{tabular}{cllllllll}
				Left data & 500&&  500 && 0.1 &&50.0& 20.0\\
				Right data & 600&&  600 && 1.0&&50.0& 20.0
			\end{tabular} \\ \hline 
		\end{tabular}
	\end{center}
	\label{trtab1}
\end{table}

 These include: (i) Test 1: pure rarefaction  and the solution is presented at time t = 1.85 are  shown in Figure \eqref{fig41}; \begin{figure}[ht!]
 	\begin{center}
 		\begin{tabular}{cc}
 			\includegraphics[width=2.5in, angle=360]{gt1}&
 			\includegraphics[width=2.5in, angle=360]{lt1}\\
 			\vspace{0.1mm}\\
 			\includegraphics[width=2.5in, angle=360]{vt1}&
 			\includegraphics[width=2.5in, angle=360]{pt1}
 		\end{tabular}
 	\end{center}
 	\caption{Test 1.  Comparison of TVD WAF scheme with the Lax-Friedrichs methods together with the exact Riemann solver at time t= 1.85.\\}
 	\label{fig41}
 \end{figure} 
\begin{figure}[ht!]
\begin{center}
	\begin{tabular}{cc}
		\includegraphics[width=2.5in, angle=360]{gt3}&
		\includegraphics[width=2.5in, angle=360]{lt3}\\
		\vspace{0.1mm}\\
		\includegraphics[width=2.5in, angle=360]{vt3}&
		\includegraphics[width=2.5in, angle=360]{pt3}
	\end{tabular}
\end{center}
\caption{Test 3. Two different numerical methods, namely, TVD WAF and Lax-Friedrichs methods are compared with the exact Riemann solver at time t= 0.015 for the velocity and mixture pressure. Where a CFL = 0.9 and 100 cells are employed for the numerical resolutions\\}
\label{fig413}
\end{figure}
\begin{figure}
\begin{center}
	\begin{tabular}{cc}
		\includegraphics[width=2.5in, angle=360]{gt4}&
		\includegraphics[width=2.5in, angle=360]{lt4}\\
		\vspace{0.1mm}\\
		\includegraphics[width=2.5in, angle=360]{vt4}&
		\includegraphics[width=2.5in, angle=360]{pt4}
	\end{tabular}
\end{center}
\caption{Test 4.Comparison of numerical methods, namely, TVD WAF and Lax–Friedrichs methods with the exact Riemann solver at time t = 0.17, CFL =0.9 and with 100 cells. \\}
\label{fig44}
\end{figure} (ii) Test 2: shock waves: The solution is composed of  a left
shock wave, a contact discontinuity and right shock wave;  where  it's displayed at  time t= 0.015 are illustrated in Figure \eqref{fig413}; and (iii) Test 3:  left shock wave, a contact discontinuity and right rarefaction wave where the solution it's displayed at  time t= 0.17 are illustrated in Figure \eqref{fig44}.  In all test cases, we compare results provided by numerical resolution with previously pvrs approximation Riemann solver for drift flow model validation.\\
In following test cases, the drift-flux model is solved in a computational domain of $[-10, 10]$ with take the initial data discontinuity  $x_{0} = 0$ and  CFL parameter  0.9 is always fixed together with transmissive boundary conditions along with the SUPERBEE limiter are employed through out the simulations.  The different symbols in all test cases represent three different methods, namely, the WAF, Lax-Friedrichs  and TVD slic  methods where the solid line refers to the exact Riemann solver of \cite{dia16}.

%
\section{Conclusion}
The WAF scheme solutions are success-fully implemented and validated with different shock-tube problems.  Comparison between WAF resolutions and other numerical solutions showed acceptable accuracy and stable enough for very large time  steps. Simulations also are compared with those obtained by the exact  Riemann solver and found an excellent agreement is achieved between  numerical and exact results

 
 
 
 
 
 
 
 
 
 
 
\begin{thebibliography}{99}
	
\bibitem{torob} 
E.F. Toro, 
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, 2009.
%
%
%


%
%
%
%  

%
%

%
%

%
%
\bibitem{dia16}
S. Kuila, T. Raja Sekhar, D. Zeidan,
On the Riemann Problem Simulation for the Drift-Flux Equations of Two-Phase Flows
International Journal of Computational Methods 13 (2016) 1650009. 
%
% 
%
\bibitem{dia15}
S. Kuila,T. Raja Sekhar, D. Zeidan ,
A Robust and accurate Riemann solver for a compressible two-phase flow model,
Applied Mathematics and Computation 265 (2015) 681-695.
%																																	
%														

%
%
%
%
%
\bibitem{uk1}
R. Torrens, L.C. Wrobel,
Weighted average flux method and flux limiters for the numerical simulation of shock waves in rigid porous media,
Int. J. Numer. Meth. Fluids 40 (2002) 1187-1207.
%
%
%
\bibitem{uk2}
R. Torrens, L.C. Wrobel, 
On the propagation of a normal shock wave through a layer of incompressible porous material,
International Journal of Numerical Methods for Heat \& Fluid Flow 13 (2002) 178-198.
%
%
%
\bibitem{kr1}
G.-S.Yeom,
A Phasic WAF-HLLC Scheme for the Compressible Two-Fluid Model on a Quadrilateral Unstructured Grid,
Numerical Heat Transfer, Part B 63 (2013) 395-424.
%
%
%
\bibitem{roe}R. Saurel, F. Petitpas, R.A. Berry, Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and
shocks in multiphase mixtures, Journal of Computational Physics 228
(2009) 1678-1712.
\bibitem{kr2}
G.-S.Yeom, K.-S. Chang, S.W. Baek, 
Robust Waf-Hll Scheme for Compressible Two-Pressure Two-Velocity Multiphase Flow Model, 
Engineering Applications of Computational Fluid Mechanics 6 (2012) 144-162.
%

%
%
\bibitem{sdj}
The WAF scheme for the isentropic drift-flux model of
compressible two-phase flows,
Computers and Fluids  (2021).






\end{thebibliography}	
\end{document}