# A MUSCL-scheme on staggered grids for the Euler equations on unstructured meshes

Charbel Ghosn, Thierry Goudon, Sebastian Minjeaud

► To cite this version:

Charbel Ghosn, Thierry Goudon, Sebastian Minjeaud. A MUSCL-scheme on staggered grids for the Euler equations on unstructured meshes. FVCA 2023 - International Conference on Finite Volumes for Complex Applications, Oct 2023, Strasbourg, France. pp.141-149, 10.1007/978-3-031-40860-1_15 hal-04294739

HAL Id: hal-04294739

https://hal.science/hal-04294739v1

Submitted on 20 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A MUSCL-scheme on staggered grids for the Euler equations on unstructured meshes

Charbel Ghosn

Thierry Goudon

Sebastian Minjeaud

Université Côte d'Azur, Inria, CNRS, LJAD

## Abstract

We set up a MUSCL approach for the 2D DDFV-like scheme introduced in [3] for the full Euler equations, discretized with unstructured meshes. It is a finite volume staggered discretization framework where numerical densities, energies and velocities are stored on different locations. For the MUSCL scheme, a combination of several approaches is needed for the reconstruction of the densities and the velocities, in order to reach the second order accuracy. This work focuses on the reconstruction of the mass density and internal energy. The improved accuracy is already sensitive on problems containing contact discontinuities. We provide a set of numerical simulations to illustrate the enhanced accuracy.

## 1 Introduction

We consider in this work the Euler system of gas dynamics

$$ egin{cases} 
\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \\
\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p = 0, \\
\partial_t (\rho E) + \nabla \cdot (\rho E \mathbf{u}) + \nabla \cdot (p \mathbf{u}) = 0,
\end{cases} 
\quad (1) $$

where $\rho$, $\mathbf{u}$, $E$ and $p$ stand for the mass density, the velocity field, the total energy and the pressure, respectively. These unknowns depend on the time and space variables $(t, x) \in [0, \infty) \times \Omega$ with $\Omega \subset \mathbb{R}^2$ a polygonal bounded domain. System (1) is supplemented by the following equation of state

$$ E = \frac{\|\mathbf{u}\|^2}{2} + e \quad \text{and} \quad p = (\gamma - 1)\rho e, $$

where $e$ is the internal energy and $\gamma > 1$.

We aim at building a second order scheme on general staggered grids. Our approach is based on the first order scheme presented in [3] and a multislope MUSCL reconstruction of the variables (see [2] for a related but different approach). Since we work on staggered grids, the densities and velocities are stored
<figure><img src="image_2.png" /><figcaption>Figure 1: Meshes and associated notations.</figcaption></figure>

at different locations of the mesh. The densities are cell-based variables and the velocities are defined on vertices. Thus, the reconstruction procedure will not be the same for these variables. We present here the MUSCL reconstruction of the mass density $\rho$ and the density of internal energy $\rho e$, the reconstruction for the velocity being an ongoing work. We use a multislope method to reconstruct the densities [4], which will be used to define the mass and internal energy fluxes. This is not enough to obtain a second order convergence but we show with the numerical simulations of 1D contact discontinuities (where the velocity is constant) that it already provides improved results.

## 2 Notation: meshes, unknowns

We construct three partitions of the domain $\Omega$: the primal mesh, the dual mesh and the diamond mesh. The steps of the construction are illustrated in Fig. 1 (see details in [3]).

The primal mesh $\mathfrak{M}$, in blue in Fig. 1, is a partition of $\Omega$ by polygonal convex subsets $K$ called “primal cells”. The centers of these cells are labeled $x_K$. The dual mesh $\mathfrak{M}^*$, in red in Fig. 1, consists of cells “built around the vertices $x_{K^*}$” of the primal mesh. The dual cell $K^*$ is formed by joining the centers $x_K$ of all cells having $K^*$ as a vertex. The diamond mesh $\mathfrak{D}$, in green in Fig. 1, is made of quadrilateral cells $D_{\sigma,\sigma^*}$ obtained by joining the endpoints of the edge $\sigma = [x_{K^*}, x_{L^*}]$ of the primal mesh to the centers $x_K$ and $x_L$ of the primal cells that share this edge. The segment $\sigma^* = [x_K, x_L]$ is an edge of the dual mesh.

We do not pay attention in this article to the description of the meshes near the boundary $\partial\Omega$ and, more generally, to boundary conditions. We refer the reader to [3] for further details. As in [3], we use the following notation: we denote $\mathfrak{s} = D_{\sigma,\sigma^*}|_{D_{\sigma',\sigma'^*}}$ the face separating two diamond cells $D_{\sigma,\sigma^*}$ and $D_{\sigma',\sigma'^*}$; for $K \in \mathfrak{M}$, we denote $\mathfrak{D}_K = \{D_{\sigma,\sigma^*} \in \mathfrak{D}, \sigma \in \partial K\}$; for $K^* \in \mathfrak{M}^*$, we similarly denote $\mathfrak{D}_{K^*} = \{D_{\sigma,\sigma^*} \in \mathfrak{D}, \sigma^* \in \partial K^*\}$; for a cell $X$ of $\mathfrak{M}$, $\mathfrak{M}^*$
or $\mathfrak{D}$ and for $\mathfrak{r} \in \partial X$, we define a unit vector $\mathbf{n}_{X,\mathfrak{r}}$ normal to the face $\mathfrak{r}$ of the cell $X$ and pointing outwards: $\mathbf{n}_{K,\sigma}$ (with $\sigma \in \partial K$ for $K \in \mathfrak{M}$), $\mathbf{n}_{K^*,\sigma^*}$ (with $\sigma^* \in \partial K^*$ for $K^* \in \mathfrak{M}^*$), and $\mathbf{n}_{D_{\sigma,\sigma^*},\mathfrak{s}}$ (with $\mathfrak{s} \in \partial D_{\sigma,\sigma^*}$ for $D_{\sigma,\sigma^*} \in \mathfrak{D}$).

The unknowns of the mesh are defined as piecewise constants over the cells:

1. The mass density $(\rho_{\sigma,\sigma^*})_{D_{\sigma,\sigma^*} \in \mathfrak{D}}$ and the internal energy $(e_{\sigma,\sigma^*})_{D_{\sigma,\sigma^*} \in \mathfrak{D}}$ are piecewise constants over the diamond cells. We set $p_{\sigma,\sigma^*} = (\gamma - 1)\rho_{\sigma,\sigma^*}e_{\sigma,\sigma^*}$.

2. The numerical velocity fields $(\mathbf{u}_K)_{K \in \mathfrak{M}}$ and $(\mathbf{u}_{K^*})_{K^* \in \mathfrak{M}^*}$ are piecewise constants over the primal and dual cells, respectively.

### 3 A MUSCL reconstruction on staggered grids

Let us first recall how the first order scheme is constructed in [3], the MUSCL scheme being described in a second part.

**The first order scheme.** The time discretization is explicit. We denote by $\delta t$ the time step and $\bar{q}$ the update at time $t + \delta t$ of any quantity $q$ at time $t$. The space discretization is a finite volume scheme based on a splitting of the mass flux inspired from the kinetic framework which involves the sound speed $c(e) = \sqrt{\gamma(\gamma - 1)e}$ of system (1). Let

$$ \mathcal{F}^{+}(\rho, c, u) = \int_{\xi>0} \xi \mathcal{M}_{[\rho, c, u]}(\xi) d\xi = \begin{cases} 0 & \text{if } u \le -c, \\ \frac{\rho}{4c}(u+c)^2 & \text{if } |u| \le c, \\ \rho u & \text{if } u \ge c, \end{cases} $$

where

$$ \xi \in \mathbb{R} \mapsto \mathcal{M}_{[\rho, c, u]}(\xi) = \frac{\rho}{2c} \mathbf{1}_{|\xi-u| \le c}. $$

The function $\mathcal{F}^{-}$ is defined by $\mathcal{F}^{-}(\rho, c, u) = -\mathcal{F}^{+}(\rho, c, -u)$, so that the functions $\mathcal{F}^{\pm}$ satisfies consistency property $\mathcal{F}^{+}(\rho, c, u) + \mathcal{F}^{-}(\rho, c, u) = \rho u$. We thus define the mass flux $\mathcal{F}_{D_{\sigma,\sigma^*},\mathfrak{s}}$ from the diamond cell $D_{\sigma,\sigma^*}$ through the interface $\mathfrak{s} = D_{\sigma,\sigma^*}|_{D_{\sigma',\sigma^*}}$ using the upwind principle as follows

$$ \mathcal{F}_{D_{\sigma,\sigma^*},\mathfrak{s}} = \mathcal{F}_{D_{\sigma,\sigma^*},\mathfrak{s}}^{+} + \mathcal{F}_{D_{\sigma,\sigma^*},\mathfrak{s}}^{-} $$

with $\mathcal{F}_{D_{\sigma,\sigma^*},\mathfrak{s}}^{+} = \mathcal{F}^{+}(\rho_{\sigma,\sigma^*}, c_{\mathfrak{s}}, u_{D_{\sigma,\sigma^*},\mathfrak{s}})$ and $\mathcal{F}_{D_{\sigma,\sigma^*},\mathfrak{s}}^{-} = \mathcal{F}^{-}(\rho_{\sigma',\sigma^*'}, c_{\mathfrak{s}}, u_{D_{\sigma,\sigma^*},\mathfrak{s}})$.

The sound speed and the normal velocity used in the expression above are averaged values defined for $\mathfrak{s} = D_{\sigma,\sigma^*}|_{D_{\sigma',\sigma^*}} = [x_K, x_{K^*}]$ by:

$$ u_{D_{\sigma,\sigma^*},\mathfrak{s}} = \frac{\mathbf{u}_K + \mathbf{u}_{K^*}}{2} \cdot \mathbf{n}_{D_{\sigma,\sigma^*},\mathfrak{s}} \quad \text{and} \quad c_{\mathfrak{s}} = c \left( \frac{e_{\sigma,\sigma^*} + e_{\sigma',\sigma^*'}}{2} \right). $$
The discrete mass equation on a cell $D_{\sigma, \sigma^*} \in \mathcal{D}$ is given by

$$ \frac{\bar{\rho}_{\sigma, \sigma^*} - \rho_{\sigma, \sigma^*}}{\delta t} + \frac{1}{|D_{\sigma, \sigma^*}|} \sum_{s \in \partial D_{\sigma, \sigma^*}} |\mathfrak{s}| \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}} = 0. $$

The averaged density on a cell $K$ of the primal mesh is defined by

$$ \rho_K = \sum_{D_{\sigma, \sigma^*} \in \mathcal{D}_K} \frac{|D_{\sigma, \sigma^*} \cap K|}{|K|} \rho_{\sigma, \sigma^*} \quad \text{for } K \in \mathfrak{M} $$

and on a cell $K^*$ of the dual mesh, we set

$$ \rho_{K^*} = \sum_{D_{\sigma, \sigma^*} \in \mathcal{D}_{K^*}} \frac{|D_{\sigma, \sigma^*} \cap K^*|}{|K^*|} \rho_{\sigma, \sigma^*} \quad \text{for } K^* \in \mathfrak{M}^*. $$

To these densities we associate averaged mass fluxes $\mathcal{F}_{K, \sigma}$ outgoing from a primal cell $K$ and $\mathcal{F}_{K^*, \sigma^*}$ outgoing from a dual cell $K^*$, with

$$ \mathcal{F}_{K, \sigma}^{\pm} = \frac{|D_{\sigma, \sigma^*} \cap K|}{|D_{\sigma, \sigma^*}|} \sum_{s \in \partial D_{\sigma, \sigma^*}} \frac{|\mathfrak{s}|}{|\sigma|} \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{\pm} - \frac{|D_{\sigma, \sigma^*} \cap L|}{|D_{\sigma, \sigma^*}|} \sum_{s \in \partial D_{\sigma, \sigma^*}} \frac{|\mathfrak{s}|}{|\sigma|} \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{\mp} $$

$$ \mathcal{F}_{K^*, \sigma^*}^{\pm} = \frac{|D_{\sigma, \sigma^*} \cap K^*|}{|D_{\sigma, \sigma^*}|} \sum_{s \in \partial D_{\sigma, \sigma^*}} \frac{|\mathfrak{s}|}{|\sigma^*|} \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{\pm} - \frac{|D_{\sigma, \sigma^*} \cap L^*|}{|D_{\sigma, \sigma^*}|} \sum_{s \in \partial D_{\sigma, \sigma^*}} \frac{|\mathfrak{s}|}{|\sigma^*|} \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{\mp} $$

Using these fluxes, we can give the definition to the momentum fluxes $\mathcal{G}_{K, \sigma}$ for the primal cells and $\mathcal{G}_{K^*, \sigma^*}$ for the dual cells

$$ \mathcal{G}_{K, \sigma} = \mathcal{F}_{K, \sigma}^{+} \mathbf{u}_K + \mathcal{F}_{K, \sigma}^{-} \mathbf{u}_L \quad \text{and} \quad \mathcal{G}_{K^*, \sigma^*} = \mathcal{F}_{K^*, \sigma^*}^{+} \mathbf{u}_{K^*} + \mathcal{F}_{K^*, \sigma^*}^{-} \mathbf{u}_{L^*}. $$

The discrete momentum equation now reads

$$ \frac{\bar{\rho}_K \bar{\mathbf{u}}_K - \rho_K \mathbf{u}_K}{\delta t} + \frac{1}{|K|} \sum_{D_{\sigma, \sigma^*} \in \mathcal{D}_K} |\sigma| \mathcal{G}_{K, \sigma} + (\nabla_d p)_K = 0, $$

$$ \frac{\bar{\rho}_{K^*} \bar{\mathbf{u}}_{K^*} - \rho_{K^*} \mathbf{u}_{K^*}}{\delta t} + \frac{1}{|K^*|} \sum_{D_{\sigma, \sigma^*} \in \mathcal{D}_{K^*}} |\sigma^*| \mathcal{G}_{K^*, \sigma^*} + (\nabla_d p)_{K^*} = 0, $$

with a suitable definition of the pressure gradients [3].

Finally, for the discretization of the internal energy equation, we define the following numerical fluxes, for all $D_{\sigma, \sigma^*} \in \mathcal{D}$ and $\mathfrak{s} = D_{\sigma, \sigma^*} |D_{\sigma', \sigma'^*}$,

$$ \mathcal{E}_{D_{\sigma, \sigma^*}, \mathfrak{s}} = e_{\sigma, \sigma^*} \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{+} + e_{\sigma', \sigma'^*} \mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{-}. $$
The discrete internal energy equation is given by

$$ \frac{\bar{\rho}_{\sigma, \sigma^*} \bar{e}_{\sigma, \sigma^*} - \rho_{\sigma, \sigma^*} e_{\sigma, \sigma^*}}{\delta t} + \frac{1}{|D_{\sigma, \sigma^*}|} \sum_{\mathfrak{s} \in \partial D_{\sigma, \sigma^*}} |\mathfrak{s}| \mathcal{E}_{D_{\sigma, \sigma^*}, \mathfrak{s}} + p_{\sigma, \sigma^*} (\nabla_d \cdot \bar{\mathbf{u}})_{\sigma, \sigma^*} = \mathcal{R}_{\sigma, \sigma^*}, \quad \forall D_{\sigma, \sigma^*} \in \mathcal{D} $$

where $\nabla_d \cdot \bar{\mathbf{u}}$ is the discrete divergence operator, see [3], and $\mathcal{R}_{\sigma, \sigma^*}$ the correction term which should keep track of possible discontinuities. This remainder is defined so that it exactly balances the kinetic energy contributions that appear when summing the internal energy equation and the kinetic energy equation. This reconstruction allows us to derive a conservative discrete equation for an averaged total energy [3]. This construction shares many similarities with the DDFV framework since it works on dual grids and the discretization aims at preserving the duality of div-grad operators, which is crucial in order to preserve local conservations.

**The MUSCL-scheme.** We discuss how we adapt the MUSCL principles to the staggered framework. We first reconstruct second order quantities at edges of the cells (primal, dual or diamond) depending on the domain where the variables are defined. Then, concerning the discretization of the mass flux, we keep unchanged the velocity defined at the interface $\mathfrak{s}$ and we shall replace the UpWind value $\rho_{\sigma, \sigma^*}$ by a MUSCL reconstruction $\rho_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML}$ of the density: it defines the upgraded mass flux $\mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML}$. For the internal energy, we combine the upgraded mass fluxes $\mathcal{F}_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML}$ with a MUSCL reconstruction of the internal energy defined from the ratio $\frac{(\rho e)_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML}}{\rho_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML}}$, where $(\rho e)_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML}$ is the MUSCL reconstruction of $\rho_{\sigma, \sigma^*} e_{\sigma, \sigma^*}$ at the interface $\mathfrak{s}$. We follow the multislope method introduced in [4]. We compute the reconstructed values at the centers of the interfaces $\mathbf{M}_{\mathfrak{s}}$ (see Fig. 2). As in the original MUSCL method, both a backward and a forward scalar slopes, respectively denoted $s_{\sigma, \mathfrak{s}}^{-}$ and $s_{\sigma, \mathfrak{s}}^{+}$, are computed for each interface $\mathfrak{s}$ of a given diamond cell $D_{\sigma, \sigma^*}$. In a classical way, we use a limiter function $\phi(s_{\sigma, \mathfrak{s}}^{-}, s_{\sigma, \mathfrak{s}}^{-})$ to ensure that no unphysical oscillation is introduced. Therefore, the reconstructed values read as follows:

$$ \rho_{D_{\sigma, \sigma^*}, \mathfrak{s}}^{ML} = \rho_{\sigma, \sigma^*} + \phi(s_{\sigma, \mathfrak{s}}^{-}, s_{\sigma, \mathfrak{s}}^{+}) \cdot \|\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_{\mathfrak{s}}\|, $$

where $\mathbf{C}_{D_{\sigma, \sigma^*}}$ is the center of the diamond cell $D_{\sigma, \sigma^*}$.

We discuss next the building of the slopes. For that, we denote by $\mathcal{W}_{\sigma, \sigma^*}$ the set of the diamond cells sharing at least a vertex with $D_{\sigma, \sigma^*}$. The idea is to determine two points $\mathbf{H}_{\mathfrak{s}}^{+}$ and $\mathbf{H}_{\mathfrak{s}}^{-}$ both located on the axis $(\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_{\mathfrak{s}})$, respectively backward and forward the point $\mathbf{C}_{D_{\sigma, \sigma^*}}$ (see Fig. 2). These points are a priori neither vertices of the grid, nor element centers. However, these points are located by construction on a line joining two element centers so that densities $\rho_{\mathbf{H}_{\mathfrak{s}}^{-}}$ and $\rho_{\mathbf{H}_{\mathfrak{s}}^{+}}$ at these points can be obtained by a linear interpolation. The backward and forward slopes are then computed as follows:
<figure><img src="image_3.png" /><figcaption>Figure 2: Forward and backward points $\mathbf{H}_s^+$ and $\mathbf{H}_s^-$</figcaption></figure>

$$s_{\sigma, s} = \frac{\rho_{\sigma, \sigma^*} - \rho_{\mathbf{H}_s^-}}{\|\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{H}_s^-\|}, \quad s_{\sigma, s}^+ = \frac{\rho_{\mathbf{H}_s^+} - \rho_{\sigma, \sigma^*}}{\|\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{H}_s^+\|}.$$

We present next the process to determine the points $\mathbf{H}_s^+$ and $\mathbf{H}_s^-$. Let $D_1^- \in W_{\sigma, \sigma^*}$ be the most backward neighboring diamond cell of $D_{\sigma, \sigma^*}$ with respect to the direction $(\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s)$ and $\mathbf{C}_{D_1^-}$ its center, in the sense that

$$\cos (\mathbf{C}_{D_1^-} \mathbf{C}_{D_{\sigma, \sigma^*}}, \mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s) = \max_{D \in W_{\sigma, \sigma^*}} \cos (\mathbf{C}_D \mathbf{C}_{D_{\sigma, \sigma^*}}, \mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s). \quad (2)$$

Let now $\mathbf{C}_{D_2^-}$ be the center of the next most backward diamond cell, provided that it is located on the other side of the axis $(\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s)$, that is

$$\cos (\mathbf{C}_{D_2^-} \mathbf{C}_{D_{\sigma, \sigma^*}}, \mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s) = \max_{D \in \bar{W}_{\sigma, \sigma^*}} \cos (\mathbf{C}_D \mathbf{C}_{D_{\sigma, \sigma^*}}, \mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s), \quad (3)$$

where $\bar{W}_{\sigma, \sigma^*}$ is the set of diamonds $D \in W_{\sigma, \sigma^*}$ different from $D_1^-$ such that

$$\sin (\mathbf{C}_D \mathbf{C}_{D_{\sigma, \sigma^*}}, \mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s) \cdot \sin (\mathbf{C}_{D_1^-} \mathbf{C}_{D_{\sigma, \sigma^*}}, \mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s) \le 0.$$

Next, we define the point $\mathbf{H}_s^-$ to be the intersection between the axis $(\mathbf{C}_{D_{\sigma, \sigma^*}} \mathbf{M}_s)$ and the line $(\mathbf{C}_{D_1^-} \mathbf{C}_{D_2^-})$. The point $\mathbf{H}_s^-$ lies in the segment $(\mathbf{C}_{D_1^-} \mathbf{C}_{D_2^-})$. Therefore, we let $(\alpha_1^-, \alpha_2^-)$ to be the barycentric coordinates of $\mathbf{H}_s^-$ with respect of $(\mathbf{C}_{D_1^-} \mathbf{C}_{D_2^-})$, that is

$$\alpha_1^- = \frac{\|\mathbf{C}_{D_2^-} \mathbf{H}_s^-\|}{\|\mathbf{C}_{D_1^-} \mathbf{C}_{D_2^-}\|} \ge 0, \quad \alpha_2^- = \frac{\|\mathbf{C}_{D_1^-} \mathbf{H}_s^-\|}{\|\mathbf{C}_{D_1^-} \mathbf{C}_{D_2^-}\|} \ge 0, \quad \alpha_1^- + \alpha_2^- = 1.$$
In a symmetric way, we determine the points $C_{D_1^+}$ and $C_{D_2^+}$ for the forward direction (we take the minimum instead of the maximum in (2) and (3)), then $\mathbf{H}_s^+$ will be the intersection between $(C_{D_{\sigma, \sigma^*}} \mathbf{M}_s)$ and $(C_{D_1^+} C_{D_2^+})$, and $(\alpha_1^+, \alpha_2^+)$ its barycentric coordinates such that

$$ \alpha_1^+ = \frac{\|C_{D_2^+} \mathbf{H}_s^+\|}{\|C_{D_1^+} C_{D_2^+}\|} \ge 0, \quad \alpha_2^+ = \frac{\|C_{D_1^+} \mathbf{H}_s^+\|}{\|C_{D_1^+} C_{D_2^+}\|} \ge 0, \quad \alpha_1^+ + \alpha_2^+ = 1. $$

Note that the computation of the coefficients $\alpha_1^\pm$ and $\alpha_2^\pm$ depends only on the mesh and it is done once for all. Finally we compute the densities at points $\mathbf{H}_s^-$ and $\mathbf{H}_s^+$ according to simple weighted means:

$$ \rho_{\mathbf{H}_s^-} = \alpha_1^- \rho_{D_1^-} + \alpha_2^- \rho_{D_2^-}, \quad \rho_{\mathbf{H}_s^+} = \alpha_1^+ \rho_{D_1^+} + \alpha_2^+ \rho_{D_2^+}. $$

Similarly, we reconstruct the density of internal energy $\rho e$ denoted by $(\rho e)_{D_{\sigma, \sigma^*}, s}^{ML}$. The reconstructed internal energy $e_{D_{\sigma, \sigma^*}, s}^{ML}$ is then defined by $(\rho e)_{D_{\sigma, \sigma^*}, s}^{ML} / \rho_{D_{\sigma, \sigma^*}, s}^{ML}$. Finally, the second order mass flux is defined by

$$ \mathcal{F}_{D_{\sigma, \sigma^*}, s}^{ML} = \mathcal{F}_{D_{\sigma, \sigma^*}, s}^{ML,+} + \mathcal{F}_{D_{\sigma, \sigma^*}, s}^{ML,-} $$

with

$$ F_{D_{\sigma, \sigma^*}, s}^{ML,+} = \mathcal{F}^+(\rho_{D_{\sigma, \sigma^*}, s}^{ML}, c_s, u_{D_{\sigma, \sigma^*}, s}) \quad \text{and} \quad \mathcal{F}_{D_{\sigma, \sigma^*}, s}^{ML,-} = \mathcal{F}^-(\rho_{D_{\sigma', \sigma'^*}, s}^{ML}, c_s, u_{D_{\sigma, \sigma^*}, s}), $$

and the internal energy flux is given by

$$ \mathcal{E}_{D_{\sigma, \sigma^*}, s}^{ML} = e_{D_{\sigma, \sigma^*}, s}^{ML} \mathcal{F}_{D_{\sigma, \sigma^*}, s}^{ML,+} + e_{D_{\sigma', \sigma'^*}, s}^{ML} \mathcal{F}_{D_{\sigma, \sigma^*}, s}^{ML,-} $$

# 4 Numerical test on a Riemann problem

We illustrate the method with a comparison between the first order scheme and the MUSCL scheme presented above. We take $\gamma = 1.4$. We consider two 1D Riemann problems. The initial data $\rho, u$ and $p$ are piecewise constant functions with a discontinuity located at $x_0 = 0.5$. The initial constants at left and right of $x_0 = 0.5$ and the final time $T$ are given in Table 4. In the two cases, the solution consists in a left rarefaction, a contact discontinuity and a right shock. The tests are performed with the 2D code on the square $[0, 1] \times [0, 0.1]$, with a triangular mesh with approximately 256 grid points in the direction $x$. Fig. 3 presents the results obtained at time $T$. The MUSCL approach does not reach the 2nd order accuracy since the velocity is not reconstructed but we can already see that the numerical diffusion at the contact discontinuity in the density and internal energy profiles is significantly reduced by the MUSCL approximation (the velocity being constant).
<table><caption>Table 1: Initial data for the Riemann problems</caption><thead><tr><th></th><th>ρ<sub>l</sub></th><th>ρ<sub>r</sub></th><th>u<sub>l</sub></th><th>u<sub>r</sub></th><th>p<sub>l</sub></th><th>p<sub>r</sub></th><th>T</th></tr></thead><tbody><tr><td>Test 1</td><td>1</td><td>0.125</td><td>0</td><td>0</td><td>1</td><td>0.1</td><td>0.25</td></tr><tr><td>Test 2</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1000</td><td>0.1</td><td>0.012</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical results (cutline along x-axis): density and internal energy at time T</figcaption></figure>
# 5 Conclusion

We present a MUSCL approach based on the first order “DDFV-like” scheme introduced in [3]. We focus here on the reconstruction of the mass density and the internal energy, inspired from [4]. This approach gives improved results for 1D contact discontinuity problems. To reach 2nd order, it is needed to equally reconstruct the velocity. The strategy, to be detailed elsewhere, is based on similar reconstructions but the point $H_5^\pm$ can be located on the edge of (primal/dual) cell and the interpolation becomes easier since discrete values are also available at cells vertices (see [1]).

# References

[1] C. Calgaro, E. Creusé, and T. Goudon. An hybrid finite volume–finite element method for variable density incompressible flows. *J. Comput. Phys.*, 227(9):4671–4696, 2008.

[2] L. Gastaldo, R. Herbin, J.-C. Latché, and N. Therme. A muscl-type segregated–explicit staggered scheme for the euler equations. *Comput. Fluids*, 175:91–110, 2018.

[3] T. Goudon, J. Llobell, and S. Minjeaud. An explicit finite volume scheme on staggered grids for the euler equations: Unstructured meshes, stability analysis, and energy conservation. *Int. J. Numer. Methods Fluids*, 94(1):76–119, 2022.

[4] C. Le Touze, A. Murrone, and H. Guillard. Multislope muscl method for general unstructured meshes. *J. Comput. Phys.*, 284:389–418, 2015.