The Cauchy problem for the standard one pressure system of two fluid flows with energy equations

TL;DR

This paper rigorously constructs approximate solutions to a two-fluid flow system with energy equations, proving convergence to Radon measures and providing a numerical method that aligns with previous solutions.

Contribution

It introduces a mathematically rigorous approach to approximate and analyze solutions of the two-fluid flow system, linking PDEs to ODEs for numerical computation.

Findings

Radon measures approximate solutions match previous numerical results

The method provides a convergent numerical scheme for the system

Numerical experiments confirm the theoretical convergence and accuracy

Abstract

We construct with full rigorous mathematical proof a family of approximate solutions to the Cauchy problem for the standard system of two fluid flows with energy equations and we pass to the limit by weak compactness to obtain Radon measures that satisfy the equations in a natural weak sense. Our method provides a convergent numerical method for the numerical calculation of these Radon measures by reducing the system of partial differential equations in the case of these approximate solutions to a system of ordinary differential equations. We observe numerically on the standard Toumi shock tube problem that the Radon measures from our method agree with the numerical solutions previously obtained by other authors with various different numerical methods. In a subsequent numerical paper, using a standard confident scheme with splittings and vanishing viscosity (independent on the above construction), we observe exactly the numerical solution given by our mathematical proof.

Figures

Equations

$$\frac{\partial }{\partial t}({\rho }_i{\alpha }_i)+\frac{\partial }{\partial x}({\rho }_i{\alpha }_i{v}_i)=0,$$

$$\frac{\partial }{\partial t}({\rho }_i{\alpha }_i{v}_i)+\frac{\partial }{\partial x}({\rho }_i{\alpha }_i(v_i{)}^2)+{\alpha }_i\frac{\partial p}{\partial x}=0,$$

$$\frac{\partial }{\partial t}({\rho }_i{\alpha }_i{e}_i)+p\frac{\partial }{\partial t}{\alpha }_i+\frac{\partial }{\partial x}({\rho }_i{\alpha }_i{e}_i{v}_i)+\frac{\partial }{\partial x}(p{\alpha }_i{v}_i)=0,$$

$$p=(K_i-1)({\rho }_i{e}_i-\frac{{\rho }_i(v_i{)}^2}{2})-K_i{p}_i^{\infty },i=1,2,$$

$$V(x+ϵ,t,ϵ)v^{-}(x+ϵ,t,ϵ)\},$$

$$\frac{d}{dt}{r}_i+[r_i,v_i{]}_x={\kappa }_1(ϵ),$$

$$\frac{d}{dt}(r_i{v}_i)+[r_i{v}_i,v_i{]}_x+{\alpha }_i{\partial }_x\overline{p}=0,$$

$$\frac{d}{dt}(r_i{e}_i)+[r_i{e}_i,v_i{]}_x+p\frac{d}{dt}{\alpha }_i+({\partial }_x\overline{p}){\alpha }_i{v}_i+p{\partial }_x(\overline{{\alpha }_i{v}_i})={\kappa }_2(ϵ),$$

$$p=(K_i-1)({\rho }_i{e}_i-\frac{{\rho }_i(v_i{)}^2}{2})-K_i{p}_i^{\infty },$$

$$v_i=\frac{(r_i{v}_i)}{r_i},$$

$$\overline{p}=p\ast (\varphi {)}_{{ϵ}^{\lambda }},\overline{{\alpha }_i{v}_i}=({\alpha }_i{v}_i)\ast (\varphi {)}_{{ϵ}^{\lambda }},i=1,2,$$

$${\rho }_i^0\in L^1(\mathbb{T}),v_i^0\in L^{\infty }(\mathbb{T}),e_i^0\in L^{\infty }(\mathbb{T}).$$

$$\Vert {\rho }_i^0(.,ϵ){\Vert }_{L^1}\le const,\Vert v_i^0(.,ϵ){\Vert }_{L^{\infty }}\le const,\Vert e_i^0(.,ϵ){\Vert }_{L^{\infty }}\le const.$$

$$K_2>K_1>1$$

$$K_2{p}_2^{\infty }>K_1{p}_1^{\infty }=0.$$

$$\exists m>0\; such\; that{r}_i(x,t,ϵ)={\alpha }_i(x,t,ϵ){\rho }_i(x,t,ϵ)\ge m>0\forall x\in \mathbb{T},\forall t\in [0,T],$$

$${\int }_{\mathbb{T}}\{\frac{\partial }{\partial t}({\rho }_i{\alpha }_i)(x,t,ϵ)\psi (x)-({\rho }_i{\alpha }_i{v}_i)(x,t,ϵ)\}{\partial }_x\psi (x)\}dx\to 0,$$

$${\int }_{\mathbb{T}}\{\frac{\partial }{\partial t}({\rho }_i{\alpha }_i{v}_i)(x,t,ϵ)\psi (x)-({\rho }_i{\alpha }_i{v}_i^2)(x,t,ϵ){\partial }_x\psi (x)+[{\alpha }_i(x,t,ϵ)\frac{\partial }{\partial x}\overline{p}(x,t,ϵ)]\psi (x)\}dx\to 0,$$

$${\int }_{\mathbb{T}}\{\frac{\partial }{\partial t}({\rho }_i{\alpha }_i{e}_i)(x,t,ϵ)\psi (x)-({\rho }_i{\alpha }_i{e}_i{v}_i)(x,t,ϵ){\partial }_x\psi (x)+p(x,t,ϵ)\frac{\partial }{\partial t}({\alpha }_i)(x,t,ϵ)\psi (x)+$$

$$[\frac{\partial }{\partial x}\overline{p}(x,t,ϵ){\alpha }_i{v}_i(x,t,ϵ)+p(x,t,ϵ)\frac{\partial }{\partial x}(\overline{{\alpha }_i{v}_i})(x,t,ϵ)]\psi (x)\}dx\to 0.$$

$${\int }_{\mathbb{T}}(\overline{p}-p)(x,t,ϵ)\psi (x)dx\to 0$$

$${\int }_{\mathbb{T}}(\overline{{\alpha }_i{v}_i}-{\alpha }_i{v}_i)(x,t,ϵ)\psi (x)dx\to 0,$$

$$(K_1-1)({\rho }_1{e}_1-\frac{{\rho }_1{v}_1^2}{2})-K_1{p}_1^{\infty }=(K_2-1)({\rho }_2{e}_2-\frac{{\rho }_2{v}_2^2}{2})-K_2{p}_2^{\infty }.$$

$${\alpha }^2\{K_1{p}_1^{\infty }-K_2{p}_2^{\infty }\}+$$

$$\alpha \{-(K_1-1)(r_1{e}_1-\frac{r_1(v_1{)}^2}{2})-K_1{p}_1^{\infty }-(K_2-1)(r_2{e}_2-\frac{r_2(v_2{)}^2}{2})+K_2{p}_2^{\infty }\}+$$

$$\{(K_1-1)(r_1{e}_1-\frac{r_1(v_1{)}^2}{2})\}=0$$

$$p\frac{d{\alpha }_j}{dt}=p\sum _{i=1}^2(\frac{d{r}_i}{dt}\frac{\partial {\alpha }_j}{\partial r_i}+\frac{d(r_i{v}_i)}{dt}\frac{\partial {\alpha }_j}{\partial (r_i{v}_i)}+\frac{d(r_i{e}_i)}{dt}\frac{\partial {\alpha }_j}{\partial (r_i{e}_i)}).$$

$$(1+p\frac{\partial {\alpha }_1}{\partial (r_1{e}_1)})\frac{d}{dt}(r_1{e}_1)+p\frac{\partial {\alpha }_1}{\partial (r_2{e}_2)}\frac{d}{dt}(r_2{e}_2)=f_1,$$

$$p\frac{\partial {\alpha }_2}{\partial (r_1{e}_1)}\frac{d}{dt}(r_1{e}_1)+(1+p\frac{\partial {\alpha }_2}{\partial (r_2{e}_2)})\frac{d}{dt}(r_2{e}_2)=f_2,$$

$$\frac{\partial \alpha }{\partial (r_1{e}_1)}[2\alpha (K_1{p}_1^{\infty }-K_2{p}_2^{\infty })+\{\{\dots \}\}]=\alpha (K_1-1)-(K_1-1)$$

Full text

AMS classification: 35D30, 35F25.

Keywords: partial differential equations, approximate solutions, fluid dynamics .

e-mail: m.colombeau@orange.fr

**this research has been partly done thanks to financial support of FAPESP, processo 2012/15780-9.

1. Introduction.

We consider the standard system used to model a mixture of two immiscible fluids from the conservation laws of mass, momentum and energy with the natural assumption that the pressures are equal inside the two fluids, as derived in [43] p. 373 and solved numerically in [39, 37, 38, 15]. We recall the system as stated in [43] p. 373

[TABLE]

[TABLE]

[TABLE]

with the state laws [39, 37, 38]

[TABLE]

where $\rho_{i},v_{i},e_{i},\alpha_{i}$ are respectively the density, the velocity, the specific total energy and the volumic proportion of fluid $i$; $p$ is the unique pressure of the mixture, so that (4) provides a compatibility relation. The constants $K_{i}$ and $p^{\infty}_{i}$ are obtained from experimental measurements on the two fluids.

This system is used in industry to model the gas kick phenomenon in offshore oil exploitation [2, 4, 33] and to model the cooling in nuclear power stations [5, 42]. This system has various peculiarities which render its numerical study far more difficult than the classical case of one single fluid. Since it has been observed [39] p. 2620 that close numerical methods can give significantly different numerical solutions, discover which is the correct one is primarily important.

The purpose of this paper is to provide with full rigorous mathematical proofs that numerical solutions obtained in [39, 37, 38] are approximate solutions of the equations. We prove existence of a solution in a weak sense. Our method provides a convergent numerical method for the numerical calculation of these approximate solutions by reducing the system of partial differential equations in the case of these approximate solutions to a system of ordinary differential equations. To this end we first construct a family of smooth approximate solutions of the initial value problem and we prove that they tend to satisfy the equations in a weak sense, i.e. when plugged into the equations one proves that the result tends to 0. Passage to the limit of approximate solutions by weak compactness gives radon measures in density, momentum and total energy that satisfy the equations in a natural weak sense. Other authors [16]-[21],[32] also introduced sequences of approximate solutions to treat nonlinear PDEs. Of course, we do not study the very difficult problem of uniqueness, which seems out of reach by the method in use here.

Then, since our approximate solutions are obtained as solutions of a system of six scalar ODEs in a Banach space, we calculate numerically these approximate solutions by means of faithful classical convergent numerical methods for ODEs. We observe numerically that one obtains same results as these authors concerning step values and locations of the discontinuities but we observe very neat peaks in liquid velocity and in gas velocity, with well defined top values and widths, that appear in our solution. Of course one could wonder the relevance of these peaks since they could be interpreted as artefacts. To clarify this point we observe numerically the stability and evolution with time of these peaks. Further, in a purely numerical paper (without proofs) [14], we obtain exactly the same result (i.e. the peaks) from a standard confident numerical scheme with splittings and vanishing viscosity. This last scheme can be adapted to both system (1-4) and to the system obtained in [39] p.2595 by adding a supplementary term to ensure hyperbolicity. One observes that this additional term ensures the disappearance of these peaks. The theoretical and numerical results in this paper and the numerical results in [14] show that these peaks are part of the solution of system (1-4) as it is, i.e. without additional terms. Further these peaks are strikingly similar to the experimental observations on the liquid flow rate $(m^{3}/s)$ when the gas kick occurs [2] figure 4 there, the survey papers [3, 6, 8, 36] and references there. It had already been noticed by B. Keyfitz et al [22]-[31] that some fully nonlinear systems which are non hyperbolic could nevertheless produce realistic physical solutions.

2. Statement of the main result. The main result is the construction of approximate solutions. Then one passes easily to the limit by weak compactness and the convergence of the numerical method follows at once from the construction. We obtain the approximate solutions as solutions of a system of 6 scalar ODEs in the Banach space $\mathcal{C}(\mathbb{T})$ of all continuous functions on the torus $\mathbb{T}=\mathbb{R}/\mathbb{Z}$. First, a simplifying notation: if $V,v$ are real valued functions of $(x,t,\epsilon)\in\mathbb{T}\times[0,T[\times]0,\eta[,\ T>0,\eta>0$, we set

$[V,v]_{x}=-\frac{1}{\epsilon}\{V(x-\epsilon,t,\epsilon)v^{+}(x-\epsilon,t,\epsilon)-V(x,t,\epsilon)|v|(x,t,\epsilon)+$

[TABLE]

where $v^{+}=max(v,0)$ and $v^{-}=max(-v,0)$ so that $v=v^{+}-v^{-}$ and $|v|=v^{+}+v^{-}$. The index $x$ in $[\dots]_{x}$ is used to avoid confusion with a simple parenthesis and because in 2-D we use $[\dots]_{y}$. We introduce two strictly positive functions $\epsilon\longmapsto\kappa_{1}(\epsilon)$ and $\epsilon\longmapsto\kappa_{2}(\epsilon)$ that tend to 0 as fast as needed; they are used in the theoretical proof. We also introduce a real number $\lambda\in]0,1[$ which is used in a convolution.

Physicists use small cubes of side $dx>0$ during times $t$ and $t+dt$ to state the laws of physics, then pass formally to the limit $dx=0$ and $dt=0$ to obtain the usual PDEs. Our method consists in letting $dt\rightarrow 0$ for fixed $dx=\epsilon>0$, then solve (theoretically from a mathematical proof and numerically from convergent schemes for ODEs) for each value of $\epsilon>0$ the system of ODEs so obtained. We pass to the limit $\epsilon=0$ on the solutions of the ODEs. For system (1-3) the method gives existence of approximate solutions, then Radon measures at their limit by weak compactness and the method gives also their numerical approximation. A physical justification of this method is given in section 5 below.

Now setting $r_{i}=\rho_{i}\alpha_{i},\ i=1,2,$ we state the ODEs we will use in the form

[TABLE]

[TABLE]

[TABLE]

together with the algebraic equations

[TABLE]

[TABLE]

where $\overline{p}$ and $\overline{\alpha_{i}v_{i}}$ are respective mollifications of $p$ and $\alpha_{i}v_{i}$:

[TABLE]

where $\phi\in\mathcal{C}_{c}^{\infty}(\mathbb{R}),\phi$ even, $\phi$ positive, $\int\phi(x)dx=1$ and $(\phi)_{\mu}(x)=\frac{1}{\mu}(\phi)(\frac{x}{\mu}).$

Remark. Via the terms $[\dots]_{x}$, which provide some kind of semi-discretization, the ODEs (6-8) model a difference between space and time: to model the space derivatives one cannot consider ”‘physical points”’ having a length smaller than the edge of a cube that contains several thousand molecules, while time can be considered as made of mathematical points. The assumption that the function $\phi$ is even models the isotropy of space away from boundaries. In absence of precisely defined solutions for nonconservative systems we believe important to stick to physical intuition as much as possible (see also section 6 concerning molecular agitation) at the same time as to provide rigorous mathematical proofs.

We assume that the initial conditions satisfy the requirements imposed by physics. More precisely $\forall x\in\mathbb{T}\ \ \rho_{i}^{0}(x)>0$ and $0<\alpha_{i}^{0}(x)<1$ (to have a real mixture) with the compatibility condition $\alpha_{1}^{0}(x)+\alpha_{2}^{0}(x)=1,e_{i}^{0}(x)>0,p^{0}(x)>0,$ with the compatibility condition imposed by the state laws (9) and

[TABLE]

When $\epsilon\rightarrow 0$ we approximate these initial conditions, for instance by convolution as in (11), by smooth functions in $x$-variable, $\rho_{i}^{0}(x,\epsilon)>0$ and $0<\alpha_{i}^{0}(x,\epsilon)<1$, with the compatibility condition $\alpha_{1}^{0}(x,\epsilon)+\alpha_{2}^{0}(x,\epsilon)=1$, $e_{i}^{0}(x,\epsilon)>0,p^{0}(x,\epsilon)>0,$ with the compatibility condition (9) for all $\epsilon$ and such that there exists $const>0$ independent on $\epsilon$ such that

[TABLE]

The strict inequalities $\rho_{i}^{0}(x,\epsilon)>0,\ 0<\alpha_{i}^{0}(x,\epsilon)<1,\ e_{i}^{0}(x,\epsilon)>0$ and $p_{i}^{0}(x,\epsilon)>0$ will be needed in the proof.

We assume that the fluids are such that

[TABLE]

and

[TABLE]

Indeed in the application to the Toumi shock tube problem in [37, 38, 39] one has $K_{1}=1.4,K_{2}=2.8,p_{1}^{\infty}=0,p_{2}^{\infty}=8.5.10^{8}.$ The fact that $K_{1}p_{1}^{\infty}=0$ plays a role in the proof of the theorem (to obtain formulas (25-27)).

We have to assume that the system of ODEs (6-11) does not lead to presence of a void region in any fluid. More precisely we have to assume that we consider solutions of (6-11) on $[0,T]$ which have the property that for $i=1,2$

[TABLE]

$\forall\epsilon>0$ small enough, which means absence of void region in any fluid.

In the sequel we assume $T$ is arbitrarily large to simplify the statement of the result. If not the result below holds on $[0,T[$ only. This limitation of absence of a void region in any fluid is needed in the proof and we have been unfortunately unable to prove the theorem from the equations and the initial data without assuming (15), which was possible in [12, 13] in the particular case of state laws $p=f(\rho)$. One has checked numerically that (15) is satisfied in the Toumi shock tube problem considered in this paper.

Theorem . Under the above assumption of absence of a void region in any fluid and if $2\lambda<1$, the system of ODEs (6-11) has a global unique solution in positive time. This solution approximates system (1-4) in the following weak sense: for all test function $\psi\in\mathcal{C}^{\infty}(\mathbb{T})$, the following limits hold $\forall t$ and for $i=1,2$ when $\epsilon\rightarrow 0$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The two state laws (4 or 9) are verified algebraically. Further $\overline{p}$ and $\overline{\alpha_{i}v_{i}}$ are mollifications of $p$ and $\alpha_{i}v_{i}$ respectively and when $\epsilon\rightarrow 0$

[TABLE]

and

[TABLE]

for all test function $\psi\in\mathcal{C}^{\infty}(\mathbb{T}).$

The terms containing these mollifications originate from the state law (9): the state laws stem from experiments done at a macroscopic order of smallness which makes a great difference with the conservation laws. The need of mollification of physical variables involved in state laws appeared in [11, 12, 13].

For fixed $\epsilon>0$ the physical variables $\overline{p}$ and $\overline{\alpha_{i}v_{i}}$ are of class $\mathcal{C}^{\infty}$ in $x,t$ variable, therefore all nonconservative products in the equations make sense for each fixed $\epsilon$. Our construction extends easily to more than two fluids (immediate) and in multidimension as in [12].

Now, before the proof, we need some preparation.

**3. Preliminary calculations.

**

$\bullet$Calculation of the volumic fractions for fixed $\epsilon$. The two state laws (9) valid for the same values of pressure imply

[TABLE]

Replacing $\rho_{i}$ by $\frac{r_{i}}{\alpha_{i}}$, multiplying by $\alpha_{1}\alpha_{2}$ and setting $\alpha_{1}=\alpha,\alpha_{2}=1-\alpha$ one obtains the following second order equation in unknown $\alpha$ in function of the 6 independent variables $r_{i},r_{i}v_{i}$ and $r_{i}e_{i}$:

[TABLE]

[TABLE]

[TABLE]

which, for each value $x,t,\epsilon$ has exactly one root $\alpha$ in the interval $]0,1[$ since we will prove that one has always $p(x,t,\epsilon)>0$: if $f(\alpha)$ denotes the left hand-side of (21) one has $f(0)=\alpha_{1}(p+K_{1}p_{1}^{\infty})>0$ and $f(1)=-\alpha_{2}(p+K_{2}p_{2}^{\infty})<0$. This value $\alpha$ from (21) is a smooth function of the independent variables $r_{i},r_{i}v_{i}$ and $r_{i}e_{i}$.

$\bullet$Resolution of the energy equations in time derivative. We develop $p\frac{d\alpha_{j}}{dt}$ in the form

[TABLE]

Using (22) the system of ODEs (6-8) can be written in the form $M\frac{dX}{dt}=N$ where $M$ is a $6\times 6$ matrix and $X=(r_{1},r_{2},r_{1}v_{1},r_{2}v_{2},r_{1}e_{1},r_{2}e_{2})^{t}$. The two energy equations (8) appear respectively in the form

[TABLE]

[TABLE]

where $f_{1}$ and $f_{2}$ are functions involving terms without time derivatives once one uses (6,7) to replace $\frac{dr_{i}}{dt}$ and $\frac{d(r_{i}v_{i})}{dt}$ from (22) by terms without time derivatives. Therefore the determinant of $M$ is $1-p(\frac{\partial\alpha}{\partial(r_{2}e_{2})}-\frac{\partial\alpha}{\partial(r_{1}e_{1})})$. The system of ODEs can be resolved in time when this determinant is nonzero. So we compute this determinant. Differentiation in the variable $(r_{1}e_{1})$ of (21) with $\alpha=\alpha(r_{1},r_{2},r_{1}v_{1},r_{2}v_{2},r_{1}e_{1},r_{2}e_{2})$ gives

[TABLE]

where we denote by $\{\{\dots\}\}$ the coefficient of $\alpha$ in (21). Similarly

[TABLE]

One obtains that

[TABLE]

We have to check that this value is always different from 1 so that the system of the two energy ODEs (8) could be solvable in time derivative. Since $K_{2}-K_{1}>0,K_{1}-1>0$ (13), $0<\alpha<1$ and since we will prove that $p(x,t,\epsilon)>0\ \forall x,t,\epsilon$ the second member of (23) is different from 1. Therefore the determinant of $M$ is always nonzero .

Therefore, after multiplication by $M^{-1}$, for each fixed $\epsilon$ one has a system of 6 ODEs in the Banach space $\mathcal{C}(\mathbb{T})$ of all continuous functions on $\mathbb{T}$ of the form

[TABLE]

with $F:\Omega\subset(\mathcal{C}(\mathbb{T}))^{6}\longmapsto(\mathcal{C}(\mathbb{T}))^{6}$ where $\Omega=\{(r_{1},r_{2})\ /\ r_{i}>0,i=1,2\}\times(\mathcal{C}(\mathbb{T}))^{4}$, because of division by $r_{i}$ in (10). For fixed $\epsilon>0$ $F$ has Lipschitz coefficients uniform in $\Omega_{M}:=\Omega\cap\{inf\{r_{i}\}\geq\frac{1}{M},\|r_{i}\|_{\infty}\leq M,\|r_{i}v_{i}\|_{\infty}\leq M,\|r_{i}e_{i}\|_{\infty}\leq M,\ i=1,2,\}$ for any $M>0$ ( the $\partial_{x}$ derivatives in (7,8) are done after convolution (11)). This remark will play a basic role in the sequel under the form that as long as a solution $(r_{i},r_{i}v_{i},r_{i}e_{i}),\ i=1,2,$ defined on some interval $[0,\delta(\epsilon)[$ takes its values in some set $\Omega_{M}$ then this solution can be extended to a larger interval $[0,\delta(\epsilon)+\eta(\epsilon)[,\ \eta(\epsilon)>0$. This will permit to prove from suitable a priori estimates existence of a global solution of (6-11) in positive time for each fixed $\epsilon$.

$\bullet$An ODE satisfied by the pressure. From the state laws (9) for the fluid 1, in which one uses that $K_{1}p_{1}^{\infty}=0$, one has

[TABLE]

and

[TABLE]

Plugging these two formulas into the energy equation (8) with $i=1$ and recalling that $r_{1}=\rho_{1}\alpha_{1}$ one obtains

[TABLE]

i.e.

$\partial_{t}(p\alpha_{1})+[p\alpha_{1},v_{1}]_{x}+\frac{K_{1}-1}{2}\{\partial_{t}(r_{1}(v_{1})^{2})+[r_{1}(v_{1})^{2},v_{1}]_{x}+2(\partial_{x}\overline{p})\alpha_{1}v_{1}+2p\partial_{x}(\overline{\alpha_{1}v_{1}})\}=$

[TABLE]

Auxiliary calculation: to the Euler equation (7) for $i=1$ we subtract the continuity equation (6) for $i=1$ multiplied by $v_{1}$. This calculation is justified because the solutions for fixed $\epsilon>0$ are smooth (one can regularize $v_{i}^{\pm}$ by setting $v_{i}^{+}-v_{i}^{-}=v_{i},v_{i}^{+}+v_{i}^{-}=\sqrt{v_{i}^{2}+a^{2}},a\in\mathbb{R},a>0$; nothing is changed in the results and proofs see [12] p. 2586-2587). One obtains

[TABLE]

Using $\partial_{t}(r_{1}(v_{1})^{2})=r_{1}v_{1}\partial_{t}v_{1}+v_{1}\partial_{t}(r_{1}v_{1})$, the above formula and (7), the term $\frac{K_{1}-1}{2}\{\dots\}$ in (25) becomes

$\frac{K_{1}-1}{2}\{[r_{1},v_{1}]_{x}(v_{1})^{2}-(v_{1})^{2}\kappa_{1}(\epsilon)-[r_{1}v_{1},v_{1}]_{x}v_{1}-\alpha_{1}(\partial_{x}\overline{p})v_{1}-$

[TABLE]

Developping from (5) the brackets $[\dots]_{x}$ that are in (26) one obtains

[TABLE]

[TABLE]

Plugging this last result into (26), then (25) one finally obtains by developing $-[p\alpha_{1},v_{1}]_{x}$ according to (5)

[TABLE]

[TABLE]

$t,\epsilon)-v_{1}(x,t,\epsilon))^{2}\}-(K_{1}-1)p\partial_{x}(\overline{\alpha_{1}v_{1}})(x,t,\epsilon)+\frac{K_{1}-1}{2}(v_{1})^{2}(x,t,\epsilon)\kappa_{1}(\epsilon)-$

[TABLE]

This formula will be used in the proof by noting that the terms in second member are positive except possibly the 3 terms involving the factor $p(x,t,\epsilon)$ since $r_{1}\geq 0,\alpha_{1}\geq 0,p\geq 0,K_{1}>1,v_{1}^{\pm}\geq 0$. In the sequel of the proof we will consider values of $(x,t)$ such that $p(x,t,\epsilon)=0$, therefore, for these values, one will have $\partial_{t}(\alpha_{1}p)(x,t,\epsilon)>0$ since $\kappa_{2}(\epsilon)>0$ and all other terms are positive or null.

$\bullet$ A simplification in an integral.

[TABLE]

proof. Dropping the notations $t$ and $\epsilon$, from (11) this integral is equal to

$\int\alpha_{i}v_{i}(x)p(x-\epsilon^{\lambda}\mu)\frac{1}{\epsilon^{\lambda}}\phi^{\prime}(\mu)d\mu dx+\int\alpha_{i}v_{i}(x-\epsilon^{\lambda}\mu)p(x)\frac{1}{\epsilon^{\lambda}}(\phi)^{\prime}(\mu)d\mu dx=$

[TABLE]

Then use the assumption that $\phi$ is even.$\Box$

Now one can start the proof of the theorem by a priori estimates.

4. Proof of the theorem. From the local existence-uniqueness theorem for ODEs in the Lipschitz case, for all $\epsilon>0$ there exists $\delta(\epsilon)>0$, depending on $\epsilon$, and a unique solution of system (6-11) on $[0,\delta(\epsilon)[$. From the choice of approximations of the initial conditions (strict positiveness) we choose $\delta(\epsilon)$ small enough so that $\forall t\in[0,\delta(\epsilon)[,\ \forall x\in\mathbb{T},\forall i=1,2$

[TABLE]

Now one assumes that the solution exists on some non necessarily small interval $[0,\delta(\epsilon)[$, for a given finite value $\delta(\epsilon)$, and that on $[0,\delta(\epsilon)[$ this solution satisfies (29). Our aim is to obtain, under assumption (15) of absence of void region in one fluid, $L^{1}$ and $L^{\infty}$ a priori estimates on the solution on $[0,\delta(\epsilon)[$, some of them uniform in $\epsilon$, so as to extend this solution to the right of $\delta(\epsilon)$ and also to use them in the proof of the approximations (16-20). The notation $const$ will denote values independent on $\epsilon$ and on $t\in[0,\delta(\epsilon)[$.

$\bullet$ First step: $L^{1}$ bounds uniform in $\epsilon$. From (6,5) and replacement of $|v|$ by $v^{+}+v^{-}$ we obtain

[TABLE]

The two integrals in $v_{i}^{+}$ simplify by translation when one integrates in $x$ on $\mathbb{T}$, as well as the two integrals in $v_{i}^{-}$. One obtains

[TABLE]

i.e. since $r_{i}=\rho_{i}\alpha_{i}>0$ (29), $\kappa_{1}(\epsilon)\rightarrow 0$ and $\delta(\epsilon)$ finite

[TABLE]

where $L^{1}$ denotes $L^{1}(\mathbb{T})$. We add the two energy equations (8) so as to eliminate the nonconservative terms $p\frac{d}{dt}\alpha_{i}$ since $\alpha_{1}+\alpha_{2}=1$. Then in the same way as for (30), using (5) and (28), one obtains a $L^{1}$ bound for $r_{1}e_{1}+r_{2}e_{2}$ uniform in $\epsilon$, from the positiveness of $r_{i}e_{i}=\rho_{i}\alpha_{i}e_{i}$ (29). Therefore

[TABLE]

which implies from the state laws (9) and $r_{i}=\rho_{i}\alpha_{i},0<\alpha_{i}<1$, that

[TABLE]

since $K_{i}p_{i}^{\infty}\geq 0$ and $\alpha_{i}p>0$ from (29). The state laws (9) and the positiveness of $\alpha_{i}p$ imply also from (31) that

[TABLE]

Since $\alpha_{1}+\alpha_{2}=1$ (33) implies

[TABLE]

Since $|v_{i}|\leq max(1,v_{i}^{2})$ (30,32) imply

[TABLE]

After these $L^{1}$ bounds uniform in $\epsilon$ that followed easily from the equations (6-11) and from the positiveness assumption (29), we are going to obtain $L^{\infty}$ bounds depending on $\epsilon$ that will permit to prove the existence of a global solution to the ODEs (6-11) for fixed $\epsilon$.

$\bullet$ Second step: $L^{\infty}$ bounds depending on $\epsilon$. The bound (34) and the convolutions (11) imply

[TABLE]

The basic point lies in bounds of $\|v_{i}\|_{\infty}$, which is more delicate. From (36) a proof based on (6,7) similar to the one in [13] pp. 204-205 to prove (31) there or to the more detailed one in [12] pp. 2581-2582, formulas (27-29,24) there, gives

[TABLE]

assuming absence of void regions (15) and small enough values $\epsilon>0$. This proof is sketched as follows. From (5,7) one develops at order one in $dt$ the quantity $(r_{i}v_{i})(x,t+dt,\epsilon)$ at time $t$ for an arbitrarily small positive increment $dt$, as well as $r_{i}(x,t+dt,\epsilon)$ from (5,6). Then one uses these developments to develop the quotient $v_{i}(x,t+dt,\epsilon)=\frac{(r_{i}v_{i})(x,t+dt,\epsilon)}{r_{i}(x,t+dt,\epsilon)}$. The first term we obtain for the quotient is bounded by a barycentric combination bounded by $\|v_{i}(.,t,\epsilon)\|_{\infty}$. The second term is bounded by $\frac{\|\partial_{x}\overline{p}\|_{\infty}}{m}$ from (15). One obtains

[TABLE]

where $R$ is a remainder obtained from the remainders in the first order developments of $(r_{i}v_{i})(x,t+dt,\epsilon)$ and $r_{i}(x,t+dt,\epsilon)$. This remainder $R$ disappears finally by dividing the interval $[0,t]$ into small intervals of length $dt=\frac{t}{n}$ and letting $n\rightarrow+\infty$, thus giving (37). In this proof we used assumption (15) of absence of void region in the fluid $i$.

We prefer to give a proof of (38) below simply based on (37) not to use assumption (15).

From (6,5) since $r_{i}$ and $v_{i}^{\pm}$ are positive

[TABLE]

From (37) and $r_{i}(x,0,\epsilon)\geq m>0$ an explicit solution of the elementary ODE $X^{\prime}(t)=-AX(t)+B(t)$ implies that

[TABLE]

for some value $m_{1}(\epsilon)>0$, see [12] pp. 2584,2585. Since $|\alpha_{i}v_{i}|\leq|v_{i}|$, (37) and the convolution (11) imply

[TABLE]

$\forall t\in[0,\delta(\epsilon)[$, for some values $M(\epsilon)$ independent on $t\in[0,\delta(\epsilon)[$.

Now we obtain a bound of $\|r_{i}(.,t,\epsilon)\|_{\infty}$ from (6,5) and the $L^{\infty}$ bound (37) of $v_{i}$ as follows:

[TABLE]

Gronwall formula implies

[TABLE]

The bounds (37) and (41) give a bound

[TABLE]

for some $M(\epsilon)$. Using the sum of the two energy equations (8), the $L^{\infty}$ bound (37) of $v_{i}$ for the second term in (8), the $L^{\infty}$ bounds (37,39) for $\alpha_{i}v_{i}$ and $\partial_{x}(\overline{\alpha_{i}v_{i}})$, the $L^{\infty}$ bounds (36) for $\partial_{x}\overline{p}$ and the bound $p\leq const\ r_{i}e_{i}$, i.e. a linear bound in $r_{i}e_{i}$, one obtains for $r_{1}e_{1}+r_{2}e_{2}$ a bound similar to the bound (40) obtained for $r_{i}$. Then, from Gronwall formula, one obtains the following $L^{\infty}$ bound as (41):

[TABLE]

for some $M(\epsilon)$, which implies from (29) that

[TABLE]

The minoration (38) and the $L^{\infty}$ bounds (41,42,43) show that on $[0,\delta(\epsilon)[$ the solution of the system of ODEs (6,11) for fixed $\epsilon$ remains in some set $\Omega_{M}\subset(\mathcal{C}(\mathbb{T}))^{6}$. Therefore as noticed above the solution extends to $[0,\delta(\epsilon)+\eta(\epsilon)[$ for some $\eta(\epsilon)>0$. Now we will prove that the strict positiveness property (29) is true on $[\delta(\epsilon),\delta(\epsilon)+\eta(\epsilon)[$, which will permit a proof by absurd that the solution exists on $[0,+\infty[$.

$\bullet$ Third step: strict positiveness (29) on $[\delta(\epsilon),\delta(\epsilon)+\eta(\epsilon)[$. We do not use directly assumption (15) because the following proof will also be needed in pressure. Let us prove that, for fixed $i$ and fixed $\epsilon$, $r_{i}(x,t,\epsilon)>0\ \forall x\in\mathbb{T}\forall t\in[\delta(\epsilon),\delta(\epsilon)+\eta(\epsilon)[$. By absurd let us assume that $\exists t_{0}\leq\delta(\epsilon)$ and $\exists x_{0}\in\mathbb{T}$ such that $r_{i}(x_{0},t_{0},\epsilon)=0$. We can choose them such that $r_{i}(x,t,\epsilon)>0\ \forall x\in\mathbb{T}\forall t<t_{0}$. Obviously $\partial_{t}r_{i}(x_{0},t_{0},\epsilon)\leq 0$. Now (6) gives

[TABLE]

which shows that $\partial_{t}r_{i}(x_{0},t_{0},\epsilon)\geq\kappa_{1}(\epsilon)>0$. Therefore we obtain a contradiction. We have proved by absurd that $r_{i}(x,t,\epsilon)>0\ \forall x\in\mathbb{T}\ \forall t\in[\delta(\epsilon),\delta(\epsilon)+\eta(\epsilon)[$. From the second degree equation (21) one has $\alpha_{i}(x,t,\epsilon)>0$ since the equation has exactly one root in $]0,1[$. From formula (27) the same proof by absurd applies to $\alpha_{1}p$ since $\alpha_{1}>0$: indeed since $p(x_{0},t_{0},\epsilon)=0$ all 3 possibly negative terms in the right hand-side of (27) disappear and the same reasoning applies. Therefore $p(x,t,\epsilon)>0\ \forall x\in\mathbb{T}\ \forall t\in[\delta(\epsilon),\delta(\epsilon)+\eta(\epsilon)[$; then the state laws (9) give the strict positiveness of $r_{i}e_{i}$.

$\bullet$ Fourth step: global solution of the ODEs. Now one can prove that for all $\epsilon>0$ the system of ODEs (6-11) has a global solution in positive time. By absurd let us assume the solution with strict positiveness property (29) ceases to exist at some time $\delta(\epsilon)>0$. We proved above from (29) on $[0,\delta(\epsilon)[$ that there exists $\eta(\epsilon)>0$ such that the solution can be continued on $[\delta(\epsilon),\delta(\epsilon)+\eta(\epsilon)[$ with again the strict positiveness property (29), which produces a contradiction.

Now it remains to prove that the solution of the system of ODEs satisfies the approximations (16-20).

$\bullet$ Fifth step: proof of the approximations (16-20). The proofs of (16),(17) and (18) are similar. Let us consider (17) and let $I$ denote the left hand-side of (17). One replaces $\partial_{t}(\rho_{i}\alpha_{i}v_{i})=\partial_{t}(r_{i}v_{i})$ by its value using (7,5). This gives

[TABLE]

$-(\alpha_{i}\partial_{x}\overline{p})(x,t,\epsilon)\psi(x)-(r_{i}(v_{i})^{2})(x,t,\epsilon)\psi^{\prime}(x)+(\alpha_{i}\partial_{x}\overline{p})(x,t,\epsilon)\psi(x)\}dx.$

The two terms $\alpha_{i}\partial_{x}\overline{p}$ simplify. One replaces $r_{i}v_{i}^{2}$ by $r_{i}v_{i}(v_{i}^{+}-v_{i}^{-})$. Then $I$ is the sum of a term involving $v_{i}^{+}$ and a similar term involving $v_{i}^{-}$. The term involving $v_{i}^{+}$ is

[TABLE]

After a change of variable this integral becomes

[TABLE]

Using $\|r_{i}v_{i}(.,t,\epsilon)\|_{L^{1}}\leq const$ (35) and $\|v_{i}^{+}\|_{\infty}\leq\frac{const}{\epsilon^{2\lambda}}$ (37) one obtains that the integral in (45) is equal to $O(\epsilon^{1-2\lambda})$. Same result for the term involving $v_{i}^{-}$ and finally $I=O(\epsilon^{1-2\lambda})$. This proof applies as well to (16) and (18) (for (18) the two nonconservative terms $p\frac{\partial}{\partial t}\alpha_{i}$ simplify). The proof of (19) is immediate from (34). For (20) one uses Holder’s formula, assumption (15) on absence of void region and (32):

[TABLE]

$\Box$

5. Passage to the limit in approximate solutions. We proved that the families $(\alpha_{i,\epsilon}\rho_{i,\epsilon}),(\alpha_{i,\epsilon}\rho_{i,\epsilon}e_{i,\epsilon}),(\alpha_{i,\epsilon}\rho_{i,\epsilon}v_{i,\epsilon}),(\alpha_{i,\epsilon}\rho_{i,\epsilon}v_{i,\epsilon}^{2})$ and $(p_{i,\epsilon})$ are $L^{1}$-stable (30-35). Therefore they are relatively compact for the *weak topology of the space $M(\mathbb{T}\times]0,T[)$ of Radon measures on $\mathbb{T}\times]0,T[$ [7]. Consequently there are convergent subsequences to Radon measures denoted respectively $(\alpha_{i}\rho_{i}),(\alpha_{i}\rho_{i}e_{i}),(\alpha_{i}\rho_{i}v_{i}),(\alpha_{i}\rho_{i}v_{i}^{2})$ and $p$.

If $\psi\in\mathcal{C}_{c}^{\infty}(\mathbb{R}\times]0,T[)$ it follows from (16) that

[TABLE]

in the sense of distributions. The momentum equations are slightly more complicated because of the nonconservative terms $\alpha_{i}\partial_{x}p$. Since $\rho_{i,\epsilon}\alpha_{i,\epsilon}v_{i,\epsilon}$ and $\rho_{i,\epsilon}\alpha_{i,\epsilon}v_{i,\epsilon}^{2}$ tend in the *weak topology (for a subsequence) to Radon measures it follows from (17) that

[TABLE]

where derivations and brackets are intended in the sense of distributions. Therefore $\alpha_{i,\epsilon}\partial_{x}\overline{p}_{\epsilon}$ tends to a distribution denoted $\alpha_{i}\partial_{x}\overline{p}=-\partial_{t}(\rho_{i}\alpha_{i}v_{i})+\partial_{x}(\rho_{i}\alpha_{i}v_{i}^{2})$ and the momentum equations are satisfied in the sense of distributions.

Now if further we assume that the velocities $v_{i,\epsilon}$ are $L^{\infty}$-stable (as physically needed but that we could not prove) the families $(\rho_{i,\epsilon}\alpha_{i,\epsilon}e_{i,\epsilon}v_{i,\epsilon}),i=1,2$ are $L^{1}$-stable. One can pass to the limit in a subsequence in the space of Radon measures for the *weak topology which gives a Radon measure denoted $\rho_{i}\alpha_{i}e_{i}v_{i}$. From (21) $\alpha_{i,\epsilon}$ is differentiable in $t$ valued in $L^{\infty}(\mathbb{T})$ and $\partial_{t}\alpha_{i}\in L^{\infty}(\mathbb{T})$ therefore the family $p_{\epsilon}\partial_{t}\alpha_{i,\epsilon}$is $L^{1}$-stable uniformly in $t$. Let $p\partial_{t}\alpha_{i}\in M(\mathbb{T})\times]0,T[$ be its limit (for the same subsequence as above for all the other limits). One proves that the term $(\partial_{x}\overline{p}_{\epsilon}\alpha_{i,\epsilon}v_{i,\epsilon}+p_{\epsilon}\partial_{x}(\overline{\alpha_{i,\epsilon}v_{i,\epsilon}})$ tends to a distribution from (18) since the first three terms in (18) tend to distributions. Then the energy equations are satisfied in the sense of distributions in the form

[TABLE]

in which we recall that each of the four terms is defined as a distribution by weak compactness.

If $\psi\in\mathcal{C}_{c}^{\infty}(\mathbb{T}\times[0,T[)$ is such that $\psi(.,0)\not=0$ one can state the exact equations taking into account the initial conditions as usual from integration by parts in time for $\epsilon>0$ then passage to the limit.

These Radon measures $(\alpha_{i}\rho_{i}),(\alpha_{i}\rho_{i}e_{i}),(\alpha_{i}\rho_{i}v_{i}),(\alpha_{i}\rho_{i}v_{i}^{2}),p,\dots$ are linked between themselves by the fact that, they are *weak limits of respective continuous functions $\alpha_{i,\epsilon}\rho_{i,\epsilon},\alpha_{i,\epsilon}\rho_{i,\epsilon}e_{i,\epsilon},\alpha_{i,\epsilon}\rho_{i,\epsilon}v_{i,\epsilon},\alpha_{i,\epsilon}\rho_{i,\epsilon}v_{i,\epsilon}^{2},p_{i,\epsilon},\dots$ when $\epsilon\rightarrow 0$. Indeed such products of Radon measures or of Radon measures and $L^{\infty}$ functions do not make sense mathematically in general. Equations (1-3) are satisfied as a sum of a few distributions obtained by *weak compactness whose sum is the null distribution: the individual products inside the notation of these terms only recall the origin of these terms from a weakly convergent sequence of approximate solutions.

This definition of these Radon measures by passage to the limit of corresponding smooth objects when $\epsilon\rightarrow 0$ appears justified from physics since the molecular structure of fluids forbids that the value $dx=\epsilon$ used to state the conservation laws could be too small: the cells $\Pi_{1\leq j\leq 3}[x_{j}-\frac{\epsilon}{2},x_{j}+\frac{\epsilon}{2}]$ should contain a certain amount of molecules so that the statement of the conservation laws could make sense physically. Therefore the ODEs (6-11) model physics for $\epsilon>0$ very small, but not arbitrarily small, and the limit case $\epsilon=0$ obtained from compactness is only an approximation of the physical situation. This justifies that the above Radon measures are only defined by limits of the corresponding smooth functions and that the apparent products inside them make no sense if not to recall their physical interpretation.

In a mathematical viewpoint our method is an inversion of order in the operations of passage to the limit $\epsilon\rightarrow 0$ (in which $\epsilon$ is the side of the small cube $\Pi_{1\leq i\leq n}[x_{i}-\frac{\epsilon}{2},x_{i}+\frac{\epsilon}{2}]$ used by physicists to obtain the equations) and resolution of the equations. Indeed to state conservation laws physicists first state (6-8) (with $\frac{d}{dt}w(x,t,\epsilon)$ replaced by $\frac{w(x,t+dt,\epsilon)-w(x,t,\epsilon)}{dt})$, then they pass formally to the limit $\epsilon=0$ to obtain (1-3). Then one attempts to solve directly (1-4). In our method we solve (6-11) for fixed $\epsilon$ then pass to the limit on the solutions to express that $\epsilon>0$ is very small. Finally this method provides Radon measures that satisfy a weak formulation of (1-3).

6. Taking into account molecular agitation. To simplify the notation let us consider the case of a single fluid and formula (1) in the familiar form

[TABLE]

Then formulas (5,6) give

[TABLE]

One considers the 3 cells $\mathcal{C}_{y}=]y-\frac{\epsilon}{2},y+\frac{\epsilon}{2}[$ with $y=x,x-\epsilon$ and $x+\epsilon$ with respective constant values $\rho(y,t,\epsilon)$ and $v(y,t,\epsilon)$ inside each cell. Formula (48) describes the transport of the matter around the cell $\mathcal{C}_{x}$: if one considers the cell interface $x-\frac{\epsilon}{2}$ then in the time interval $[t,t+dt]$ the amount of matter $\rho v^{+}(x-\epsilon,t,\epsilon)dt$ crosses this interface to the right and the amount of matter $\rho v^{-}(x,t,\epsilon)dt$ crosses this interface to the left. Therefore (48) models the transport of matter according to its macroscopic velocity $v$ as described in [9] p. 1906 provided the condition $|v|dt<\epsilon$ to eliminate the influence of more remote cells during time $dt$.

In various circumstances the mathematical proofs imposed or suggested the replacement of $v^{\pm}$ by $v^{\pm}+\mu$ with $\mu>0$ large enough. This was done in [12] pp. 2586-2588 and in [1] pp. 1208-1211. The physical meaning of this replacement is that between times $t$ and $t+dt$ an additional amount of matter $\mu\rho(x-\epsilon,t,\epsilon)dt$ crosses the interface $x-\frac{\epsilon}{2}$ to the right and an additional amount of matter $\mu\rho(x,t,\epsilon)dt$ crosses the interface $x-\frac{\epsilon}{2}$ to the left. This clearly can be interpreted as a model of molecular agitation with mean velocity value (in all senses: to the right and same to the left) of molecules equal to $\mu$, on which the macroscopic velocity $v$ is superposed.

This change of $v^{\pm}$ into $v^{\pm}+\mu$ does not modify the proof in this paper. One can even change $v^{\pm}$ into $v^{\pm}+\frac{\mu}{\epsilon^{\beta}},\ 0<\beta<1$: then one simply notices that in (45) $\|v_{i}^{\pm}\|_{\infty}$ becomes $\|v_{i}^{\pm}\|_{\infty}+\frac{\mu}{\epsilon^{\beta}}$ which gives a bound $O(\epsilon^{1-2\lambda})+O(\epsilon^{1-\beta})$ for the right hand side of (45). This change could play a role in a search of a physically admissible solution since it permits to prove that the method in this paper gives the Kruzhkov entropy solution [1]. This change also plays a role numerically by adding some vanishing viscosity, see [12] pp. 2586,2587. Indeed it has been noticed that the presence of $\mu>0$ large enough is needed in the numerical tests below: in figures 1 and 2 one uses $\mu=300$.

7. Numerical calculation of the solutions. In this section we calculate numerically the solution in the case of the Toumi shock tube problem [41, 39, 37, 38]. In figure 1 we compare the approximate solutions obtained from our theoretical construction with the result of a standard ”‘transport-averaging-pressure correction”’ splitted scheme which is a direct extension of the scheme in [13] section 7 and is described in [14]. We obtain an exact superposition. In figure 2 we compare the approximate solutions from our theoretical construction with the numerical solution obtained from standard numerical methods of scientific computing in [39] pp. 2615,2617, [37] p. 437, [38] p. 497, in which these authors insert into system (1-4) an additional term to render it hyperbolic to improve its mathematical and numerical properties, [39] p. 2595. We observe close results in pressure, gas temperature and liquid temperature (same step values and same location of the discontinuities). In gas volume fraction, gas velocity and liquid velocity we still observe same step values and same location of discontinuities, but in these three cases we observe a significant difference in form of peaks: a ”‘down-peak”’ in gas volume fraction and two very neat ”‘up-peaks”’ of well defined height in gas and liquid velocities. The features of these peaks are well defined: they are invariant under changes of the time steps and the CFL number. They are identical for the two different schemes tested in figure 1. Therefore these peaks are not artefacts of calculation since they are obtained from two completely different numerical numerical methods. Further quite similar peaks in shape and values are observed in engineering codes and observations on the gas kick [2, 3, 6, 8, 36].

By testing both absence and presence of the additional term introduced by many authors to render the system hyperbolic by the same standard ”‘transport-averaging-pressure correction”’ order 1 numerical scheme adapted from [9] and [13] section 7, see [14], one shows that the additional term is responsible of the disappearance of these peaks: without the additional term the scheme produces the result in figure 1, with the additional term the scheme produces the result in [37, 38, 39]. The additional term is motivated to ensure hyperbolicity of the system. However in some cases of fully nonlinear systems of physics it has been observed in [22]-[31] that hyperbolicity is not always indispensable to produce a well-posed solution.

From the experimental result reproduced in [2], figure 4, the peak obtained from the asymptotic solution in this paper for system (1-4) and, also in absence of the additional term, by the ”‘transport-averaging-pressure correction”’ scheme of [14] appears to be related to the gas kick phenomenon. In figure 3 we observe the evolution of the peak in liquid flow rate according to time. Its velocity, top value and width appear of an order of magnitude completely compatible with observations on the gas kick, see [2, 3, 6, 8, 36] and articles quoted there. These authors report experimental results and results from numerical codes describing the gas and liquid flow rates ($m^{3}/s$ or $kg/s$) at the top of the well as a function of time. The results reported in this paper concern solutions of the Riemann problem i.e. they represent the gas and liquid physical variables (in particular velocities) inside the tube at a given time. From them one can at once obtain the flow rates at the end of the tube as a function of time: constant velocity before the peaks arrive at the end of the tube, then peaks for the flow rates as a function of time when the”’ kick”’ goes out of the tube, completely similar to the results reported by engineers [2, 3, 6, 8, 36] and articles quoted there. But the experimental results concern real situations with two coaxial tubes and the gas kick comes from the annular domain between the tubes etc, therefore a detailed comparison is not possible within the scope of this mathematical paper.

INSERT FIGURE 1

Figure 1. *Comparison of the asymptotic solution constructed in this paper (black +) and the result from the transport-averaging-correction scheme of [14] with $\delta=0$ (red, continuous line). One observes a perfect coincidence.

INSERT FIGURE 2

Figure 2. *Comparison of the asymptotic solution constructed in this paper (black,+) with the numerical solution constructed by various authors (red, continuous line) by adding an additional term in the equations. This last solution is obtained here with the schemes of [14] after having checked it gives same results.

INSERT FIGURE 3

Figure 3. *The liquid flow rate ($kg/s$) at time t=0.03 (top left panel), t=0.06 (top right panel), t=0.12 (bottom panel, in a twice longer tube). One observes that the top value increases very slowly with time, and that its width increases proportionally to time: one interval=2.5 centimeter in each panel.

In the three figures the system of ODEs (6-8) is solved by the explicit Euler order 1 scheme. Space which is 100 meters long is divided into 4000 cells (8000 cells in the bottom panel of figure 3 where we consider a 200 meters long tube); the solution is given at time $t=0.06,r=\frac{\Delta t}{\Delta x}=0.0005$ in all tests, $\mu=300$ in figures 1 and 2 and $\mu=$200, 300 and 400 in figure 3 from top left panel to bottom panel with various values of time (t=0.03,0.06,0.12)$,\lambda=1$ in (11) with $\epsilon$ equal to the space step $\Delta x$.

Remark: the case of one fluid. One has checked that the adaptation of the method in this paper to the case of one fluid has exactly given numerically the known solutions on the four 1-D Toro tests in [40] and (easy extension of the proof to $n$-D, $n=2,3,\dots$ as explained in [12]) the six 2-D Lax tests in [34, 35].

8. Conclusion. Since various numerical schemes giving same result for conservative systems can give really different results for nonconservative systems such as (1-4) the search of mathematically well defined solutions is particularly important. In this paper we have proved that one can obtain approximate and Radon measure solutions in a weak sense. By reducing the system of partial differential equations in the case of these approximate solutions to a system of ordinary differential equations we observe numerically on the standard Toumi shock tube problem that the Radon measures from our method agree with the numerical solutions previously obtained by other authors with various different numerical methods, modulo very neat peaks of well defined limited height and width in liquid and gas velocities, which suggest the gas kick phenomenon that appear in our method, engineering codes and experimental observations. Indeed in a subsequent numerical paper, we observe exactly the numerical solution given by the theoretical mathematical proof presented in the present paper.