Viscous Regularization for the Non-equilibrium Seven-Equation Two-Phase Flow Model

Abstract

In this paper, a viscous regularization is derived for the non-equilibrium seven-equation two-phase flow model (SEM). This regularization, based on an entropy condition, is an artificial viscosity stabilization technique that selects a weak solution satisfying an entropy-minimum principle. The viscous regularization ensures nonnegativity of the entropy residual, is consistent with the viscous regularization for Euler equations when one phase disappears, does not depend on the spatial discretization scheme chosen, and is compatible with the generalized Harten entropies. We investigate the behavior of the proposed viscous regularization for two important limit-cases. First, a Chapman–Enskog expansion is performed for the regularized SEM and we show that the five-equation flow model of Kapila is recovered with a well-scaled viscous regularization. Second, a low-Mach asymptotic limit of the regularized seven-equation flow model is carried out whereby the scaling of the non-dimensional numbers associated with the viscous terms is determined such that an incompressible two-phase flow model, with a properly scaled regularization, is recovered. Both limit-cases are illustrated with one-dimensional numerical results, including two-phase flow shock tube tests and steady-state two-phase flows in converging-diverging nozzles. A continuous finite element discretization is employed for all numerical simulations.

Appendices

Appendix 1: Entropy Equation for the Multi-dimensional Seven Equation Model Without Viscous Regularization

This appendix provides the steps that lead to the derivation of the phasic entropy equation of the seven-equation two-phase flow Model [6]. For the purpose of this appendix, two phases are considered with no interphase mass or heat transfer and denoted by the indexes j and k. In the seven-equation two-phase flow Model, each phase obeys to the following set of equations (Eq. (55)):

$$\begin{aligned}&\partial _t \left( \alpha _k A\right) + A {u}_{int} \cdot {\nabla }\alpha _k = A \mu _P \left( P_k - P_j \right) \end{aligned}$$

(55a)

$$\begin{aligned}&\partial _t \left( \alpha _k \rho _k A \right) + {\nabla }\! \cdot \!\left( \alpha _k \rho _k {u}_k A \right) = 0\end{aligned}$$

(55b)

$$\begin{aligned}&\partial _t \left( \alpha _k \rho _k {u}_k A \right) + {\nabla }\! \cdot \!\left[ \alpha _k A \left( \rho _k {u}_k \otimes {u}_k + P_k {\mathbb {I}} \right) \right] \nonumber \\&\quad = \alpha _k P_k {\nabla }A + P_{int} A {\nabla }\alpha _k + A \lambda _u \left( {u}_j - {u}_k \right) \end{aligned}$$

(55c)

$$\begin{aligned}&\partial _t \left( \alpha _k \rho _k E_k A \right) + {\nabla }\! \cdot \!\left[ \alpha _k A {u}_k \left( \rho _k E_k + P_k \right) \right] \nonumber \\&\quad = P_{int} A {u}_{int} \cdot {\nabla }\alpha _k - A\mu _P {\bar{P}}_{int} \left( P_k-P_j \right) \nonumber \\&\qquad +\, \bar{{u}}_{int} A \lambda _u \left( {u}_j - {u}_k \right) \end{aligned}$$

(55d)

where $\rho _k$, ${u}_k$, $E_k$ and $P_k$ denote the density, velocity, specific total energy, and pressure of phase k, respectively. $\mu _P$ and $\lambda _u$ and the pressure and velocity relaxation parameters, respectively. We recall that we assume that the cross section A is only function of space: $\partial _t A = 0$ (a value of $A \ne 1$ is mostly of practical important for 1D nozzle problems). Variables with subscript int correspond to the interfacial variables; their definitions are given in Eq. (56).

$$\begin{aligned} \left\{ \begin{array}{lll} P_{int} = {\bar{P}}_{int} - \frac{{\nabla }\alpha _k}{|| {\nabla }\alpha _k ||} \frac{Z_k Z_j}{Z_k + Z_j} \left( {u}_k-{u}_j \right) \\ {\bar{P}}_{int} = \frac{Z_k P_j + Z_j P_k}{Z_k + Z_j} \\ {u}_{int} = \bar{{u}}_{int} - \frac{{\nabla }\alpha _k}{|| {\nabla }\alpha _k ||} \frac{P_k - P_j}{Z_k + Z_j} \\ \bar{{u}}_{int} = \frac{Z_k {u} _k + Z_j {u}_j}{Z_k + Z_j} \end{array} \right. \end{aligned}$$

(56)

where $Z_k = \rho _k c_k$ and $Z_j = \rho _j c_j$ are the impedances of phases k and j, respectively. The speed of sound is denoted by the symbol c.

The first step in proving the entropy minimum principle for Eq. (55) consists of recasting these equations using the primitive variables $(\alpha _k, \rho _k, {u}_k, e_k)$, where $e_k$ is the specific internal energy of phase k. We introduce the material derivative $\frac{D (\cdot )}{Dt} = \partial _t (\cdot ) + {u}_k \cdot {\nabla }(\cdot )$ for simplicity.

The continuity equation can be expressed as follows:

$$\begin{aligned} \alpha _k A \frac{D \rho _k}{Dt} + \rho _k A \mu _P \left( P_k-P_j \right) + \rho _k A \left( {u}_k-{u}_{int} \right) \cdot {\nabla }\alpha _k + \rho _k \alpha _k {\nabla }\! \cdot \!\left( A {u}_k \right) = 0 . \end{aligned}$$

(57)

The momentum and continuity equations are combined to yield an equation for the velocity:

$$\begin{aligned} \alpha _k \rho _k A \frac{D{u}_k}{Dt} + {\nabla }\left( \alpha _k A P_k \right) = \alpha _k P_k {\nabla }A + P_{int} A {\nabla }\alpha _k + A \lambda _u \left( {u}_j-{u}_k \right) . \end{aligned}$$

(58)

A kinetic energy equation is obtained by taking the vector scalar product of the previous results with ${u}_k$ to yield:

$$\begin{aligned} \alpha _k \rho _k A \frac{D\left( {u}_k^2/2\right) }{Dt} + {u}_k {\nabla }\left( \alpha _k A P_k \right) = {u}_k \Big ( \alpha _k P_k {\nabla }A + P_{int} A {\nabla }\alpha _k + A \lambda _u \left( {u}_j-{u}_k \right) \Big ) . \end{aligned}$$

(59)

The internal energy equation is obtained by subtracting the above kinetic energy equation from the total energy equation:

$$\begin{aligned} \alpha _k \rho _k A \frac{D e_k}{Dt} + \alpha _k P_k {\nabla }\! \cdot \!\left( A {u}_k \right)= & {} P_{int} A \left( {u}_{int}-{u}_k \right) \cdot {\nabla }\alpha _k - {\bar{P}}_{int} A \mu _P \left( P_k-P_j \right) \nonumber \\&+\, A \lambda _u \left( {u}_j-{u}_k \right) \cdot \left( \bar{{u}}_{int}-{u}_k \right) . \end{aligned}$$

(60)

In the next step, we assume the existence of a phasic entropy $s_k$ that is function of the density $\rho _k$ and the internal energy $e_k$. Using the chain rule,

$$\begin{aligned} \frac{Ds_k}{Dt} = (s_\rho )_k \frac{D \rho _k}{Dt} + (s_e)_k \frac{De_k}{Dt}, \end{aligned}$$

(61)

we combine the density and internal energy equations ($\rho _k (s_\rho )_k \times Eq.~(57) + (s_e)_k \times Eq.~(60))$ to obtain the following entropy equation:

$$\begin{aligned}&\alpha _k \rho _k A \frac{Ds_k}{Dt} + \underbrace{\alpha _k \left( P_k (s_e)_k + \rho _k^2 (s_\rho )_k \right) {\nabla }\! \cdot \!\left( A {u}_k \right) }_\text {(a)}\nonumber \\&\quad = (s_e)_k A \left[ P_{int}({u}_{int}-{u}_k)\cdot {\nabla }\alpha _k - {\bar{P}}_{int} A \mu _P (P_k-P_j) + A \lambda _u (\bar{{u}}_{int}-{u}_k) \cdot ({u}_j-{u}_k)\right] \nonumber \\&\qquad -\, \rho _k^2 (s_\rho )_k \left[ \mu _P A (P_k-P_j) + A({u}_k-{u}_{int}) \cdot {\nabla }\alpha _k\right] \end{aligned}$$

(62)

where $(s_e)_k$ and $(s_\rho )_k$ denote the partial derivatives of entropy $s_k$ with respect to the internal energy $e_k$ and the density $\rho _k$, respectively. The term denoted by (a) on the left-hand side of Eq. (62) can be set to zero by invoking the Gibbs relation from the second law of thermodynamics:

$$\begin{aligned} T_k ds_k = de_k - \frac{P_k}{\rho _k^2} d \rho _k \text { with } (s_e)_k = \frac{1}{T_k} \text { and } (s_\rho )_k = - \frac{P_k}{\rho _k^2} (s_e)_k \end{aligned}$$

(63)

which yields

$$\begin{aligned} P_k (s_e)_k + \rho _k^2 (s_\rho )_k = 0 . \end{aligned}$$

(64)

Finally, Eq. (62) is as follows:

$$\begin{aligned} ((s_e)_k)^{-1} \alpha _k \rho _k \frac{Ds_k}{Dt}= & {} \underbrace{\left[ P_{int} ({u}_{int}-{u}_k) + P_k ({u}_k-{u}_{int}) \right] \cdot {\nabla }\alpha _k}_\text {(b)} + \underbrace{\mu _P (P_k-P_j)(P_k-{\bar{P}}_{int})}_{(\mathrm{c})} \nonumber \\&+\, \underbrace{\lambda _u({u}_j-{u}_k)\cdot (\bar{{u}}_{int}-{u}_k)}_{(\mathrm{d})} \end{aligned}$$

(65)

The right-hand side of Eq. (65) has been split into three terms, (b), (c), and (d); next we analyze each of these terms separately. The terms (c) and (d) can be easily recast by using the definitions of $\bar{{u}}_{int}$ and ${\bar{P}}_{int}$ given in Eq. (56):

$$\begin{aligned} \mu _P (P_k-P_j)(P_k-{\bar{P}}_{int}) = \mu _P \frac{Z_k}{Z_k+Z_j} (P_j - P_k)^2 ,\nonumber \\ \lambda _u({u}_j-{u}_k)\cdot (\bar{{u}}_{int}-{u}_k) = \lambda _u \frac{Z_j}{Z_k+Z_j} ({u}_j - {u}_k)^2 . \end{aligned}$$

(66)

By definition, $\mu _P$, $\lambda _u$, and $Z_k$ are all positive. Thus, the above terms (c) and (d) are unconditionally positive.

We now inspect term (b). Once again, we use the definitions of $P_{int}$ and ${u}_{int}$ and the following relations:

$$\begin{aligned} {u}_{int}-{u}_k= & {} \frac{Z_j}{Z_k+Z_j}({u}_j-{u}_k) - \frac{{\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert } \frac{Pk-P_j}{Z_k+Z_j} ,\nonumber \\ P_{int}-P_k= & {} \frac{Z_k}{Z_k+Z_j} (P_j-P_k) - \frac{{\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert } \frac{Z_k Z_j}{Z_k+Z_j} ({u}_k-{u}_j) . \nonumber \end{aligned}$$

Then, term (b) becomes:

$$\begin{aligned}&\left[ P_{int} ({u}_{int}-{u}_k) + P_k ({u}_k-{u}_{int}) \right] \cdot {\nabla }\alpha _k = (P_{int}-P_k)({u}_{int}-{u}_k)\cdot {\nabla }\alpha _k\nonumber \\&\quad = \frac{Z_k}{\left( Z_k+Z_j \right) ^2} {\nabla }\alpha _k \cdot \left[ Z_j ({u}_j-{u}_k)(P_j-P_k)+\frac{{\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert } Z_j^2 ({u}_j-{u}_k)^2 \right. \nonumber \\&\qquad +\left. \frac{{\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert }(P_k-P_j)^2 + \frac{{\nabla }\alpha _k \cdot {\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert ^2}(P_k-P_j)Z_j ({u}_k-{u}_j) \right] \end{aligned}$$

(67)

The above equation is factorized by $\Vert {\nabla }\alpha _k \Vert$ and then recast under a quadratic form using $\frac{{\nabla }\alpha _k \cdot {\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert ^2} = 1$. This yields:

$$\begin{aligned}&\left[ ({u}_{int}-{u}_k)P_{int} + ({u}_k-{u}_{int})P_k \right] {\nabla }\alpha _k \nonumber \\&\quad = \Vert {\nabla }\alpha _k \Vert \frac{Z_k }{\left( Z_k+Z_j \right) ^2} \left[ Z_j ({u}_j-{u}_k) + \frac{{\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert }(P_k-P_j)\right] ^2 \end{aligned}$$

(68)

Thus, using Eq. (65), Eq. (66), Eq. (67) and Eq. (68), the entropy equation obtained in [6] holds and is recalled here for convenience:

$$\begin{aligned} (s_{e})_k^{-1} \alpha _k \rho _k A \frac{Ds_k}{Dt}= & {} \mu _P \frac{Z_k}{Z_k+Z_j} (P_j - P_k)^2 + \lambda _u \frac{Z_j}{Z_k+Z_j} ({u}_j -{u}_k)^2 \nonumber \\&+\, \Vert {\nabla }\alpha _k \Vert \frac{Z_k}{\left( Z_k+Z_j \right) ^2} \left[ Z_j ({u}_j-{u}_k)+\frac{{\nabla }\alpha _k}{\Vert {\nabla }\alpha _k \Vert }(P_k-P_j)\right] ^2. \end{aligned}$$

Appendix 2: Compatibility of the Viscous Regularization for the Seven-Equation Two-Phase Model with the Generalized Harten Entropies

We investigate in this appendix whether the viscous regularization of the seven-equation two-phase model derived in Sect. 3 is compatible with some or all generalized entropies identified in Harten et al. [27]. Considering the single-phase Euler equations, Harten et al. [27] demonstrated that a function $\rho {\mathscr {H}}(s)$ is called a generalized entropy and is strictly concave if ${\mathscr {H}}$ is twice differentiable and

$$\begin{aligned} {\mathscr {H}}' (s) \ge 0, \ \ \ \ {\mathscr {H}}'(s)c_p^{-1} - {\mathscr {H}}'' \ge 0, \ \forall \left( \rho , e \right) \in {\mathbb {R}}_+^2 , \end{aligned}$$

(69)

where $c_p \left( \rho , e \right) = T \partial _T s \left( \rho , e \right)$ is the specific heat at constant pressure (T is a function of e and $\rho$ through the equation of state). Because the seven-equation two-phase model was initially derived by assuming that each phase obeys the single-phase Euler equation, we want to investigate whether the above property still holds when considering the seven-equation model with viscous regularization included. To do so, we consider a phasic generalized entropy, ${\mathscr {H}}_k(s_k)$ and a phasic specific heat at constant pressure, $c_{p,k} \left( \rho _k, e_k \right) = T_k \partial _{T_k} s_k \left( \rho _k, T_k \right)$ characterized by Eq. (69). The objective is to find an entropy inequality verified by $\rho _k {\mathscr {H}}_k(s_k)$.

We start from the entropy inequality verified by $s_k$,

$$\begin{aligned} \alpha _k \rho _k A \frac{Ds_k}{Dt} = {f}_k \cdot {\nabla }s_k + {\nabla }\! \cdot \!\left( \alpha _k A \rho _k \kappa _k {\nabla }s_k \right) - \alpha _k \rho _k A \kappa _k Q_k + (s_e)_k \alpha _k A \rho _k \mu _k {\nabla }^s {u}_k : {\nabla }{u}_k.\nonumber \\ \end{aligned}$$

(70)

Eq. (70) is multiplied by ${\mathscr {H}}_k'(s_k)$ to yield:

$$\begin{aligned} \alpha _k \rho _k A \frac{D{\mathscr {H}}_k(s_k)}{Dt}= & {} {\nabla }\! \cdot \!\left( \alpha _k A \rho _k \kappa _k {\nabla }{\mathscr {H}}_k (s_k) \right) - {\mathscr {H}}_k''(s_k) \alpha _k A \kappa _k \rho _k \Vert {\nabla }s_k \Vert ^2 \nonumber \\&+\, {\mathscr {H}}_k'(s_k) {f}_k \cdot {\nabla }s_k - {\mathscr {H}}_k'(s_k)\alpha _k \rho _k A \kappa _k Q_k \nonumber \\&+\, {\mathscr {H}}_k'(s_k)(s_e)_k \alpha _k A \rho _k \mu _k {\nabla }^s {u}_k : {\nabla }{u}_k \end{aligned}$$

(71)

Let us now multiply the continuity equation of phase k by ${\mathscr {H}}_k (s_k)$ and add the result to the above equation to obtain:

$$\begin{aligned}&\partial _t \left( \alpha _k \rho _k A {\mathscr {H}}_k(s_k)\right) + {\nabla }\! \cdot \!\left( \alpha _k \rho _k {u}_k A {\mathscr {H}}_k(s_k) \right) \nonumber \\&\qquad -\,{\nabla }\! \cdot \!\left[ \alpha _k A \rho _k \kappa _k {\nabla }{\mathscr {H}}_k (s_k) + \alpha _k A \kappa _k {\mathscr {H}}_k (s_k) {\nabla }\rho _k \right. \left. + A \kappa _k \rho _k {\mathscr {H}}_k (s_k) {\nabla }\alpha _k\right] \nonumber \\&\quad = \underbrace{-{\mathscr {H}}_k''(s_k) \alpha _k A \kappa _k \rho _k \Vert {\nabla }s_k \Vert ^2 - {\mathscr {H}}_k'(s_k) \alpha _k A \kappa _k \rho _k Q_k}_{{\mathbb {T}}_0}\nonumber \\&\qquad +\, \underbrace{ {\mathscr {H}}_k'(s_k)(s_e)_k \alpha _k A \rho _k \mu _k {\nabla }^s {u}_k : {\nabla }{u}_k}_{{\mathbb {T}}_1} . \end{aligned}$$

(72)

As in Sect. 3, the left-hand side of Eq. (72) is split into two residuals denoted by ${\mathbb {T}}_0$ and ${\mathbb {T}}_1$ in order to study the sign of each of them. Obviously the sign of ${\mathbb {T}}_1$ is positive since it is assumed that ${\mathscr {H}}_k'(s_k) \ge 0$. To investigate the sign of ${\mathbb {T}}_0$, we use Eq. (69) to get:

$$\begin{aligned} - {\mathbb {T}}_0 \le {\mathscr {H}}_k'(s_k) \alpha _k A \kappa _k \rho _k \left( c_{p,k}^{-1} \Vert {\nabla }s_k\Vert ^2 + Q_k\right) . \end{aligned}$$

(73)

The right-hand side of Eq. (73) is a quadratic form that was already defined in Appendix 5 of [24] and can be recast in the matrix form $X^t_k {\mathbb {S}} X_k$ where ${\mathbb {S}}$ is a $2 \times 2$ matrix and the vector $X_k$ was previously defined in Sect. 3. In [24], matrix ${\mathbb {S}}$ is shown to be negative semi-definite which allows us to conclude that ${\mathbb {T}}_0$ is unconditionally positive using Eq. (73). Then, knowing the sign of the two residuals ${\mathbb {T}}_0$ and ${\mathbb {T}}_1$, we conclude that:

$$\begin{aligned}&\partial _t \left( \alpha _k \rho _k A {\mathscr {H}}_k(s_k)\right) + {\nabla }\! \cdot \!\left( \alpha _k \rho _k {u}_k A {\mathscr {H}}_k(s_k) \right) \nonumber \\&\quad -\, {\nabla }\! \cdot \!\left[ \alpha _k A \rho _k \kappa _k {\nabla }{\mathscr {H}}_k (s_k) + \alpha _k A \kappa _k {\mathscr {H}}_k (s_k) {\nabla }\rho _k \right. \left. +\, A \kappa _k \rho _k {\mathscr {H}}_k (s_k) {\nabla }\alpha _k\right] \ge 0 . \end{aligned}$$

Subsequently, we conclude that an entropy inequality is satisfied for all generalized entropies $\rho _k {\mathscr {H}}_k (s_k)$ when using the viscous regularization derived in Sect. 3 for the seven-equation two-phase model. Note that the above inequality holds as well for the total entropy of the system (i.e., summation over the phasic entropy statements).