follow Hérard’s original work [25] for barotropic fluids and choose for the
mathematical entropy function the total specific energy
$E_{k}=e_{k}(\rho_{k})+0.5{\mbox{\boldmath$v$}}_{k}^{2}$ with internal energy
$e_{k}$ determined by

$\displaystyle
e_{k}^{\prime}(\rho_{k})=\frac{p_{k}(\rho_{k})}{\rho_{k}^{2}}\quad\mbox{or,
equiv., }\quad e_{k}^{\prime}(\tau_{k})=-p_{k}(\tau_{k})$ (4.57)

motivated by (2.14). For smooth solutions of the model (2.2), (2.4) we obtain
the evolution equations for the fractional total energies

$\displaystyle{\partial_{t}(\alpha_{k}\,\rho_{k}\,E_{k})}+{\nabla\cdot(\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k}\,(E_{k}+p_{k}/\rho_{k}))}=\overline{S}_{\alpha\rho
E,k}+{\overline{\Pi}}_{k},$ (4.58)

with relaxation term $\overline{S}_{\alpha\rho E,k}$ and production term
${\overline{\Pi}}_{k}$

$\displaystyle\overline{S}_{\alpha\rho E,k}$
$\displaystyle=(e_{k}+p_{k}/\rho_{k}-0.5{\mbox{\boldmath$v$}}_{k}^{2})S_{\alpha\rho,k}-p_{k}S_{\alpha,k}+{\mbox{\boldmath$v$}}_{k}\cdot{\mbox{\boldmath$S$}}_{\alpha\rho\mbox{\scriptsize\boldmath${\mbox{\boldmath$v$}}$},k},$
(4.59) $\displaystyle{\overline{\Pi}}_{k}$
$\displaystyle=p_{k}({\mbox{\boldmath$V$}}_{I}-{\mbox{\boldmath$v$}}_{k})\cdot\nabla\alpha_{k}-\sum_{l=1,\neq
k}^{K}P_{k,l}{\mbox{\boldmath$v$}}_{k}\cdot\nabla\alpha_{l}$
$\displaystyle=-\sum_{l=1,\neq
k}^{K}\left(p_{k}({\mbox{\boldmath$V$}}_{I}-{\mbox{\boldmath$v$}}_{k})+P_{k,l}{\mbox{\boldmath$v$}}_{k}\
\right)\cdot\nabla\alpha_{l}.$ (4.60)

We emphasize that the relaxation term $\overline{S}_{\alpha\rho E,k}$ is
different from the corresponding term determined by the relaxation models
(2.45), (2.46), (2.48) and (2.50). It is interesting to note that in contrast
to the single-phase case the energy equation is not redundant in the multi-
component case due to the non-conservative products and the relaxation terms.
Similar to the isothermal case, see Eqn. (4.19), we determine the heat flux as

$\displaystyle{\nabla\cdot(\alpha_{k}\,{\mbox{\boldmath$q$}}_{k})}=S_{\alpha\rho
E,k}-\overline{S}_{\alpha\rho E,k}-\sum_{l=1,\neq
k}^{K}P_{k,l}\,{{\mbox{\boldmath$V$}}_{I}}\cdot{\nabla\alpha_{l}}-{\overline{\Pi}}_{k}$
(4.61)

comparing the difference of (2.8) and (4.58).

Summing these evolution equations over all components we obtain the evolution
equation of the total energy

$\displaystyle{\partial_{t}\left(\sum_{k=1}^{K}\alpha_{k}\,\rho_{k}\,E_{k}\right)}+{\nabla\cdot\left(\sum_{k=1}^{K}(\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k}\,(E_{k}+p_{k}/\rho_{k}))\right)}=\sum_{k=1}^{K}\overline{S}_{\alpha\rho
E,k}+\sum_{k=1}^{K}{\overline{\Pi}}_{k},$ (4.62)

Comparing (4.62) with (4.4) we may choose for the entropy-entropy flux pair

$\displaystyle
U({\mbox{\boldmath$u$}}):=\sum_{k=1}^{K}\alpha_{k}\rho_{k}E_{k},$ (4.63a)
$\displaystyle
F_{i}({\mbox{\boldmath$u$}}):=\sum_{k=1}^{K}(\alpha_{k}\,\rho_{k}\,v_{k,i}\,(E_{k}+p_{k}/\rho_{k})),\
i=1,\ldots,d.$ (4.63b)

The compatibility conditions (4.3) hold by derivation of (4.62). In the
isentropic case non-strict convexity of $U$ has been proven in [54].
Alternatively, the proof of Theorem 4.1 for the isothermal case also applies
in case of isentropic fluids where

$U({\mbox{\boldmath$u$}})=\sum_{k=1}^{K}\alpha_{k}\rho_{k}(e_{k}+{\mbox{\boldmath$v$}}_{k}^{2}/2)=\sum_{k=1}^{K}\alpha_{k}U_{k}({\mbox{\boldmath$u$}}_{k})$
(4.64)

with $e_{k}=e_{k}(\tau_{k},s_{k})$ and $s_{k}=const$. For the Hessian
${\mbox{\boldmath$U$}}_{k}^{\prime\prime}$ we again obtain (4.34). Then
convexity of $e_{k}$ with respect to $\tau_{k}$ according to (2.16) ensures
that the Hessian is positive semi-definite.

We emphasize that in case of an isothermal fluid the mixture total energy is
not an admissible entropy function because (4.58) does not account for heat
exchange. To ensure the entropy inequality (4.6) we have to verify that the
right-hand side of (4.62) is not positive. This will provide us with
admissibility criteria for the interfacial pressures as well as the relaxation
terms.

####  4.2.1 Entropy production due to interfacial velocity and pressure

Obviously, we cannot control the sign of $\Pi_{k}$. Therefore we determine the
interfacial pressures $P_{k,l}$ and the interfacial velocity
${\mbox{\boldmath$V$}}_{I}$ such that the sum
${\overline{\Pi}}:=\sum_{k=1}^{K}{\overline{\Pi}}_{k}$ vanishes and the
conservation constraints (2.10) or, equivalently, (2.11) hold. Assuming that
the interfacial velocity can be written as a linear combination of the
velocities of the components according to (3.42) there exists a unique closure
for the interfacial pressures.

######  Theorem 4.6.

(Entropy production due to interfacial states) For any convex combination
(3.42) for the interfacial velocity ${\mbox{\boldmath$V$}}_{I}$ there uniquely
exist interfacial pressures

$P_{k,l}=\beta_{k}p_{l}+p_{k}(1-\beta_{k})$ (4.65)

such that the production term
${\overline{\Pi}}=\sum_{k=1}^{K}{\overline{\Pi}}_{k}$ vanishes and the
conservation constraints (2.10) or, equivalently, (2.11) hold. In particular,
the interfacial pressures are all positive and

$P_{I}=\sum_{k=1}^{K}p_{k}(1-\beta_{k}).$ (4.66)

###### Proof.

First of all, we note that by the conservation constraint (2.10) and the
definition (4.2) of ${\overline{\Pi}}_{k}$ it holds

$\displaystyle{\overline{\Pi}}=\sum_{k=1}^{K}{\overline{\Pi}}_{k}=\sum_{k=1}^{K}\sum_{l=1,\neq
k}^{K}(P_{k,l}-p_{k})({\mbox{\boldmath$V$}}_{I}-{\mbox{\boldmath$v$}}_{k})\cdot{\nabla\alpha_{l}}.$
(4.67)

Obviously, this term coincides with $T_{ref}\Pi=\sum_{k=1}^{K}\Pi_{k}$ with
$\Pi_{k}$ determined by (4.16) in the isothermal case. Since by the saturation
constraint (2.1) we may substitute ${\nabla\alpha_{K}}$ by the linearly
independent gradients ${\nabla\alpha_{k}}$, $k=1,\ldots,K-1$, we may rewrite
${\overline{\Pi}}$ in terms of the linearly independent gradients. Then
${\overline{\Pi}}$ vanishes if and only if the cofactors vanish, i.e.,

$\sum_{k=1,k\neq
l}^{K-1}\left((P_{k,l}-p_{k})-(P_{k,K}-p_{k})\right)({\mbox{\boldmath$v$}}_{k}-{\mbox{\boldmath$V$}}_{I})+(P_{K,l}-p_{K})({\mbox{\boldmath$v$}}_{K}-{\mbox{\boldmath$V$}}_{I})={\mbox{\boldmath$0$}},\quad
l=1,\ldots,K-1.$ (4.68)

Next we use that the interfacial velocity ${\mbox{\boldmath$V$}}_{I}$ is a
convex combination of the single component velocities
${\mbox{\boldmath$v$}}_{k}$, i.e., (3.42) holds. Then we may rewrite the
velocity differences ${\mbox{\boldmath$v$}}_{k}-{\mbox{\boldmath$V$}}_{I}$ in
terms of the differences
${\mbox{\boldmath$v$}}_{i}-{\mbox{\boldmath$v$}}_{i+1}$, $i=1,\ldots,K-1$.
Rearranging (4.68) in terms of these yields $(K-1)^{2}$ equations for the
$(K-1)K$ unknowns $P_{k,l}$, $k,l=1,\ldots,K$, $k\neq l$. Additional $K-1$
equations are provided by the conservation constraints (2.11). Altogether we
obtain a linear system of equations for the interfacial pressures. This system
is identical to the linear system derived in [41], Eqns. (6.20) and (6.22),
where the temperatures are formally replaced by 1. This system has a unique
solution. The proof is identical to the one of the analogue theorem in [41] in
case of the full Baer-Nunziato model.

$\square$

We would like to mention that in [25] the interfacial pressures in case of
barotropic fluids are given for $K=2,3,4$. The proof is shortly sketched
without investigating the regularity of the resulting linear system.
Furthermore, as in the isothermal case the interfacial pressures (4.65)
coincide with those presented in [41] for the full multi-component Baer-
Nunziato model at local thermal equilibrium. Thus, using for the convex
combination (3.42) the coefficients (3.43), then the interfacial pressure and
velocity again are given by (4.37), (4.38) and (4.39).

####  4.2.2 Entropy production due to relaxation

To be admissible with the entropy inequality (4.6) the production
$\sum_{k=1}^{K}\overline{S}_{\alpha\rho E,k}$ has to be non-positive. Note
that in abuse of physical notation we call this term entropy production rather
than energy production with respect to the mathematical concept of entropy-
entropy flux pairs because energy is known to be conserved. In the following
we consider one-by-one the contributions of pressure and velocity relaxation.
The analysis is similar to [41] in case of the physical entropy production.
Note that all the relaxation processes are in agreement with the constraints
due to conservation (2.9).

######  Theorem 4.7.

(Entropy production due to pressure relaxation) Let the relaxation parameter
be non-negative, i.e., $\theta_{p}\geq 0$. Furthermore, let the coefficients
$d_{k,l}$ be positive and symmetric, i.e.

$\displaystyle 0\leq d_{k,l}=d_{l,k},\quad l,k=1,\ldots,K,\ l\neq k.$ (4.69)

Then the entropy production due to pressure relaxation (2.47a) is non-
positive, i.e.,

$S_{U}^{p}=\sum_{k=1}^{K}S_{U,k}^{p}\leq 0$ (4.70)

with entropy production of component $k=1,\ldots,K$

$\displaystyle S_{U,k}^{p}=\overline{S}_{\alpha\rho
E,k}^{p}=-\theta_{p}p_{k}\sum_{l=1}^{K}d_{k,l}(p_{k}-p_{l}).$ (4.71)

In particular, this holds true for the model (2.45).

###### Proof.

First of all, we note that by symmetry of arbitrary coefficients
$a_{k,l}=a_{l,k}$ it holds

$\displaystyle\sum_{k=1}^{K-1}\sum_{l=k+1}^{K}a_{k,l}=\sum_{l=2}^{K}\sum_{k=1}^{l-1}a_{k,l}=\sum_{k=2}^{K}\sum_{l=1}^{k-1}a_{l,k}$

From this we conclude with $a_{k,l}=p_{k}^{2}d_{k,l}$ and
$a_{k,l}=p_{k}d_{k,l}p_{l}$

$\displaystyle\sum_{l,k=1,l\neq
k}^{K}p_{k}^{2}d_{k,l}=\sum_{k=2}^{K}\sum_{l=1}^{k-1}p_{k}^{2}d_{k,l}+\sum_{k=1}^{K-1}\sum_{l=k+1}^{K}p_{k}^{2}d_{k,l}=\sum_{k=2}^{K}\sum_{l=1}^{k-1}(p_{k}^{2}+p_{l}^{2})d_{k,l},$
$\displaystyle\sum_{l,k=1,l\neq
k}^{K}p_{k}p_{l}d_{k,l}=\sum_{k=2}^{K}\sum_{l=1}^{k-1}p_{k}p_{l}d_{k,l}+\sum_{k=1}^{K-1}\sum_{l=k+1}^{K}p_{k}p_{l}d_{k,l}=2\sum_{k=2}^{K}\sum_{l=1}^{k-1}p_{k}p_{l}d_{k,l}.$

Combining these results we obtain

$\displaystyle\sum_{l,k=1,l\neq
k}^{K}p_{k}d_{k,l}(p_{k}-p_{l})=\sum_{k=2}^{K}\sum_{l=1}^{k-1}d_{k,l}(p_{k}^{2}-2p_{k}p_{l}+p_{l}^{2})=\sum_{k=2}^{K}\sum_{l=1}^{k-1}d_{k,l}(p_{k}-p_{l})^{2}.$

Obviously, this term is non-negative and vanishes only if all pressures
coincide. Thus, by assumption the energy production $\overline{S}_{\alpha\rho
E,k}^{p}$ is non-positive. In particular, this holds true for the special
model (2.45) because $d_{k,l}=\alpha_{k}\alpha_{l}$ satisfies the constraint
(4.69)

$\square$

Since the proof holds for an arbitrary number of components $K$ the theorem
generalizes the result in [25] for $K=4$. For $K=2$ the result is well-known.

For the velocity relaxation we obtain a similar result.

######  Theorem 4.8.

(Entropy production due to velocity relaxation) Let the relaxation parameter
be non-negative, i.e., $\theta_{v}\geq 0$. Furthermore, let the coefficients
$e_{k,l}$ be non-negative and symmetric, i.e.

$\displaystyle 0\leq e_{k,l}=e_{l,k},\quad l,k=1,\ldots,K,\ l\neq k.$ (4.72)

Then the entropy production due to velocity relaxation (2.47b) is non-
positive, i.e.,

$S_{U}^{v}=\sum_{k=1}^{K}S_{U,k}^{v}\leq 0$ (4.73)

with entropy production of component $k=1,\ldots,K$

$\displaystyle S_{U,k}^{v}=\overline{S}_{\alpha\rho
E,k}^{v}=-\theta_{v}\sum_{l=1}^{K}e_{k,l}\,{\mbox{\boldmath$v$}}_{k}\cdot({\mbox{\boldmath$v$}}_{k}-{\mbox{\boldmath$v$}}_{l}).$
(4.74)

In particular, this holds true for the model (2.46).

###### Proof.

The proof is in full analogy to the pressure relaxation where we replace
$d_{k,l}$, $p_{k}$ and $p_{l}$ by $e_{k,l}$, $v_{k,i}$ and $v_{l,i}$,
respectively. Then we obtain

$\displaystyle\sum_{l,k=1,l\neq
k}^{K}e_{k,l}{\mbox{\boldmath$v$}}_{k}\cdot({\mbox{\boldmath$v$}}_{k}-{\mbox{\boldmath$v$}}_{l})=\sum_{i=1}^{d}\sum_{k=2}^{K}\sum_{l=1}^{k-1}e_{k,l}(v_{k,i}-v_{l,i})^{2}=\sum_{k=2}^{K}\sum_{l=1}^{k-1}e_{k,l}({\mbox{\boldmath$v$}}_{k}-{\mbox{\boldmath$v$}}_{l})^{2}.$

Obviously, this term is non-negative and vanishes only if all velocities
coincide. Thus, by assumption the energy production $\overline{S}_{\alpha\rho
E,k}^{v}$ is non-positive. In particular, this holds true for the special
model (2.46) because
$e_{k,l}=\alpha_{k}\rho_{k}\alpha_{l}\rho_{l}/\sum_{j=1}^{K}\alpha_{j}\rho_{j}$
satisfies the constraint (4.72).

$\square$

Since there is no temperature relaxation and no relaxation of chemical
potentials, there is no corresponding energy production.

##  5 Relaxation to equilibrium

The main objective of our work is to develop new, efficient and robust
relaxation procedures. For this purpose we split the relaxation process from
the fluid motion by formally performing an operator splitting. Since the
relaxation times differ for the different relaxation types, we solve the
initial value problem

${\frac{d\,{\mbox{\boldmath$w$}}(t)}{dt}}={\mbox{\boldmath$S$}}^{\xi}({\mbox{\boldmath$w$}}(t)),\quad
t\in[t_{n},t_{n+1}],\qquad{\mbox{\boldmath$w$}}(t_{n})={\mbox{\boldmath$w$}}^{0}$
(5.1)

with
${\mbox{\boldmath$S$}}^{\xi}=((S^{\xi}_{\alpha_{k}})_{k=1,\ldots,K-1},(S^{\xi}_{\rho,k},{\mbox{\boldmath$S$}}^{\xi}_{\rho{\mbox{\boldmath$v$}},k})_{k=1,\ldots,K})^{T}$
separately for each relaxation type $\xi$. To avoid the explicit computation
of the relaxation times $\theta_{\xi}$ in (2.45), (2.46), (2.48) and (2.50),
we may perform the change of variables ${\overline{t}}:=(t_{n}-t)\theta_{\xi}$
and
${\overline{{\mbox{\boldmath$w$}}}}({\overline{t}}):={\mbox{\boldmath$w$}}(t)$
in (5.1), i.e.,

$\frac{d\,{\overline{{\mbox{\boldmath$w$}}}}({\overline{t}})}{d{\overline{t}}}=\frac{1}{\theta^{\xi}}{\mbox{\boldmath$S$}}^{\xi}({\overline{{\mbox{\boldmath$w$}}}}({\overline{t}})),\quad{\overline{t}}\in[0,\Delta
t\theta_{\xi}],\qquad{\overline{{\mbox{\boldmath$w$}}}}(0)={\overline{{\mbox{\boldmath$w$}}}}^{0}.$
(5.2)

Then by definition of the relaxation terms, the relaxation parameter cancels
on the right-hand side of (5.2). Since for all relaxation processes the
conservation constraints (2.9) are satisfied we conclude from (5.1) that

$\sum_{k=1}^{K}\frac{d\,{\overline{\alpha}}_{k}({\overline{t}})}{d{\overline{t}}}=0,\
\frac{d\,{\overline{\rho}}({\overline{t}})}{d{\overline{t}}}=0,\
\frac{d\,({\overline{\rho}}\,{\overline{{\mbox{\boldmath$v$}}}})({\overline{t}})}{d{\overline{t}}}={\mbox{\boldmath$0$}}.$
(5.3)

Hence, the bulk quantities for mixture density, momentum and energy as well as
the saturation condition (2.1) remain constant during the relaxation process
and it holds

$\sum_{k=1}^{K}\alpha_{k}^{\infty}=\sum_{k=1}^{K}\alpha_{k}^{0}=1,\
\rho^{\infty}=\rho^{0},\
(\rho\,{\mbox{\boldmath$v$}})^{\infty}=(\rho\,{\mbox{\boldmath$v$}})^{0}.$
(5.4)

Furthermore, the conservation of bulk mass and momentum imply that the bulk
velocity remains constant

$\frac{d\,{\overline{{\mbox{\boldmath$v$}}}}}{d{\overline{t}}}={\mbox{\boldmath$0$}}$
(5.5)

and, thus, it holds

${\mbox{\boldmath$v$}}^{\infty}={\mbox{\boldmath$v$}}^{0}.$ (5.6)

We are only interested in the equilibrium state rather the transient
relaxation behavior as considered, for instance, in [27, 28]. Therefore we
assume that the relaxation process is infinitely fast, i.e.,
$\theta_{\xi}\to\infty$ and we do not need to determine the relaxation
parameters. To compute the equilibrium state, where the source terms (2.45),
(2.46), (2.48) and (2.50) vanish, we perform integration of the ODE system to
infinity. This results in a system of algebraic equations for the equilibrium
state that will be derived and solved in the subsequent sections. For ease of
notation we will use $w$ instead of ${\overline{{\mbox{\boldmath$w$}}}}$ in
(5.2). In the following we consider one by one the different relaxation
processes to equilibrium. Note that the equilibrium state is independent of
the order of the relaxation procedures.

Velocity relaxation. First of all, we investigate velocity relaxation. This
relaxation process is the same for isentropic and isothermal fluids. The
procedure to determine the velocity equilibrium we proceed similar to Saurel
and Le Métayer in [56] for the full Baer-Nunziato model. The equilibrium state
of the velocity relaxation process is determined by solving the system of ODEs

$\displaystyle\frac{d\,\alpha_{k}}{d{\overline{t}}}=0,$ (5.7a)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}}{d{\overline{t}}}=0,$ (5.7b)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k}}{d{\overline{t}}}=\alpha_{k}\,\rho_{k}({\mbox{\boldmath$v$}}-{\mbox{\boldmath$v$}}_{k})$
(5.7c)

resulting from (5.2) with the source terms (2.46). Since we assume that the
relaxation process is infinitely fast, the solution of the system of ODEs
converges towards the equilibrium state where the right-hand side vanishes.
This holds true for

${\mbox{\boldmath$v$}}_{k}={\mbox{\boldmath$v$}}^{\infty},\quad{\overline{t}}\to\infty.$
(5.8)

Integration of (5.7a), (5.7b) over $[0,\infty]$ then results in the algebraic
equations

$\displaystyle\alpha_{k}^{\infty}=\alpha_{k}^{0},$ (5.9a)
$\displaystyle\alpha_{k}^{\infty}\,\rho_{k}^{\infty}=\alpha_{k}^{0}\,\rho_{k}^{0}\quad\mbox{equiv.}\quad\rho_{k}^{\infty}=\rho_{k}^{0}.$
(5.9b)

From the equilibrium condition (5.8) and (5.6) we conclude that

${\mbox{\boldmath$v$}}_{k}^{\infty}={\mbox{\boldmath$v$}}^{\infty}={\mbox{\boldmath$v$}}^{0}.$
(5.10)

Finally, we end up with the algebraic system (5.6), (5.9a), (5.9b) and (5.10)
by which we determine the velocity equilibrium state. The procedure is
summarized in the following algorithm.

######  Algorithm 5.1 (Velocity relaxation procedure).

  1. 1.

Compute the equilibrium mixture velocity according to (2.12) as

${\mbox{\boldmath$v$}}^{\infty}=\left.\sum_{k=1}^{K}\alpha_{k}^{0}\rho_{k}^{0}{\mbox{\boldmath$v$}}_{k}^{0}\right/\sum_{k=1}^{K}\alpha_{k}^{0}\rho_{k}^{0};$
(5.11)

  2. 2.

Update the velocities of the components using (5.10);

For pressure relaxation, thermal relaxation and chemical relaxation we have to
distinguish between isentropic fluids and isothermal fluids.

###  5.1 Isothermal case

In this section we always assume for each component an isothermal stiffened
gas, see Example 2.2.

####  5.1.1 Pressure relaxation

For the pressure relaxation process we have to solve the system of ODEs

$\displaystyle\frac{d\,\alpha_{k}}{d{\overline{t}}}=\alpha_{k}\,(p_{k}-p),$
(5.12a) $\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}}{d{\overline{t}}}=0,$
(5.12b)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k}}{d{\overline{t}}}=0$
(5.12c)

resulting from (5.2) with the source terms (2.45). Since we assume that the
relaxation process is infinitely fast, the solution of the system of ODEs
converges towards the equilibrium state where the right-hand side vanishes.
This holds true for

$p_{k}^{\infty}=p^{\infty},\qquad{\overline{t}}\to\infty.$ (5.13)

Integration of (5.12b), (5.12c) over $[0,\infty]$ then results in the
algebraic equations

$\displaystyle\alpha_{k}^{\infty}\rho_{k}^{\infty}=\alpha_{k}^{0}\rho_{k}^{0},$
(5.14a)
$\displaystyle\mathbf{v}_{k}^{\infty}=\mathbf{v}_{k}^{0},\,k=1,\dots,K.$
(5.14b)

In the equilibrium state it must hold

$\displaystyle p_{k}^{\infty}=p^{\infty},\,k=1,\dots,K.$ (5.15)

In addition, for an isothermal fluid all temperatures are in equilibrium,
i.e., $T_{k}=T_{ref}$, $k=1,\dots,K$. Assume all components to be stiffened
gases than it holds

$\displaystyle p_{k}=C_{k}\rho_{k}-\pi_{k},$ (5.16)

where $C_{k}=\tilde{c}_{k}^{2}=Tc_{v,k}(\gamma_{k}-1)>0$. From the saturation
condition

$\displaystyle\sum_{k=1}^{K}\alpha_{k}=1$ (5.17)

and mass conservation (5.14a) we obtain

$\displaystyle\prod_{k=1}^{K}\rho_{k}^{\infty}=\sum_{k=1}^{K}(\alpha_{k}^{0}\rho_{k}^{0})\prod_{j=1,j\neq
k}^{K}\rho_{j}^{\infty}.$ (5.18)

Let us fix an arbitrary component $k_{0}$. Then we obtain from (5.16)

$\displaystyle\rho_{k}^{\infty}=\frac{C_{k_{0}}}{C_{k}}\rho_{k_{0}}^{\infty}+\frac{\pi_{k}-\pi_{k_{0}}}{C_{k}},\,k=1,\dots,K,\,k\neq
k_{0}.$ (5.19)

Putting (5.19) into (5.18) we get a polynomial for $\rho_{k_{0}}^{\infty}$

$\displaystyle\prod_{k=1}^{K}\left(\frac{C_{k_{0}}}{C_{k}}\rho_{k_{0}}^{\infty}+\frac{\pi_{k}-\pi_{k_{0}}}{C_{k}}\right)=\sum_{k=1}^{K}(\alpha_{k}^{0}\rho_{k}^{0})\prod_{j=1,j\neq
k}^{K}\left(\frac{C_{k_{0}}}{C_{j}}\rho_{k_{0}}^{\infty}+\frac{\pi_{j}-\pi_{k_{0}}}{C_{j}}\right).$
(5.20)

In order to obtain a polynomial of smallest possible degree we introduce the
following sets.

$\displaystyle{\cal{E}}_{k_{0}}:=\\{k\,:\,\pi_{k}=\pi_{k_{0}},\,k\in\\{1,\dots,K\\}\\},\quad{\cal{N}}_{k_{0}}=\\{k\,:\,\pi_{k}\neq\pi_{k_{0}},\,k\in\\{1,\dots,K\\}\\}.$
(5.21)

Depending on the choice of $k_{0}$ the polynomial to determine the equilibrium
density $\rho_{k_{0}}^{\infty}$ is of degree $\\#{\cal{N}}_{k_{0}}+1$.
Accordingly we recommend $k_{0}$ to be chosen such that the cardinality of
${\cal{N}}_{k_{0}}$ is smallest.

Special case 1. Let us assume that $\pi_{1}=\pi_{2}\dots=\pi_{K}=:\pi$. Then
there exists only one root of the polynomial (5.20). Applying (5.19) the
equilibrium phasic densities $\rho_{j}^{\infty}$ are determined by

$\displaystyle\rho_{j}^{\infty}=\frac{1}{C_{j}}\sum_{k=1}^{K}\alpha_{k}^{0}\rho_{k}^{0}C_{k}=\frac{\sum_{k=1}^{K}\alpha_{k}^{0}(p_{k}^{0}+\pi)}{C_{j}}\quad\mbox{for
}j=1,\dots,K.$ (5.22)

The equilibrium volume fractions can be calculated employing mass conservation
(5.14a).

Special case 2. Let us assume that $K=2$ with $\pi_{1}=0$ and $\pi_{2}\neq 0$.
Then it holds

$\displaystyle
X_{1}=\frac{-\pi_{2}+\sum_{k=1}^{2}\alpha_{k}^{0}\rho_{k}^{0}C_{k}}{C_{1}},$
(5.23a) $\displaystyle
X_{2}=-\frac{\pi_{2}}{C_{1}}(\alpha_{1}^{0}\rho_{1}^{0}),$ (5.23b)
$\displaystyle\rho_{1,\pm}^{\infty}=\frac{X_{1}}{2}\pm\sqrt{\frac{X_{1}^{2}}{4}-X_{2}}.$
(5.23c)

Because $-\pi_{2}<0$ one has to choose the upper sign, i.e.,
$\rho_{1}^{\infty}=\rho_{1,+}^{\infty}$. Using the saturation condition (5.17)
and mass conservation (5.14a) one can find the equilibrium volume fractions
and finally $\rho_{2}^{\infty}$.

Special case 3. Let us assume that $K=3$ with $\pi_{1}=\pi_{3}=0$ and
$\pi_{2}\neq 0$. Then it holds

$\displaystyle
X_{1}=\frac{-\pi_{2}+\sum_{k=1}^{3}\alpha_{k}^{0}\rho_{k}^{0}C_{k}}{C_{1}},$
(5.24a) $\displaystyle
X_{2}=-\frac{\pi_{2}}{C_{1}^{2}}(\alpha_{1}^{0}\rho_{1}^{0}C_{1}+\alpha_{3}^{0}\rho_{3}^{0}C_{3}),$
(5.24b)
$\displaystyle\rho_{1,\pm}^{\infty}=\frac{X_{1}}{2}\pm\sqrt{\frac{X_{1}^{2}}{4}-X_{2}}.$
(5.24c)

Because $-\pi_{2}<0$ one has to choose the upper sign. Analogously one can
find $\rho_{3}^{*}$. Again, using the saturation condition (5.17) and mass
conservation (5.14a) one can find the equilibrium volume fractions and finally
$\rho_{2}^{\infty}$.

Finally, after having determined the equilibrium phasic masses the equilibrium
pressure can be computed by

$\displaystyle
p^{\infty}=\sum_{k=1}^{K}\alpha_{k}^{\infty}p_{k}^{\infty}=\sum_{k=1}^{K}\alpha_{k}^{\infty}\left(C_{k}\rho_{k}^{\infty}-\pi_{k}\right)=\sum_{k=1}^{K}\left(C_{k}\alpha_{k}^{0}\rho_{k}^{0}-\alpha_{k}^{\infty}\pi_{k}\right).$
(5.25)

For convenience of the reader we summarize the relaxation procedure in the
following algorithm.

######  Algorithm 5.2 (Pressure relaxation procedure in case of isothermal
stiffened gases).

  1. 1.

Determine the index $k_{0}$ such that the cardinality of the set ${\cal
N}_{k_{0}}$ determined by (5.21) is smallest.

  2. 2.

Determine the roots $\rho_{k_{0}}^{\infty,i}$, $i=1,\ldots,\\#{\cal
N}_{k_{0}}+1$, of the polynomial (5.20);

  3. 3.

Single out the unique physically admissible solution by verifying the
conditions $\alpha_{k}^{\infty,i}\in(0,1)$ and $\rho_{k}^{\infty,i}>0$ for all
$k=1,\ldots,K$ where the volume fractions and phasic densities are determined
by (5.19) and (5.14a);

  4. 4.

By means of this solution update the phasic densities $\rho_{k}^{\infty}$ and
volume fractions $\alpha_{k}^{\infty}$ using (5.19) and (5.14a), respectively.

  5. 5.

Compute the equilibrium pressure using (5.25);

For the special cases discussed above we can prove the existence and
uniqueness of admissible equilibrium states, see Step 4 in the above
algorithm.

######  Theorem 5.1 (Properties of equilibrium state).

For all of the three cases (1) $K$ components with $\pi_{1}=\ldots=\pi_{K}$
(2) two components with $\pi_{2}\neq\pi_{1}=0$, (3) three components with
$\pi_{1}=\pi_{2}\neq\pi_{3}$ there exist unique equilibrium pressure,
equilibrium volume fractions and equilibrium phasic densities with
$p^{\infty}>-\pi$, $\alpha_{k}^{\infty}\in[0,1]$ and $\rho_{k}^{\infty}>0$
with minimal pressure $\pi:=\max\\{\pi_{k}\,;\,k=1,\ldots,K\\}$ provided that
$p_{k}^{0}>-\pi_{k}$, $k=1,\ldots,K$.

Proof:
Case 1: In this case there only exists one root $\rho_{k_{0}}^{\infty}$ of
(5.20) for arbitrary but fixed $k_{0}\in\\{1,\ldots,L\\}$. By means of (5.19)
the remaining equilibrium phasic densities are determined by (5.22). All of
them are positive. Furthermore, the equilibrium volume fractions are
determined by the mass conservation (5.14a) as

$\displaystyle\alpha_{j}^{\infty}=\frac{\alpha_{j}^{0}\rho_{j}^{0}}{\rho_{j}^{\infty}}=\frac{\alpha_{j}^{0}\rho_{j}^{0}}{\frac{1}{C_{j}}\sum_{k=1}^{K}\alpha_{k}^{0}\rho_{k}^{0}C_{k}}\in[0,1].$

Case 2: In this case $\pi_{2}\neq\pi_{1}=0$. Since the discriminant in (5.23c)
is positive, there exist two roots of the polynomial (5.20) from which we
choose $\rho_{1}^{\infty}=\rho_{1,+}^{\infty}$. By (5.16) we can rewrite the
expression (5.23a) as

$\displaystyle
X_{1}=\frac{\alpha_{1}^{0}\rho_{1}^{0}C_{1}+\alpha_{2}^{0}\,(p_{2}^{0}+\pi_{2})}{C_{1}}>0$

that is positive if $p_{2}$ is larger than the minimal pressure $-\pi_{2}$.
This implies $\rho_{1}^{\infty}\geq X_{1}>0$ a and, thus, because of (5.19)

$\displaystyle\rho_{2}^{\infty}=\frac{C_{1}}{C_{2}}\rho_{1}^{\infty}+\frac{\pi_{2}}{C_{2}}>0.$

For the equilibrium volume fractions we conclude using mass conservation
(5.14a) and saturation (5.17)

$\displaystyle
0\leq\alpha_{1}^{\infty}=\frac{\alpha_{1}^{0}\rho_{1}^{0}}{\rho_{1}^{\infty}}<\frac{\alpha_{1}^{0}\rho_{1}^{0}C_{1}}{\alpha_{1}^{0}\rho_{1}^{0}C_{1}+\alpha_{2}^{0}\,(p_{2}^{0}+\pi_{2})}\leq
1,$ $\displaystyle 0\leq\alpha_{2}^{\infty}=1-\alpha_{1}^{\infty}\leq 1.$

Case 3: In this case we assume that $\pi_{1}=\pi_{2}\neq\pi_{3}$. Since the
discriminant in (5.24c) is positive, there exist two roots of the polynomial
(5.20) from which we choose $\rho_{1}^{\infty}=\rho_{1,+}^{\infty}$. By (5.16)
we can rewrite the expression (5.24a) as

$\displaystyle
X_{1}=\frac{\alpha_{1}^{0}\rho_{1}^{0}C_{1}+\alpha_{3}^{0}\rho_{3}^{0}C_{3}+\alpha_{2}^{0}\,(p_{2}^{0}+\pi_{2})}{C_{1}}>0$

that is positive if $p_{2}$ is larger than the minimal pressure $-\pi_{2}$.
This implies $\rho_{1}^{\infty}\geq X_{1}>0$ and, thus, because of (5.19)

$\displaystyle\rho_{2}^{\infty}=\frac{C_{1}}{C_{2}}\rho_{1}^{\infty}+\frac{\pi_{2}}{C_{2}}>0,\quad\rho_{3}^{\infty}=\frac{C_{1}}{C_{3}}\rho_{1}^{\infty}>0.$

For the equilibrium volume fractions we conclude using mass conservation
(5.14a) and saturation (5.17)

$\displaystyle
0\leq\alpha_{1}^{\infty}=\frac{\alpha_{1}^{0}\rho_{1}^{0}}{\rho_{1}^{\infty}}<\frac{\alpha_{1}^{0}\rho_{1}^{0}C_{1}}{\alpha_{1}^{0}\rho_{1}^{0}C_{1}+\alpha_{3}^{0}\rho_{3}^{0}C_{3}+\alpha_{2}^{0}\,(p_{2}^{0}+\pi_{2})}\leq
1,$ $\displaystyle
0\leq\alpha_{3}^{\infty}=\frac{\alpha_{3}^{0}\rho_{3}^{0}}{\rho_{3}^{\infty}}=\frac{\alpha_{3}^{0}\rho_{3}^{0}C_{3}}{C_{1}\rho_{1}^{\infty}}<\frac{\alpha_{3}^{0}\rho_{3}^{0}C_{3}}{\alpha_{1}^{0}\rho_{1}^{0}C_{1}+\alpha_{3}^{0}\rho_{3}^{0}C_{3}+\alpha_{2}^{0}\,(p_{2}^{0}+\pi_{2})}\leq
1,$ $\displaystyle
0\leq\alpha_{2}^{\infty}=1-\alpha_{1}^{\infty}-\alpha_{3}^{\infty}\leq 1.$

Finally, for all cases we conclude for the equilibrium pressure (5.25) with
(5.16) and mass conservation (5.14a)

$\displaystyle
p^{\infty}=\sum_{k=1}^{K}\left(\alpha_{k}^{0}p_{k}^{0}-(\alpha_{k}^{\infty}-\alpha_{k}^{0})\pi_{k}\right).$

Assuming that the phasic pressures are larger than the minimal pressures,
i.e., $p_{k}^{0}>-\pi_{k}$, then we can estimate the equilibrium pressure by
$-\pi$ using the saturation (5.17).

####  5.1.2 Relaxation of chemical potentials

For two components the equilibrium state of the chemical potential relaxation
process is determined by solving the system of ODEs

$\displaystyle\frac{d\,\alpha_{k}}{d{\overline{t}}}=(-1)^{k+1}\frac{\dot{m}}{\varrho},\
k=1,2,$ (5.26)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}}{d{\overline{t}}}=(-1)^{k+1}\dot{m},\
k=1,2,$ (5.27)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k}}{d{\overline{t}}}=(-1)^{k+1}\dot{m}\,{\mbox{\boldmath$V$}}_{I},\
k=1,2.$ (5.28)

resulting from (5.2) with the source terms (2.48). Again we assume that the
relaxation process is infinitely fast such that the solution of the system of
ODEs converges towards the equilibrium state, where the right-hand side
vanishes. According to [62, 64] this holds true for

$\displaystyle p_{1}^{\infty}=p_{2}^{\infty}$ $\displaystyle=$ $\displaystyle
p^{\infty},\quad{\overline{t}}\to\infty,$ (5.29) $\displaystyle
g_{1}^{\infty}$ $\displaystyle=$ $\displaystyle
g_{2}^{\infty},\quad{\overline{t}}\to\infty.$ (5.30)

This implies $\dot{m}=\dot{m}(g_{2}-g_{1})=0$ and $p^{\infty}=p(T_{ref})$ as
well as $\rho_{k}^{\infty}=\rho_{k,ref}$. Here $p(T_{ref})$ denotes the
pressure at saturation temperature and $\rho_{k,ref}$ the corresponding
densities. Since we consider isothermal fluids, the temperature uniquely
defines the equilibrium state on the saturation line where both phases
coexist. Note that an equilibrium state does not necessarily exist. This can
easily be checked as follows.

The saturation condition and conservation of mass imply

$M:=\alpha_{1}^{0}\rho_{1}^{0}+(1-\alpha_{1}^{0})\rho_{2}^{0}=\alpha_{1}^{\infty}\rho_{1}^{\infty}+(1-\alpha_{1}^{\infty})\rho_{2}^{\infty}\,.$
(5.31)

Then one of the three cases applies: (i) if (5.31) has a solution
$\alpha_{1}^{\infty}\in(0,1)$, then an equilibrium state exists; (ii) for non-
positive solutions the vapor phase completely vanishes; (iii) otherwise the
liquid phase will be extinct. In the latter two cases we set the volume
fraction of the vanishing phase to a small positive value $\varepsilon$ and
determine the densities using pressure equilibrium
$p_{1}^{\infty}=p_{2}^{\infty}$.

The full relaxation procedure is given in the following algorithm.

######  Algorithm 5.3 (Chemical Relaxation Procedure for $K=2$).

  1. 1.

Solve equation (5.31) to identify (i) equilibrium, (ii) total condensation or
(iii) total evaporation case

  2. 2.

Determine the missing quantities

     * •

if (i): $\alpha_{1}^{\infty}$ is the solution of (5.31),
$p^{\infty}=p(T_{ref})$, $\rho_{k}^{\infty}=\rho_{k,ref}$, $k=1,2$

     * •

if (ii): $\alpha_{1}^{\infty}:=\varepsilon$,
$\alpha_{2}^{\infty}=1-\alpha_{1}^{\infty}$,
$\rho_{1}^{\infty}=\frac{\tilde{c}_{2}^{2}M-\alpha_{2}^{\infty}(\pi_{2}-\pi_{1})}{\alpha_{2}^{\infty}\tilde{c}_{1}^{2}+\alpha_{1}^{\infty}\tilde{c}_{2}^{2}}$,
$\rho_{2}^{\infty}=\frac{M-\alpha_{1}^{\infty}\rho_{1}^{\infty}}{\alpha_{2}^{\infty}}$

     * •

if (iii): $\alpha_{2}^{\infty}:=\varepsilon$,
$\alpha_{1}^{\infty}=1-\alpha_{2}^{\infty}$,
$\rho_{1}^{\infty}=\frac{\tilde{c}_{2}^{2}M-\alpha_{2}^{\infty}(\pi_{2}-\pi_{1})}{\alpha_{2}^{\infty}\tilde{c}_{1}^{2}+\alpha_{1}^{\infty}\tilde{c}_{2}^{2}}$,
$\rho_{2}^{\infty}=\frac{M-\alpha_{1}^{\infty}\rho_{1}^{\infty}}{\alpha_{2}^{\infty}}$

Note that in case (ii) and (iii) we determine $\rho_{2}^{\infty}$ by the mass
conservation (5.31). From this we deduce together with (2.30) for an
isothermal stiffened gas the equilibrium density $\rho_{1}^{\infty}$. Since
here $k=1$ represents water vapor modeled by an ideal gas, it holds
$\pi_{1}=0$.

A direct consequence of the algorithm is the following theorem.

######  Theorem 5.2.

The chemical relaxation procedure given by Alg. 5.3 provides a unique
solution. In particular, the resulting equilibrium pressure is positive.

The equilibrium state of the chemical potential relaxation process for three
components is determined by solving the system of ODEs

$\displaystyle\frac{d\,\alpha_{k}}{d{\overline{t}}}=\frac{\dot{m}}{\varrho_{k}},\
k=1,2,\quad\frac{d\,\alpha_{3}}{d{\overline{t}}}=-\dot{m}\,\left(\frac{1}{\varrho_{1}}+\frac{1}{\varrho_{2}}\right),$
(5.32)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}}{d{\overline{t}}}=(-1)^{k+1}\dot{m},\
k=1,2,\quad\frac{d\,\alpha_{3}\,\rho_{3}}{d{\overline{t}}}=0,$ (5.33)
$\displaystyle\frac{d\,\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k}}{d{\overline{t}}}=(-1)^{k+1}\dot{m}\,{\mbox{\boldmath$V$}}_{I},\
k=1,2,\quad\frac{d\,\alpha_{3}\,\rho_{3}\,{\mbox{\boldmath$v$}}_{3}}{d{\overline{t}}}={\mbox{\boldmath$0$}}.$
(5.34)

resulting from (5.2) with the source terms (2.50). According to [62, 64] at
equilibrium it holds

$\displaystyle p_{1}^{\infty}=p_{2}^{\infty}=p_{3}^{\infty}$ $\displaystyle=$
$\displaystyle p^{\infty},\quad{\overline{t}}\to\infty,$ (5.35)
$\displaystyle\mu_{1}^{\infty}$ $\displaystyle=$ $\displaystyle
g_{2}^{\infty},\quad{\overline{t}}\to\infty.$ (5.36)

Note that due to these equilibrium conditions the mass flux $\dot{m}$ vanishes
at equilibrium and, thus, the right-hand sides in (5.32), (5.33) and (5.34)
become zero. Obviously, the equations (5.32) are in agreement with the
saturation condition (2.1), i.e.,

$\alpha_{1}^{\infty}+\alpha_{2}^{\infty}+\alpha_{3}^{\infty}=1.$ (5.37)

In the following we start with relaxed values for velocity and pressure
obtained by the previous mechanical relaxation procedures and, thus, it holds

$p_{1}^{0}=p_{2}^{0}=p_{3}^{0},\quad{\mbox{\boldmath$v$}}_{1}^{0}={\mbox{\boldmath$v$}}_{2}^{0}={\mbox{\boldmath$v$}}_{3}^{0}.$
(5.38)

During the chemical relaxation procedure we now assume that mechanical
equilibrium is maintained but the corresponding equilibrium values may change.
Due to conservation of total mass and total momentum it is obvious that the
equilibrium velocity does not change, i.e.,

${\mbox{\boldmath$v$}}^{\infty}={\mbox{\boldmath$v$}}_{k}^{\infty}={\mbox{\boldmath$v$}}_{k}^{0}={\mbox{\boldmath$v$}}^{0},\quad
k=1,2,3.$ (5.39)

The mass equations (5.33) induce

$\displaystyle\alpha_{3}^{\infty}\rho_{3}^{\infty}$ $\displaystyle=$
$\displaystyle\alpha_{3}^{0}\rho_{3}^{0},$ (5.40)
$\displaystyle\alpha_{1}^{\infty}\rho_{1}^{\infty}+\alpha_{2}^{\infty}\rho_{2}^{\infty}$
$\displaystyle=$
$\displaystyle\alpha_{1}^{0}\rho_{1}^{0}+\alpha_{2}^{0}\rho_{2}^{0}=:M\,.$
(5.41)

Motivated by the equilibrium condition (5.36) the aim is now to derive a
function

$f_{\mu}(\alpha_{1}\rho_{1}):=\mu_{1}^{\infty}(\alpha_{1}\rho_{1})-g_{2}^{\infty}(\alpha_{1}\rho_{1})$
(5.42)

depending only on the product $\alpha_{1}\rho_{1}$ such that the root
$\alpha_{1}^{\infty}\rho_{1}^{\infty}$ is the solution for the relaxed mass
density for the first component. For this purpose we have to express
$\alpha_{1}^{\infty}$, $\alpha_{2}^{\infty}$ and $p^{\infty}$ in terms of
$\alpha_{1}^{\infty}\rho_{1}^{\infty}$, see below. Then the chemical potential
for water vapor $\mu_{1}$ and the Gibbs free energies $g_{k}$, $k=1,2$, can be
written as

$\displaystyle\mu_{1}^{\infty}(\alpha_{1}\rho_{1})$ $\displaystyle=$
$\displaystyle
g_{1}^{\infty}(\alpha_{1}\rho_{1})+c_{1}^{2}\ln\left(\frac{\alpha_{1}(\alpha_{1}\rho_{1})}{1-\alpha_{2}(\alpha_{1}\rho_{1})}\right),$
(5.43) $\displaystyle g_{k}^{\infty}(\alpha_{1}\rho_{1})$ $\displaystyle=$
$\displaystyle
c_{k}^{2}\ln\left(\frac{p^{\infty}(\alpha_{1}\rho_{1})+\pi_{k}}{p_{k,ref}+\pi_{k}}\right).$
(5.44)

Here $p_{k,ref}$ denotes the saturation pressure to a given reference
temperature. For the particular configuration at hand it holds
$p_{1,ref}=p_{2,ref}=p(T_{ref})$. If an inert component is present a root of
(5.42) always exists. Nevertheless, the equilibrium state can be very close to
total condensation or total evaporation. To avoid numerical instabilities a
small amount of the particular phase should be kept. In a first step one has
to figure out, which of the following four cases occurs:

  * (i-a)

condensation process with equilibrium solution,

  * (i-b)

nearly total condensation,

  * (ii-a)

evaporation process with equilibrium solution,

  * (ii-b)

nearly total evaporation.

In case (i-b) and (ii-b) the result can be directly obtained. In the cases
(i-a) and (ii-a) a bisection method is provided to find the equilibrium state.
Due to the fact that the pressure relaxation method is simple and always gives
a unique, physical solution, we base the bisection method to find the
equilibrium state on the pressure relaxation procedure.

According to thermodynamics a condensation process (i) and an evaporation
process (ii) are characterized by a positive or negative sign of
$\mu_{1}-g_{2}$, respectively. Thus, using the data from the pressure
relaxation procedure we may identify condensation and evaporation processes.

Condensation. If a condensation process is identified by
$f_{\mu}(\alpha_{1}^{0}\rho_{1}^{0})=\mu_{1}(\alpha_{1}^{0}\rho_{1}^{0})-g_{2}(\alpha_{1}^{0}\rho_{1}^{0})>0$
the expression $\alpha_{1}^{0}\rho_{1}^{0}$ is too large and
$\alpha_{1}\rho_{1}$ has to decrease. The smallest admissible value for this
expression is $\alpha_{1}^{*}\rho_{1}^{*}=tol>0.$ Using the pressure
relaxation procedure with $\alpha_{1}^{*}\rho_{1}^{*}$ instead of
$\alpha_{1}^{0}\rho_{1}^{0}$ we determine the corresponding values for all
variables of the phases. Using these data one has to check the sign of the
difference of the chemical potentials. If still $\mu_{1}-g_{2}\geq 0$ holds,
then nearly total condensation will occur. We keep
$\alpha_{1}^{*}\rho_{1}^{*}=tol>0$ and the corresponding data. Otherwise the
interval $[\alpha_{1}^{*}\rho_{1}^{*}=tol,\alpha_{1}^{0}\rho_{1}^{0}]$ is
admissible for the bisection method.

Evaporation. If an evaporation process is identified by
$f_{\mu}(\alpha_{1}^{0}\rho_{1}^{0})=\mu_{1}(\alpha_{1}^{0}\rho_{1}^{0})-g_{2}(\alpha_{1}^{0}\rho_{1}^{0})<0$
the expression $\alpha_{1}^{0}\rho_{1}^{0}$ is too small and
$\alpha_{1}\rho_{1}$ has to increase. The largest admissible value for this
expression is $\alpha_{1}^{*}\rho_{1}^{*}=M-tol.$ Analogously to the previous
case we determine the corresponding values for all variables of the phases and
check the sign of the difference of the chemical potentials. Either we
identify nearly total evaporation and keep the data or we find the admissible
interval $[\alpha_{1}^{0}\rho_{1}^{0},M-tol]$ for the bisection method.

We summarize the relaxation procedure in the following algorithm.

######  Algorithm 5.4 (Chemical Relaxation Procedure for $K=3$).

  1. 1.

Identify condensation or evaporation processes (i) or (ii) by the sign of
$f_{\mu}(\alpha_{1}^{0}\rho_{1}^{0})=\mu_{1}(\alpha_{1}^{0}\rho_{1}^{0})-g_{2}(\alpha_{1}^{0}\rho_{1}^{0})$
and determine the initial interval for the bisection method

  2. 2.

     * •

if (i-a): Apply the bisection method to find the root
$\overline{\alpha_{1}\rho_{1}}$ of (5.42) in the interval
$[\alpha_{1}^{*}\rho_{1}^{*}=tol,\alpha_{1}^{0}\rho_{1}^{0}]$ and perform
pressure relaxation with $\overline{\alpha_{1}\rho_{1}}$;

     * •

if (i-b): Apply the pressure relaxation procedure to
$\alpha_{1}^{*}\rho_{1}^{*}=tol$;

     * •

if (ii-a): Apply the bisection method to find the root
$\overline{\alpha_{1}\rho_{1}}$ of (5.42) in the interval
$[\alpha_{1}^{0}\rho_{1}^{0},\alpha_{1}^{*}\rho_{1}^{*}]$ with
$\alpha_{1}^{*}\rho_{1}^{*}=M-tol$ and perform pressure relaxation with
$\overline{\alpha_{1}\rho_{1}}$;

     * •

if (ii-b): Apply the pressure relaxation procedure to
$\alpha_{1}^{*}\rho_{1}^{*}=M-tol$.

A direct consequence of the algorithm is the following theorem.

######  Theorem 5.3.

The chemical relaxation procedure given by Alg. 5.4 provides a unique
solution. In particular, the resulting equilibrium pressure is positive.

The above relaxation procedures Algs. 5.3 and 5.4 can be considered as an
essential improvement of the original ones presented in [62, 64] due to the
following aspects: In the original procedure the Gibbs free energies are
relaxed, i.e., the influence of the mixture entropy, which is described by the
difference of the chemical potential and the Gibbs free energy of the vapor
phase, is neglected. This extra term cancels in pure phases. Moreover, in the
original procedure a nested iteration method was used. We now reduce the
numerical costs significantly by simplifying the equilibrium system to a
_scalar_ equation depending only on $\alpha_{1}\rho_{1}$ that has to be solved
by a _single_ iteration process. Moreover, for $K=2$ the solution can be
determined directly. Thus, the computational costs are significantly reduced.

Finally we point out that in complete analogy the above procedure can also be
applied to mixtures with an arbitrary number of additional inert ideal gases.

###  5.2 Isentropic case

Similar to isothermal fluids we have to solve the system of ODEs (5.12) for
the pressure relaxation process implying conservation of mass and momentum
(5.14a).

For the full multi-component Baer-Nunziato model pressure relaxation has been
investigated in [35] for two- and multi-component models. Lallemand et al.
[35] give eight different procedures based on a closure of the interfacial
pressure that become physically non-admissible for more than two components,
see [41]. These procedures predict different values for the equilibrium
pressure.

Pressure-temperature relaxation. Alternatively, uniqueness can be enforced by
the additional condition of temperature equilibrium. For this purpose it was
suggested in [17] to simultaneously relax pressure and temperature and, thus,
determine a unique equilibrium pressure by avoiding iteration processes. This
results in an efficient and stable procedure. However, for isentropic fluids
this approach does not work as can be verified for ideal gases: in the
equilibrium state it must hold

$\displaystyle
p_{k}^{\infty}=p^{\infty},\,T_{k}^{\infty}=T^{\infty},\,k=1,\dots,K.$ (5.45)

Assume $K=2$ and both components to be ideal gases, see Example 2.3 with
$\pi_{k}=0$. Then from the equations of state we have the relations

$\displaystyle
p_{k}=A_{k}\rho_{k}^{\gamma_{k}},\,T_{k}=B_{k}\rho_{k}^{\gamma_{k}-1}$ (5.46)

with constants $A_{k}$, $B_{k}$ determined by (2.38). From (5.45) and (5.46)
we obtain

$\displaystyle(\rho_{1}^{\infty})^{\gamma_{1}-\gamma_{2}}=\left(\frac{A_{2}}{A_{1}}\right)^{\gamma_{2}-1}\left(\frac{B_{1}}{B_{2}}\right)^{\gamma_{2}}\,.$
(5.47)

Obviously, for $\gamma_{1}\neq\gamma_{2}$ the parameters of the equations of
state define a unique equilibrium density, that does not depend on the initial
state. Analogously one can find $\rho_{2}^{\infty}$. Using (5.14a) one can
obtain the corresponding volume fractions $\alpha_{1}^{\infty}$ and
$\alpha_{2}^{\infty}$, that usually will not satisfy the saturation condition

$\displaystyle\alpha_{1}^{\infty}+\alpha_{2}^{\infty}=1\,.$ (5.48)

Pressure relaxation. Therefore for isentropic fluids we relax the pressure
without simultaneously relaxing the temperature, i.e., in the equilibrium
state it must hold

$\displaystyle p_{k}^{\infty}=p^{\infty},\,k=1,\dots,K.$ (5.49)

Assume all $K$ components to be stiffened gases. Then from the equations of
state we have the relations

$\displaystyle p_{k}+\pi_{k}=A_{k}\rho_{k}^{\gamma_{k}}$ (5.50)

with constant $A_{k}$ determined by (2.38a). Then conservation of mass (5.14a)
and the equilibrium conditions (5.49) lead to

$\displaystyle
A_{1}\left(\frac{\alpha_{1}^{0}\rho_{1}^{0}}{\alpha_{1}^{\infty}}\right)^{\gamma_{1}}-\pi_{1}=A_{k}\left(\frac{\alpha_{k}^{0}\rho_{k}^{0}}{\alpha_{k}^{\infty}}\right)^{\gamma_{k}}-\pi_{k},\,k=2,\dots,K\,.$
(5.51)

Solving these equations for $\alpha_{k}^{\infty}$ gives

$\displaystyle\alpha_{k}^{\infty}=\left[\frac{A_{k}(\alpha_{k}^{0}\rho_{k}^{0})^{\gamma_{k}}(\alpha_{1}^{\infty})^{\gamma_{1}}}{A_{1}(\alpha_{1}^{0}\rho_{1}^{0})^{\gamma_{1}}+(\pi_{k}-\pi_{1})(\alpha_{1}^{\infty})^{\gamma_{1}}}\right]^{\frac{1}{\gamma_{k}}}$
(5.52)

depending on the only unknown $\alpha_{1}^{\infty}$. The numbering of the
components should be chosen such that

$\min_{k}\\{\pi_{k},k=1,\dots,K\\}=\pi_{1}.$

Using (5.52) and the saturation condition we obtain

$\displaystyle
F(\alpha_{1}^{\infty})=\alpha_{1}^{\infty}+\sum_{k=2}^{K}\left[\frac{A_{k}(\alpha_{k}^{0}\rho_{k}^{0})^{\gamma_{k}}(\alpha_{1}^{\infty})^{\gamma_{1}}}{A_{1}(\alpha_{1}^{0}\rho_{1}^{0})^{\gamma_{1}}+(\pi_{k}-\pi_{1})(\alpha_{1}^{\infty})^{\gamma_{1}}}\right]^{\frac{1}{\gamma_{k}}}-1=0.$
(5.53)

It holds $F(0)<0$, $F(1)>0$ and the smooth function $F$ is strictly increasing
for $\alpha_{1}^{\infty}\in[0,1]$. Accordingly equation (5.53) has a unique
solution for $\alpha_{1}^{\infty}\in[0,1]$ which gives us a unique positive
equilibrium pressure $p^{\infty}$ satisfying

$\displaystyle\min_{k=1,\ldots,K}p_{0,k}\leq
p^{\infty}\leq\max_{k=1,\ldots,K}p_{0,k}.$

This is not true in general for the full Baer-Nunziato model as discussed
above.

The relaxation procedure is summarized in the following algorithm.

######  Algorithm 5.5 (Pressure relaxation procedure in case of isentropic
stiffened gases).

  1. 1.

Determine $\alpha_{1}^{\infty}\in[0,1]$ as the unique root of the function
(5.53) in $[0,1]$;

  2. 2.

Determine $\alpha_{k}^{\infty,i}$, $k=2,\ldots,K$, from (5.52);

  3. 3.

Compute the equilibrium pressure $p^{\infty}$ using (5.51) as

$\displaystyle
p^{\infty}=A_{1}\left(\frac{\alpha_{1}^{0}\rho_{1}^{0}}{\alpha_{1}^{\infty}}\right)^{\gamma_{1}}-\pi_{1}=A_{k}\left(\frac{\alpha_{k}^{0}\rho_{k}^{0}}{\alpha_{k}^{\infty}}\right)^{\gamma_{k}}-\pi_{k},\,k=2,\dots,K\,.$
(5.54)

  4. 4.

Update the phasic densities by $\rho_{k}^{\infty}=\rho_{k}(p^{\infty})$.

Since in general pressure and temperature cannot be in equilibrium at the same
time for isentropic fluids, we will not discuss temperature relaxation and
chemical relaxation because pressures relax faster than temperatures and
chemical potentials.

##  6 Numerical results

###  6.1 Discretization

For a condensed presentation of the discretization it is convenient to rewrite
the system (2.4a), (2.4b), (4.58) and (2.2) in matrix-vector representation

$\displaystyle{\partial_{t}{\mbox{\boldmath$u$}}}+{\nabla\cdot{\mbox{\boldmath$f$}}^{i}({\mbox{\boldmath$u$}})}+\sum_{i=1}^{d}{\mbox{\boldmath$G$}}^{i}({\mbox{\boldmath$u$}})\partial_{{x_{i}}}{\mbox{\boldmath$u$}}={\mbox{\boldmath$s$}}({\mbox{\boldmath$u$}})$
(6.1)

where we make use of the convention (2.3). This forms a coupled system for the
volume fractions
${\mbox{\boldmath$\alpha$}}:=(\alpha_{1},\ldots,\alpha_{K-1})^{T}$ and the
vectors
${\mbox{\boldmath$u$}}_{k}:=(\alpha_{k}\,\rho_{k},\alpha_{k}\,\rho_{k}\,{\mbox{\boldmath$v$}}_{k})^{T}$
of the conserved quantities of phase $k=1,\ldots,K$ condensed in the vector
${\mbox{\boldmath$u$}}:=({\mbox{\boldmath$\alpha$}}^{T},{\mbox{\boldmath$u$}}_{1}^{T},\ldots,{\mbox{\boldmath$u$}}_{K}^{T})^{T}$.
The flux ${\mbox{\boldmath$f$}}^{i}$ in the $i$th coordinate direction and the
source term $s$ are determined by

$\displaystyle{\mbox{\boldmath$f$}}^{i}({\mbox{\boldmath$u$}})^{T}:=\left({\mbox{\boldmath$0$}}_{K-1}^{T},{\mbox{\boldmath$f$}}^{i}_{1}(\alpha_{1},{\mbox{\boldmath$u$}}_{1})^{T},\ldots,{\mbox{\boldmath$f$}}^{i}_{K}(\alpha_{K},{\mbox{\boldmath$u$}}_{K})^{T}\right),$
(6.2)
$\displaystyle{\mbox{\boldmath$s$}}({\mbox{\boldmath$u$}}):=\left({\mbox{\boldmath$S$}}_{\alpha}({\mbox{\boldmath$u$}})^{T},{\mbox{\boldmath$S$}}_{1}({\mbox{\boldmath$u$}})^{T},\ldots,{\mbox{\boldmath$S$}}_{K}({\mbox{\boldmath$u$}})^{T}\right)$
(6.3)

with phasic fluxes and sources

$\displaystyle{\mbox{\boldmath$f$}}^{i}_{k}(\alpha_{k},{\mbox{\boldmath$u$}}_{k}):=(\alpha_{k}\rho_{k}v_{k,i},\alpha_{k}\rho_{k}{\mbox{\boldmath$v$}}_{1}^{T}v_{k,i}+\alpha_{k}p_{k}{\mbox{\boldmath$e$}}_{d,i}^{T},$
(6.4)
$\displaystyle{\mbox{\boldmath$S$}}_{k}({\mbox{\boldmath$u$}}):=(S_{\alpha\rho,k}({\mbox{\boldmath$u$}}),{\mbox{\boldmath$S$}}_{\alpha\rho\mbox{\scriptsize\boldmath${\mbox{\boldmath$v$}}$},k}^{T}({\mbox{\boldmath$u$}}))^{T},$
(6.5)
$\displaystyle{\mbox{\boldmath$S$}}_{\alpha}({\mbox{\boldmath$u$}}):=(S_{\alpha,1}({\mbox{\boldmath$u$}}),\ldots,S_{\alpha,K-1}({\mbox{\boldmath$u$}}))^{T}.$
(6.6)

Here ${\mbox{\boldmath$e$}}_{d,i}\in\mathbb{R}^{d}$ denotes the $ith$ unit
vector in $\mathbb{R}^{d}$ and ${\mbox{\boldmath$0$}}_{K-1}$ the zero vector
in $\mathbb{R}^{K-1}$. The nonlinear product in (6.1) is characterized by the
matrices

$\displaystyle{\mbox{\boldmath$G$}}^{i}({\mbox{\boldmath$u$}})=\left(\begin{matrix}V_{I,i}{\mbox{\boldmath$I$}}_{K-1}&{\mbox{\boldmath$0$}}_{(K-1)\times(K(d+1))}\\\
\left(\begin{matrix}{\mbox{\boldmath$0$}}_{K-1}^{T}\\\
{\mbox{\boldmath$e$}}_{d,i}\otimes{\mbox{\boldmath$d$}}_{1}\end{matrix}\right)&{\mbox{\boldmath$0$}}_{(d+1)\times(K(d+1))}\\\
\vdots&\vdots\\\ \left(\begin{matrix}{\mbox{\boldmath$0$}}_{K-1}^{T}\\\
{\mbox{\boldmath$e$}}_{d,i}\otimes{\mbox{\boldmath$d$}}_{K}\end{matrix}\right)&{\mbox{\boldmath$0$}}_{(d+1)\times(K(d+1))}\end{matrix}\right),$
(6.7)

where ${\mbox{\boldmath$I$}}_{K-1}$ denotes the identity matrix in
$\mathbb{R}^{(K-1)\times(K-1)}$ and the components of the vectors
${\mbox{\boldmath$d$}}_{k}$, $k=1,\ldots,K$ are defined by
$d_{k,l}=P_{k,l}-P_{k,K}$, $l=1,\ldots,K-1$. Note that $P_{k,k}$ by convention
ist set to 0.

Since we are interested in flows near to equilibrium, we do not discretize the
coupled system (6.1) but perform an operator splitting according to Godunov or
Strang, see [60], resulting in the system of equations of fluid motion

$\displaystyle{\partial_{t}{\mbox{\boldmath$u$}}}+{\nabla\cdot{\mbox{\boldmath$f$}}^{i}({\mbox{\boldmath$u$}})}+\sum_{i=1}^{d}{\mbox{\boldmath$G$}}^{i}({\mbox{\boldmath$u$}})\partial_{{x_{i}}}{\mbox{\boldmath$u$}}={\mbox{\boldmath$0$}}$
(6.8)

and the system of relaxation

$\displaystyle{\frac{d\,{\mbox{\boldmath$u$}}}{dt}}={\mbox{\boldmath$s$}}({\mbox{\boldmath$u$}}).$
(6.9)

This procedure is in analogy to [55]. Note that the differential equation
(6.8) describes a variation in time and space, whereas the system (6.9) is an
evolution in time only although the state $u$ depends on both time and space.
Therefore we use different symbols for the differentiation.

In each time step we thus perform alternately the evolution of the fluid and
the relaxation process:

For the fluid equations (6.8) we apply a path-conservative DG scheme. For this
purpose, we extended our in-house Runge-Kutta DG solver [31, 10] by a path-
conservative discretization of the nonlinear products in (6.8) following [45].
For details on the derivation of the spatial path-conservative DG scheme we
refer to [2].

The efficiency of the scheme is improved by local grid adaption where we
employ the multi-resolution concept based on multiwavelets. The key idea is to
perform a multi-resolution analysis on a sequence of nested grids providing a
decomposition of the data on a coarse scale and a sequence of details that
encode the difference of approximations on subsequent resolution levels. The
detail coefficients become small when the underlying data are locally smooth
and, hence, can be discarded when dropping below a threshold value
$\varepsilon_{thresh}$. By means of the thresholded sequence a new, locally
refined grid is determined. Details on this concept can be found in [31, 10,
11].

Since our interest is in equilibrium processes we replace the transient
relaxation procedure (6.9) by integration to infinity, i.e., we determine the
equilibrium states by algebraic equations according to Sect. 5.2 and Sect. 5.1
in case of an isentropic and isothermal fluid, respectively. The relaxation to
equilibrium is performed one-by-one. First we determine the equilibrium
velocity by the velocity relaxation according to Algorithm 5.1. Then we
perform pressure relaxation applying Algorithm 5.5 or 5.2. In addition, in
case of an isothermal fluid we determine the state at chemical equilibrium
according to Algorithm 5.3 and 5.4 for a two-component and three-component
mixture, respectively. In case of a higher order DG scheme the solution on
each element is given by a polynomial of higher degree. Since the equilibrium
states are determined locally, for each element we have to perform a
projection of the equilibrium states to the polynomial space.

For this purpose, we first determine for each element the equilibrium states
in all quadrature points. From these we compute an interpolation polynomial
that interpolates these equilibrium states. The polynomial degree of the
interpolation polynomial is chosen smaller than the degree of the DG
polynomial. Thus, projecting the interpolation polynomial onto the DG space
will preserve the equilibrium states in the quadrature points. Note that in
the DG scheme we only access the DG polynomials in the quadrature points.

We would like to remark that the relaxation process is applied whenever the
states are out of equilibrium, i.e., after applying the limiter, performing
grid adaptation or a Runge-Kutta stage in the time evolution.

###  6.2 Validation

For validation we consider two configurations. First, we investigate a two-
phase configuration where each component is complemented by the isothermal
stiffened gas model, see Example 2.2. The results are compared to a sharp-
interface isothermal Euler model, see [18, 21]. For another configuration we
consider a four-component mixture of isothermal ideal gases. The results are
compared to an isothermal phase-field model discussed in [19].

####  6.2.1 How to choose material parameters

In the literature one can find numerous numerical examples for two- and multi-
phase flows with phase transition at similar states with different sets of
material parameters. Often these parameters are chosen arbitrarily. This in
many cases leads to unphysical phase states, unphysical equilibrium states and
unphysical processes. This can be avoided when appropriate tables are used,
see for example [61]. For water vapor and liquid water we proceed as follows,
where water vapor is modeled as an ideal gas and liquid water is modeled as a
stiffened gas.

######  Algorithm 6.1 (Procedure to determine material parameters).

  * •

Vapor:

    * –

Choose some reference temperature $T_{ref}$.

    * –

Identify the corresponding saturation pressure $p_{ref}=p(T_{ref})$ from the
table.

    * –

Identify the corresponding specific heat capacity at constant volume
$c_{v,1}$, the specific volume $1/\rho_{ref,1}$ for the vapor phase and the
specific entropy $s_{0,1}$ for the vapor phase from the table.

    * –

$\pi_{1}=0$

    * –

Determine $\gamma_{1}$ using (2.30):
$\gamma_{1}=1+\frac{p_{ref}}{\rho_{ref,1}c_{v,1}T_{ref}}$.

  * •

Liquid:

    * –

Determine $q_{1}$ using (2.31): $q_{1}=(s_{0,1}-c_{v,1})T_{ref}$.

    * –

Identify the corresponding specific heat capacity at constant volume
$c_{v,2}$, the specific volume $1/\rho_{ref,2}$, the specific entropy
$s_{0,2}$ and the speed of sound $c_{2}$ for the liquid phase from the table.

    * –

Determine $\gamma_{2}$ using (2.23b) and (2.24):
$\gamma_{2}=\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{c_{2}^{2}}{c_{v,2}T_{ref}}}$.

    * –

Determine $\pi_{2}$ using (2.30):
$\pi_{2}=\rho_{ref,2}(\gamma_{2}-1)c_{v,2}T_{ref}-p_{ref}$.

    * –

Determine $q_{2}$ using (2.31): $q_{2}=(s_{0,2}-c_{v,2})T_{ref}$.

Choosing different values for $T_{ref}$ we find the following set of
parameters recorded in Table 1. Note that the above procedure applies for
stiffened gas in general because we start from a reference state. Thus, we may
apply it for both the isothermal and the isentropic case, see Examples 2.2 and
2.3. Note that for water vapor we do not need to specify the sound speed
because for an ideal gas we set $\pi=0$.

component |  $1/\rho_{ref}$[m3/kg] |  $c_{v}$ [J/kg/K] |  $c$ [m/s] |  $s_{0}$ [J/kg]
---|---|---|---|---
water vapor (1) | 2.35915 | 1530.4 | - | 7478.1
liquid water (2) | 0.00103594 | 3818.8 | 1552.7 | 1192.7
Table 1: Material parameters taken from [61] for $T_{ref}=363.15$ K with
$p_{ref}=70182.4$ Pa.

####  6.2.2 Numerical example: vapor–vapor compression problem

To investigate nucleation phenomena we consider the following 1d-Riemann
problem determined by a compression problem for an initially single component
vapor. The initial data are listed in Table 2 with material parameters given
in Table 1. Two sets of computations are performed where we either perform
mechanical relaxation, see Fig. 1, or mechanical and chemical relaxation, see
Fig. 2, i.e., the computations are performed with and without mass transfer.
In the figures we present the volume of water vapor (top, left), the pressure
of the mixture (top, right), the velocity of the mixture (bottom, left) and
the density of the mixture (bottom, right).

The computational domain is $\Omega=(-2,2)$. The simulations are run to time
$t=0.001$. We perform second order $(p=2)$ computations using 3 cells on the
base level and 10 levels of refinement. The CFL number is 0.1 resulting in
about 61,000 time steps to reach the final time. The constant $c_{thresh}$ for
the thresholding performed in the grid adaptation step is chosen as $10^{-5}$.
This small value was needed to ensure that all discontinuities are resolved on
the highest refinement level where limiting is performed. At the final time
the adaptive grid is about 82% of the uniformly refined grid.

In this configuration we consider “pure” water vapor that is slightly polluted
with liquid water because the model does not allow for pure states. Pressure
and velocity are chosen in mechanical equilibrium. The initial pressure is
slightly below the saturation pressure. At the center $x=0$ the fluid is
compressed due to the opposed velocities causing two compression/shock waves
running to the left and right of the center. This results in a significant
pressure increase near the center where the pressure exceeds the saturation
pressure. According to the pressure relaxation we observe a significant drop
in the water vapor volume that is caused by a rise in the water vapor density.

There is a significant difference between the computations without and with
mass transfer. When no chemical relaxation is performed no mass is transferred
from water vapor to liquid water. This is reflected in the results for the
density of the mixture. In case of no mass transfer the value is still close
to the initial value. However, when mass transfer occurs the density of the
mixture significantly increases in the vicinity of the center. Note that this
behavior is not fully resolved by the Baer-Nunziato model because it requires
a much finer grid to fully capture the condensation of the water vapor. Note
that the wave structure consists of two outer shock waves and two phase
boundaries in between, which enclose the nucleated liquid. However, due to the
very small volume of the created liquid this cell is not fully resolved and we
refer to [37] for further insights on the numerical difficulties.

To validate the above results we compare the results with computations
performed for an isothermal Euler model where two phases are separated by a
sharp interface. Phase transition is accounted for by the jump condition at
the interface. For details we refer to [18, 21]. Numerical computations of
this model are performed independently applying a finite volume scheme in both
of the two domains to either side of the interface. At the interface a two-
phase Riemann problem is solved to compute the fluxes where the jump
conditions account for mass transfer or neglect mass transfer.
We observe a perfect agreement of the two models in case of no mass transfer,
see Fig. 1. When accounting for mass transfer then wave speeds of the two
models perfectly match but there is an obvious severe discrepancy in the
values of the intermediate states, see the red lines in Fig. 2. By numerical
experiments we observe that in this case a scaling of the kinetic relation in
the sharp-interface model by a factor of $300$ results in a perfect agreement
of the states. Note that in the Baer-Nunziato model we relax to the
equilibrium state, i.e., the relaxation velocity is infinite, whereas in the
isothermal Euler model the relaxation velocity is finite. Thus, increasing the
relaxation velocity in the isothermal Euler model it comes closer to the
equilibrium state.

$\alpha_{1,L}$ |  $p_{L}$ [Pa] |  $v_{L}$ [m/s] | $\alpha_{1,R}$ |  $p_{R}$ [Pa] |  $v_{R}$ [m/s]
---|---|---|---|---|---
$1-10^{-8}$ | 69985.63 | 2.7 | $1-10^{-8}$ | 69985.63 | -2.7
Table 2: Initial data for Riemann problem: vapor–vapor compression problem.
Figure 1: Vapor-vapor compression: without phase transition. Figure 2: Vapor-
vapor compression: with phase transition.

####  6.2.3 Numerical example: two components — two phases

Another interesting test case for further validation of the model is given by
the flow of two components with two phases, i.e., $k=1,2$ (component 1,
liquid/gas) and $k=3,4$ (component 2, liquid/gas), in one dimension. The
initial data are listed in Table 3. All components are modeled by isothermal
ideal gases, i.e., $\pi_{1}=\pi_{2}=\pi_{3}=\pi_{4}=0$ in (2.30). The
isothermal speeds of sounds are chosen as $c_{1}=600$ m/s, $c_{2}=1000$ m/s,
$c_{3}=650$ m/s, $c_{4}=1100$ m/s. The heat capacities are $c_{v,1}=1100$ J/kg
K, $c_{v,2}=1200$ J/kg K, $c_{v,3}=2100$ J/kg K, $c_{v,4}=2400$ J/kg K. The
reference temperature is $T=293.15$ K. From these values the parameters
$\gamma_{k}$ can be computed by (2.32). Note that the liquid phases are also
modeled by ideal gases instead of stiffened gases. This allows to rewrite the
data of the Baer-Nunziato model into the data of the phase field model [19]
and the same material parameters can be used. For stiffened gases this is no
longer obvious because of the minimal pressure $\pi$.

Computations for this configuration are only performed using velocity and
pressure relaxation. The computational domain is $\Omega=(-2,2)$. The
simulations are run to time $t=0.001$. We perform second order $(p=2)$
computations using 3 cells on the base level and 10 levels of refinement. The
CFL number is 0.1 resulting in about 30,000 time steps to reach the final
time. The constant $c_{thresh}$ for the thresholding performed in the grid
adaptation step is chosen as $10^{-5}$. This small value was needed to ensure
that all discontinuities are resolved on the highest refinement level where
limiting is performed. At the final time the adaptive grid is about 94% of the
uniformly refined grid.

In this configuration the pressure and velocity are chosen in mechanical
equilibrium. At the center $x=0$ the fluid is compressed due to the opposed
velocities causing two compression/shock waves running to the left and right
of the center, see Fig. 3. These waves can be observed in the mixture
velocity, the mixture pressure and the mixture densities of component 1 and
component 2. The volume fractions only jump at the material interface near the
center. We note that the mixture velocities and the mixture pressures are
continuous across the material wave but the mixture densities suffer from a
jump. This indicates that the material wave is a contact wave.

For validation of the Baer-Nunziato model we compare the results with
computations using the phase-field model [19]. Here each component is
described by a single density while the phase-field parameter indicates the
present phase. Accordingly the density $\rho_{1}$ of the phase-field model is
related to $\alpha_{1}\rho_{1}+\alpha_{2}\rho_{2}$ of the Baer-Nunziato model
and $\rho_{2}$ of the phase-field model is related to
$\alpha_{3}\rho_{3}+\alpha_{4}\rho_{4}$ of the Baer-Nunziato model. Further we
compare the mixture velocity and the mixture pressure of both models. The
results of the two models show a very good agreement except at the jumps were
the phase-field model shows some smearing of the discontinuity.

$\alpha_{1,L}$ | $\alpha_{2,L}$ | $\alpha_{3,L}$ |  $p_{L}$ [Pa] |  $v_{L}$ [m/s] | $\alpha_{1,R}$ | $\alpha_{2,R}$ | $\alpha_{3,R}$ |  $p_{R}$ [Pa] |  $v_{R}$ [m/s]
---|---|---|---|---|---|---|---|---|---
$0.2$ | $0.3$ | $0.1$ | $100000$ | $60$ | $0.5$ | $0.1$ | $0.2$ | $110000$ | $-50$
Table 3: Initial data for Riemann problem: two components — two phases.
Figure 3: two components — two phases.

##  7 Conclusion

In this work we have investigated Baer-Nunziato-type models for barotropic
fluids with an arbitrary number of components. In the literature frequently
barotropic fluids are confined to isentropic fluids, cf. [25, 3, 54] whereas
here we also consider isothermal fluids.

First of all, we have determined eigenvalues and eigenvectors for the quasi-1D
model in case of general barotropic fluids, i.e., the phasic pressures only
depend on the phasic densities. For non-resonance states a full set of
linearly independent eigenvectors has been derived. The different fields have
been characterized as either genuinely nonlinear or linearly degenerated
assuming that the phasic fundamental derivatives of thermodynamics are
positive and choosing a particular convex combination for the interfacial
velocity. For all fields we have given a full set of Riemann invariants whose
gradients are linearly independent. Furthermore we have verified that the
quasi-1D equations is symmetrizable away from non-resonance states
constructing a symmetrizer.

In case of isothermal fluids and isentropic fluids we have constructed
appropriate entropy-entropy flux pairs and have verified that (i) the
(mathematical) entropy functions are convex but not strictly convex even if
the phasic entropy functions are strictly convex and (ii) the (mathematical)
entropy inequality holds, i.e., the sum of the entropy production due to the
non-conservative products and the relaxation processes is non-increasing. It
is worthwhile mentioning that for isothermal fluids and isentropic fluids we
have used different entropy functions, namely, the mixture of the phasic
entropies and the mixture of phasic total energies, respectively. In both
cases we have to account for the heat flux in the balance laws for the phasic
total energies to render these equations to be redundant in case of barotropic
fluids. From a physical point of view this is mandatory. In the isothermal
case heat production due to the work performed in the system has to be removed
from the system to keep the temperatures constant. On the other hand, for
isentropic fluids there is an energy production due to the relaxation
processes. For both cases corresponding heat fluxes have been identified.
Furthermore, since we cannot control the sign of the energy production due to
the non-conservative products we have determined the interfacial velocity and
the interfacial pressures such that this production term vanishes. Starting
from an interfacial velocity determined by a convex combination of phasic
velocities a unique set of interfacial pressures has been identified that
coincides for both isothermal fluids and isentropic fluids.

The relaxation processes to the equilibrium state has been investigated for
both the isothermal case and the isentropic case. Mechanical relaxation has
been considered in both cases whereas relaxation of the chemical potentials is
only meaningful for isothermal fluids. Algorithms for the computation of the
corresponding equilibrium states are given.

Finally, we perform numerical computations in one space dimension by means of
an adaptive path-conservative Runge-Kutta DG scheme. Computations for
isothermal fluids have been performed for the configuration of a vapor-vapor
expansion problem near pressure saturation. The results have been validated by
means of a sharp-interface model separating two pure phases modelled by the
isothermal Euler equations. The results showed a perfect agreement when
neglecting mass transfer. When accounting for phase change then the kinetic
relation describing the mass transfer in the sharp-interface model had to be
scaled to obtain a perfect agreement. For another configuration we considered
a mixture of two components with two phases. Computations were performed for
mechanical equilibrium but without relaxation of the chemical potentials. The
results were compared to a phase field model. Again, the results showed a very
good agreement between the models.

## References

  * [1] M. Baer and J. Nunziato.  A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials.  International Journal of Multiphase Flow, 12(6):861–889, 1986.
  * [2] B. Blaut.  A Path-conservative DG Scheme for Conservation Laws, 2021.  Master thesis (in German).
  * [3] Boukili, Hamza and Hérard, Jean-Marc.  Relaxation and simulation of a barotropic three-phase flow model.  ESAIM: M2AN, 53(3):1031–1059, 2019.
  * [4] F. Coquel, J.-M. Hérard, K. Saleh, and N. Seguin.  Two properties of two-velocity two-pressure models for two-phase flows.  Int. J. Multiphase Flows, 12(3):593–600, 2014.
  * [5] D. Drew.  Mathematical modeling of two-phase flow.  Ann. Rev. Fluid Mech., 15:261–291, 1983.
  * [6] D. Drew and S. Passman.  Theory of Multicomponent Fluids, volume 135 of Applied Mathematical Sciences.  Springer, 1999.
  * [7] W. Dreyer, F. Duderstadt, M. Hantke, and G. Warnecke.  Bubbles in liquids with phase transition. part 1: On phase change of a single vapor bubble in liquid water.  Contin. Mech. Thermodyn., 24:461–483, 2012.
  * [8] T. Flåtten and H. Lund.  Relaxation two-phase flow models and the subcharacteristic condition.  Mathematical Models and Methods in Applied Sciences, 21(12):2379–2407, 2011.
  * [9] T. Gallouët, J.-M. Hérard, and N. Seguin.  Numerical modelling of two-phase flows using the two-fluid two-pressure approach.  Mathematical Models and Methods in Applied Sciences, 14(5):663–700, 2004.
  * [10] N. Gerhard and S. Müller.  Adaptive multiresolution discontinuous Galerkin schemes for conservation laws: multi-dimensional case.  Computational and Applied Mathematics, 35(2):321–349, 2016.
  * [11] N. Gerhard, S. Müller, and A. Sikstel.  A wavelet-free approach for multiresolution-based grid adaptation for conservation laws.  Communications on Applied Mathematics and Computation, 4:108–142, 2021.
  * [12] P. Goatin and P. LeFloch.  The Riemann problem for a class of resonant hyperbolic systems of balance laws.  Annales Inst. Henri Poincare, 21(6):881–902, 2004.
  * [13] E. Godlewski and P.-A. Raviart.  Hyperbolic Systems of Conservation Laws, volume 3-4 of Mathématiques & Applications (Paris).  Springer, 1991.
  * [14] S. K. Godunov, T. Y. Mikhaîlova, and E. I. Romenskiî.  Systems of thermodynamically coordinated laws of conservation invariant under rotations.  Siberian Mathematical Journal, 37(4):690–705, jul 1996.
  * [15] S. K. Godunov and E. Romenski.  Elements of Continuum Mechanics and Conservation Laws.  Springer US, 2003.
  * [16] S. K. Godunov and E. I. Romenski.  Thermodynamics, conservation laws, and symmetric forms of differential equations in mechanics of continuous media.  In Computational Fluid Dynamics Review 95, pages 19–31. John Wiley, NY, 1995.
  * [17] E. Han, M. Hantke, and S. Müller.  Efficient and robust relaxation procedures for multi-component mixtures including phase transition.  Journal of Computational Physics, 338:217–239, 2017.
  * [18] M. Hantke, W. Dreyer, and G. Warnecke.  Exact solutions to the Riemann problem for compressible isothermal Euler equations for two phase flows with and without phase transition.  Quarterly of Applied Mathematics, 71(3):509 – 540, 2013.
  * [19] M. Hantke, C. Matern, G. Warnecke, and H. Yaghi.  The Riemann problem for a two-phase mixture hyperbolic system with phase function and multi-component equation of state.  Quarterly of Applied Mathematics, 82:451–466, June 2024.
  * [20] M. Hantke, S. Müller, and L. Grabowsky.  News on Baer-Nunziato-type model at pressure equilibrium.  Continuum Mechanics and Thermodynamics, 2020.
  * [21] M. Hantke and F. Thein.  A general existence result for isothermal two-phase flows with phase transition.  Journal of Hyperbolic Differential Equations, 16(04):595–637, 2019\.
  * [22] Hantke, M. and Müller, S.  Analysis and simulation of a new multi-component two-phase flow model with phase transitions and chemical reactions.  Quart. Appl. Math., 76(2):253–287, 2018.
  * [23] Hantke, M. and Müller, S.  Closure conditions for a one temperature non-equilibrium multi-component model of Baer-Nunziato type.  ESAIM: ProcS, 66:42–60, 2019.
  * [24] J.-M. Hérard.  A three-phase flow model.  Mathematical and Computer Modelling, 45(5):732–755, 2007.
  * [25] J.-M. Hérard.  A class of compressible multiphase flow models.  Comptes Rendus Mathematique, 354(9):954–959, 2016.
  * [26] J.-M. Hérard and O. Hurisse.  A fractional step method to compute gas-liquid flows.  Computers and Fluids, 55:57–69, 2012.
  * [27] Hérard, J.-M. and Jomée, G.  Pressure relaxation in some multiphase flow models.  ESAIM: ProcS, 72:19–40, 2023.
  * [28] Hérard, J.-M. and Jomée, G.  Relaxation process in a hybrid two-phase flow model.  2023\.
  * [29] Hérard, J.-M. and Mathis, H.  A three-phase flow model with two miscible phases.  ESAIM: M2AN, 53(4):1373–1389, 2019.
  * [30] Hérard, Jean-Marc, Hurisse, Olivier, and Quibel, Lucie.  A four-field three-phase flow model with both miscible and immiscible components.  ESAIM: M2AN, 55:S251–S278, 2021.
  * [31] N. Hovhannisyan, S. Müller, and R. Schäfer.  Adaptive multiresolution discontinuous galerkin schemes for conservation laws.  Mathematics of Computation, 83:113–151, 01 2014.
  * [32] I. Müller.  Thermodynamics.  Pitman, London, 1985.
  * [33] A. Kapila, R. Menikoff, J. Bdzil, S. Son, and D. Stewart.  Two-phase modelling of DDT in granular materials: Reduced equations.  Phys. Fluid, 13:3002–3024, 2001.
  * [34] T. Kato.  The cauchy problem for quasi-linear symmetric hyperbolic systems.  Arch. Rational Mech. Anal., 58(3):181–205, 1975.
  * [35] M.-H. Lallemand, A. Chinnayya, and O. Le Metayer.  Pressure relaxation procedures for multiphase compressible flows.  International Journal for Numerical Methods in Fluids, 49(1):1–56, 2005.
  * [36] Linga, G. and Flåtten, T.  A hierarchy of non-equilibrium two-phase flow models.  ESAIM: ProcS, 66:109–143, 2019.
  * [37] S. May and F. Thein.  Explicit implicit domain splitting for two phase flows with phase transition.  Physics of Fluids, 35(1):016108, 2023.
  * [38] I. Müller and W. Müller.  Fundamentals of Thermodynamics and Applications.  Springer-Verlag, Berlin, 2009.
  * [39] S. Müller.  Multiresolution schemes for conservation laws.  In R. DeVore and A. Kunoth, editors, Multiscale, nonlinear and adaptive approximation, pages 379–408. Berlin: Springer, 2009.  Dedicated to Wolfgang Dahmen on the occasion of his 60th birthday.
  * [40] S. Müller, M. Hantke, and P. Richter.  Closure conditions for non-equilibrium multi-component models.  IGPM Preprint 414, RWTH Aachen University, 2014.
  * [41] S. Müller, M. Hantke, and P. Richter.  Closure conditions for non-equilibrium multi-component models.  Continuum Mechanics and Thermodynamics, 28(4):1157–1189, Jul 2016\.
  * [42] I. Peshkov, M. Pavelka, E. Romenski, and M. Grmela.  Continuum mechanics and thermodynamics in the hamilton and the Godunov-type formulations.  Continuum Mechanics and Thermodynamics, 30(6):1343–1378, Nov 2018\.
  * [43] F. Petitpas, R. Saurel, E. Franquet, and A. Chinnayya.  Modeling detonation waves in condensed materials: multiphase cj conditions and multidimensional computations.  Shock Waves, 19:377–401, 2009.
  * [44] B. Re and R. Abgrall.  Non-equilibrium model for weakly compressible multi-component flows: The hyperbolic operator.  In Lecture Notes in Mechanical Engineering, pages 33–45. Springer International Publishing, 2020.
  * [45] S. Rhebergen, O. Bokhove, and J. van der Vegt.  Discontinuous galerkin finite element methods for hyperbolic nonconservative partial differential equations.  Journal of Computational Physics, 227(3):1887–1922, 2008.
  * [46] L. Río-Martín and M. Dumbser.  High-order ader discontinuous galerkin schemes for a symmetric hyperbolic model of compressible barotropic two-fluid flows.  Communications on Applied Mathematics and Computation, Nov 2023\.
  * [47] M. G. Rodio and R. Abgrall.  An innovative phase transition modeling for reproducing cavitation through a five-equation model and theoretical generalization to six and seven-equation models.  International Journal of Heat and Mass Transfer, pages –, 05 2015\.
  * [48] E. Romenski.  Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechanics.  Math. Comput. Modell., 28(10):115–130, 1998.
  * [49] E. Romenski.  Thermodynamics and hyperbolic systems of balance laws in continuum mechanics.  In: Toro E.F. (eds) Godunov Methods. Springer, Boston, MA, pages 745–761, 2001.
  * [50] E. Romenski, D. Drikakis, and E. Toro.  Conservative models and numerical methods for compressible two-phase flow.  Journal of Scientific Computing, 42(1):68, Jul 2009.
  * [51] E. Romenski, A. D. Resnyansky, and E. F. Toro.  Conservative hyperbolic formulation for compressible two-phase flow with different phase pressures and temperatures.  Quart. Appl. Math., 65(2):259–279, 2007.
  * [52] T. Ruggeri and M. Sugiyama.  Classical and Relativistic Rational Extended Thermodynamics of Gases.  Springer, Cham, 2021.
  * [53] K. Saleh.  Analyse et Simulation Numérique par Relaxation d’Écoulements Diphasiques Compressibles.  PhD thesis, Université Pierre et Marie Curie, Paris, 2012.
  * [54] K. Saleh and N. Seguin.  Some mathematical properties of a barotropic multiphase flow model.  ESAIM: ProcS, 69:70–78, 2020.
  * [55] R. Saurel and R. Abgrall.  A multiphase Godunov method for compressible multifluid and multiphase flows.  Journal of Computational Physics, 150(2):425–467, 1999.
  * [56] R. Saurel and O. LeMetayer.  A multiphase flow model for compressible flows with interfaces, shocks, detonation waves and cavitation.  J. Fluid Mech., 43:239–271, 2001.
  * [57] R. Saurel, F. Petitpas, and R. Abgrall.  Modelling phase transition in metastable liquids: Application to cavitating and flashing flows.  Journal of Fluid Mechanics, 607:313–350, 07 2008.
  * [58] F. Thein.  Results for Two Phase Flows with Phase Transition.  PhD thesis, Otto-von-Guericke Universität, Magdeburg, 2018.
  * [59] F. Thein, E. Romenski, and M. Dumbser.  Exact and numerical solutions of the Riemann problem for a conservative model of compressible two-phase flows.  Journal of Scientific Computing, 93(3):83, Nov 2022.
  * [60] E. F. Toro.  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction.  Berlin: Springer, 1997.
  * [61] W. Wagner and H.-J. Kretzschmar.  International Steam Tables.  Berlin - Heidelberg: Springer, 2008.
  * [62] A. Zein.  Numerical methods for multiphase mixture conservation laws with phase transition.  PhD thesis, Otto-von-Guericke University, Magdeburg, 2010.
  * [63] A. Zein, M. Hantke, and G. Warnecke.  Modeling phase transition for compressible two-phase flows applied to metastable liquids.  J. Comput. Phys., 229(8):2964–2998, 2008.
  * [64] A. Zein, M. Hantke, and G. Warnecke.  On the modelling and simulation of a laser-induced cavitation bubble.  Int. J. Num. Meth. Fluids., 73(2):172–203, 2013.