A high-order nonconservative approach for hyperbolic equations in fluid dynamics

TL;DR

This paper introduces a high-order nonconservative numerical approach for hyperbolic fluid dynamics equations, enabling oscillation-free solutions and effective handling of nonlinear equations of state, with extensions to multiphase and multidimensional flows.

Contribution

It demonstrates how to obtain relevant weak solutions using a pressure-based formulation, addressing long-standing issues with nonconservative schemes in fluid dynamics.

Findings

Relevance of pressure-based formulation for weak solutions

Oscillation-free solutions in nonconservative schemes

Extensions to multiphase and multidimensional flows

Abstract

It is well known, thanks to Lax-Wendroff theorem, that the local conservation of a numerical scheme for a conservative hyperbolic system is a simple and systematic way to guarantee that, if stable, a scheme will provide a sequence of solutions that will converge to a weak solution of the continuous problem. In [1], it is shown that a nonconservative scheme will not provide a good solution. The question of using, nevertheless, a nonconservative formulation of the system and getting the correct solution has been a long-standing debate. In this paper, we show how get a relevant weak solution from a pressure-based formulation of the Euler equations of fluid mechanics. This is useful when dealing with nonlinear equations of state because it is easier to compute the internal energy from the pressure than the opposite. This makes it possible to get oscillation free solutions, contrarily to classical conservative methods. An extension to multiphase flows is also discussed, as well as a multidimensional extension.

Tables

Table 1. Table 1: Parameters for the Cochran & Chan EOS ( 27 ).

${\rho }_0$	$A_1$	$E_1$	$A_2$	$E_2$	$\Gamma$
(kg/m3)	(GPa)		(GPa)
$1134$	$0.819181$	$4.52969$	$1.50835$	$1.42144$	$1.19$

Table 2. Table 2: Initial fluid properties and parameters for the test case.

Phase	Fluid	${\alpha }_1$	$\rho \left[\frac{kg}{m^3}\right]$	$\gamma$	$P_{\infty }\left[Pa\right]$
1	Epoxy	$0.5954$	$1185$	$2.43$	$5.3\cdot {10}^9$
2	Spinel	$0.4046$	$3622$	$1.62$	$141\cdot {10}^9$

Equations

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ \rho \mathbf{u}\\ E\end{array}\right)+\text{ div}\left(\begin{array}{c}\rho \mathbf{u}\\ \rho \mathbf{u}\otimes \mathbf{u}+p\text{Id}\\ (E+p)\mathbf{u}\end{array}\right)=0,$$

$$E=e+\frac{1}{2}\rho {\mathbf{u}}^2.$$

$$p=p(\rho ,e)=p(\rho ,E-\frac{1}{2}\rho {\mathbf{u}}^2).$$

$$c^2=\kappa h,h=\frac{e+p}{\rho }.$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ \rho \mathbf{u}\\ e\end{array}\right)+\left(\begin{array}{c}\text{ div }\rho \mathbf{u}\\ \text{div }(\rho \mathbf{u}\otimes \mathbf{u}+p\text{Id})\\ \mathbf{u}\cdot \nabla e+(e+p)\text{ div }\mathbf{u}\end{array}\right)=0$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ \rho \mathbf{u}\\ p\end{array}\right)+\left(\begin{array}{c}\text{ div }\rho \mathbf{u}\\ \text{div }(\rho \mathbf{u}\otimes \mathbf{u}+p\text{Id})\\ \mathbf{u}\cdot \nabla p+\rho c^2\text{ div }\mathbf{u}\end{array}\right)=0.$$

$$\frac{\partial U}{\partial t}+\frac{\partial f(U)}{\partial x}=0.$$

$$U_{j}^{n+1}=U_{j}^{n}-\lambda\big{(}\hat{f}_{j+1/2}-\hat{f}_{j-1/2})$$

$$U_j^{n+1}=U_j^n-\lambda ({\Phi }_j^{j+1/2}+{\Phi }_j^{j-1/2}).$$

$${\Phi }_j^{j+1/2}+{\Phi }_{j+1}^{j+1/2}=f(U_{j+1})-f(U_j)$$

$${\Phi }_j^{j+1/2}={\hat{f}}_{j+1/2}-f(U_j)\text{and}{\Phi }_{j+1}^{j+1/2}=f(U_{j+1})-{\hat{f}}_{j+1/2}.$$

$$\frac{\partial U}{\partial t}+\text{div }\mathbf{f}(U)=0.$$

$$\begin{split}U^{(0)}&=U^{n},\\ \frac{U^{(1)}-U^{{(0)}}}{\Delta t}&+\text{div }\mathbf{f}(U^{(0)})=0,\\ \frac{U^{(2)}-U^{{(0)}}}{\Delta t}&+\frac{1}{2}\bigg{(}\text{div }\mathbf{f}(U^{(0)})+\text{div }\mathbf{f}(U^{(1)})\bigg{)}=0,\\ U^{n+1}&=U^{(2)}.\end{split}$$

$$\frac{U^{(l+1)}-U^{(0)}}{\Delta t}+\text{ DIV }\mathcal{F}(U^{(l)},U^{(0)})=0$$

$$\text{ DIV }\mathcal{F}(U^{(l)},U^{{(0)}})=\frac{1}{2}\bigg{(}\text{div }\mathbf{f}(U^{(0)})+\text{div }\mathbf{f}(U^{(l)})\bigg{)}$$

$$\text{ DIV }\mathcal{F}(U^{(l)},U^{{(0)}})=\text{div }\mathbf{f}\bigg{(}\frac{U^{(0)}+U^{(l)}}{2}\bigg{)}$$

$$V_h=\{u\in L^2(\Omega ):\text{ for any }{K}_j,j\in \mathcal{J},u_{|K_j}\in {\mathbb{P}}^1(K_j)\}.$$

$$\begin{split}|C_{i}|\dfrac{U_{\sigma}^{(l+1)}-U_{\sigma}^{{({0})}}}{\Delta t}&+\sum_{K\ni\sigma}\Phi_{\sigma}^{K}(U^{(l)}_{h},U^{{(0)}}_{h})=0\\ {\text{with}\quad}\Phi_{\sigma}^{K}(U^{(l)}_{h},U^{{(0)}}_{h})&=\beta_{\sigma}^{K}(U^{(l)}_{h},U^{{(0)}}_{h})\bigg{(}\int_{K}\dfrac{U_{h}^{(l)}-U_{h}^{{(0)}}}{\Delta t}dx+\int_{\partial K}\mathcal{F}(U^{(l)}_{h},U^{{(0)}}_{h})\cdot{\mathbf{n}}{\,d\Gamma}\bigg{)}\end{split}$$

$$\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})=\frac{1}{N_{K}}\bigg{(}\int_{K}\dfrac{U_{h}^{(l)}-U_{h}^{{(0)}}}{\Delta t}dx+\int_{\partial K}\mathcal{F}(U^{(l)}_{h},U^{{(0)}}_{h})\cdot{\mathbf{n}}{\,d\Gamma}\bigg{)}+\alpha_{K}\big{(}\frac{U_{h}^{(l)}+U_{h}^{{(0)}}}{2}-\widehat{U}_{K}\big{)}$$

$${\hat{U}}_K=\sum _{\sigma \in K}\frac{U_{\sigma }^{(l)}+U_{\sigma }^{(0)}}{2}.$$

$$X=\sum _{\xi }{\mathbf{\ell }}_{\xi }(X){\mathbf{r}}_{\xi }.$$

$${\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})=\sum_{\xi}\mathbf{\ell}_{\xi}\big{[}\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{]}\mathbf{r}_{\xi}}$$

$$\mathbf{\ell}_{\xi}\big{[}\Phi^{K}\big{]}=\sum_{\sigma\in K}\mathbf{\ell}_{\xi}\big{[}\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{]},$$

$$\Phi_{\sigma}^{K,\star}(U^{(l)}_{h},U^{{(0)}}_{h})=\sum_{\xi}\beta_{\sigma}^{K,\star}\mathbf{\ell}_{\xi}\big{(}\Phi^{K}\big{)}\mathbf{r}_{\xi}$$

$$\beta_{\sigma}^{K,\star}=(1-\Theta_{\xi})\dfrac{\max\Big{(}\frac{\mathbf{\ell}_{\xi}\big{(}\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{)}}{\mathbf{\ell}_{\xi}\big{[}\Phi^{K}\big{]}},0\Big{)}}{\sum\limits_{\sigma^{\prime}\in K}\max\Big{(}\frac{\mathbf{\ell}_{\xi}\big{(}\Phi_{\sigma^{\prime}}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{)}}{\mathbf{\ell}_{\xi}\big{[}\Phi^{K}\big{]}},0\Big{)}}+\Theta_{\xi}\mathbf{\ell}_{\xi}\big{(}\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{)}$$

$$\Theta_{\xi}=\dfrac{\big{|}\mathbf{\ell}_{\xi}\big{(}\Phi_{\sigma}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{)}\big{|}}{\sum\limits_{\sigma^{\prime}\in K}\big{|}\mathbf{\ell}_{\xi}\big{(}\Phi_{\sigma^{\prime}}^{K,Rus}(U^{(l)}_{h},U^{{(0)}}_{h})\big{)}\big{|}}.$$

$$dE=de+\mathbf{u}\cdot d\mathbf{m}-\frac{1}{2}||\mathbf{u}|{|}^{2}d\rho ,$$

$$|C_{i}|\big{(}U_{i}^{n+1}-U_{i}^{n})+\sum_{j\in\mathcal{V}_{i}}\hat{\mathbf{f}}_{ij}=0$$

$$\delta \hat{\mathbf{f}}:=\sum _{j\in {\mathcal{V}}_i}{\hat{\mathbf{f}}}_{ij}.$$

$$|C_i|({\rho }_i^{n+1}-{\rho }_i^n)+\delta {\hat{f}}_{\rho }$$

$$|C_i|({\mathbf{m}}_i^{n+1}-{\mathbf{m}}_i^n)+\delta {\hat{f}}_{\mathbf{m}}$$

$$|C_i|(E_i^{n+1}-E_i^n)+\delta {\hat{f}}_E$$

Full text

1 Introduction

According to the Lax-Wendroff theorem, it is well known that, when considering the numerical approximation of a system of hyperbolic PDEs written in conservative form, the numerical scheme must be written in conservation form, too. It is also known that, for a sequence of meshes with characteristic sizes tending to zero, if a sequence of solutions remains bounded and if its subsequence converges in some norm in $L^{p}$, $p\geq 1$, then the limit solution is the weak solution of the original PDE. Moreover, if the scheme satisfies a discrete entropy inequality, then the limit solution will automatically satisfy an entropy inequality. If conservation is lost, then there is no hope to get any meaningful solution, see [1] for the analysis.

However, for engineering purposes, the conservative formulation of the behavior of a mechanical system is not necessarily the best one. Consider for example the Euler equations of fluid dynamics. The system of PDEs is

[TABLE]

supplemented by initial and boundary conditions. As usual, $\rho$ stands for the density, $\mathbf{u}$ for the velocity, and the total energy is

[TABLE]

The pressure $p$ is related to these variables via an equation of state (EOS):

[TABLE]

The system (1) is hyperbolic if $\kappa:=\dfrac{\partial p}{\partial e}>0$ since the speed of sound $c$ is defined by

[TABLE]

However, for engineering purposes, the relevant variables are not the conserved ones but rather the primitive ones, namely density, velocity and internal energy or pressure. When the solution is smooth, system (1) can be equivalently written as:

[TABLE]

or

[TABLE]

These equations are valid for smooth flows and cannot be considered for discontinuities. Nevertheless, there have been several attempts to solve the Euler equations with either formulation (3) or (4), including solutions with shocks. One example of such methods is Karni’s hybrid scheme [2, 3], where formulation (4) is used along slip lines only thanks to a switch in the scheme. In any case, this method violates strict conservation.

This has been a long ongoing debate on how to, nevertheless, use formulations (3) or (4) for the numerical approximation of Euler equations valid for all kinds of flows, e. g. with complex equations of state. In the case of nonlinear equations of state, i. e. when the pressure explicitly depends on the density, the pressure obtained by the numerical scheme cannot be uniform across contact discontinuities. The reason for that behaviour is that, on one hand, the density is evaluated from the mass conservation, and, on the other hand, one evaluates the pressure via energy and density. If, in addition, we want the pressure to be constant across contact discontinuities, this puts a constraint that is in general not compatible with the updated densities, momentum and total energy, see [4] for a short discussion.

In our knowledge, in the Eulerian framework, there exists only one approach to this problem which is described in [5]. In this paper, we propose a simpler and more general framework for dealing with nonconservative formulations. We solve equations (3) or (4) in a way which is compatible with local conservation and the continuity of pressure and velocity across contact discontinuities. To achieve this, we rely on a finite volume formulation that uses residuals instead of fluxes. In a flux formulation, the unknowns are approximations of the average values of the conserved variables, and they are balanced by the sum of normal fluxes across the boundary of the control volume. This assumes that the control volume has a polygonal shape. In general, these control volumes are interpreted as cells of a dual mesh which is made of simplices. In the residual formulation, one starts by a mesh whose elements are simplices, and interprets the unknowns as approximations of the point values of the conserved variables. These unknowns, for any given degree of freedom, are then updated by a sum of the local residuals over all the elements that share this degree of freedom. Given any element $K$, the local conservation is recovered by requiring that the sum of the local residuals for that element is the normal flux over the boundary of $K$ of some consistent approximation of the flux. It is easy to show that any flux formulation leads to a residual form, and the opposite is also true. However, the fluxes that are computed depend not only on the solution on both sides of the face of the control volume, see [6] for details.

The format of this paper is the following. In Section 2, we recall how one can get a residual distribution formulation for the system (3) that is equivalent to a flux formulation of (1). This enables us to get a relation on the increment of the energy that can be generalised for the residual formulation. In Section 3, we show how to use this principle, first on the energy-based formulation of the Euler equations (3) and then on the pressure-based formulation (4) for several kinds of equations of state. In Section 4, this is further generalised to multiphase flows where the phases may have very complex equations of state. Finally, we give some concluding remarks in Section 5.

2 From a conservative to a nonconservative formulation

2.1 A residual formulation of a finite volume scheme

The main advantage of the residual formulation can be understood from the one dimensional setting. Consider the problem:

[TABLE]

With standard notations, a generic finite volume scheme writes:

[TABLE]

where $\lambda=\Delta t/\Delta x$ and $\hat{f}_{j+1/2}$ is the flux between the states $U_{j}^{n}$ and $U_{j+1}^{n}$. High order accuracy in space amounts to tune the arguments of the flux, while high order in time can be reached via a strong stability preserving (SSP) scheme. To fix the conservation problem one must fix the scheme to recover a flux form, i.e to work directly with the fluxes. This is not easy from an algebraic point of view.

It is known, see for example [7] that any finite volume scheme can be rewritten in terms of a distribution of the residual. Consider for example, and for simplicity, the one dimensional case, its generalisation to any kind of control volume is straightforward, see again [7].

The residual formulation writes (in its simplest form) as

[TABLE]

The conservation is recovered if for any element $[x_{j},x_{j+1}]$ one gets

[TABLE]

for any order of accuracy. One can go from a flux formulation to a residual formulation by defining, for example:

[TABLE]

The local conservation is a consequence of relation (5). If we start from a nonconservative formulation in residual form, one can check the conservation if one can provide linear transformations of this residual to obtain a form satisfying (5).

2.2 Unsteady residual distribution formulation for the conservative case

Consider a multidimensional hyperbolic system in the form

[TABLE]

Recall the residual distribution approach from [8]. We start by a Runge-Kutta (RK) formulation: knowing $U^{n}$, we define:

[TABLE]

Next for $l=0,1$ we rewrite each sub-step as:

[TABLE]

where

[TABLE]

This will provide an approximation that is second order in time. Without loss of accuracy, we can also write it as

[TABLE]

and the adapted modifications for the RK scheme.

We assume that the computational domain $\Omega$ is covered by non-overlapping simplices $\{K_{j}\}_{j\in\mathcal{J}}$, $\Omega=\cup_{j\in\mathcal{J}}\{K_{j}\}$. The elements $K_{j}$ are segments in 1D, triangles/quadrilateral in 2D and tetrahedrons/hexahedrons in 3D. In order to simplify the notations, we denote by $\mathbb{P}^{1}(K)$ the set of polynomials of degree $1$ on triangles/tetrahedrons or by $\mathbb{Q}^{1}{(K)}$ the set of polynomials of degree $1$ on quadrilaterals/hexahedrons. Both guarantee a second order accurate approximation of any $C^{1}$ function. We introduce the approximation space

[TABLE]

We denote by $\sigma$ the degrees of freedom, i.e. the vertices of the element $K_{j}$. Then the residual distribution approximation of (6) reads: for any degree of freedom $\sigma$, define first $U^{{(0)}}=U^{n}$ and for $l=0,1$, do

[TABLE]

In (9), the parameter $\beta_{\sigma}^{K}(U^{(l)}_{h},U^{{(0)}}_{h})$ is a matrix (for Euler system of equations $3\times 3$ in 1D, $4\times 4$ in 2D, etc.), and it is constructed as a formal extension of the same construction for scalar problems that is guaranteed to have an $L^{\infty}$ stability bound, see [8] for details.

To compute $\beta_{\sigma}$, we can proceed as follows.

In (9), instead of the $\beta$-formulation, we first consider the Rusanov residual:

[TABLE]

with

[TABLE]

In (10), $N_{K}$ corresponds to the number of degrees of freedom in $K$.

• 2.

Then we consider the eigen-decomposition of $A_{\mathbf{n}}=A(\widehat{U}_{K})n_{x}+B(\widehat{U}_{K})n_{y}$ (in 2D, while in 1D it would simply be $A(\widehat{U}_{K})$. The matrices $A$ and $B$ are the Jacobians of the $x$- and $y$- component of the flux $\mathbf{f}$ with respect to the state $\widehat{U}_{K}$. Here, the vector ${\mathbf{n}}$ is a unit vector in the direction of the velocity field when it is nonzero, or any arbitrary direction otherwise. Of course, this direction is not relevant in one dimension. The right eigenvectors of $A_{\mathbf{n}}$ are denoted by $\{\mathbf{r}_{\xi}\}$. We denote by $\{\mathbf{\ell}_{\xi}\}$ the left eigenvectors of $A_{\mathbf{n}}$, so that any state $X$ can be written as:

[TABLE]

• 3.

Then we decompose, for any degree of freedom $\sigma$:

[TABLE]

We note that for any $\xi$,

[TABLE]

where $\Phi^{K}=\sum_{\sigma\in K}\Phi_{\sigma}^{K,Rus}$ is the total residual.

• 4.

The next step is to define $\Phi_{\sigma}^{K,\star}(U^{(l)}_{h},U^{{(0)}}_{h})$ as:

[TABLE]

Note that $\Theta_{\xi}\in[0,1]$. This guarantees that the scheme is second order in time and space and (formally) non-oscillatory, see [8, 9] for more details.

2.3 Discrete conservation for a nonconservative formulation

As a first exercise, let us show that any explicit scheme for (1) can be rewritten as an explicit scheme for (3) or (4). It is well known that

[TABLE]

where $\mathbf{m}=\rho\mathbf{u}$ is the momentum, and thus the question is to see how one can use this relation for discrete problems.

To simplify the notations, we start from a simple Euler time stepping method, but it is important to note, that more general algorithms can be treated similarly. Setting $U=(\rho,\mathbf{m},E)^{T}$, where $\mathbf{m}=\rho\mathbf{u}$, we start from the finite volume scheme

[TABLE]

where $C_{i}$ is the area/volume of the control volume, $\mathcal{V}_{i}$ is the set of neighbouring cells to cell $i$, and $\hat{\mathbf{f}}_{ij}$ is the numerical flux between the cells $C_{i}$ and $C_{j}$. In the following, since $\hat{\mathbf{f}}_{ij}$ doesn’t play any role, we write:

[TABLE]

This finite volume scheme in the component-wise form reads:

[TABLE]

Next, we introduce the $\Delta$ operator

[TABLE]

From (14) we see that

[TABLE]

and since $\overline{\rho\mathbf{u}}=\overline{\rho}\;\overline{\mathbf{u}}$ and $\Delta\mathbf{m}=\overline{\rho}\Delta\mathbf{u}+\underline{\mathbf{u}}\Delta\rho$, we can write

[TABLE]

From this last equation we see that

[TABLE]

The main problem is that (13a)-(13b)-(16) is identical to (13) with the same properties and same drawbacks, so no progress has yet been made.

Let us again consider the scheme (12), and in particular its flux term. It is known, see for example [7], that any finite volume scheme can be rewritten in terms of distribution of a residual. Consider for simplicity, the one dimensional case, its generalisation to any kind of control volume is straightforward, see again [7]. Relation (12), using standard notations, writes:

[TABLE]

i.e.

[TABLE]

where we have set

[TABLE]

Hence, on the element $[x_{j},x_{j+1}]$ we can define two residuals

[TABLE]

and we notice the conservation relation:

[TABLE]

Coming back to the system (13), written as a residual distribution scheme in the form (17), the residuals $\Phi_{\sigma}^{j+1/2}$, $\sigma=j$, $j+1$, have 3 components, namely $\Phi_{\sigma,\rho}^{j+1/2}$, $\Phi_{\sigma,\mathbf{m}}^{j+1/2}$ and $\Phi_{\sigma,E}^{j+1/2}$ for the mass, momentum and total energy, respectively. From (16a) and (16b) we see that we can define a residual for the internal energy as follows:

[TABLE]

We obviously have

[TABLE]

Defining the total residuals as

[TABLE]

for $\xi=e$ and $E$, we see that (20) rewrites as

[TABLE]

Relation (21) represents our target relation.

3 Conservative approximation of the Euler equations in primitive variables

If we are able to define the residuals ${\Phi_{j,e}^{j+1/2}}$ and ${\Phi_{j+1,e}^{j+1/2}}$ that satisfy (20)-(21), then using the same technique as in [10], we can show that if the conditions of the Lax-Wendroff theorem for the sequence of approximation hold, then the limit solution is a weak solution of (1). A sketch of the proof is recalled in A. In this section, we show how to achieve this requirement on system (3) using the residual formulation described in Section 2.2. To describe the principle of the method, we start with system (3) and we deal then later with system (4) in the case of non linear equations of state. The discussion will be general, however, the numerical experiments will be one dimensional in order to compare the results with exact solutions.

3.1 Residual distribution formulation for unsteady problems in non-conservative form

To describe our method we first consider system (3) and start by considering the first order case, and then extend it to the second order in time and space.

3.1.1 First order scheme

We start by a first order case in order to illustrate the method. The temporal scheme is a simple one step Euler scheme, and the space residual will also be first order accurate. We set $U=(\rho,\mathbf{m},e)^{T}$, and knowing $U^{n}$ at every degree of freedom, $U^{n+1}$ is obtained by:

[TABLE]

This scheme has the same format as the one of Section 2 that has been used to derive (21). To be specific, we calculate the Rusanov residual for the first two components of $U$ by applying formula (10), which for the first order gives

[TABLE]

For the internal energy, we define the Rusanov residual as:

[TABLE]

The choice of the approximation (\mathbf{u}\cdot\nabla e+\rho h\text{ div }\mathbf{u}\big{)}^{\star\star} in (23) is dictated only by the accuracy constraint, namely that for a smooth solution we must have

[TABLE]

where $k$ is the degree of the interpolation and $d$ is the number of spacial dimensions, see [11] for more details.

Experimentally, one can see that under a CFL like constraint, the numerical solution converges, but the limit is not a weak solution of the problem, as expected. The reason is that the conservation relation (20) is not satisfied. Hence, from the density and the momentum equation, we can compute the velocity at time $t_{n+1}$, and then we modify the residual on the internal energy by setting:

[TABLE]

The perturbations $r_{\sigma}^{K}$ are chosen such that the relation (20) is satisfied,

[TABLE]

There is no reason to favour one degree of freedom over another, therefore we set:

[TABLE]

3.1.2 Second order scheme

We shall use scheme (6) to get second order of accuracy. One of its properties is that again each stage of the scheme looks like a first order Runge-Kutta scheme but with a modified residual since one needs to take into account both the time and space increments in relation (21). More specifically, we set $U=(\rho,\mathbf{m},e)^{T}$ and proceed as follows.

For the first step: we first compute the velocity $\mathbf{u}^{(1)}$ and relation (21) is written with

[TABLE]

• 2.

For the second step, knowing the states $U^{{(0)}}=U^{n}$ and $U^{(1)}$, we first compute $\mathbf{u}^{(2)}=\mathbf{u}^{n+1}$ and then

[TABLE]

These modifications are mandatory to get proper conservation (see (14a)-(14b)).

This means that ${\Phi_{E}^{K}}$ is now:

[TABLE]

where the total energy flux $\mathcal{F}_{E}$ is evaluated according to (7) or (8). The approximation (\mathbf{u}\cdot\nabla e+\rho h\text{ div }\mathbf{u}\big{)}^{\star\star} in (23) satisfies second order accuracy e.g. by means of arithmetic averages of the variables $U^{(0)}$ and $U^{(l)}$ with $l=0,1$. The correction is:

[TABLE]

Remark: Once again it is important to emphasize, that it is possible to use the velocity quantities at $l+1$ to evaluate (26), since one first solves the mass and momentum equations for $l+1$ and only then focuses on the nonconservative internal energy (or pressure) equation.

3.2 Ensuring conservation in the case of non-linear equations of state

In the introduction, we have mentioned the problem of evaluating the pressure across contact discontinuities for non linear equations of state. Consider an equation of state (EOS) to be

[TABLE]

and assume it to be such that

[TABLE]

Note that for a perfect gas $\gamma=\kappa+1$, see [12] for details. Let us recall an example from [5]. We consider Cochran and Chan EOS given by

[TABLE]

We set up a shock tube problem with the EOS parameters listed in Table 1 and consider piecewise-constant initial data given by

$u=1000$ [m/s], $p=20\cdot 10^{9}$ [Pa]

• 2.

$\rho_{L}=1134$, $\rho_{R}=500$ [kg/m3]

Using the conservative scheme described in Section 3 we get the results displayed in Figure 1. This results are compared to those of [5], which have been obtained with a HLLC scheme with a MUSCL extrapolation on the physical variables.

In both cases, we see that the pressure and the velocity do not stay uniform as they should be. It is well known that even with the Godunov scheme, the evaluation of the pressure across the contact discontinuity can be problematic, even for single fluids, see [4, 13] for the analysis. The reason is that the equation of state that relates the internal energy to the density and the pressure can be highly nonlinear. The internal energy is obtained from the total energy and the kinetic energy, and thus the pressure via:

[TABLE]

where we have put conserved variables emphasised by using the notation $\langle\cdot\rangle$. The problem is that across a contact, the pressure is uniform, and in the case of a highly nonlinear EOS, there is no reason that the relation

[TABLE]

will guarantee a uniform pressure when $\langle\rho\rangle$ changes across a contact discontinuity. Up to our knowledge, the only Eulerian method that provides correct values is described in [5]. It is however quite complicated and tuned for Cartesian meshes.

A way to solve this issue is to start from (4) and use the same ideas as before. At the continous level, we have

[TABLE]

therefore

[TABLE]

and this is the relation to mimic in the numerical scheme. Hence, we wish to satisfy

[TABLE]

where $\widetilde{\dfrac{\partial e}{\partial\rho}}_{\sigma}$, $\widetilde{\dfrac{\partial e}{\partial p}}_{\sigma}$ are approximations of $\dfrac{\partial e}{\partial\rho}$ and $\dfrac{\partial e}{\partial p}$ at the degree of freedom $\sigma$ which need to be determined and $r_{\sigma}^{p}$ corresponds to the corrections on the pressure residuals. As before, we will assume that $r_{\sigma}^{p}=r^{p}$ is independent of $\sigma$. Note that we only perturb the pressure residual and not the density residual. The reason is that we wish to keep an explicit scheme: first we update the density, then the velocity and finally the pressure. If we can omit this constraint, more freedom can be obtained.

Relation (29) can be rewritten as:

[TABLE]

The question is how to define the terms $\widetilde{\dfrac{\partial e}{\partial\rho}}_{\sigma}$ and $\widetilde{\dfrac{\partial e}{\partial p}}_{\sigma}$. Once this is done we can get $r^{p}$. One constraint is that, if initially the pressure and the velocity are uniform, we keep this property at the next time step. We will start to see how relation (29) behaves in case of a 1D flow with uniform velocities and pressures. We start by the first order scheme, and then go to the second order.

The goal is to find ’good’ approximations of partial derivatives of the internal energy (28) so that the correction for the pressure vanishes for uniform pressures, without violating conservation. For this, we consider the behaviour of the internal energy increment. For any $\sigma$, we can write, for any $\lambda\in\mathbb{R}$

[TABLE]

with

[TABLE]

When $p$ and $u$ are uniform, then the right hand side of (29) reduces to

[TABLE]

If for any $\lambda$ we have $\widetilde{\dfrac{\partial e}{\partial p}}_{\sigma,\lambda}=0$, then $\Delta_{e}=0$. Otherwise, one can find $\lambda\in\mathbb{R}$ such that if

[TABLE]

then

[TABLE]

Clearly, we have to set $r^{p}=0$, if $p$ is uniform so that uniformity is preserved at the next time step.

It is important to note, that (31) is explicit in pressure in the case that the chosen equation of state guarantees that

[TABLE]

This is for example always achieved when using the Mie-Grüneisen EOS

[TABLE]

as done in the current work with the Cochran Chan EOS. In case the reader might be interested in an EOS not fulfilling (33), it is possible to simplify (31) and consider

[TABLE]

To detect contact discontinuities and preserve the uniformity of pressure and velocity, in the one-dimensional case we proceed as follows:

if

[TABLE]

then we set $r^{p}=0$, as this would mean that we are going across a contact wave.

• 2.

In any other case, $r^{p}$ is evaluated from (32b).

3.3 Numerical results

System (1) has been tested on a set of very demanding benchmark problems with two different typologies of equations of state.

3.3.1 Perfect gas EOS

The first test consists in a shock tube on a domain $x=[0,1]m$ with a diaphragm located at $x_{d}=0.5m$ at the initial conditions left and right of $x_{d}$

$\rho_{l}=100\;[\frac{kg}{m^{3}}]$, $u_{l}=0\;[\frac{m}{s}]$, $p_{l}=10^{9}\;[Pa]$,

• 2.

$\rho_{r}=1\;[\frac{kg}{m^{3}}]$, $u_{r}=0\;[\frac{m}{s}]$, $p_{r}=10^{5}\;[Pa]$.

The closing equations of state for system (1) are for a perfect gas and read $p(\rho,e)=(\gamma-1)\rho e$, with $\gamma=1.4$.

The results are shown at $t=45\;\mu s$ and have been obtained with a grid of $N=5000$ nodes and a $\text{CFL}=0.5$. The choice of a high number of cells for this test case is intended to show how the proposed numerical approximation converges to the exact solution for a very strong ”Sod”-like problem.

Figures 2 and 3 show the comparison of the results given by the exact solution and the conservative and nonconservative approximations presented in this work. The behaviour of these solutions show a good overlap. Both the conservative and nonconservative results are characterized by a glitch at the sonic point, which is less pronounced in the nonconservative case. The glitch itself can be easily corrected by adding some entropy correction, see e. g. [14] for a possible fix, but what matters in these results is that, this difference is also due to the fact, that the nonconservative approximation results in a slightly more diffused solution. This can be particularly seen in the zooms of the pressure (Fig. 2(c)), the density and of the velocity (Fig. 3(b) and (d)).

3.3.2 Nonlinear EOS

In order to check the quality of the approximation of the nonconservative scheme both for the pressure and energy formulations, a second testcase on a Riemann problem with a strong discontinuity has been evaluated with the choice of the Cochran-Chan EOS, as described in Section 3.2. The values for the EOS are those of Table 1 and we set111Note that the results are not sensitive to the choice of $\varepsilon$. $\varepsilon_{1}=\varepsilon=10^{-6}$ in (35). The considered domain $x=[0,1]m$ is split at $x_{d}=0.5m$ in a left and right state, where initially the values are set to

$\rho_{l}=1134\;[\frac{kg}{m^{3}}]$, $u_{l}=0\;[\frac{m}{s}]$, $p_{l}=2\cdot 10^{10}[Pa]$

• 2.

$\rho_{r}=120\;[\frac{kg}{m^{3}}]$, $u_{r}=0\;[\frac{m}{s}]$, $p_{r}=2\cdot 10^{5}[Pa]\,.$

The solution displayed in Figure 4 at a final time $50\cdot 10^{-6}{s}$ show an excellent approximation of the contact discontinuity wave for both nonconservative approximations. It is interesting to notice how the pressure formulation is less oscillatory then the energy one. In general, for both approximations, the shock propagates at the same speed.

4 Extension to multiphase flows

4.1 Kapila’s five equation model

Let us consider a different set of equations in the framework of compressible multiphase flows given by the five equation model of Kapila et al. [15] and shown in [16] to be the formal limit of the Baer and Nunziato model [17] when the relaxation parameters simultaneously tend to infinity though being proportional.

[TABLE]

In this two-phase system, $\alpha_{1}$ is the volume fraction of phase $\Sigma_{1}$, while the volume of the second phase $\Sigma_{2}$ is given by $\alpha_{2}=1-\alpha_{1}$. The density of the phase $\Sigma_{1}$ (respectively $\Sigma_{2}$) is $\rho_{1}$ (respectively $\rho_{2}$). The mixture density is given by $\rho=\alpha_{1}\rho_{1}+\alpha_{2}\rho_{2}$. We assume a single velocity $u$ and a single pressure $p$. This allows to consider mixture quantities for the energy and momentum conservation laws, while the mass conservation law is described for each phase separately. The internal energy of each phase $\Sigma_{j}$ is given by $e_{j}=e_{j}(p,\rho_{j})$ with $j=1,2$, and the mixture internal energy reads $e=\alpha_{1}e_{1}+\alpha_{2}e_{2}$. The total energy $E$ is the sum of the internal energy and the kinetic energy. In (36a), $c_{j}$ represents the speed of sound of phase $\Sigma_{j}$ and, in general, this transport equation is nonconservative. System (36) is a hyperbolic model and the mixture speed of sound $c$ is defined via

[TABLE]

with $c_{j}$ being the speed of sound of a phase $j$.

The Baer and Nunziato seven equation model [17], from which (36) has been derived, considers two phases which are described by a set of a mass, momentum and energy conservation laws for each phase and an additional transport equation which links the two sets of equations in terms of the volume fractions. In case of mechanical relaxation, which means that we assume a very large interface between the two phases, we can consider the pressures of each phase to be identical. The same also holds for the velocities. Following [16], since the pressures of the two phases are equal, their Lagrangian derivatives are equal, too. Therefore, it is possible to write that the entropies $s_{1}$,$s_{2}$ are constant and that $p_{1}(\rho_{1},s_{1})=p_{2}(\rho_{2},s_{2})$, leading to

[TABLE]

This allows to reformulate the transport equation as

[TABLE]

and to reduce the original seven equations model of [17] to the one of [15].

In the following, in order to be able to rewrite the system (36) in terms of primitive variables, we need the differential relations linking the pressure and the internal energy to the densities, $\alpha_{j}\rho_{j}$ and the volume fraction $\alpha_{1}$, since these are independent parameters. To achieve this, we start from:

[TABLE]

Since $e=\alpha_{1}e_{1}+\alpha_{2}e_{2}$, $p_{1}=p_{2}=p$ and $d(\alpha_{j}\rho_{j})=\alpha_{j}d\rho_{j}+\rho_{j}d\alpha_{j}$, we get:

[TABLE]

We rewrite

[TABLE]

[TABLE]

with

[TABLE]

4.2 Numerical results

Similarly as in the previous section, we test the system (36) on two different test cases with different equations of state and physical properties.

4.2.1 Stiffened Gas EOS

To test the capabilities of the nonconservative approximation in the context of compressible multiphase flows, a very severe benchmark problem with high differences in the pressure along a shock tube for epoxy and spinel has been taken from literature [18, 19, 20, 21].

The considered Riemann problem has a domain of $1m$ length and a discontinuity at $x_{d}=0.6m$. The parameters for each phase are shown in Table 2 and we use the stiffened gas EOS for both phases, which is part of the Mie Grüneisen EOS family and reads $p=(\gamma-1)\rho_{j}e_{j}-\gamma p_{\infty}$.

A mixture of epoxy and spinel is set up, with $u=u_{1}=u_{2}=0\;[\frac{m}{s}]$ both on the left and on the right of the discontinuity, while on the left the pressure is set to $p_{L}=p_{1L}=p_{2L}=2\cdot 10^{11}{[Pa]}$ and on the right to $p_{R}=p_{1R}=p_{2R}=1\cdot 10^{5}[Pa]$.

The final time is $t=29\;\mu s$ and the CFL number is set to $0.5$. The results are displayed in Figures 5 and 6.

The obtained solutions show an excellent contact discontinuity approximation. The difference in the plateaus is caused by the nonconservative form of the systems, which affects the numerical dissipation, while the rarefaction waves are identical. Note however that the complexity of the pressure formulation is lower: the pressure is a primary parameter, and it doesn’t need to be computed from the internal energy. In other words, there is no need to get $p$ knowing the mass fraction and the densities via the formula

[TABLE]

The inversion of the latter relation can be costly for nonlinear EOS. Our goal is certainly not to get the best possible solution in this context, and there exist methods that provide better results on this kind of problems, see e. g. [22, 23]. Our goal is only to show the versatility of our approach.

4.2.2 Mixed stiffened gas and Cochran-Chan EOS

The following test case is a more challenging variation of the previous test case of section 4.2.1. Let us assume the first phase to be described by a stiffened gas EOS with parameters $\gamma_{1}=2.43$, $p_{\infty}=5.3\;10^{9}$, and the other one by the Cochran-Chan EOS with the parameters given in Table 1. We initially set the same conditions as in section 4.2.1 and change only the volume fractions of the involved phases: $\alpha_{l}=0.5954$ and $\alpha_{r}=0.2$. The results are displayed in Figures 7 and 8.

The obtained solutions show again an excellent contact discontinuity approximation without any oscillations. The fact that the plateaus are on different levels is not unexpected, as pointed out before in 4.2.1.

5 Conclusions

In this paper, we have described a technique that enables the usage of a nonconservative formulation which nevertheless guarantees the conservation of the involved quantities. This method uses a residual distribution discretization with simple additional conditions which are able to guarantee the convergence to the correct weak solution. The emphasis is put on nonlinear equations of states where the pressure depends nonlinearly on the density. The presented approximation is then generalized to a multiphase system. The numerical tests have been done in one dimension in order to compare the results to the exact solutions. Extension to multidimensional problems can be done following the lines of [8] for the second order. Extension to higher order of accuracy could easily be obtained using the ideas of [24, 25], where the Runge-Kutta type timestepping scheme that we use here can be reinterpreted as a particular version of a deferred correction method. In this case, the relation (21) stays the same provided the residuals are defined according to [24] and the residual on the internal energy behaves like $O(h^{k+d})$, where $k$ is the expected order and $d$ the spatial dimension. We emphasize that the algebra remains identical. This will be demonstrated in a forthcoming paper.

Acknowledgements.

P. B. and S. T. have been funded by SNSF project # 200021_153604. R. A. has been funded in part by the same grant.

Appendix A A remark on the convergence to a weak solution

There are two ways of showing this. We provide the main idea for the system (3). Assuming that we have a scheme for this system that satisfies (21), and to simplify the derivations, we shall use the first order time scheme. Then one can define a residual for the total energy by simply setting

[TABLE]

By construction we have

[TABLE]

This shows that the sequence of solutions will converge to a weak solution. Using the results of [6], one can compute numerical fluxes for the density, momentum and total energy, so that local conservation is guaranteed. It is also simple to extend the proof of the Lax-Wendroff type theorem of [10] using the same conditions.