A Projection Method for Compressible Generic Two-Fluid Model

We denote by the subscript $k=g,l$ for phase $g$ and phase $l$, respectively. Let $\phi_{k}$ be the volume fractions, $\rho_{k}$ be the densities, $\bm{u}_{k}$ be the velocities, $p_{k}$ be the pressures. The governing equations are given by

The stress tensor $\tau_{k}(\bm{u}_{k})$ are of the form

where $D(\bm{u}_{k})$ denotes the symmetric part of the velocity gradient $\nabla\bm{u}_{k}$, $\mu_{k}$ are the dynamic viscosities and $\lambda_{k}$ are the second viscosities. The drag force $F_{D,k}$ are assumed to be of the form:

where $C_{D}$ is assumed to be a positive constant.

Let $\alpha_{k}=\phi_{k}\rho_{k}$, (2.1) becomes

In the following, we will show that with given $\alpha_{g}$ and $\alpha_{l}$, the quantities $\phi_{g}$, $\phi_{l}$, $\rho_{g}$, $\rho_{l}$, and $p$ can be uniquely determined.

In this work, the stability issues are particularly payed attention and a discrete version associated to the numerical scheme is desired. Let us recall the a priori estimates associated to problem (2.8a)-(2.8c) with zero forcing term and zero velocities on $\Gamma$ for all time[15]. The following estimates present the positiveness of the mass, mass conservation, and the energy stability, respectively. For $k=g,l$, we have

• •

  Positiveness of the mass

  $$\alpha_{k}(x,t)>0,~{}~{}~{}\forall(x,t)\in\Omega_{T}$$

  (2.16)

• •

  Mass conservation

  $$\int_{\Omega}\alpha_{k}(x,t)dx=\int_{\Omega}\alpha_{k}(x,0)dx,~{}~{}~{}\forall t\in(0,T)$$

  (2.17)

• •

  Energy stability

  		$\displaystyle\sum_{k=g,l}\left[\frac{1}{2}\frac{d}{dt}\int_{\Omega}\alpha_{k}(x,t)\bm{u}_{k}(x,t)^{2}dx+\frac{d}{dt}\int_{\Omega}\alpha_{k}(x,t)e_{k}(\rho_{k}(x,t))dx\right.$		(2.18)
  		$\displaystyle\left.+\int_{\Omega}\phi_{k}(x,t)\tau_{k}(\bm{u}_{k}(x,t)):\nabla\bm{u}_{k}(x,t)dx\right]+\int_{\Omega}C_{D}\phi_{g}(x,t)\phi_{l}(x,t)|\bm{u}_{g}(x,t)-\bm{u}_{l}(x,t)|^{3}dx=0$

In the above, $e_{k}(\cdot)$ is the potential energy dervied from the equation of state such that

We may choose proper constants $\rho_{k,ref}$ so that the following expression makes sense:

The estimate (2.18) is a little bit different from that proposed by [15] because of extra capillary terms in their work. Nevertheless, they can be obtained by a completely same way.

We proceed the following to get the time-discrete solution:

• •

  Step 1. Prediction of the mass fraction $\widetilde{\alpha}_{k}^{m+1}$:

  $$\frac{\widetilde{\alpha}_{k}^{m+1}-\alpha_{k}^{m}}{\delta t}+\nabla\cdot(\widetilde{\alpha}_{k}^{m+1}\bm{u}_{k}^{m})=0$$

  (3.1)

• •

  Step 2. Solve the intermediate volume fraction $\widetilde{\phi}_{k}^{m+1}$ and the density $\widetilde{\rho}_{k}^{m+1}$ by

  $$\begin{cases}\widetilde{\phi}_{g}^{m+1}+\widetilde{\phi}_{l}^{m+1}=1\\ \widetilde{\phi}_{k}^{m+1}\widetilde{\rho}_{k}^{m+1}=\widetilde{\alpha}_{k}^{m+1},~{}~{}k=g,l\\ \zeta_{g}(\widetilde{\rho}_{g}^{m+1})=\zeta_{l}(\widetilde{\rho}_{l}^{m+1})\end{cases}$$

  (3.2)

• •

  Step 3. Renormalization of the intermediate pressure $\widetilde{p}_{k}^{m+1}$:

  $$\nabla\cdot\left(\frac{\widetilde{\phi}_{k}^{m+1}}{\widetilde{\rho}_{k}^{m+1}}\nabla\widetilde{p}_{k}^{m+1}\right)=\nabla\cdot\left(\sqrt{\frac{\widetilde{\phi}_{k}^{m+1}\widetilde{\phi}_{k}^{m}}{\widetilde{\rho}_{k}^{m+1}\widetilde{\rho}_{k}^{m}}}\nabla p^{m}\right)$$

  (3.3)

  with the boundary constraint:

  $$\frac{\widetilde{\phi}_{k}^{m+1}}{\widetilde{\rho}_{k}^{m+1}}\frac{\partial\widetilde{p}_{k}^{m+1}}{\partial n}=\sqrt{\frac{\widetilde{\phi}_{k}^{m+1}\widetilde{\phi}_{k}^{m}}{\widetilde{\rho}_{k}^{m+1}\widetilde{\rho}_{k}^{m}}}\frac{\partial p^{m}}{\partial n}~{}~{}~{}\text{on}~{}\Gamma$$

  (3.4)

• •

  Step 4. Solve the intermediate velocities $\widetilde{\bm{u}}_{k}^{m+1}$

  		$\displaystyle\frac{\widetilde{\alpha}_{k}^{m+1}\widetilde{\bm{u}}_{k}^{m+1}-\alpha_{k}^{m}\bm{u}_{k}^{m}}{\delta t}+\nabla\cdot(\widetilde{\alpha}_{k}^{m+1}\bm{u}_{k}^{m}\otimes\widetilde{\bm{u}}_{k}^{m+1})+\widetilde{\phi}_{k}^{m+1}\nabla\widetilde{p}_{k}^{m+1}$		(3.5)
  	$\displaystyle-$	$\displaystyle\nabla\cdot(\widetilde{\phi}_{k}^{m+1}\tau_{k}(\widetilde{\bm{u}}_{k}^{m+1}))+C_{D}\widetilde{\phi}_{g}^{m+1}\widetilde{\phi}_{l}^{m+1}|\widetilde{\bm{u}}_{g}^{m+1}-\widetilde{\bm{u}}_{l}^{m+1}|(\widetilde{\bm{u}}_{k}^{m+1}-\widetilde{\bm{u}}_{\widetilde{k}}^{m+1})=0$

• •

  Step 5. (Projection) Solve $\overline{\bm{u}}_{k}^{m+1}$, $p^{m+1}$, $\phi_{k}^{m+1}$, $\rho_{k}^{m+1}$:

  $$\begin{cases}\displaystyle\widetilde{\alpha}_{k}^{m+1}\frac{\overline{\bm{u}}_{k}^{m+1}-\widetilde{\bm{u}}_{k}^{m+1}}{\delta t}+\widetilde{\phi}_{k}^{m+1}\nabla(p^{m+1}-\widetilde{p}_{k}^{m+1})=0\\ \displaystyle\frac{\phi_{k}^{m+1}\rho_{k}^{m+1}-\phi_{k}^{m}\rho_{k}^{m}}{\delta t}+\nabla\cdot(\widetilde{\phi}_{k}^{m+1}\rho_{k}^{m+1}\overline{\bm{u}}_{k}^{m+1})=0\\ \phi_{g}^{m+1}+\phi_{l}^{m+1}=1\\ \rho_{k}^{m+1}=\zeta_{k}^{-1}(p^{m+1})\end{cases}$$

  (3.6)

  with the boundary constraint

  $$\frac{\partial p^{m+1}}{\partial n}=\frac{\partial\widetilde{p}_{k}^{m+1}}{\partial n}$$

  (3.7)

• •

  Step 6. Renormalization of the end-of-step velocities $\bm{u}_{k}^{m+1}$:

  $$\sqrt{\alpha_{k}^{m+1}}\bm{u}_{k}^{m+1}=\sqrt{\widetilde{\alpha}_{k}^{m+1}}\overline{\bm{u}}_{k}^{m+1},~{}~{}\alpha_{k}^{m+1}=\phi_{k}^{m+1}\rho_{k}^{m+1},~{}~{}k=g,l$$

  (3.8)

For simplicity, we consider the problem with $\bm{u}_{k}=0$ on $\Gamma$, $k=g,l$. The case with nonhomogeneous boundary conditions or outlet can be modified accordingly by a standard way.

Let us introduce finite element approximations $Y_{h}\subset H_{1}$ for densities $\rho_{k,h}$ (also for intermediate step $\widetilde{\rho}_{k,h}$), volume fractions $\phi_{k,h}$ (also for intermediate step $\widetilde{\phi}_{k,h}$), their products $\alpha_{k,h}$ (also for intermediate step $\widetilde{\alpha}_{k,h}$), and for pressure $p_{h}$ (also for intermediate step $\widetilde{p}_{k,h}$); $\bm{X}_{0,h}\subset\bm{H}_{0}^{1}$ for intermediate step velocities $\widetilde{\bm{u}}_{k,h}$.

To maintain the positiveness of $\alpha_{k}$, the weak formulation for Step 1 - (3.1) reads: For $m\geq 0$, find $\widetilde{\alpha}_{k,h}^{m+1}\in Y_{h}$ such that for $\widetilde{\alpha}_{k,h}^{m+1}|_{\Gamma_{in}}=\phi_{k,in}\rho_{k,in}$ and such that for all $q_{h}\in Y_{h}$ with $q_{h}|_{\Gamma_{in}}=0$,

The positiveness of the above scheme is shown in Remark 2.

At Step 2, the intermediate $\widetilde{\phi}_{k,h}^{m+1}$ and $\widetilde{\rho}_{k,h}^{m+1}$ can be obtain pointwisely by the given $\widetilde{\alpha}_{k,h}^{m+1}$ without error generation by spatial discretization since all of them belong to the same space $Y_{h}$. One may conduct a root-finding technique (e.g. Ridder’s method[27]) to solve $\widetilde{\phi}_{k,h}^{m+1}$ and $\widetilde{\rho}_{k,h}^{m+1}$ pointwisely.

The renormalization (Step 3) for intermediate pressure $\widetilde{p}_{k}^{m+1}$ can be proceeded by the following problem: For $m\geq 0$, find $\widetilde{p}_{k}^{m+1}\in Y_{h}$ such that for all $w_{h}\in Y_{h}$,

We note that the boundary terms are eliminated because of (3.4).

The weak formulation for Step 4 (advection-diffusion step) reads: For $m\geq 0$, find $\widetilde{\bm{u}}_{k,h}^{m+1}\in\bm{X}_{0,h}$ such that for all $\bm{v}_{h}\in\bm{X}_{0,h}$,

The above equation can be solved separately for $g$ and $l$ and the nonlinear term can be tackled by iteration (see for example [28]).

The projection step (Step 5) is handled by (3.10) and the according weak formulation reads: For $m\geq 0$, find $\alpha_{k}^{m+1}\in Y_{h}$ such that for all $q_{h}\in Y_{h}$,

with

The above system can be solved iteratively by the way provided in Remark 1. In (3.16), functions of $\alpha_{g,h}^{m+1}$ and $\alpha_{l,h}^{m+1}$: $C^{2}$, $\rho_{l}$, and $\rho_{g}$ are taken to be their projection on $Y_{h}$. Similar to (3.13), the boundary terms of (3.15) are eliminated because of the boundary constraint (3.7).