# Simulation of sharp interface multi-material flows involving an arbitrary number of component through an extended five-equation model

Marie Billaud-Friess, S. Kokh

► To cite this version:

Marie Billaud-Friess, S. Kokh. Simulation of sharp interface multi-material flows involving an arbitrary number of component through an extended five-equation model. Journal of Computational Physics, 2014, 1307.8375, pp.1-52. 10.1016/j.jcp.2014.05.012. hal-01007455

HAL Id: hal-01007455

https://hal.science/hal-01007455

Submitted on 6 Jun 2024

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 Attribution 4.0 International License
# Simulation of sharp interface multi-material flows involving an arbitrary number of component through an extended five-equation model

Marie Billaud-Friess, S. Kokh

► To cite this version:

Marie Billaud-Friess, S. Kokh. Simulation of sharp interface multi-material flows involving an arbitrary number of component through an extended five-equation model. Journal of Computational Physics, 2014, 1307.8375, pp.1-52. 10.1016/j.jcp.2014.05.012. hal-01007455

HAL Id: hal-01007455

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

Submitted on 6 Jun 2024

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 Attribution 4.0 International License
Simulation of sharp interface multi-material flows
involving an arbitrary number of components through
an extended five-equation model

Marie Billaud Friess<sup>a,*</sup>, Samuel Kokh<sup>b,c</sup>

<sup>a</sup>LUNAM Université, GeM UMR CNRS 6183, Ecole Centrale Nantes, Université de Nantes, 1 rue de la Noë, BP 92101,44321 Nantes Cedex 3, France

<sup>b</sup>Maison de la Simulation USR 3441, Digiteo Labs - bât. 565 - PC 190, CEA Saclay, 91191 Gif-sur-Yvette, France

<sup>c</sup>DEN/DANS/DM2S/STMF, CEA Saclay, 91191 Gif-sur-Yvette, France

Abstract

In this paper, we present an anti-diffusive method dedicated to the simulation of interface flows on Cartesian grids involving an arbitrary number *m* of compressible components. Our work is two-fold: first, we introduce a *m*-component flow model that generalizes a classic two material five-equation model. In that way, interfaces are localized thanks to color function discontinuities and a pressure equilibrium closure law is used to complete this new model. The resulting model is demonstrated to be hyperbolic under simple assumptions and consistent. Second, we present a discretization strategy for this model relying on an Lagrange-Remap scheme. Here, the projection step involves an anti-dissipative mechanism allowing to prevent numerical diffusion of the material interfaces. The proposed solver is built ensuring consistency and stability properties but also that the sum of the color functions remains equal to one. The resulting scheme is first order accurate and conservative for the mass, momentum, energy and partial masses. Furthermore, the obtained discretization preserves Riemann invariants like pressure and velocity at the interfaces. Finally, validation computations of this numerical method are performed on several tests in one and two dimensions. The accuracy of the method is also compared to results obtained with the upwind Lagrange-Remap scheme.

Keywords: Multi-component flows, Compressible flows, Lagrange-Remap
Anti-diffusive scheme

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: marie.billaud-friess@ec-nantes.fr (Marie Billaud Friess),
samuel.kokh@cea.fr (Samuel Kokh)</footnote>
Simulation of sharp interface multi-material flows
involving an arbitrary number of components through
an extended five-equation model

Marie Billaud Friess$^{a, *}$, Samuel Kokh$^{b, c}$

<sup>a</sup>LUNAM Université, GeM UMR CNRS 6183, Ecole Centrale Nantes, Université de Nantes, 1 rue de la Noë, BP 92101,44321 Nantes Cedex 3, France

<sup>b</sup>Maison de la Simulation USR 3441, Digiteo Labs - bât. 565 - PC 190, CEA Saclay, 91191 Gif-sur-Yvette, France

<sup>c</sup>DEN/DANS/DM2S/STMF, CEA Saclay, 91191 Gif-sur-Yvette, France

Abstract

In this paper, we present an anti-diffusive method dedicated to the simulation of interface flows on Cartesian grids involving an arbitrary number *m* of compressible components. Our work is two-fold: first, we introduce a *m*-component flow model that generalizes a classic two material five-equation model. In that way, interfaces are localized thanks to color function discontinuities and a pressure equilibrium closure law is used to complete this new model. The resulting model is demonstrated to be hyperbolic under simple assumptions and consistent. Second, we present a discretization strategy for this model relying on an Lagrange-Remap scheme. Here, the projection step involves an anti-dissipative mechanism allowing to prevent numerical diffusion of the material interfaces. The proposed solver is built ensuring consistency and stability properties but also that the sum of the color functions remains equal to one. The resulting scheme is first order accurate and conservative for the mass, momentum, energy and partial masses. Furthermore, the obtained discretization preserves Riemann invariants like pressure and velocity at the interfaces. Finally, validation computations of this numerical method are performed on several tests in one and two dimensions. The accuracy of the method is also compared to results obtained with the upwind Lagrange-Remap scheme.

Keywords: Multi-component flows, Compressible flows, Lagrange-Remap
Anti-diffusive scheme

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: marie.billaud-friess@ec-nantes.fr (Marie Billaud Friess),
samuel.kokh@cea.fr (Samuel Kokh)</footnote>
# 1. Introduction

The present paper deals with the simulation of compressible flows that involve *m* distinct materials separated by sharp interfaces on a fixed Cartesian mesh. In our physical framework, there is no velocity jump across the material front and the interfaces are passively advected by the local velocity. We suppose that all diffusive processes are negligible and that each component is equipped with a specific Equation of State (EOS).

The most straightforward simulation strategy consists in considering *m* subdomains with free boundaries in such way that each subdomain is occupied by a single component throughout the computation. Such method is usually called front tracking [14, 29, 35, 63]. It implies implementing a tracking procedure for the interfaces and coping with the boundary conditions for each subdomain thanks to appropriate jump conditions.

We consider here another popular approach that relies on a single-fluid representation of the whole multi-material medium. The material interfaces are represented by loci of discontinuity of the medium physical properties. In our case these discontinuities will produce a switch between the various EOS of each component. Tracking these discontinuities is commonly achieved by introducing additional parameters. For example, in the well-known level set method widely studied in the case $m = 2$ [24, 26, 45, 48, 49, 57, 67, 69, 70] the additional parameter is the signed distance function to the material interface. This function is then evolved thanks to an additional Partial Differential Equation (PDE). Few works propose an extension of the level set method to treat physical situations with $m > 2$ components [69, 70] and ensuring conservation properties like partial masses or momentum conservation can be a complex task.

Another approach relies on introducing discontinuous parameters $Z_k$, $k = 1, \dots, m$, often referred to as *color functions*, such that $Z_k = 1$ (resp. $Z_k = 0$) in regions occupied (resp. not occupied) by the sole fluid $k$. It is possible to use physical parameters like the mass fractions of the volume fractions as color functions. Then, the flow and the interface locations are governed by a system of PDEs formed by the evolution equation of the physical unknown parameters of the components and the *m* additional color functions. Unfortunately, standard discretization techniques like Finite Difference or Finite Volume methods tend to spread the discontinuities that represent the interfaces into several-grid cell wide regions due to numerical diffusion as depicted in figure 1. This raises two issues: first, these transition zones may not be physically relevant as our primary model was designed for modelling sharp interfaces. Second, the diffused interface regions may expand over an important portion of the computational domain. Both issues can be circumvented by implementing interfaces reconstruction techniques. This strategy has been widely used in the case $m = 2$ with the Volume of Fluid (VOF) method [35, 42, 56] and recent works like the Moment of Fluid (MOF) method have successfully addressed the case of
# 1. Introduction

The present paper deals with the simulation of compressible flows that involve $m$ distinct materials separated by sharp interfaces on a fixed Cartesian mesh. In our physical framework, there is no velocity jump across the material front and the interfaces are passively advected by the local velocity. We suppose that all diffusive processes are negligible and that each component is equipped with a specific Equation of State (EOS).

The most straightforward simulation strategy consists in considering $m$ subdomains with free boundaries in such way that each subdomain is occupied by a single component throughout the computation. Such method is usually called front tracking [14, 29, 35, 63]. It implies implementing a tracking procedure for the interfaces and coping with the boundary conditions for each subdomain thanks to appropriate jump conditions.

We consider here another popular approach that relies on a single-fluid representation of the whole multi-material medium. The material interfaces are represented by loci of discontinuity of the medium physical properties. In our case these discontinuities will produce a switch between the various EOS of each component. Tracking these discontinuities is commonly achieved by introducing additional parameters. For example, in the well-known level set method widely studied in the case $m = 2$ [24, 26, 45, 48, 49, 57, 67, 69, 70] the additional parameter is the signed distance function to the material interface. This function is then evolved thanks to an additional Partial Differential Equation (PDE). Few works propose an extension of the level set method to treat physical situations with $m > 2$ components [69, 70] and ensuring conservation properties like partial masses or momentum conservation can be a complex task.

Another approach relies on introducing discontinuous parameters $Z_k$, $k = 1, \dots, m$, often referred to as color functions, such that $Z_k = 1$ (resp. $Z_k = 0$) in regions occupied (resp. not occupied) by the sole fluid $k$. It is possible to use physical parameters like the mass fractions of the volume fractions as color functions. Then, the flow and the interface locations are governed by a system of PDEs formed by the evolution equation of the physical unknown parameters of the components and the $m$ additional color functions. Unfortunately, standard discretization techniques like Finite Difference or Finite Volume methods tend to spread the discontinuities that represent the interfaces into several-grid cell wide regions due to numerical diffusion as depicted in figure 1. This raises two issues: first, these transition zones may not be physically relevant as our primary model was designed for modelling sharp interfaces. Second, the diffused interface regions may expand over an important portion of the computational domain. Both issues can be circumvented by implementing interfaces reconstruction techniques. This strategy has been widely used in the case $m = 2$ with the Volume of Fluid (VOF) method [35, 42, 56] and recent works like the Moment of Fluid (MOF) method have successfully addressed the case of
$m > 2$ [28, 40]. While these approaches provide a true sharp description of the interfaces in the discrete setting, they often remain complex to develop and implement.

An alternate strategy have been used in [3, 4] for the case $m = 2$. In this framework, although the diffused interfaces remain a priori not physically relevant, they are consistent (up to a numerical truncation error) with the target sharp discontinuities. This approach belong to a family of methods that was popularized through several publications in the past years in [1, 2, 3, 4, 38, 39, 46, 53, 54, 59, 60, 61, 63, 68]. Concerning the control of the numerical diffusion produced at the material fronts, a special Lagrange-Remap method was proposed in [41] that spares the difficulty of interfaces reconstruction. This numerical scheme encompasses an anti-diffusive discretization for updating the color functions following the lines of [18, 19, 20, 43].

The present work addresses two questions: first, it proposes an extension for $m \ge 2$ of the isobaric five-equation interface capture model used in [3, 4]. Second, it presents an extension of the anti-diffusive Lagrange-Remap strategy of [41] for this model. Preliminary results have been announced in [9], and we intend to provide here a detailed study of the model within a larger choice of EOS for the material components.

Before closing this introduction, we wish to refer the reader to several other significant efforts in the modelling and the simulation of multi-component flows like [11, 12, 13] by means of a phase-field description of the material interfaces and [16, 43] for a study modelling possible mixture laws for of $m$-component flows. Moreover, we would like to mention that a similar extension of the five-equation model of [3, 4] has been considered by other authors independantly in [29].

The paper is structured as follows. In section 2, we present our $m$-component interface capture model. Under simple hypotheses, we show that the isobaric closure is a consistent definition of a generalized EOS for the whole medium and that, granted simple thermodynamical assumptions, the overall PDE system is hyperbolic. We then detail the construction of the numerical scheme that relies on a Lagrange-Remap splitting as described in section 3. In the section 4, we specify the anti-diffusive strategy used for the remap step, following the work of [36]. Finally we show numerical results in 1D and 2D involving up to $m = 5$ materials in section 5.

## 2. The $m$-component interface flow model

In this part, we propose to extend the two-material five-equation model with isobaric closure [3, 4] to multi-material flows with an arbitrary number of components.
$m > 2$ [28, 40]. While these approaches provide a true sharp description of the interfaces in the discrete setting, they often remain complex to develop and implement.

An alternate strategy have been used in [3, 4] for the case $m = 2$. In this framework, although the diffused interfaces remain a priori not physically relevant, they are consistent (up to a numerical truncation error) with the target sharp discontinuities. This approach belong to a family of methods that was popularized through several publications in the past years in [1, 2, 3, 4, 38, 39, 46, 53, 54, 59, 60, 61, 63, 68]. Concerning the control of the numerical diffusion produced at the material fronts, a special Lagrange-Remap method was proposed in [41] that spares the difficulty of interfaces reconstruction. This numerical scheme encompasses an anti-diffusive discretization for updating the color functions following the lines of [18, 19, 20, 43].

The present work addresses two questions: first, it proposes an extension for $m \ge 2$ of the isobaric five-equation interface capture model used in [3, 4]. Second, it presents an extension of the anti-diffusive Lagrange-Remap strategy of [41] for this model. Preliminary results have been announced in [9], and we intend to provide here a detailed study of the model within a larger choice of EOS for the material components.

Before closing this introduction, we wish to refer the reader to several other significant efforts in the modelling and the simulation of multi-component flows like [11, 12, 13] by means of a phase-field description of the material interfaces and [16, 43] for a study modelling possible mixture laws for of $m$-component flows. Moreover, we would like to mention that a similar extension of the five-equation model of [3, 4] has been considered by other authors independantly in [29].

The paper is structured as follows. In section 2, we present our $m$-component interface capture model. Under simple hypotheses, we show that the isobaric closure is a consistent definition of a generalized EOS for the whole medium and that, granted simple thermodynamical assumptions, the overall PDE system is hyperbolic. We then detail the construction of the numerical scheme that relies on a Lagrange-Remap splitting as described in section 3. In the section 4, we specify the anti-diffusive strategy used for the remap step, following the work of [36]. Finally we show numerical results in 1D and 2D involving up to $m = 5$ materials in section 5.

## 2. The $m$-component interface flow model

In this part, we propose to extend the two-material five-equation model with isobaric closure [3, 4] to multi-material flows with an arbitrary number of components.
<figure><img src="image_2.png" /><figcaption>Figure 1: Medium with multiple materials separated by interfaces. Left: the fluids are sharply separated. Right: Numerical diffusion spreads the interface over few cells. A specific EOS as to be defined in this region.</figcaption></figure>

## 2.1. Evolution equations

Let us consider a medium composed of $m \ge 2$ compressible materials. Each component is equipped with an EOS that is given as a smooth mapping $(\rho_k, p_k) \mapsto e_k(\rho_k, p_k)$, where $\rho_k, p_k, e_k$ denote respectively the density, partial pressure and specific internal energy of the fluid $k = 1, \dots, m$. Let us note $\xi_k = (\partial\rho_k e_k / \partial p_k)_{\rho_k}$, we assume each pure fluid EOS to satisfy the following hypotheses

$$ \left(\frac{\partial p_k}{\partial \rho_k}\right)_{e_k} + \frac{p_k}{\rho_k^2} \left(\frac{\partial p_k}{\partial e_k}\right)_{\rho_k} > 0, \qquad (1) $$

$$ \xi_k = \left( \frac{\partial \rho_k e_k}{\partial p_k} \right)_{\rho_k} > 0. \qquad (2) $$

Relation (1) enables the classic definition of the sound velocity $c_k > 0$ for the fluid $k$ by setting

$$ c_k^2 = \left(\frac{\partial p_k}{\partial \rho_k}\right)_{e_k} + \frac{p_k}{\rho_k^2} \left(\frac{\partial p_k}{\partial e_k}\right)_{\rho_k}. \qquad (3) $$

Relation (2) implies that for a fixed given $\rho_k > 0$, the mapping $p_k \mapsto e_k(\rho_k, p_k)$ is strictly increasing, thus one-to-one.

So as to localize each fluid $k$, we introduce a color function noted $Z_k$ that takes the value $Z_k = 1$ in the regions of the computational domain solely occupied by the fluid $k$ and the value $Z_k = 0$ otherwise. Finally, this quantity is supposed to remain bounded i.e.

$$ Z_k \in [0, 1] \qquad (4) $$

and to satisfy the constraint

$$ \sum_{k=1}^{m} Z_k = 1. \qquad (5) $$

All the flow components share the same velocity $u$. The density $\rho$ and the
<figure><img src="image_2.png" /><figcaption>Figure 1: Medium with multiple materials separated by interfaces. Left: the fluids are sharply separated. Right: Numerical diffusion spreads the interface over few cells. A specific EOS as to be defined in this region.</figcaption></figure>

## 2.1. Evolution equations

Let us consider a medium composed of $m \ge 2$ compressible materials. Each component is equipped with an EOS that is given as a smooth mapping $(\rho_k, p_k) \mapsto e_k(\rho_k, p_k)$, where $\rho_k, p_k, e_k$ denote respectively the density, partial pressure and specific internal energy of the fluid $k = 1, \dots, m$. Let us note $\xi_k = (\partial\rho_k e_k / \partial p_k)_{\rho_k}$, we assume each pure fluid EOS to satisfy the following hypotheses

$$ \left( \frac{\partial p_k}{\partial \rho_k} \right)_{e_k} + \frac{p_k}{\rho_k^2} \left( \frac{\partial p_k}{\partial e_k} \right)_{\rho_k} > 0, \quad (1) $$

$$ \xi_k = \left( \frac{\partial \rho_k e_k}{\partial p_k} \right)_{\rho_k} > 0. \quad (2) $$

Relation (1) enables the classic definition of the sound velocity $c_k > 0$ for the fluid $k$ by setting

$$ c_k^2 = \left( \frac{\partial p_k}{\partial \rho_k} \right)_{e_k} + \frac{p_k}{\rho_k^2} \left( \frac{\partial p_k}{\partial e_k} \right)_{\rho_k}. \quad (3) $$

Relation (2) implies that for a fixed given $\rho_k > 0$, the mapping $p_k \mapsto e_k(\rho_k, p_k)$ is strictly increasing, thus one-to-one.

So as to localize each fluid $k$, we introduce a color function noted $Z_k$ that takes the value $Z_k = 1$ in the regions of the computational domain solely occupied by the fluid $k$ and the value $Z_k = 0$ otherwise. Finally, this quantity is supposed to remain bounded i.e.

$$ Z_k \in [0, 1] \quad (4) $$

and to satisfy the constraint

$$ \sum_{k=1}^{m} Z_k = 1. \quad (5) $$

All the flow components share the same velocity $u$. The density $\rho$ and the
specific internal energy $\rho e$ of the m-component medium are defined by

$$ \rho = \sum_{k=1}^{m} Z_k \rho_k \quad \text{and} \quad \rho e = \sum_{k=1}^{m} Z_k \rho_k e_k. \qquad (6) $$

We note $E = e + ||u||^2/2$ the specific total energy of the medium and we also set

$$ \xi = \sum_{k=1}^{m} Z_k \xi_k. \qquad (7) $$

If we neglect all dissipative effects and volumic source terms, we consider that the flow is governed by the following system of $2m+2$ equations:

$$ \partial_t(\rho u) + \nabla \cdot (\rho u \otimes u) + \nabla p = 0, \qquad (8a) $$

$$ \partial_t(\rho E) + \nabla \cdot (u(\rho E + p)) = 0, \qquad (8b) $$

$$ \partial_t(\rho_k Z_k) + \nabla \cdot (\rho_k Z_k u) = 0, \quad k = 1, \dots, m, \qquad (8c) $$

$$ \partial_t Z_k + u \cdot \nabla Z_k = 0, \quad k = 1, \dots, m. \qquad (8d) $$

The system of equations (8) will be referred to as the m-component model. This system is quasi-conservative since the equation (8d) is not conservative. We shall see in Appendix A that this is not an issue: the non-conservative product $u \cdot \nabla Z_k$ is well-defined and one can recast system (8) into an equivalent fully-conservative formulation.

**Remark 1.** System (8) only involves conservation of partial masses $Z_k \rho_k$. Conservation of global mass $\rho$ can be retrieved by summing over k equations (8c). In addition, using (5) system (8) can be recast using the evolution equation for $\rho$, $\rho u$, $\rho E$ and $\rho_k Z_k$, $Z_k$ for $k = 1, \dots, m-1$. In particular for $m = 2$, we recover the so-called five-equation model of [3, 4].

## 2.2. Multicomponent Pressure Law: Isobaric Closure

Several closure relations that define EOSs for multi-component system have already been examined in the literature for the case of two-component or multi-component medium (see for example [3, 4, 16, 27, 43, 46]). We propose here a direct extension of the so-called isobaric closure law used in [3, 4]. Given the fluid parameters $Z_k$, $\rho_k$ and $\rho e$, we define $p$ by

$$ \rho e = \sum_{k=1}^{m} Z_k \rho_k e_k(\rho_k, p). \qquad (9) $$

The above pressure definition is consistent for a wide range of EOSs in the sense that it is possible to find a unique solution $p$ satisfying (9). We have the following proposition.
specific internal energy ρe of the m-component medium are defined by

$$ \rho = \sum_{k=1}^{m} Z_k \rho_k \quad \text{and} \quad \rho e = \sum_{k=1}^{m} Z_k \rho_k e_k. \qquad (6) $$

We note $E = e + ||\mathbf{u}||^2/2$ the specific total energy of the medium and we also set

$$ \xi = \sum_{k=1}^{m} Z_k \xi_k. \qquad (7) $$

If we neglect all dissipative effects and volumic source terms, we consider that the flow is governed by the following system of $2m+2$ equations:

$$ \partial_t(\rho\mathbf{u}) + \nabla \cdot (\rho\mathbf{u} \otimes \mathbf{u}) + \nabla p = 0, \qquad (8a) $$

$$ \partial_t(\rho E) + \nabla \cdot (\mathbf{u}(\rho E + p)) = 0, \qquad (8b) $$

$$ \partial_t(\rho_k Z_k) + \nabla \cdot (\rho_k Z_k \mathbf{u}) = 0, \quad k = 1, \dots, m, \qquad (8c) $$

$$ \partial_t Z_k + \mathbf{u} \cdot \nabla Z_k = 0, \quad k = 1, \dots, m. \qquad (8d) $$

The system of equations (8) will be referred to as the m-component model. This system is quasi-conservative since the equation (8d) is not conservative. We shall see in Appendix A that this is not an issue: the non-conservative product $\mathbf{u} \cdot \nabla Z_k$ is well-defined and one can recast system (8) into an equivalent fully-conservative formulation.

**Remark 1.** System (8) only involves conservation of partial masses $Z_k \rho_k$. Conservation of global mass $\rho$ can be retrieved by summing over k equations (8c). In addition, using (5) system (8) can be recast using the evolution equation for $\rho$, $\rho\mathbf{u}$, $\rho E$ and $\rho_k Z_k$, $Z_k$ for $k = 1, \dots, m-1$. In particular for $m = 2$, we recover the so-called five-equation model of [3, 4].

## 2.2. Multicomponent Pressure Law: Isobaric Closure

Several closure relations that define EOSs for multi-component system have already been examined in the literature for the case of two-component or multi-component medium (see for example [3, 4, 16, 27, 43, 46]). We propose here a direct extension of the so-called isobaric closure law used in [3, 4]. Given the fluid parameters $Z_k$, $\rho_k$ and $\rho e$, we define $p$ by

$$ \rho e = \sum_{k=1}^{m} Z_k \rho_k e_k(\rho_k, p). \qquad (9) $$

The above pressure definition is consistent for a wide range of EOSs in the sense that it is possible to find a unique solution $p$ satisfying (9). We have the following proposition.
**Proposition 1.** Suppose that for $k = 1, \dots, m$ and any fixed $\rho_k > 0$

$$ \lim_{p_k \to 0} \rho_k e_k(\rho_k, p_k) = 0, \qquad (10) $$

$$ \lim_{p_k \to +\infty} \rho_k e_k(\rho_k, p_k) = +\infty, \qquad (11) $$

then under the assumption (2), there exists a single $p$ verifying (8).

*Proof.* For the given values $\mathcal{Z}_k, \rho_k > 0$, $k = 1, \dots, m$ and $\rho e > 0$, we consider the mapping $\Xi : p \mapsto \sum_{k=1}^m \mathcal{Z}_k \rho_k e_k(\rho_k, p) - \rho e$. The assumptions (10) and (11) imply respectively that $\lim_{p \to 0} \Xi = -\rho e < 0$ and $\lim_{p \to +\infty} \Xi = +\infty$. The theorem of intermediate value then states that there exists $\bar{p}$ such that $\Xi(\bar{p}) = 0$. According to (2), $d\Xi/dp = \sum_{k=1}^m \mathcal{Z}_k \xi_k > 0$ and thus $\Xi$ is one-to-one. This implies that $\bar{p}$ is unique. $\square$

The proposition 1 is valid for a wide range of EOSs since it does not require any specific analytical form for the EOSs and enables the use of tabulated data. Assumptions (10)-(11) are sufficient but not necessary conditions for obtaining a proper definition of the pressure $p$. Indeed, for a $m$-component flow involving $m$ Mie-Gruneisen materials (10) and (11) are not necessarily verified, but simple calculations allow to properly define $p$ as in proposition 2.

**Proposition 2.** Suppose that for each $k = 1, \dots, m$ the component $k$ is a Mie-Gruneisen material whose EOS reads

$$ (\rho_k, p_k) \mapsto \rho_k e_k = \rho_k e_k^{\text{ref}}(\rho_k) + \frac{p_k - p_k^{\text{ref}}(\rho_k)}{\Gamma_k(\rho_k)}, $$

$$ (\rho_k, e_k) \mapsto p_k = p_k^{\text{ref}}(\rho_k) + \Gamma_k(\rho_k) \rho_k [e_k - e_k^{\text{ref}}(\rho_k)], $$

where $\rho_k \mapsto \Gamma_k(\rho_k)$, $\rho_k \mapsto e_k^{\text{ref}}(\rho_k)$ and $\rho_k \mapsto p_k^{\text{ref}}(\rho_k)$ are real valued functions. Furthermore, we have

$$ \xi_k = \frac{1}{\Gamma_k}, $$

$$ c_k^2 = \frac{dp_k^{\text{ref}}}{d\rho_k} + \frac{p_k - p_k^{\text{ref}}}{\Gamma_k} \frac{d\Gamma_k}{d\rho_k} - \rho_k \Gamma_k \frac{de_k^{\text{ref}}}{d\rho_k} + (\Gamma_k + 1) \frac{p_k - p_k^{\text{ref}}}{\rho_k}. $$

For the $m$-component system (8), the unique solution pressure $p$ defined by the isobaric closure (9) is

$$ p = \left[ \rho e - \sum_{k=1}^{m} \mathcal{Z}_k \rho_k e_k^{\text{ref}} + \sum_{k=1}^{m} \mathcal{Z}_k \frac{p_k^{\text{ref}}}{\Gamma_k} \right] \times \left[ \sum_{k=1}^{m} \frac{\mathcal{Z}_k}{\Gamma_k} \right]^{-1}. $$

Two very common choices of EOS are often found in the literature: the Stiffened Gas and the Perfect Gas which belong to the class of Mie-Gruneisen materials with

$$ \text{Perfect Gas:} \quad \Gamma_k = \gamma_k - 1, \quad p_k^{\text{ref}} = 0, \quad e_k^{\text{ref}} = 0, \qquad (12) $$

$$ \text{Stiffened Gas:} \quad \Gamma_k = \gamma_k - 1, \quad p_k^{\text{ref}} = -\gamma_k \pi_k, \quad e_k^{\text{ref}} = 0, \qquad (13) $$
**Proposition 1.** Suppose that for $k = 1, \dots, m$ and any fixed $\rho_k > 0$

$$ \lim_{p_k \to 0} \rho_k e_k(\rho_k, p_k) = 0, \qquad (10) $$

$$ \lim_{p_k \to +\infty} \rho_k e_k(\rho_k, p_k) = +\infty, \qquad (11) $$

then under the assumption (2), there exists a single p verifying (8).

*Proof.* For the given values $\mathcal{Z}_k, \rho_k > 0$, $k = 1, \dots, m$ and $\rho e > 0$, we consider the mapping $\Xi : p \mapsto \sum_{k=1}^m \mathcal{Z}_k \rho_k e_k(\rho_k, p) - \rho e$. The assumptions (10) and (11) imply respectively that $\lim_{p \to 0} \Xi = -\rho e < 0$ and $\lim_{p \to +\infty} \Xi = +\infty$. The theorem of intermediate value then states that there exists $\bar{p}$ such that $\Xi(\bar{p}) = 0$. According to (2), $d\Xi/dp = \sum_{k=1}^m \mathcal{Z}_k \xi_k > 0$ and thus $\Xi$ is one-to-one. This implies that $\bar{p}$ is unique. $\square$

The proposition 1 is valid for a wide range of EOSs since it does not require any specific analytical form for the EOSs and enables the use of tabulated data. Assumptions (10)-(11) are sufficient but not necessary conditions for obtaining a proper definition of the pressure $p$. Indeed, for a $m$-component flow involving $m$ Mie-Gruneisen materials (10) and (11) are not necessarily verified, but simple calculations allow to properly define $p$ as in proposition 2.

**Proposition 2.** Suppose that for each $k = 1, \dots, m$ the component $k$ is a Mie-Gruneisen material whose EOS reads

$$ (\rho_k, p_k) \mapsto \rho_k e_k = \rho_k e_k^{\text{ref}}(\rho_k) + \frac{p_k - p_k^{\text{ref}}(\rho_k)}{\Gamma_k(\rho_k)}, $$

$$ (\rho_k, e_k) \mapsto p_k = p_k^{\text{ref}}(\rho_k) + \Gamma_k(\rho_k) \rho_k [e_k - e_k^{\text{ref}}(\rho_k)], $$

where $\rho_k \mapsto \Gamma_k(\rho_k)$, $\rho_k \mapsto e_k^{\text{ref}}(\rho_k)$ and $\rho_k \mapsto p_k^{\text{ref}}(\rho_k)$ are real valued functions. Furthermore, we have

$$ \xi_k = \frac{1}{\Gamma_k}, $$

$$ c_k^2 = \frac{dp_k^{\text{ref}}}{d\rho_k} + \frac{p_k - p_k^{\text{ref}}}{\Gamma_k} \frac{d\Gamma_k}{d\rho_k} - \rho_k \Gamma_k \frac{de_k^{\text{ref}}}{d\rho_k} + (\Gamma_k + 1) \frac{p_k - p_k^{\text{ref}}}{\rho_k}. $$

For the $m$-component system (8), the unique solution pressure $p$ defined by the isobaric closure (9) is

$$ p = \left[ \rho e - \sum_{k=1}^{m} \mathcal{Z}_k \rho_k e_k^{\text{ref}} + \sum_{k=1}^{m} \mathcal{Z}_k \frac{p_k^{\text{ref}}}{\Gamma_k} \right] \times \left[ \sum_{k=1}^{m} \frac{\mathcal{Z}_k}{\Gamma_k} \right]^{-1}. $$

Two very common choices of EOS are often found in the literature: the Stiffened Gas and the Perfect Gas which belong to the class of Mie-Gruneisen materials with

$$ \text{Perfect Gas:} \quad \Gamma_k = \gamma_k - 1, \quad p_k^{\text{ref}} = 0, \quad e_k^{\text{ref}} = 0, \qquad (12) $$

$$ \text{Stiffened Gas:} \quad \Gamma_k = \gamma_k - 1, \quad p_k^{\text{ref}} = -\gamma_k \pi_k, \quad e_k^{\text{ref}} = 0, \qquad (13) $$
where $\gamma_k$ and $\pi_k$ are constants.

### 2.3. Hyperbolicity and eigenstructure

We now briefly study the hyperbolicity of the m-component model with isobaric closure (more details are given in Appendix A). System (8) in one space dimension for smooth solutions formulated for the primitive variables $\mathbf{V} = (\rho, u, p, Z_1, \dots, Z_m, Y_1, \dots, Y_{m-1})^T$ reads

$$ \partial_t \mathbf{V} + \mathbf{A}(\mathbf{V})\partial_x \mathbf{V} = 0, \quad (14) $$

where $\mathbf{A}(\mathbf{V})$ is a $(2m+2) \times (2m+2)$ matrix defined by

$$ \mathbf{A}(\mathbf{V}) = \begin{pmatrix} \mathbf{A}_1 & 0 \\ 0 & u\mathbf{I}_{2m-1} \end{pmatrix}, \quad \mathbf{A}_1(\mathbf{V}) = \begin{pmatrix} u & \rho & 0 \\ 0 & u & 1/\rho \\ 0 & \rho c^2 & u \end{pmatrix} $$

and $\mathbf{I}_{2m-1}$ is the $(2m-1) \times (2m-1)$ identity matrix.

**Proposition 3.** The sound velocity $c$ of the m-component model defined by relation

$$ \rho\xi c^2 = \sum_k Z_k \rho_k \xi_k c_k^2. \quad (15) $$

is a positive real number and the matrix $\mathbf{A}(\mathbf{V})$ is diagonalizable. The eigenvalues of $A$ are $\lambda_1 = u - c, \lambda_2 = u + c, \lambda_3 = \dots = \lambda_{2m+2} = u$. The fields associated with the eigenvalues $\lambda_1$ and $\lambda_2$ are genuinely nonlinear, those associated with $\lambda_k$ are linearly degenerate for $k = 3, \dots, 2m+2$. Therefore, the m-component system (8) is hyperbolic.

## 3. Lagrange-Remap scheme

In this section, we present a general quasi-conservative finite-volume scheme for approximating the solution of system (8). Here, we use a two-step Lagrange-Remap strategy [21, 31]: the first step accounts for acoustic effects while the second step deals with material transport of the components. For the sake of simplicity, we present our numerical scheme in the case of one-dimensional problems. The extension of the scheme to multi-dimensional problems when using a Cartesian discretization of the computational domain is achieved in this work thanks to a simple dimensional splitting detailed in section 3.4. Here, we consider a regular discretization of the real line into a set of cells $([x_{i-1/2}, x_{i+1/2}])_{i \in \mathbb{Z}}$, where $x_{i+1/2} = i\Delta x$ for $i \in \mathbb{Z}$ and $\Delta x > 0$ is the space step. Let $x_i = \frac{x_{i+1/2}+x_{i-1/2}}{2}$ be the center of the cell $i$, for $i \in \mathbb{Z}$. Let $(t^n)_{n \in \mathbb{N}}$ be the regular sequence of instants $t^n = n\Delta t$, where $\Delta t$ is the time step. If $(x, t) \mapsto a(x, t)$ is a fluid parameter, we consider the discrete variable $a_i^n$ to be an approximation of

$$ \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} a(x, t^n) dx, \quad i \in \mathbb{Z}, n \in \mathbb{N}. $$
where $\gamma_k$ and $\pi_k$ are constants.

### 2.3. Hyperbolicity and eigenstructure

We now briefly study the hyperbolicity of the m-component model with isobaric closure (more details are given in Appendix A). System (8) in one space dimension for smooth solutions formulated for the primitive variables $\mathbf{V} = (\rho, u, p, Z_1, \dots, Z_m, Y_1, \dots, Y_{m-1})^T$ reads

$$ \partial_t \mathbf{V} + \mathbf{A}(\mathbf{V})\partial_x \mathbf{V} = 0, \qquad (14) $$

where $\mathbf{A}(\mathbf{V})$ is a $(2m+2) \times (2m+2)$ matrix defined by

$$ \mathbf{A}(\mathbf{V}) = \begin{pmatrix} \mathbf{A}_1 & 0 \\ 0 & u\mathbf{I}_{2m-1} \end{pmatrix}, \quad \mathbf{A}_1(\mathbf{V}) = \begin{pmatrix} u & \rho & 0 \\ 0 & u & 1/\rho \\ 0 & \rho c^2 & u \end{pmatrix} $$

and $\mathbf{I}_{2m-1}$ is the $(2m-1) \times (2m-1)$ identity matrix.

**Proposition 3.** The sound velocity $c$ of the m-component model defined by relation

$$ \rho\xi c^2 = \sum_k Z_k \rho_k \xi_k c_k^2. \qquad (15) $$

is a positive real number and the matrix $\mathbf{A}(\mathbf{V})$ is diagonalizable. The eigenvalues of $A$ are $\lambda_1 = u - c, \lambda_2 = u + c, \lambda_3 = \dots = \lambda_{2m+2} = u$. The fields associated with the eigenvalues $\lambda_1$ and $\lambda_2$ are genuinely nonlinear, those associated with $\lambda_k$ are linearly degenerate for $k = 3, \dots, 2m+2$. Therefore, the m-component system (8) is hyperbolic.

## 3. Lagrange-Remap scheme

In this section, we present a general quasi-conservative finite-volume scheme for approximating the solution of system (8). Here, we use a two-step Lagrange-Remap strategy [21, 31]: the first step accounts for acoustic effects while the second step deals with material transport of the components. For the sake of simplicity, we present our numerical scheme in the case of one-dimensional problems. The extension of the scheme to multi-dimensional problems when using a Cartesian discretization of the computational domain is achieved in this work thanks to a simple dimensional splitting detailed in section 3.4. Here, we consider a regular discretization of the real line into a set of cells $([x_{i-1/2}, x_{i+1/2}])_{i \in \mathbb{Z}}$, where $x_{i+1/2} = i\Delta x$ for $i \in \mathbb{Z}$ and $\Delta x > 0$ is the space step. Let $x_i = \frac{x_{i+1/2}+x_{i-1/2}}{2}$ be the center of the cell $i$, for $i \in \mathbb{Z}$. Let $(t^n)_{n \in \mathbb{N}}$ be the regular sequence of instants $t^n = n\Delta t$, where $\Delta t$ is the time step. If $(x, t) \mapsto a(x, t)$ is a fluid parameter, we consider the discrete variable $a_i^n$ to be an approximation of

$$ \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} a(x, t^n) dx, \quad i \in \mathbb{Z}, n \in \mathbb{N}. $$
Discrete values located at the cell boundary $x = x_{i+1/2}$ are denoted using the subscript $(i + 1/2)$; for the sake of readability discrete variables after the Lagrangian step are designed by $\tilde{a}$, such as $\tilde{a}_{i+1/2}$ or $\tilde{a}_i$. Finally, we denote by $[[a]]_i = a_{i+1/2} - a_{i-1/2}$ the flux difference in a cell $i$.

### 3.1. Lagrangian step

The Lagrangian step consists in solving the Euler equations in Lagrangian coordinates (see Appendix A. or for more details [21, 31]). Following [20, 41], we perform this task using the acoustic scheme [21]. The acoustic scheme can be obtained using a Suliciu-type relaxation approach [10, 14, 55] and it can be considered as a particular HLLC [6, 65] solver for the system written in Lagrangian coordinates. Let us mention that other approximate Riemann solvers may be used to achieve the approximation of the Lagrangian step like for example a Roe-type linearization [51] or a Rusanov scheme [52]. If we note $L_i = 1 + \frac{\Delta t}{\Delta x}[[u^n]]_i$, this numerical scheme provides the update relations

$$ 
\left\{
\begin{aligned}
L_i(\tilde{\rho}u)_i &= (\rho u)_i^n - \frac{\Delta t}{\Delta x}[[p^n]]_i, \\
L_i(\tilde{\rho}E)_i &= (\rho E)_i^n - \frac{\Delta t}{\Delta x}[[p^n u^n]]_i, \\
L_i(\tilde{\rho}_k \tilde{Z}_k)_i &= (\rho_k Z_k)_i^n, \\
\tilde{Z}_{k,i} &= Z_{k,i}^n,
\end{aligned}
\right.
\qquad (16)
$$

where the fluxes $p_{i+1/2}^n$ and $u_{i+1/2}^n$ are defined by

$$ 
\left\{
\begin{aligned}
p_{i+1/2}^n &= \frac{p_i^n + p_{i+1}^n}{2} - \frac{1}{2}(\rho c)_i^n (u_{i+1}^n - u_i^n), \\
u_{i+1/2}^n &= \frac{u_i^n + u_{i+1}^n}{2} - \frac{1}{2}(\rho c)_i^n (p_{i+1}^n - p_i^n), \\
(\rho c)_i^n &= \sqrt{\max((\rho c^2)_i^n, (\rho c^2)_{i+1}^n) \min(\rho_i^n, \rho_{i+1}^n)},
\end{aligned}
\right.
\qquad (17)
$$

and $(\rho c^2)_i^n$ is computed according to (A.7).

As in [20, 41], $\Delta t$ is chosen so that it verifies the following Courant-Friedrichs-Lewy (CFL) condition

$$ \frac{\Delta t}{\Delta x} \max_{i \in \mathbb{Z}} \left( |u_{i+1/2}^n|, (\rho c)_{i+1/2}^n / \min(\rho_i^n, \rho_{i+1}^n) \right) \le C_{\text{CFL}}, \qquad (18) $$

where $0 \le C_{\text{CFL}} \le 1$. The condition (18) ensures stability for the Lagrange step (16).

### 3.2. Remap step

The remap step accounts for the material transport of the fluid [21, 31]. During this phase, the updated Lagrangian variables issued from (16) are projected
Discrete values located at the cell boundary $x = x_{i+1/2}$ are denoted using the subscript $(i + 1/2)$; for the sake of readability discrete variables after the Lagrangian step are designed by $\tilde{\cdot}$, such as $\tilde{a}_{i+1/2}$ or $\tilde{a}_i$. Finally, we denote by $[[a]]_i = a_{i+1/2} - a_{i-1/2}$ the flux difference in a cell $i$.

### 3.1. Lagrangian step

The Lagrangian step consists in solving the Euler equations in Lagrangian coordinates (see Appendix A. or for more details [21, 31]). Following [20, 41], we perform this task using the acoustic scheme [21]. The acoustic scheme can be obtained using a Suliciu-type relaxation approach [10, 14, 55] and it can be considered as a particular HLLC [6, 65] solver for the system written in Lagrangian coordinates. Let us mention that other approximate Riemann solvers may be used to achieve the approximation of the Lagrangian step like for example a Roe-type linearization [51] or a Rusanov scheme [52]. If we note $L_i = 1 + \frac{\Delta t}{\Delta x}[[u^n]]_i$, this numerical scheme provides the update relations

$$ 
\left\{
\begin{aligned}
L_i(\tilde{\rho}u)_i &= (\rho u)_i^n - \frac{\Delta t}{\Delta x}[[p^n]]_i, \\
L_i(\tilde{\rho}E)_i &= (\rho E)_i^n - \frac{\Delta t}{\Delta x}[[p^n u^n]]_i, \\
L_i(\tilde{\rho}_k \tilde{Z}_k)_i &= (\rho_k Z_k)_i^n, \\
\tilde{Z}_{k,i} &= Z_{k,i}^n,
\end{aligned}
\right.
\qquad (16)
$$

where the fluxes $p_{i+1/2}^n$ and $u_{i+1/2}^n$ are defined by

$$ 
\left\{
\begin{aligned}
p_{i+1/2}^n &= \frac{p_i^n + p_{i+1}^n}{2} - \frac{1}{2}(\rho c)_{i+1/2}^n (u_{i+1}^n - u_i^n), \\
u_{i+1/2}^n &= \frac{u_i^n + u_{i+1}^n}{2} - \frac{1}{2}(\rho c)_{i+1/2}^n (p_{i+1}^n - p_i^n), \\
(\rho c)_{i+1/2}^n &= \sqrt{\max((\rho c^2)_i^n, (\rho c^2)_{i+1}^n) \min(\rho_i^n, \rho_{i+1}^n)},
\end{aligned}
\right.
\qquad (17)
$$

and $(\rho c^2)_i^n$ is computed according to (A.7).

As in [20, 41], $\Delta t$ is chosen so that it verifies the following Courant-Friedrichs-Lewy (CFL) condition

$$ \frac{\Delta t}{\Delta x} \max_{i \in \mathbb{Z}} \left( |u_{i+1/2}^n|, (\rho c)_{i+1/2}^n / \min(\rho_i^n, \rho_{i+1}^n) \right) \le C_{\text{CFL}}, \qquad (18) $$

where $0 \le C_{\text{CFL}} \le 1$. The condition (18) ensures stability for the Lagrange step (16).

### 3.2. Remap step

The remap step accounts for the material transport of the fluid [21, 31]. During this phase, the updated Lagrangian variables issued from (16) are projected
onto the Eulerian mesh. Following the classic lines of [31] and as suggested by [8, 20, 41] we achieve this task by setting

$$ \mathbf{W}_i^{n+1} = \tilde{\mathbf{W}}_i - \frac{\Delta t}{\Delta x} [u^n \tilde{\mathbf{W}}_i] + \frac{\Delta t}{\Delta x} [u^n]_i \tilde{\mathbf{W}}_i, \quad (19) $$

where $\mathbf{W} = [\rho u, \rho E, \rho_1 Z_1, \dots, \rho_m Z_m, Z_1, \dots, Z_m]^T$. Here, the fluxes $[u^n]_i$ are computed using $u_{i+1/2}^n$ and $u_{i-1/2}^n$ given by (17) and the fluxes $\tilde{\rho}_{i+1/2}$, $(\tilde{\rho}u)_{i+1/2}$, $(\tilde{\rho}E)_{i+1/2}$, $(\tilde{\rho}_k \tilde{Z}_k)_{i+1/2}$ and $\tilde{Z}_{k,i+1/2}$ are defined by

$$ \left\{ \begin{aligned} \tilde{\rho}_{i+1/2} &= \sum_{k=1}^{m} \tilde{Z}_{k,i+1/2} \tilde{\rho}_{k,i+1/2}, \\ (\tilde{\rho}e)_{i+1/2} &= \sum_{k=1}^{m} \tilde{Z}_{k,i+1/2} (\tilde{\rho}_k e_k)_{i+1/2}, \\ (\tilde{\rho}u)_{i+1/2} &= \tilde{\rho}_{i+1/2} \tilde{u}_{i+1/2}, \\ (\tilde{\rho}E)_{i+1/2} &= (\tilde{\rho}e)_{i+1/2} + \tilde{\rho}_{i+1/2} (\tilde{u}_{i+1/2})^2/2, \\ (\tilde{\rho}_k \tilde{Z}_k)_{i+1/2} &= \tilde{Z}_{k,i+1/2} \tilde{\rho}_{k,i+1/2}. \end{aligned} \right. \quad (20) $$

Then, we choose to define all the phasic values $\tilde{\rho}_{k,i+1/2}$, $\tilde{e}_{k,i+1/2}$ and the velocity $\tilde{u}_{i+1/2}$ by taking their upwind value with respect to the velocity $u_{i+1/2}^n$, namely

$$ (\tilde{\rho}_k, \tilde{\rho}_k e_k, \tilde{u})_{i+1/2} = \begin{cases} (\tilde{\rho}_k, \tilde{\rho}_k e_k, \tilde{u})_i, & \text{if } u_{i+1/2}^n > 0, \\ (\tilde{\rho}_k, \tilde{\rho}_k e_k, \tilde{u})_{i+1}, & \text{if } u_{i+1/2}^n \le 0. \end{cases} \quad (21) $$

At this step the color function fluxes $\tilde{Z}_{k,i+1/2}$ can be computed in several manners. A very common way is to use an upwind flux

$$ \tilde{Z}_{k,i+1/2} = Z_{k,i+1/2}^{up} = \begin{cases} Z_{k,i}^n & \text{if } u_{i+1/2}^n > 0, \\ Z_{k,i+1}^n & \text{otherwise.} \end{cases} \quad (22) $$

Nevertheless, despite its simplicity and satisfying good properties: unit constraint (5), stability, it is well known that the resulting scheme is very diffusive in particular for the transport of $Z_k$ (see section 5). To overcome this curse, high order schemes as Minmod limiter, Superbee limiter or Ultrabee limiter schemes can be envisaged. Nevertheless, it was demonstrated [36] that such schemes present defects for the passive transport of $m > 2$ scalar functions since they can not respect at same time stability and unit constraint. Here, we propose and detail an anti-diffusive procedure to appropriately define the fluxes $\tilde{Z}_{k,i+1/2}$ in section 4 insuring stability and unit constraint (5), and that allows accurate computation of $Z_k$ during remap step.

**Remark 2.** The cell-centered quantities $\tilde{\rho}_{k,i}$ and $\tilde{e}_{k,i}$ are defined by $\tilde{\rho}_{k,i} = (\tilde{\rho}_k \tilde{Z}_k)_i / \tilde{Z}_{k,i}$, $\tilde{e}_{k,i} = e_k(\tilde{\rho}_{k,i}, \tilde{p}_i)$. If $|\tilde{Z}_{k,i}|$ is lower than a user-defined threshold value, we set $\tilde{\rho}_{k,i} = \tilde{e}_{k,i} = 0$ while evaluating (20).
onto the Eulerian mesh. Following the classic lines of [31] and as suggested by [8, 20, 41] we achieve this task by setting

$$ \mathbf{W}_i^{n+1} = \tilde{\mathbf{W}}_i - \frac{\Delta t}{\Delta x} [u^n \tilde{\mathbf{W}}_i] + \frac{\Delta t}{\Delta x} [u^n]_i \tilde{\mathbf{W}}_i, \quad (19) $$

where $\mathbf{W} = [\rho u, \rho E, \rho_1 Z_1, \dots, \rho_m Z_m, Z_1, \dots, Z_m]^T$. Here, the fluxes $[u^n]_i$ are computed using $u_{i+1/2}^n$ and $u_{i-1/2}^n$ given by (17) and the fluxes $\tilde{\rho}_{i+1/2}, (\tilde{\rho}u)_{i+1/2}, (\tilde{\rho}E)_{i+1/2}, (\tilde{\rho}_k \tilde{Z}_k)_{i+1/2}$ and $\tilde{Z}_{k,i+1/2}$ are defined by

$$ \left\{ \begin{aligned} \tilde{\rho}_{i+1/2} &= \sum_{k=1}^{m} \tilde{Z}_{k,i+1/2} \tilde{\rho}_{k,i+1/2}, \\ (\tilde{\rho}e)_{i+1/2} &= \sum_{k=1}^{m} \tilde{Z}_{k,i+1/2} (\tilde{\rho}_k e_k)_{i+1/2}, \\ (\tilde{\rho}u)_{i+1/2} &= \tilde{\rho}_{i+1/2} \tilde{u}_{i+1/2}, \\ (\tilde{\rho}E)_{i+1/2} &= (\tilde{\rho}e)_{i+1/2} + \tilde{\rho}_{i+1/2} (\tilde{u}_{i+1/2})^2/2, \\ (\tilde{\rho}_k \tilde{Z}_k)_{i+1/2} &= \tilde{Z}_{k,i+1/2} \tilde{\rho}_{k,i+1/2}. \end{aligned} \right. \quad (20) $$

Then, we choose to define all the phasic values $\tilde{\rho}_{k,i+1/2}, \tilde{e}_{k,i+1/2}$ and the velocity $\tilde{u}_{i+1/2}$ by taking their upwind value with respect to the velocity $u_{i+1/2}^n$, namely

$$ (\tilde{\rho}_k, \tilde{\rho}_k e_k, \tilde{u})_{i+1/2} = \begin{cases} (\tilde{\rho}_k, \tilde{\rho}_k e_k, \tilde{u})_i, & \text{if } u_{i+1/2}^n > 0, \\ (\tilde{\rho}_k, \tilde{\rho}_k e_k, \tilde{u})_{i+1}, & \text{if } u_{i+1/2}^n \le 0. \end{cases} \quad (21) $$

At this step the color function fluxes $\tilde{Z}_{k,i+1/2}$ can be computed in several manners. A very common way is to use an upwind flux

$$ \tilde{Z}_{k,i+1/2} = Z_{k,i+1/2}^{up} = \begin{cases} Z_{k,i}^n & \text{if } u_{i+1/2}^n > 0, \\ Z_{k,i+1}^n & \text{otherwise.} \end{cases} \quad (22) $$

Nevertheless, despite its simplicity and satisfying good properties: unit constraint (5), stability, it is well known that the resulting scheme is very diffusive in particular for the transport of $Z_k$ (see section 5). To overcome this curse, high order schemes as Minmod limiter, Superbee limiter or Ultrabee limiter schemes can be envisaged. Nevertheless, it was demonstrated [36] that such schemes present defects for the passive transport of $m > 2$ scalar functions since they can not respect at same time stability and unit constraint. Here, we propose and detail an anti-diffusive procedure to appropriately define the fluxes $\tilde{Z}_{k,i+1/2}$ in section 4 insuring stability and unit constraint (5), and that allows accurate computation of $Z_k$ during remap step.

**Remark 2.** The cell-centered quantities $\tilde{\rho}_{k,i}$ and $\tilde{e}_{k,i}$ are defined by $\tilde{\rho}_{k,i} = (\tilde{\rho}_k \tilde{Z}_k)_i / \tilde{Z}_{k,i}$, $\tilde{e}_{k,i} = e_k(\tilde{\rho}_{k,i}, \tilde{p}_i)$. If $|\tilde{Z}_{k,i}|$ is lower than a user-defined threshold value, we set $\tilde{\rho}_{k,i} = \tilde{e}_{k,i} = 0$ while evaluating (20).
### 3.3. Overall algorithm

Let us summarize the overall algorithm detailed in the previous sections. For updating the variables at each instant $t^n$ to their values at $t^{n+1}$, we proceed following the steps of the algorithm 1.

<caption>**Algorithm 1** Lagrange-Remap procedure</caption>

**I - Lagrangian step:**
1: Compute the acoustic fluxes $p_{i+1/2}^n, u_{i+1/2}^n$ thanks to (17).
2: Compute the time step $\Delta t$ according to the CFL constraint (18).
3: Update $\tilde{\rho}u, \tilde{\rho}E, \tilde{\rho}_k\tilde{Z}_k, \tilde{Z}_k$ with (16).

**II - Remap step:**
1: Compute the fluxes $\tilde{Z}_{k,i+1/2}$ for each interface $i + 1/2$.
2: Compute the remap fluxes of the conservative quantities $\tilde{\rho}_{i+1/2}, \tilde{\rho}u_{i+1/2}, \tilde{\rho}E_{i+1/2}, (\tilde{\rho}_k\tilde{Z}_k)_{i+1/2}$ with (20).
3: Update $(\rho u)^{n+1}, (\rho E)^{n+1}, (\rho_k Z_k)^{n+1}, Z_k^{n+1}$ with (19).

Combining the Lagrange (16) and Remap (19) steps, we obtain the following quasi-conservative scheme

$$ \left\{ \begin{aligned} (\rho u)_i^{n+1} &= (\rho u)_i^n & - \frac{\Delta t}{\Delta x} [u^n \tilde{\rho} u]_i & - \frac{\Delta t}{\Delta x} [p^n]_i, \\ (\rho E)_i^{n+1} &= (\rho E)_i^n & - \frac{\Delta t}{\Delta x} [u^n \tilde{\rho} E]_i & - \frac{\Delta t}{\Delta x} [p^n u^n]_i, \\ (\rho_k Z_k)_i^{n+1} &= (\rho_k Z_k)_i^n & - \frac{\Delta t}{\Delta x} [u^n \tilde{\rho}_k Z_k]_i, \\ Z_{k,i}^{n+1} &= Z_{k,i}^n & - \frac{\Delta t}{\Delta x} [u^n \tilde{Z}_k]_i & + \frac{\Delta t}{\Delta x} [u^n]_i \tilde{Z}_{k,i}. \end{aligned} \right. \quad (23) $$

As in [41], the overall discretization strategy is conservative with respect to the momentum, total energy and partial masses.

### 3.4. Extension to the multi-dimensional case

Without loss of generality, we consider a two-dimensional problem that is discretized over a Cartesian grid. Noting $\mathbf{u} = (u_1, u_2)^T$ the velocity, the m-component system (8) reads

$$ \left\{ \begin{aligned} \partial_t \mathbf{U} + \partial_{x_1} \mathbf{F}_1(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) + \partial_{x_2} \mathbf{F}_2(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) &= 0, \\ \partial_t \mathcal{Z}_k + u_1 \partial_{x_1} \mathcal{Z}_k + u_2 \partial_{x_2} \mathcal{Z}_k &= 0, \quad k = 1, \dots, m, \end{aligned} \right. \quad (24) $$

where $\mathbf{U} = [\rho u_1, \rho u_2, \rho E, \rho_1 \mathcal{Z}_1, \dots, \rho_m \mathcal{Z}_m]^T$, and

$$ \begin{aligned} \mathbf{F}_1(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) &= [\rho u_1^2 + p, \rho u_1 u_2, (\rho E + p)u_1, \rho_1 \mathcal{Z}_1 u_1, \dots, \rho_m \mathcal{Z}_m u_1]^T, \\ \mathbf{F}_2(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) &= [\rho u_1 u_2, \rho u_2^2 + p, (\rho E + p)u_2, \rho_1 \mathcal{Z}_1 u_2, \dots, \rho_m \mathcal{Z}_m u_2]^T. \end{aligned} $$
### 3.3. Overall algorithm

Let us summarize the overall algorithm detailed in the previous sections. For updating the variables at each instant $t^n$ to their values at $t^{n+1}$, we proceed following the steps of the algorithm 1.

<caption>**Algorithm 1** Lagrange-Remap procedure</caption>

**I - Lagrangian step:**
1: Compute the acoustic fluxes $p_{i+1/2}^n, u_{i+1/2}^n$ thanks to (17).
2: Compute the time step $\Delta t$ according to the CFL constraint (18).
3: Update $\tilde{\rho}u, \tilde{\rho}E, \tilde{\rho}_k\tilde{Z}_k, \tilde{Z}_k$ with (16).

**II - Remap step:**
1: Compute the fluxes $\tilde{Z}_{k,i+1/2}$ for each interface $i + 1/2$.
2: Compute the remap fluxes of the conservative quantities $\tilde{\rho}_{i+1/2}, \tilde{\rho}u_{i+1/2}, \tilde{\rho}E_{i+1/2}, (\tilde{\rho}_k\tilde{Z}_k)_{i+1/2}$ with (20).
3: Update $(\rho u)^{n+1}, (\rho E)^{n+1}, (\rho_k Z_k)^{n+1}, Z_k^{n+1}$ with (19).

Combining the Lagrange (16) and Remap (19) steps, we obtain the following quasi-conservative scheme

$$ \left\{ \begin{array}{l} (\rho u)_i^{n+1} = (\rho u)_i^n - \frac{\Delta t}{\Delta x} [u^n \tilde{\rho} u]_i - \frac{\Delta t}{\Delta x} [p^n]_i, \\ (\rho E)_i^{n+1} = (\rho E)_i^n - \frac{\Delta t}{\Delta x} [u^n \tilde{\rho} E]_i - \frac{\Delta t}{\Delta x} [p^n u^n]_i, \\ (\rho_k Z_k)_i^{n+1} = (\rho_k Z_k)_i^n - \frac{\Delta t}{\Delta x} [u^n \tilde{\rho}_k Z_k]_i, \\ Z_{k,i}^{n+1} = Z_{k,i}^n - \frac{\Delta t}{\Delta x} [u^n \tilde{Z}_k]_i + \frac{\Delta t}{\Delta x} [u^n]_i \tilde{Z}_{k,i}. \end{array} \right. \quad (23) $$

As in [41], the overall discretization strategy is conservative with respect to the momentum, total energy and partial masses.

### 3.4. Extension to the multi-dimensional case

Without loss of generality, we consider a two-dimensional problem that is discretized over a Cartesian grid. Noting $\mathbf{u} = (u_1, u_2)^T$ the velocity, the m-component system (8) reads

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_{x_1} \mathbf{F}_1(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) + \partial_{x_2} \mathbf{F}_2(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) = 0, \\ \partial_t \mathcal{Z}_k + u_1 \partial_{x_1} \mathcal{Z}_k + u_2 \partial_{x_2} \mathcal{Z}_k = 0, \quad k = 1, \dots, m, \end{array} \right. \quad (24) $$

where $\mathbf{U} = [\rho u_1, \rho u_2, \rho E, \rho_1 \mathcal{Z}_1, \dots, \rho_m \mathcal{Z}_m]^T$, and

$$ \begin{aligned} \mathbf{F}_1(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) &= [\rho u_1^2 + p, \rho u_1 u_2, (\rho E + p)u_1, \rho_1 \mathcal{Z}_1 u_1, \dots, \rho_m \mathcal{Z}_m u_1]^T, \\ \mathbf{F}_2(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) &= [\rho u_1 u_2, \rho u_2^2 + p, (\rho E + p)u_2, \rho_1 \mathcal{Z}_1 u_2, \dots, \rho_m \mathcal{Z}_m u_2]^T. \end{aligned} $$
We approximate the solution of (24) by means of a dimensional splitting strategy. This boils down to solve successively in each space direction

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_{x_1} \mathbf{F}_1(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) = 0, \\ \partial_t \mathcal{Z}_k + u_1 \partial_{x_1} \mathcal{Z}_k = 0, \quad k = 1, \dots, m, \end{array} \right. $$

then

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_{x_2} \mathbf{F}_2(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) = 0, \\ \partial_t \mathcal{Z}_k + u_2 \partial_{x_2} \mathcal{Z}_k = 0, \quad k = 1, \dots, m. \end{array} \right. $$

In the first (resp. second) step, the velocity component $u_2$ (resp. $u_1$) is a scalar passively advected by the flow. The solution of the first (resp. second) step is simply approximated using the algorithm 1 detailed in section 3.3 supplemented with the update of the variable $\rho u_2$ (resp. $\rho u_1$) achieved thanks to an upwind scheme.

## 4. Anti-diffusive strategy for computing the color functions

In this part we focus on the construction of the numerical fluxes $\tilde{Z}_{k,i+1/2}$ required during the projection step of the Lagrange-Remap scheme presented in the previous section. As in [8, 41], we perform a specific discretization for this transport stage that involves an anti-dissipative algorithm for the advection of the $m$ color functions $\mathcal{Z}_k$. This transport algorithm has been proposed in [36]: it relies on a recursive construction for the flux of each color function $\mathcal{Z}_k$ and it guarantees that both relations (4)-(5) are satisfied. Here, we use the same notation as the previous section, and focus on the case of one-dimensional problems.

The choice of the fluxes $\tilde{Z}_{k,i+1/2}$ relies on the alternative

* if $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n < 0$ (resp. $u_{i+3/2}^n < 0$ and $u_{i+1/2}^n > 0$), then we set $\tilde{Z}_{k,i+1/2}$ to the upwind value of $\mathcal{Z}_k$, with respect to the sign of $u_{i+1/2}^n$,

* if $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (resp. $u_{i+3/2}^n < 0$ and $u_{i+1/2}^n < 0$), then we choose for $\tilde{Z}_{k,i+1/2}$ the closest value to the downwind value of $\mathcal{Z}_k$, with respect to the sign of $u_{i+1/2}^n$ within a trust interval $I_{k,i+1/2}^n$ that provides important features to our numerical scheme.

In the next sections, we present a procedure to define a trust interval $I_{k,i+1/2}^n$ that ensures:

1. consistency for $\tilde{Z}_{k,i+1/2}$,

2. stability for $\mathcal{Z}_k$,

3. verifies relation (5) for the discrete unknowns.
We approximate the solution of (24) by means of a dimensional splitting strategy. This boils down to solve successively in each space direction

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_{x_1} \mathbf{F}_1(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) = 0, \\ \partial_t \mathcal{Z}_k + u_1 \partial_{x_1} \mathcal{Z}_k = 0, \quad k = 1, \dots, m, \end{array} \right. $$

then

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_{x_2} \mathbf{F}_2(\mathbf{U}, \mathcal{Z}_1, \dots, \mathcal{Z}_m) = 0, \\ \partial_t \mathcal{Z}_k + u_2 \partial_{x_2} \mathcal{Z}_k = 0, \quad k = 1, \dots, m. \end{array} \right. $$

In the first (resp. second) step, the velocity component $u_2$ (resp. $u_1$) is a scalar passively advected by the flow. The solution of the first (resp. second) step is simply approximated using the algorithm 1 detailed in section 3.3 supplemented with the update of the variable $\rho u_2$ (resp. $\rho u_1$) achieved thanks to an upwind scheme.

## 4. Anti-diffusive strategy for computing the color functions

In this part we focus on the construction of the numerical fluxes $\tilde{Z}_{k,i+1/2}$ required during the projection step of the Lagrange-Remap scheme presented in the previous section. As in [8, 41], we perform a specific discretization for this transport stage that involves an anti-dissipative algorithm for the advection of the $m$ color functions $\mathcal{Z}_k$. This transport algorithm has been proposed in [36]: it relies on a recursive construction for the flux of each color function $\mathcal{Z}_k$ and it guarantees that both relations (4)-(5) are satisfied. Here, we use the same notation as the previous section, and focus on the case of one-dimensional problems.

The choice of the fluxes $\tilde{Z}_{k,i+1/2}$ relies on the alternative

* if $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n < 0$ (resp. $u_{i+3/2}^n < 0$ and $u_{i+1/2}^n > 0$), then we set $\tilde{Z}_{k,i+1/2}$ to the upwind value of $\mathcal{Z}_k$, with respect to the sign of $u_{i+1/2}^n$,

* if $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (resp. $u_{i+3/2}^n < 0$ and $u_{i+1/2}^n < 0$), then we choose for $\tilde{Z}_{k,i+1/2}$ the closest value to the downwind value of $\mathcal{Z}_k$, with respect to the sign of $u_{i+1/2}^n$ within a trust interval $I_{k,i+1/2}^n$ that provides important features to our numerical scheme.

In the next sections, we present a procedure to define a trust interval $I_{k,i+1/2}^n$ that ensures:

1. consistency for $\tilde{Z}_{k,i+1/2}$,

2. stability for $\mathcal{Z}_k$,

3. verifies relation (5) for the discrete unknowns.
### 4.1. Consistency and stability criterions

First we exhibit a sufficient criterion ensuring a consistency constraint for $\tilde{Z}_{k,i+1/2}$. Let us note

$$m_{k,i+1/2}^n = \min(Z_{k,i}^n, Z_{k,i+1}^n), \quad M_{k,i+1/2}^n = \max(Z_{k,i}^n, Z_{k,i+1}^n).$$

As in [20, 18, 41], consistency for the flux $\tilde{Z}_{k,i+1/2}$ is imposed by requiring

$$\tilde{Z}_{k,i+1/2} \in [m_{k,i+1/2}^n, M_{k,i+1/2}^n]. \qquad (25)$$

We now seek a sufficient condition on $\tilde{Z}_{k,i+1/2}$ that guarantees the stability for the color function $Z_k$ in a neighborhood cell of the interface $i+1/2$. We assume $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (resp. $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$): as $Z_k$ is evolved by a pure transport equation, stability within cell $i$ (resp. $i+1$) is ensured by a Harten-Leroux criterion [33, 44]

$$Z_{k,i}^{n+1} \in [m_{k,i-1/2}^n, M_{k,i-1/2}^n] \quad (\text{resp. } Z_{k,i+1}^{n+1} \in [m_{k,i+3/2}^n, M_{k,i+3/2}^n]) \qquad (26)$$

One should note this enforces a local maximum principle in the cell $i$ (resp. $i+1$). From relation (26) and following same lines as in the proof of [41, § 3.3.3, proposition 3.1], we are able to deduce a sufficient condition on $\tilde{Z}_{k,i+1/2}$ for ensuring stability for $Z_k$ when $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (resp. $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$) that is resumed by the following proposition.

**Proposition 4.** *Assuming that both conditions for consistency (25) and CFL (18) hold.*

(a) If $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$, set

$$\begin{align*} a_{k,i+1/2}^n &= Z_{k,i}^n + (M_{k,i-1/2}^n - Z_{k,i}^n) \left( \frac{u_{i-1/2}^n}{u_{i+1/2}^n} - \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right) \\ A_{k,i+1/2}^n &= Z_{k,i}^n + (m_{k,i-1/2}^n - Z_{k,i}^n) \left( \frac{u_{i-1/2}^n}{u_{i+1/2}^n} - \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right), \end{align*}$$

then

$$\begin{align*} \tilde{Z}_{k,i+1/2} \in [a_{i+1/2}^n, A_{i+1/2}^n] &\Rightarrow Z_{k,i}^{n+1} \in [m_{k,i+1/2}^n, M_{k,i+1/2}^n] \\ &\Rightarrow \text{stability for } Z_k \text{ in the cell } i, \end{align*}$$

(b) If $u_{i+3/2}^n < 0$ et $u_{i+1/2}^n < 0$, set

$$\begin{align*} a_{k,i+1/2}^n &= Z_{k,i+1}^n + (M_{k,i+3/2}^n - Z_{k,i+1}^n) \left( \frac{u_{i+3/2}^n}{u_{i+1/2}^n} + \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right) \\ A_{k,i+1/2}^n &= Z_{k,i+1}^n + (m_{k,i+3/2}^n - Z_{k,i+1}^n) \left( \frac{u_{i+3/2}^n}{u_{i+1/2}^n} + \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right), \end{align*}$$
### 4.1. Consistency and stability criterions

First we exhibit a sufficient criterion ensuring a consistency constraint for $\tilde{Z}_{k,i+1/2}$. Let us note

$$m_{k,i+1/2}^n = \min(Z_{k,i}^n, Z_{k,i+1}^n), \quad M_{k,i+1/2}^n = \max(Z_{k,i}^n, Z_{k,i+1}^n).$$

As in [20, 18, 41], consistency for the flux $\tilde{Z}_{k,i+1/2}$ is imposed by requiring

$$\tilde{Z}_{k,i+1/2} \in [m_{k,i+1/2}^n, M_{k,i+1/2}^n]. \qquad (25)$$

We now seek a sufficient condition on $\tilde{Z}_{k,i+1/2}$ that guarantees the stability for the color function $Z_k$ in a neighborhood cell of the interface $i+1/2$. We assume $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (resp. $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$): as $Z_k$ is evolved by a pure transport equation, stability within cell $i$ (resp. $i+1$) is ensured by a Harten-Leroux criterion [33, 44]

$$Z_{k,i}^{n+1} \in [m_{k,i-1/2}^n, M_{k,i-1/2}^n] \quad (\text{resp. } Z_{k,i+1}^{n+1} \in [m_{k,i+3/2}^n, M_{k,i+3/2}^n]) \qquad (26)$$

One should note this enforces a local maximum principle in the cell $i$ (resp. $i+1$). From relation (26) and following same lines as in the proof of [41, § 3.3.3, proposition 3.1], we are able to deduce a sufficient condition on $\tilde{Z}_{k,i+1/2}$ for ensuring stability for $Z_k$ when $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (resp. $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$) that is resumed by the following proposition.

**Proposition 4.** *Assuming that both conditions for consistency (25) and CFL (18) hold.*

(a) If $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$, set

$$\begin{align*} a_{k,i+1/2}^n &= Z_{k,i}^n + (M_{k,i-1/2}^n - Z_{k,i}^n) \left( \frac{u_{i-1/2}^n}{u_{i+1/2}^n} - \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right) \\ A_{k,i+1/2}^n &= Z_{k,i}^n + (m_{k,i-1/2}^n - Z_{k,i}^n) \left( \frac{u_{i-1/2}^n}{u_{i+1/2}^n} - \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right), \end{align*}$$

then

$$\begin{align*} \tilde{Z}_{k,i+1/2} \in [a_{i+1/2}^n, A_{i+1/2}^n] &\Rightarrow Z_{k,i}^{n+1} \in [m_{k,i+1/2}^n, M_{k,i+1/2}^n] \\ &\Rightarrow \text{stability for } Z_k \text{ in the cell } i, \end{align*}$$

(b) If $u_{i+3/2}^n < 0$ et $u_{i+1/2}^n < 0$, set

$$\begin{align*} a_{k,i+1/2}^n &= Z_{k,i+1}^n + (M_{k,i+3/2}^n - Z_{k,i+1}^n) \left( \frac{u_{i+3/2}^n}{u_{i+1/2}^n} + \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right) \\ A_{k,i+1/2}^n &= Z_{k,i+1}^n + (m_{k,i+3/2}^n - Z_{k,i+1}^n) \left( \frac{u_{i+3/2}^n}{u_{i+1/2}^n} + \frac{\Delta x}{\Delta t} \frac{1}{u_{i+1/2}^n} \right), \end{align*}$$
then

$$ \tilde{Z}_{k,i+1/2} \in [a_{i+1/2}^n, A_{i+1/2}^n] \Rightarrow \mathcal{Z}_{k,i+1}^{n+1} \in [m_{k,i+3/2}^n, M_{k,i+3/2}^n] \\ \Rightarrow \text{stability for } \mathcal{Z}_k \text{ in the cell } i+1. $$

Combining both criterion (25) and bounds given in proposition 4 we can state as in [8, 18, 20, 41] the following property: for $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (or. $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$) both consistency and stability are satisfied if

$$ \tilde{Z}_{k,i+1/2} \in [\omega_{k,i+1/2}^n, \Omega_{k,i+1/2}^n] \quad (27) $$

where

$$ [\omega_{k,i+1/2}^n, \Omega_{k,i+1/2}^n] = [m_{k,i+1/2}^n, M_{k,i+1/2}^n] \cap [a_{k,i+1/2}^n, A_{k,i+1/2}^n] \quad (28) $$

### 4.2. Recursive construction of the trust interval $I_{k,i+1/2}^n$

A first simple choice for constructing the trust interval is to take $I_{k,i+1/2}^n = [\omega_{k,i+1/2}^n, \Omega_{k,i+1/2}^n]$. Although this choice provides stability and consistency for $\mathcal{Z}_k$, it may unfortunately fail to comply with (5). This flaw can indeed be revealed with a pure transport test of an arbitrary number of materials (see [36]). To overcome this difficulty, a solution should consists in computing only the first $m-1$ color functions and deduce the last one from the unity constraint (5). In that case, this is the constraint (4) that is violated (see again [36]).

In consequence, we propose to construct a sequence of intervals $[d_{k,i+1/2}^n, D_{k,i+1/2}^n] \subset [\omega_{i+1/2}^n, \Omega_{i+1/2}^n]$, $k = 1, \dots, m$ thanks to a recursive process like [36] as follows.

(a) For $k=1$, set

$$ d_{1,i+1/2}^n = \max \left( \omega_{1,i+1/2}^n, 1 - \sum_{l=2}^{m-1} \Omega_{l,i+1/2}^n \right), $$

$$ D_{1,i+1/2}^n = \min \left( \Omega_{1,i+1/2}^n, 1 - \sum_{l=2}^{m-1} \omega_{l,i+1/2}^n \right). $$

(b) For $k=2, \dots, m-1$, suppose that $\tilde{Z}_{l,i+1/2}$ are already known for $l \le k$ and set

$$ d_{k,i+1/2}^n = \max \left( \omega_{k,i+1/2}^n, 1 - \sum_{l=1}^{k-1} \tilde{Z}_{l,i+1/2}^n - \sum_{l=k+1}^{m-1} \Omega_{l,i+1/2}^n \right), $$

$$ D_{k,i+1/2}^n = \min \left( \Omega_{k,i+1/2}^n, 1 - \sum_{l=1}^{k-1} \tilde{Z}_{l,i+1/2}^n - \sum_{l=k+1}^{m-1} \omega_{l,i+1/2}^n \right). $$
(c) For $k = m$, suppose that $\tilde{Z}_{l,i+1/2}$ are already known for $l < m$ and set for the last flux

$$d_{k,i+1/2}^n = \max\left(\omega_{k,i+1/2}^n, 1 - \sum_{l=1}^{m-1} \tilde{Z}_{l,i+1/2}\right),$$

$$D_{k,i+1/2}^n = \min\left(\Omega_{k,i+1/2}^n, 1 - \sum_{l=1}^{m-1} \tilde{Z}_{l,i+1/2}\right).$$

As pointed out in [36], we are ensured that $I_{k,i+1/2}^n \neq \emptyset$ since the fluxes are chosen in their admissibility interval $[\omega_{i+1/2}^n, \Omega_{i+1/2}^n]$. Moreover $[\omega_{i+1/2}^n, \Omega_{i+1/2}^n] \neq \emptyset$, since the upwind choice for $\tilde{Z}_{k,i+1/2}^n$ belongs to $[\omega_{i+1/2}^n, \Omega_{i+1/2}^n]$.

Then, the following result proved in [36] can be stated.

**Theorem 1.** Suppose that the CFL constraint (18) and $Z_{k,i}^n \in [0, 1]$, $\sum_{k=1}^m Z_{k,i}^n = 1$ are satisfied. If the fluxes are chosen such that $\tilde{Z}_{k,i+1/2} \in I_{k,i+1/2}^n = [d_{k,i+1/2}^n, D_{k,i+1/2}^n]$ then for all $k \in \{1, \dots, m\}$ and for all $i \in \mathbb{Z}$

$$1. \quad Z_{k,i}^{n+1} \in [0, 1],$$

$$2. \quad \sum_{k=1}^{m} Z_{k,i}^{n+1} = 1.$$

In that way, we have defined a trust interval $I_{k,i+1/2}^n$ such that the numerical approximation of $Z_k$ provides the three following features: stability, consistency and a discrete equivalent of (5) if

$$\tilde{Z}_{k,i+1/2} \in I_{k,i+1/2}^n = [d_{k,i+1/2}^n, D_{k,i+1/2}^n] \quad (29)$$

when $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (or $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$).

**Remark 3.** Here, the sequence of intervals $[d_{k,i+1/2}^n, D_{k,i+1/2}^n]$ are computed in the order of the increasing index $k$. Of course, this is not the only possible choice and the intervals construction order does not significantly impact the numerical result. This point will be illustrated in section 5.1. We also refer to [36, § 2] in the context of passive scalar convection.

### 4.3. Choice of $\tilde{Z}_{k,i+1/2}$

We detail here the definition of $\tilde{Z}_{k,i+1/2}$ that was drafted at the beginning of section 4 using similar lines as [8, 18, 20, 36, 41].

First, when $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (or $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$), we construct the trust interval $I_{k,i+1/2}^n$ given by (29).
(c) For $k = m$, suppose that $\tilde{Z}_{l,i+1/2}$ are already known for $l < m$ and set for the last flux

$$d_{k,i+1/2}^n = \max\left(\omega_{k,i+1/2}^n, 1 - \sum_{l=1}^{m-1} \tilde{Z}_{l,i+1/2}\right),$$

$$D_{k,i+1/2}^n = \min\left(\Omega_{k,i+1/2}^n, 1 - \sum_{l=1}^{m-1} \tilde{Z}_{l,i+1/2}\right).$$

As pointed out in [36], we are ensured that $I_{k,i+1/2}^n \neq \emptyset$ since the fluxes are chosen in their admissibility interval $[\omega_{i+1/2}^n, \Omega_{i+1/2}^n]$. Moreover $[\omega_{i+1/2}^n, \Omega_{i+1/2}^n] \neq \emptyset$, since the upwind choice for $\tilde{Z}_{k,i+1/2}^n$ belongs to $[\omega_{i+1/2}^n, \Omega_{i+1/2}^n]$.

Then, the following result proved in [36] can be stated.

**Theorem 1.** Suppose that the CFL constraint (18) and $Z_{k,i}^n \in [0, 1]$, $\sum_{k=1}^m Z_{k,i}^n = 1$ are satisfied. If the fluxes are chosen such that $\tilde{Z}_{k,i+1/2} \in I_{k,i+1/2}^n = [d_{k,i+1/2}^n, D_{k,i+1/2}^n]$ then for all $k \in \{1, \dots, m\}$ and for all $i \in \mathbb{Z}$

$$1. \quad Z_{k,i}^{n+1} \in [0, 1],$$

$$2. \quad \sum_{k=1}^{m} Z_{k,i}^{n+1} = 1.$$

In that way, we have defined a trust interval $I_{k,i+1/2}^n$ such that the numerical approximation of $Z_k$ provides the three following features: stability, consistency and a discrete equivalent of (5) if

$$\tilde{Z}_{k,i+1/2} \in I_{k,i+1/2}^n = [d_{k,i+1/2}^n, D_{k,i+1/2}^n] \quad (29)$$

when $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (or $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$).

**Remark 3.** Here, the sequence of intervals $[d_{k,i+1/2}^n, D_{k,i+1/2}^n]$ are computed in the order of the increasing index $k$. Of course, this is not the only possible choice and the intervals construction order does not significantly impact the numerical result. This point will be illustrated in section 5.1. We also refer to [36, § 2] in the context of passive scalar convection.

### 4.3. Choice of $\tilde{Z}_{k,i+1/2}$

We detail here the definition of $\tilde{Z}_{k,i+1/2}$ that was drafted at the beginning of section 4 using similar lines as [8, 18, 20, 36, 41].

First, when $u_{i+1/2}^n > 0$ and $u_{i-1/2}^n > 0$ (or $u_{i+1/2}^n < 0$ and $u_{i+3/2}^n < 0$), we construct the trust interval $I_{k,i+1/2}^n$ given by (29).
Then $\mathcal{Z}_{k,i+1/2}$ is set to the closest value to the downwind value $\mathcal{Z}_{k,i+1/2}^{do}$ (according to sign of $u_{i+1/2}^n$) within $I_{k,i+1/2}^n$ (see figure 2) in order to limit the numerical diffusion as in [18, 20, 36, 41]. Namely

$$ \tilde{Z}_{k,i+1/2} = \min_{Z \in I_{k,i+1/2}^n} |\mathcal{Z}_{k,i+1/2}^{do} - Z| \quad (30) $$

$$ \text{with } \mathcal{Z}_{k,i+1/2}^{do} = \begin{cases} \mathcal{Z}_{k,i+1}^n & \text{if } u_{i+1/2} > 0, \\ \mathcal{Z}_{k,i}^n & \text{otherwise.} \end{cases} $$

**Remark 4.** The trust interval $I_{k,i+1/2}^n$ is clearly defined for only two cases ($u_{i+1/2} > 0$ and $u_{i-1/2} > 0$) or ($u_{i+1/2} < 0$ and $u_{i+3/2} < 0$). In other cases, the upwind value is selected as default value for $\tilde{Z}_{k,i+1/2}$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Choice of the flux $\tilde{Z}_{i+1/2}$, three possible cases: 1) $\mathcal{Z}_{k,i+1/2}^{do} > D_{k,i+1/2}^n$, 2) $d_{k,i+1/2}^n \le \mathcal{Z}_{k,i+1/2}^{do} \le D_{k,i+1/2}^n$, 3) $\mathcal{Z}_{k,i+1/2}^{do} < d_{k,i+1/2}^n$.</figcaption></figure>

# 5. Numerical results

In this section, we present numerical results with one-dimensional and two-dimensional tests cases for the m-component interface flow model (8). We will compare simulations performed with the proposed Lagrange-Remap scheme (see algorithm 1) with anti-diffusive (30) and upwind (22) construction of the color function fluxes during the Remap step.

## 5.1. Test 1: One-dimensional five-material passive transport

We consider a one-dimensional periodic domain of length 1 m. At the initial instant six interfaces located at the positions

$$ X_1 = 0, \quad X_2 = 0.1, \quad X_3 = 0.25, \quad X_4 = 0.7, \quad X_5 = 0.9, \quad X_6 = 1, $$

separate five constant pure fluid states $k = 1, \dots, 5$. Due to periodicity, the discontinuity located at $x = X_6$, matches the discontinuity at $x = X_1$. Each region $k = 1, \dots, 5$ delimited by $X_k < x < X_{k+1}$ is occupied by a different material. The fluid in region $k \in \{1, 5\}$ is a perfect gas whose EOS is given by
Then $\mathcal{Z}_{k,i+1/2}$ is set to the closest value to the downwind value $\mathcal{Z}_{k,i+1/2}^{do}$ (according to sign of $u_{i+1/2}^n$) within $I_{k,i+1/2}^n$ (see figure 2) in order to limit the numerical diffusion as in [18, 20, 36, 41]. Namely

$$ \tilde{Z}_{k,i+1/2} = \min_{Z \in I_{k,i+1/2}^n} |\mathcal{Z}_{k,i+1/2}^{do} - Z| \quad (30) $$

with $\mathcal{Z}_{k,i+1/2}^{do} = \begin{cases} \mathcal{Z}_{k,i+1}^n & \text{if } u_{i+1/2} > 0, \\ \mathcal{Z}_{k,i}^n & \text{otherwise.} \end{cases}$

**Remark 4.** The trust interval $I_{k,i+1/2}^n$ is clearly defined for only two cases ($u_{i+1/2} > 0$ and $u_{i-1/2} > 0$) or ($u_{i+1/2} < 0$ and $u_{i+3/2} < 0$). In other cases, the upwind value is selected as default value for $\tilde{Z}_{k,i+1/2}$.

<table><caption>Figure 2: Choice of the flux $\tilde{Z}_{i+1/2}$, three possible cases: 1) $Z_{k,i+1/2}^{do} > D_{k,i+1/2}^n$, 2) $d_{k,i+1/2}^n \le Z_{k,i+1/2}^{do} \le D_{k,i+1/2}^n$, 3) $Z_{k,i+1/2}^{do} < d_{k,i+1/2}^n$.</caption><tr><td>Case 1</td><td>$d_{k,i+1/2}^n$</td><td>$\tilde{Z}_{k,i+1/2}^n = D_{k,i+1/2}^n$</td><td>$Z_{k,i+1/2}^{do}$</td></tr><tr><td>Case 2</td><td>$d_{k,i+1/2}^n$</td><td>$\tilde{Z}_{k,i+1/2}^n = Z_{k,i+1/2}^{do}$</td><td>$D_{k,i+1/2}^n$</td></tr><tr><td>Case 3</td><td>$Z_{k,i+1/2}^{do}$</td><td>$\tilde{Z}_{k,i+1/2}^n = d_{k,i+1/2}^n$</td><td>$D_{k,i+1/2}^n$</td></tr></table>

# 5. Numerical results

In this section, we present numerical results with one-dimensional and two-dimensional tests cases for the m-component interface flow model (8). We will compare simulations performed with the proposed Lagrange-Remap scheme (see algorithm 1) with anti-diffusive (30) and upwind (22) construction of the color function fluxes during the Remap step.

## 5.1. Test 1: One-dimensional five-material passive transport

We consider a one-dimensional periodic domain of length 1 m. At the initial instant six interfaces located at the positions

$$ X_1 = 0, \quad X_2 = 0.1, \quad X_3 = 0.25, \quad X_4 = 0.7, \quad X_5 = 0.9, \quad X_6 = 1, $$

separate five constant pure fluid states $k = 1, \dots, 5$. Due to periodicity, the discontinuity located at $x = X_6$, matches the discontinuity at $x = X_1$. Each region $k = 1, \dots, 5$ delimited by $X_k < x < X_{k+1}$ is occupied by a different material. The fluid in region $k \in \{1, 5\}$ is a perfect gas whose EOS is given by
<table><caption>Table 1: test 1, one-dimensional five-material passive transport. EOS parameters.</caption><thead><tr><th>material</th><th>γ<sub>k</sub></th><th>π<sub>k</sub> (Pa)</th><th>a<sub>k</sub> (Pa.m<sup>6</sup>/kg<sup>2</sup>)</th><th>b<sub>k</sub> (m<sup>3</sup>/kg)</th></tr></thead><tbody><tr><td>k = 1</td><td>1.6</td><td>×</td><td>×</td><td>×</td></tr><tr><td>k = 2</td><td>4.4</td><td>6.10<sup>8</sup></td><td>×</td><td>×</td></tr><tr><td>k = 3</td><td>1.4</td><td>×</td><td>5</td><td>1.10<sup>-3</sup></td></tr><tr><td>k = 4</td><td>2.4</td><td>2.10<sup>8</sup></td><td>×</td><td>×</td></tr><tr><td>k = 5</td><td>1.6</td><td>×</td><td>×</td><td>×</td></tr></tbody></table>

<table><caption>Table 2: test 1, one-dimensional five-material passive transport. Initial location of each material and initial density values.</caption><thead><tr><th>location</th><th>k</th><th>ρ<sub>k</sub> (kg.m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>X<sub>1</sub> = 0 ≤ x &lt; X<sub>2</sub></td><td>1</td><td>50.0</td></tr><tr><td>X<sub>2</sub> ≤ x &lt; X<sub>3</sub></td><td>2</td><td>1000.0</td></tr><tr><td>X<sub>3</sub> ≤ x &lt; X<sub>4</sub></td><td>3</td><td>500.0</td></tr><tr><td>X<sub>4</sub> ≤ x &lt; X<sub>5</sub></td><td>4</td><td>1200.0</td></tr><tr><td>X<sub>5</sub> ≤ x &lt; X<sub>6</sub> = 1</td><td>5</td><td>150.0</td></tr></tbody></table>

(12), while region $k \in \{2, 4\}$ contains a stiffened gas whose EOS is (13). A Van der Waals fluid that is governed by the EOS

$$p_k = \left( \frac{\gamma_k - 1}{1 - b_k \rho_k} \right) (\rho_k e_k + a_k \rho_k^2) - a_k \rho_k^2.$$

fills region $k = 3$. The initial value of both pressure $p$ and velocity $u$ are uniform and set to $p = 10^5$ Pa and $u = 100$ m.s$^{-1}$.

The EOS parameters used for each component are displayed in table 1. Let us underline that these initial data did not produce any phasic value for the Van der Waals fluid out of the region where the is hyperbolic. More precisely $c_k^2$ for $k = 1, \dots, 5$ and $c^2$ remained positive values whenever they need to be evaluated by the algorithm throughout the whole computation. The initial density within each regions are given in table 2. The computational domain is discretized over a 100-cell grid and the simulation is performed with a CFL coefficient $C_{CFL} = 0.9$.

Figures 3, 4 and 5 respectively display the variable $Z_k$, $Y_k$, $k = 1, \dots, 5$ and $\rho$ for two instants $t_{end} = 0.01$ s corresponding to one turn through the initial position, and $t_{end} = 1.5$ s after 150 turns. We can see that the upwind scheme performs very poorly from the first moments as the approximate solution has been flatten out by numerical diffusion. Despite the use of a coarse mesh, the solution obtained with the anti-diffusive scheme remains sharp and shows a very good agreement with the exact solution, especially for both color functions and mass fractions. To quantify the evolution of the numerical diffusion of the color functions during the computation, we count the number of cells where the variables $Z_k$ are numerically diffused. For a given time $t^n$, a cell $i$ is considered to be a diffusion cell for the material $k$ if $\varepsilon \le Z_{k,i} \le 1 - \varepsilon$, where $\varepsilon = 10^{-6}$. The
<table><caption>Table 1: test 1, one-dimensional five-material passive transport. EOS parameters.</caption><thead><tr><th>material</th><th>γ<sub>k</sub></th><th>π<sub>k</sub> (Pa)</th><th>a<sub>k</sub> (Pa.m<sup>6</sup>/kg<sup>2</sup>)</th><th>b<sub>k</sub> (m<sup>3</sup>/kg)</th></tr></thead><tbody><tr><td>k = 1</td><td>1.6</td><td>×</td><td>×</td><td>×</td></tr><tr><td>k = 2</td><td>4.4</td><td>6.10<sup>8</sup></td><td>×</td><td>×</td></tr><tr><td>k = 3</td><td>1.4</td><td>×</td><td>5</td><td>1.10<sup>-3</sup></td></tr><tr><td>k = 4</td><td>2.4</td><td>2.10<sup>8</sup></td><td>×</td><td>×</td></tr><tr><td>k = 5</td><td>1.6</td><td>×</td><td>×</td><td>×</td></tr></tbody></table>

<table><caption>Table 2: test 1, one-dimensional five-material passive transport. Initial location of each material and initial density values.</caption><thead><tr><th>location</th><th>k</th><th>ρ<sub>k</sub> (kg.m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>X<sub>1</sub> = 0 ≤ x &lt; X<sub>2</sub></td><td>1</td><td>50.0</td></tr><tr><td>X<sub>2</sub> ≤ x &lt; X<sub>3</sub></td><td>2</td><td>1000.0</td></tr><tr><td>X<sub>3</sub> ≤ x &lt; X<sub>4</sub></td><td>3</td><td>500.0</td></tr><tr><td>X<sub>4</sub> ≤ x &lt; X<sub>5</sub></td><td>4</td><td>1200.0</td></tr><tr><td>X<sub>5</sub> ≤ x &lt; X<sub>6</sub> = 1</td><td>5</td><td>150.0</td></tr></tbody></table>

(12), while region $k \in \{2, 4\}$ contains a stiffened gas whose EOS is (13). A Van der Waals fluid that is governed by the EOS

$$p_k = \left( \frac{\gamma_k - 1}{1 - b_k \rho_k} \right) (\rho_k e_k + a_k \rho_k^2) - a_k \rho_k^2.$$

fills region $k = 3$. The initial value of both pressure $p$ and velocity $u$ are uniform and set to $p = 10^5$ Pa and $u = 100$ m.s$^{-1}$.

The EOS parameters used for each component are displayed in table 1. Let us underline that these initial data did not produce any phasic value for the Van der Waals fluid out of the region where the is hyperbolic. More precisely $c_k^2$ for $k = 1, \dots, 5$ and $c^2$ remained positive values whenever they need to be evaluated by the algorithm throughout the whole computation. The initial density within each regions are given in table 2. The computational domain is discretized over a 100-cell grid and the simulation is performed with a CFL coefficient $C_{CFL} = 0.9$.

Figures 3, 4 and 5 respectively display the variable $Z_k$, $Y_k$, $k = 1, \dots, 5$ and $\rho$ for two instants $t_{end} = 0.01$ s corresponding to one turn through the initial position, and $t_{end} = 1.5$ s after 150 turns. We can see that the upwind scheme performs very poorly from the first moments as the approximate solution has been flatten out by numerical diffusion. Despite the use of a coarse mesh, the solution obtained with the anti-diffusive scheme remains sharp and shows a very good agreement with the exact solution, especially for both color functions and mass fractions. To quantify the evolution of the numerical diffusion of the color functions during the computation, we count the number of cells where the variables $Z_k$ are numerically diffused. For a given time $t^n$, a cell $i$ is considered to be a diffusion cell for the material $k$ if $\varepsilon \le Z_{k,i} \le 1 - \varepsilon$, where $\varepsilon = 10^{-6}$. The
<table><caption>Table 3: test 1, one-dimensional five-material passive transport. Difference $e_1$ for $t_{end} = 0.01$ s</caption><thead><tr><th>a</th><th>ρ</th><th>p</th><th>u</th></tr></thead><tbody><tr><td>e<sub>1</sub>(a)</td><td>8.85 10<sup>-12</sup></td><td>1.68 10<sup>-11</sup></td><td>6.27 10<sup>-14</sup></td></tr></tbody></table>

percent of diffusion cells for both anti-diffusive and upwind scheme thorough the computation is represented in figure 7. We observe that the anti-diffusive scheme performs computation with at most 2% of diffusion cells contrary to the upwind scheme that completely diffuses the profile of the color functions for $t = 0.01$s.

In addition, we note that the maximum principle $0 \le Z_k \le 1$, $k = 1, \dots, 5$ and the discrete equivalent of (5) are verified throughout the computation for $Z_k$ as expected. We can observe that the mass fractions $Y_k$ comply with $Y_k \in [0, 1]$ and $\sum_k Y_k = 1$. Finally, we also verify that the pressure and velocity profiles initially constant are preserved during the computation for each scheme as depicted in figure 6 as for the five-equation two-component model (see [41]). A proof of this property is given in Appendix C.

At last, we test the influence of the material ordering (then the order of the color function fluxes construction) on the numerical solution computed with the anti-diffusive scheme. To this end, we propose to solve again the one-dimensional five-material passive transport test with the anti-diffusive scheme using the following material numbering permutation

fluid 1 → fluid 2, fluid 2 → fluid 1, fluid 3 → fluid 5,

fluid 4 → fluid 3, fluid 5 → fluid 4.

This permutation can be summarized thanks to the application $\sigma$ over the set of indices $k$ defined such that $\sigma(1) = 2, \sigma(2) = 1, \sigma(3) = 5, \sigma(4) = 3, \sigma(5) = 4$. In order to compare the numerical solutions obtained with the anti-diffusive scheme for both numbering of the materials, we have measured the maximal relative difference over the time interval $[0, t_{end}]$ i.e.

$$e_1(a) = \max_{t \in [0, t_{end}]} \max_{x \in [0,1]} \left| \frac{a(x,t) - \hat{a}(x,t)}{a(x,t) + \hat{a}(x,t)} \right|$$

for $a \in \{\rho, p, u\}$ and the absolute difference

$$e_2(a_k) = \max_{t \in [0, t_{end}]} \max_{x \in [0,1]} \left| a_k(x,t) - \hat{a}_{\sigma(k)}(x,t) \right|$$

for $a_k \in \{Z_k, Y_k\}$. Here, the quantity $a$ computed for the test with modified material order is noted with $\hat{a}$. The tables 3-4 summarizes the $e_1$ and $e_2$ values obtained for the anti-diffusive scheme with $t_{end} = 0.01$s. As we can see, the order of material influence seems to be negligible since both the errors $e_1$ and $e_2$ are at most of magnitude $10^{-11}$.
<table><caption>Table 3: test 1, one-dimensional five-material passive transport. Difference $e_1$ for $t_{end} = 0.01$ s</caption><thead><tr><th>a</th><th>ρ</th><th>p</th><th>u</th></tr></thead><tbody><tr><td>e<sub>1</sub>(a)</td><td>8.85 10<sup>-12</sup></td><td>1.68 10<sup>-11</sup></td><td>6.27 10<sup>-14</sup></td></tr></tbody></table>

percent of diffusion cells for both anti-diffusive and upwind scheme thorough the computation is represented in figure 7. We observe that the anti-diffusive scheme performs computation with at most 2% of diffusion cells contrary to the upwind scheme that completely diffuses the profile of the color functions for $t = 0.01$s.

In addition, we note that the maximum principle $0 \le Z_k \le 1$, $k = 1, \dots, 5$ and the discrete equivalent of (5) are verified throughout the computation for $Z_k$ as expected. We can observe that the mass fractions $Y_k$ comply with $Y_k \in [0, 1]$ and $\sum_k Y_k = 1$. Finally, we also verify that the pressure and velocity profiles initially constant are preserved during the computation for each scheme as depicted in figure 6 as for the five-equation two-component model (see [41]). A proof of this property is given in Appendix C.

At last, we test the influence of the material ordering (then the order of the color function fluxes construction) on the numerical solution computed with the anti-diffusive scheme. To this end, we propose to solve again the one-dimensional five-material passive transport test with the anti-diffusive scheme using the following material numbering permutation

fluid 1 → fluid 2, fluid 2 → fluid 1, fluid 3 → fluid 5,

fluid 4 → fluid 3, fluid 5 → fluid 4.

This permutation can be summarized thanks to the application $\sigma$ over the set of indices $k$ defined such that $\sigma(1) = 2$, $\sigma(2) = 1$, $\sigma(3) = 5$, $\sigma(4) = 3$, $\sigma(5) = 4$. In order to compare the numerical solutions obtained with the anti-diffusive scheme for both numbering of the materials, we have measured the maximal relative difference over the time interval $[0, t_{end}]$ i.e.

$$e_1(a) = \max_{t \in [0, t_{end}]} \max_{x \in [0,1]} \left| \frac{a(x,t) - \hat{a}(x,t)}{a(x,t) + \hat{a}(x,t)} \right|$$

for $a \in \{\rho, p, u\}$ and the absolute difference

$$e_2(a_k) = \max_{t \in [0, t_{end}]} \max_{x \in [0,1]} \left| a_k(x,t) - \hat{a}_{\sigma(k)}(x,t) \right|$$

for $a_k \in \{Z_k, Y_k\}$. Here, the quantity $a$ computed for the test with modified material order is noted with $\hat{a}$. The tables 3-4 summarizes the $e_1$ and $e_2$ values obtained for the anti-diffusive scheme with $t_{end} = 0.01$s. As we can see, the order of material influence seems to be negligible since both the errors $e_1$ and $e_2$ are at most of magnitude $10^{-11}$.
<figure><img src="image_4.png" /><figcaption>Figure 3: test 1, one-dimensional five-material passive transport. Profile of the color functions $Z_k$ and of $\sum_{k=1}^5 Z_k$ for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme and the anti-diffusive solver.</figcaption></figure>

<table><caption>Table 4: test 1, one-dimensional five-material passive transport. Difference $e_2$ for $t_{end} = 0.01 \text{ s}$</caption><thead><tr><th></th><th>Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>Z<sub>3</sub></th><th>Z<sub>4</sub></th><th>Z<sub>5</sub></th></tr></thead><tbody><tr><th>a<sub>k</sub></th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>e<sub>2</sub>(a<sub>k</sub>)</th><td>6.12 10<sup>-12</sup></td><td>3.43 10<sup>-12</sup></td><td>3.45 10<sup>-12</sup></td><td>3.12 10<sup>-12</sup></td><td>6.15 10<sup>-12</sup></td></tr></tbody><thead><tr><th></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th><th>Y<sub>4</sub></th><th>Y<sub>5</sub></th></tr></thead><tbody><tr><th>a<sub>k</sub></th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>e<sub>2</sub>(a<sub>k</sub>)</th><td>1.79 10<sup>-11</sup></td><td>6.84 10<sup>-12</sup></td><td>7.49 10<sup>-12</sup></td><td>2.03 10<sup>-11</sup></td><td>2.02 10<sup>-11</sup></td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 3: test 1, one-dimensional five-material passive transport. Profile of the color functions $\mathcal{Z}_k$ and of $\sum_{k=1}^5 \mathcal{Z}_k$ for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme and the anti-diffusive solver.</figcaption></figure>

<table><caption>Table 4: test 1, one-dimensional five-material passive transport. Difference $e_2$ for $t_{end} = 0.01 \text{ s}$</caption><thead><tr><th></th><th>Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>Z<sub>3</sub></th><th>Z<sub>4</sub></th><th>Z<sub>5</sub></th></tr></thead><tbody><tr><th>a<sub>k</sub></th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>e<sub>2</sub>(a<sub>k</sub>)</th><td>6.12 10<sup>-12</sup></td><td>3.43 10<sup>-12</sup></td><td>3.45 10<sup>-12</sup></td><td>3.12 10<sup>-12</sup></td><td>6.15 10<sup>-12</sup></td></tr></tbody><thead><tr><th></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th><th>Y<sub>4</sub></th><th>Y<sub>5</sub></th></tr></thead><tbody><tr><th>a<sub>k</sub></th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>e<sub>2</sub>(a<sub>k</sub>)</th><td>1.79 10<sup>-11</sup></td><td>6.84 10<sup>-12</sup></td><td>7.49 10<sup>-12</sup></td><td>2.03 10<sup>-11</sup></td><td>2.02 10<sup>-11</sup></td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 4: test 1, one-dimensional five-material passive transport. Profile of $\mathcal{Y}_k$ and of $\sum_{i=1}^{5} \mathcal{Y}_k$ for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme and the anti-diffusive solver.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4: test 1, one-dimensional five-material passive transport. Profile of $\mathcal{Y}_k$ and of $\sum_{k=1}^{5} \mathcal{Y}_k$ for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme and the anti-diffusive solver.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: test 1, one-dimensional five-material passive transport. Profile of density for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme, the anti-diffusive solver and the exact solution.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: test 1, one-dimensional five-material passive transport. Profile of the pressure (top) and velocity (bottom) for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme, the anti-diffusive solver and the exact solution.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: test 1, one-dimensional five-material passive transport. Profile of density for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme, the anti-diffusive solver and the exact solution.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: test 1, one-dimensional five-material passive transport. Profile of the pressure (top) and velocity (bottom) for $t \in \{0.01 \text{ s}, 1.5 \text{ s}\}$. Comparison between the upwind scheme, the anti-diffusive solver and the exact solution.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: test 1, one-dimensional five-material passive transport. Percent of cells in the domain where the color functions $Z_k$ are numerically diffused for $t \in \{0s, 0.01s\}$ for both anti-diffusive (left) and upwind (right) schemes.</figcaption></figure>

<table><caption>Table 5: test 2, one-dimensional three-material Riemann problems juxtaposition. Initial values for pressure and density in the domain.</caption><thead><tr><th>location</th><th>k</th><th>ρ<sub>k</sub> (kg.m<sup>-3</sup>)</th><th>p<sub>k</sub> (Pa)</th><th>γ<sub>k</sub></th></tr></thead><tbody><tr><td>0 &lt; x &lt; 0.4</td><td>k = 1</td><td>1.0</td><td>1.0</td><td>1.6</td></tr><tr><td>0.4 ≤ x ≤ 0.6</td><td>k = 2</td><td>0.125</td><td>0.1</td><td>2.4</td></tr><tr><td>0.6 &lt; x &lt; 1</td><td>k = 3</td><td>0.1</td><td>0.1</td><td>1.4</td></tr></tbody></table>

## 5.2. Test 2: One-dimensional three-material Riemann problems juxtaposition

We consider now a test that consists of two side-by-side Riemann problems within a one-dimensional 1 m long domain. At $t = 0$ three perfect gases are separated by two interfaces located at $x = 0.4$ m and $x = 0.6$ m. Initially the fluid is at rest in the whole domain and both pressure and density are constant in each region $0 &lt; x &lt; 0.4$, $0.4 \le x \le 0.6$ and $0.6 &lt; x &lt; 1$. These values are displayed in table 5 along with the specific heat ratio value of each fluid. The initial discontinuity at $x = 0.6$ is a stationary material discontinuity. The exact solution of this problem is composed by a set of waves that are depicted in figure 8: the Riemann problem centered at $x = 0.4$ produces a leftward going rarefaction wave and two waves that travels towards the right end of the domain. The fastest of these waves is a shock that propagates at speed $D_1 \approx 2.2780 \text{ m.s}^{-1}$. This shock hits at $t_{\text{shock}} \approx 0.0878 \text{ s}$ the stationary discontinuity located at $x = 0.6$ and triggers another two-component Riemann problem centered at $x = 0.6$. This Riemann problem is also composed by a rarefaction wave that travels leftwards and two waves moving towards the right end of the domain: a material discontinuity and a shock wave of speed $D_2 \approx 2.034 \text{ m.s}^{-1}$. The computational domain is discretized over a 500-cell mesh and we perform simulations with both the anti-diffusive and the Lagrange-
<figure><img src="image_7.png" /><figcaption>Figure 7: test 1, one-dimensional five-material passive transport. Percent of cells in the domain where the color functions $Z_k$ are numerically diffused for $t \in \{0s, 0.01s\}$ for both anti-diffusive (left) and upwind (right) schemes.</figcaption></figure>

<table><caption>Table 5: test 2, one-dimensional three-material Riemann problems juxtaposition. Initial values for pressure and density in the domain.</caption><thead><tr><th>location</th><th>k</th><th>ρ<sub>k</sub> (kg.m<sup>-3</sup>)</th><th>p<sub>k</sub> (Pa)</th><th>γ<sub>k</sub></th></tr></thead><tbody><tr><td>0 &lt; x &lt; 0.4</td><td>k = 1</td><td>1.0</td><td>1.0</td><td>1.6</td></tr><tr><td>0.4 ≤ x ≤ 0.6</td><td>k = 2</td><td>0.125</td><td>0.1</td><td>2.4</td></tr><tr><td>0.6 &lt; x &lt; 1</td><td>k = 3</td><td>0.1</td><td>0.1</td><td>1.4</td></tr></tbody></table>

## 5.2. Test 2: One-dimensional three-material Riemann problems juxtaposition

We consider now a test that consists of two side-by-side Riemann problems within a one-dimensional 1 m long domain. At $t = 0$ three perfect gases are separated by two interfaces located at $x = 0.4$ m and $x = 0.6$ m. Initially the fluid is at rest in the whole domain and both pressure and density are constant in each region $0 &lt; x &lt; 0.4$, $0.4 \le x \le 0.6$ and $0.6 &lt; x &lt; 1$. These values are displayed in table 5 along with the specific heat ratio value of each fluid. The initial discontinuity at $x = 0.6$ is a stationary material discontinuity. The exact solution of this problem is composed by a set of waves that are depicted in figure 8: the Riemann problem centered at $x = 0.4$ produces a leftward going rarefaction wave and two waves that travels towards the right end of the domain. The fastest of these waves is a shock that propagates at speed $D_1 \approx 2.2780 \text{ m.s}^{-1}$. This shock hits at $t_{\text{shock}} \approx 0.0878 \text{ s}$ the stationary discontinuity located at $x = 0.6$ and triggers another two-component Riemann problem centered at $x = 0.6$. This Riemann problem is also composed by a rarefaction wave that travels leftwards and two waves moving towards the right end of the domain: a material discontinuity and a shock wave of speed $D_2 \approx 2.034 \text{ m.s}^{-1}$. The computational domain is discretized over a 500-cell mesh and we perform simulations with both the anti-diffusive and the Lagrange-
<figure><img src="image_9.png" /><figcaption>Figure 8: test 2, one-dimensional three-material Riemann problems juxtaposition. Wave structure of the solution in the (x,t)-plane.</figcaption></figure>

Remap upwind scheme using $C_{CFL} = 0.8$. Transparent boundary conditions are set at each end of the domain.

Figure 9 displays the profile of the density, pressure and velocity and figures 10 and 11 show respectively the color functions and mass fractions at instant $t_{end} = 0.12$ s. As we can see, the agreement between the solution obtained with both solvers and the exact solution is good. As for the two-component solver of [41], the behavior of the upwind solver and the anti-diffusive solvers are quite similar for the approximation of the shock waves and the rarefaction waves. On the contrary the anti-diffusive mechanism is very efficient for the resolution of the contact waves that coincide with the material interfaces. Indeed, the density jump at the interface between fluid 1 and 2 (resp. 2 and 3) is captured with very few cells of numerical diffusion: as depicted in figure 12 the anti-diffusive scheme generates at most two cells of numerical diffusion for each variable $Z_k$ (i.e. 0.4% of cells). Furthermore, we again satisfy the maximum principle (4) and discrete unit constraint (5).

A slight undershoot across both interfaces between the different materials is present on the density profiles for the anti-diffusive scheme. This is typical of the proposed method, indeed similar observations have been made two-component case [8, 41]. Numerical experiment shows that this defect disappears with mesh refinement demonstrating that the numerical solution converges to the exact one. In order to illustrate this point, we run the same three-material shock tube with a finer mesh of 50000 cells. Results for the density, the pressure and the velocity are displayed in figure 13 and as expected the numerical solution converge to the exact one, in particular the overshoots at the interfaces for the density profile have significantly shrunk. The results for the mass fractions and color functions agree with the exact solutions without any surprising phenomena.

As for the one dimensional five-material passive transport problem, we examine now the influence of the material ordering on the computation of the
<figure><img src="image_8.png" /><figcaption>Figure 8: test 2, one-dimensional three-material Riemann problems juxtaposition. Wave structure of the solution in the (x,t)-plane.</figcaption></figure>

Remap upwind scheme using $C_{CFL} = 0.8$. Transparent boundary conditions are set at each end of the domain.

Figure 9 displays the profile of the density, pressure and velocity and figures 10 and 11 show respectively the color functions and mass fractions at instant $t_{end} = 0.12$ s. As we can see, the agreement between the solution obtained with both solvers and the exact solution is good. As for the two-component solver of [41], the behavior of the upwind solver and the anti-diffusive solvers are quite similar for the approximation of the shock waves and the rarefaction waves. On the contrary the anti-diffusive mechanism is very efficient for the resolution of the contact waves that coincide with the material interfaces. Indeed, the density jump at the interface between fluid 1 and 2 (resp. 2 and 3) is captured with very few cells of numerical diffusion: as depicted in figure 12 the anti-diffusive scheme generates at most two cells of numerical diffusion for each variable $Z_k$ (i.e. 0.4% of cells). Furthermore, we again satisfy the maximum principle (4) and discrete unit constraint (5).

A slight undershoot across both interfaces between the different materials is present on the density profiles for the anti-diffusive scheme. This is typical of the proposed method, indeed similar observations have been made two-component case [8, 41]. Numerical experiment shows that this defect disappears with mesh refinement demonstrating that the numerical solution converges to the exact one. In order to illustrate this point, we run the same three-material shock tube with a finer mesh of 50000 cells. Results for the density, the pressure and the velocity are displayed in figure 13 and as expected the numerical solution converge to the exact one, in particular the overshoots at the interfaces for the density profile have significantly shrunk. The results for the mass fractions and color functions agree with the exact solutions without any surprising phenomena.

As for the one dimensional five-material passive transport problem, we examine now the influence of the material ordering on the computation of the
<figure><img src="image_10.png" /><figcaption>Figure 9: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the density (top left) and zoom on the second interface (top right), pressure and velocity in the domain at instant $t_{end} = 0.12$ s for the 500-cell mesh.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 9: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the density (top left) and zoom on the second interface (top right), pressure and velocity in the domain at instant $t_{end} = 0.12$ s for the 500-cell mesh.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the color functions $Z_k$, $k = 1, 2, 3$ and of $Z_1 + Z_2 + Z_3$ at instant $t_{end} = 0.12$ s for the 500-cell mesh.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 10: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the color functions $Z_k$, $k = 1, 2, 3$ and of $Z_1 + Z_2 + Z_3$ at instant $t_{end} = 0.12$ s for the 500-cell mesh.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the mass fractions $Y_k$, $k = 1, 2, 3$ and of $Y_1 + Y_2 + Y_3$ at instant $t_{end} = 0.12$ s for the 500-cell mesh.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: test 2, one-dimensional three-material Riemann problems juxtaposition. Percent of cells in the domain where the color functions $Z_k$ are numerically diffused for $t \in \{0, 0.12\}$ s for both anti-diffusive (left) and upwind (right) schemes.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the mass fractions $Y_k$, $k = 1, 2, 3$ and of $Y_1 + Y_2 + Y_3$ at instant $t_{end} = 0.12$ s for the 500-cell mesh.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: test 2, one-dimensional three-material Riemann problems juxtaposition. Percent of cells in the domain where the color functions $Z_k$ are numerically diffused for $t \in \{0, 0.12\}$ s for both anti-diffusive (left) and upwind (right) schemes.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the density (top left) and zoom on the second interface (top right), pressure and velocity in the domain at instant $t_{end} = 0.12$ s for the 50000-cell mesh.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 13: test 2, one-dimensional three-material Riemann problems juxtaposition. Profile of the density (top left) and zoom on the second interface (top right), pressure and velocity in the domain at instant $t_{end} = 0.12$ s for the 50000-cell mesh.</figcaption></figure>
<table><caption>Table 6: test 2, one-dimensional three-material Riemann problems juxtaposition. Differences $e_1$ for $t_{end} = 0.12s$</caption><thead><tr><th>a</th><th>ρ</th><th>p</th><th>u</th></tr></thead><tbody><tr><td>e<sub>1</sub>(a)</td><td>2.19 10<sup>-14</sup></td><td>9.70<sup>-15</sup></td><td>3.96 10<sup>-12</sup></td></tr></tbody></table>

<table><caption>Table 7: test 2, one-dimensional three-material Riemann problems juxtaposition. Differences $e_2$ for $t_{end} = 0.12s$</caption><thead><tr><th>a<sub>k</sub></th><th>Z&#773;<sub>1</sub></th><th>Z&#773;<sub>2</sub></th><th>Z&#773;<sub>3</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th></tr></thead><tbody><tr><td>e<sub>2</sub>(a<sub>k</sub>)</td><td>2.79 10<sup>-14</sup></td><td>2.98 10<sup>-14</sup></td><td>2.99 10<sup>-14</sup></td><td>6.47 10<sup>-14</sup></td><td>6.47 10<sup>-14</sup></td><td>3.09 10<sup>-14</sup></td></tr></tbody></table>

numerical solution with the anti-diffusive scheme. This time, the following per-mutation has been considered $\sigma(1) = 2, \sigma(2) = 3, \sigma(3) = 1$. To quantify the impact of the fluids ordering change, the absolute and relative difference $e_1, e_2$ for $t_{end} = 0.12s$ are computed and their numerical values are given in the tables 6-7. Once again, both $e_1, e_2$ remain very small (at most $10^{-12}$), demonstrating that the fluid order (and then the computation order of the fluxes) does not impact the numerical solution in a significant manner.

At last, in order to examine the accuracy and convergence rate of the proposed numerical schemes, we perform several simulations for different meshes with a number of cells varying between 100 and 10000. We compute the $L^1$ relative error for each variable $\rho, u, p, Z_k, Y_k$ to deduce the convergence rates of both upwind and anti-diffusive scheme see table 8. Like for the two-material case [41], we observe that the accuracy is improved by the anti-diffusive scheme (up to first order) for the variables that are sensitive to contact discontinuities that carry material interfaces.

### 5.3. Test 3: One-dimensional three-material Riemann problems juxtaposition involving high pressure ratios

We now test the proposed scheme against a three-material Riemann problem with high pressure ratios. We consider again a one-dimensional 1 m long domain occupied by three materials separated by two interfaces localized at $x = 0.75$ m and $x = 0.95$ m. At the initial time, the flow is at rest and the first region $0 < x < 0.75$ is occupied by a stiffened gas while the two remaining regions

<table><caption>Table 8: test 2, one-dimensional three-material Riemann problems juxtaposition. Convergence rate estimates for all the variables computed with both upwind scheme and anti-diffusive scheme.</caption><thead><tr><th>Variable</th><th>ρ</th><th>p</th><th>u</th><th>Z&#773;<sub>1</sub></th><th>Z&#773;<sub>2</sub></th><th>Z&#773;<sub>3</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th></tr></thead><tbody><tr><td>Upwind</td><td>0.646</td><td>0.776</td><td>0.796</td><td>0.545</td><td>0.535</td><td>0.519</td><td>0.483</td><td>0.493</td><td>0.517</td></tr><tr><td>Anti-Diff.</td><td>0.786</td><td>0.783</td><td>0.807</td><td>1.002</td><td>1.031</td><td>1.112</td><td>1.003</td><td>1.032</td><td>1.101</td></tr></tbody></table>
<table><caption>Table 6: test 2, one-dimensional three-material Riemann problems juxtaposition. Differences $e_1$ for $t_{end} = 0.12s$</caption><thead><tr><th>a</th><th>ρ</th><th>p</th><th>u</th></tr></thead><tbody><tr><td>e<sub>1</sub>(a)</td><td>2.19 10<sup>-14</sup></td><td>9.70<sup>-15</sup></td><td>3.96 10<sup>-12</sup></td></tr></tbody></table>

<table><caption>Table 7: test 2, one-dimensional three-material Riemann problems juxtaposition. Differences $e_2$ for $t_{end} = 0.12s$</caption><thead><tr><th>a<sub>k</sub></th><th>Z&#773;<sub>1</sub></th><th>Z&#773;<sub>2</sub></th><th>Z&#773;<sub>3</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th></tr></thead><tbody><tr><td>e<sub>2</sub>(a<sub>k</sub>)</td><td>2.79 10<sup>-14</sup></td><td>2.98 10<sup>-14</sup></td><td>2.99 10<sup>-14</sup></td><td>6.47 10<sup>-14</sup></td><td>6.47 10<sup>-14</sup></td><td>3.09 10<sup>-14</sup></td></tr></tbody></table>

numerical solution with the anti-diffusive scheme. This time, the following per-mutation has been considered $\sigma(1) = 2, \sigma(2) = 3, \sigma(3) = 1$. To quantify the impact of the fluids ordering change, the absolute and relative difference $e_1, e_2$ for $t_{end} = 0.12s$ are computed and their numerical values are given in the tables 6-7. Once again, both $e_1, e_2$ remain very small (at most $10^{-12}$), demonstrating that the fluid order (and then the computation order of the fluxes) does not impact the numerical solution in a significant manner.

At last, in order to examine the accuracy and convergence rate of the proposed numerical schemes, we perform several simulations for different meshes with a number of cells varying between 100 and 10000. We compute the $L^1$ relative error for each variable $\rho, u, p, Z_k, Y_k$ to deduce the convergence rates of both upwind and anti-diffusive scheme see table 8. Like for the two-material case [41], we observe that the accuracy is improved by the anti-diffusive scheme (up to first order) for the variables that are sensitive to contact discontinuities that carry material interfaces.

### 5.3. Test 3: One-dimensional three-material Riemann problems juxtaposition involving high pressure ratios

We now test the proposed scheme against a three-material Riemann problem with high pressure ratios. We consider again a one-dimensional 1 m long domain occupied by three materials separated by two interfaces localized at $x = 0.75$ m and $x = 0.95$ m. At the initial time, the flow is at rest and the first region $0 < x < 0.75$ is occupied by a stiffened gas while the two remaining regions

<table><caption>Table 8: test 2, one-dimensional three-material Riemann problems juxtaposition. Convergence rate estimates for all the variables computed with both upwind scheme and anti-diffusive scheme.</caption><thead><tr><th>Variable</th><th>ρ</th><th>p</th><th>u</th><th>Z&#773;<sub>1</sub></th><th>Z&#773;<sub>2</sub></th><th>Z&#773;<sub>3</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>Y<sub>3</sub></th></tr></thead><tbody><tr><td>Upwind</td><td>0.646</td><td>0.776</td><td>0.796</td><td>0.545</td><td>0.535</td><td>0.519</td><td>0.483</td><td>0.493</td><td>0.517</td></tr><tr><td>Anti-Diff.</td><td>0.786</td><td>0.783</td><td>0.807</td><td>1.002</td><td>1.031</td><td>1.112</td><td>1.003</td><td>1.032</td><td>1.101</td></tr></tbody></table>
<table><caption>Table 9: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Initial values for pressure and density in the domain.</caption><thead><tr><th>location</th><th>k</th><th>ρ (kg.m<sup>-3</sup>)</th><th>p (Pa)</th><th>γ<sub>k</sub></th><th>π<sub>k</sub></th></tr></thead><tbody><tr><td>0 ≤ x &lt; 0.75</td><td>1</td><td>1000</td><td>1 10<sup>9</sup></td><td>4.4</td><td>6 10<sup>8</sup></td></tr><tr><td>0.75 ≤ x ≤ 0.95</td><td>2</td><td>50</td><td>1 10<sup>5</sup></td><td>2.4</td><td>×</td></tr><tr><td>0.95 &lt; x ≤ 1</td><td>3</td><td>1</td><td>1 10<sup>5</sup></td><td>1.4</td><td>×</td></tr></tbody></table>

respectively $0.75 \le x \le 0.95$ and $0.95 < x < 1$ are occupied by two different perfect gases. The EOSs parameters as well as initial states for this test are provided in table 9. The computational domain is discretized over a 2000-cell mesh. As the acoustic solver is only first order accurate the mesh is chosen sufficiently fine to capture the waves dynamic especially in the smaller regions of the domain occupied by the perfect gases. We perform simulations with both anti-diffusive and upwind schemes using $C_{CFL} = 0.8$ until $t_{end} = 270$ µs. Transparent boundary conditions are used at each boundary of the domain.

The wave pattern for this test is the same as the one depicted in figure 8. Indeed, the initial pressure ratio between the two first material is $10^4$ which generates a strong shock that propagates towards the right at velocity $D_1 = 819.92$ m.s$^{-1}$, and a rarefaction moving to the left. At time $t_{shock} = 26.07$ µs this shock reaches the second interface (as depicted in figure 8) and generates a new shock moving two the right at the speed $D_2 = 1271$ m.s$^{-1}$ and a rarefaction wave traveling leftwards between both interfaces. The profiles of the pressure, velocity and density computed with both schemes are depicted in figure 14 along with the exact solution. Both numerical solutions are in good agreement with the exact solution and illustrate the ability of the proposed scheme to deal with high pressure ratio. For this test anti-diffusive scheme provide approximations for the color functions and mass fractions that remain very sharp (see figure 15 and 16 respectively) and that verify the maximum principle and the discrete unit constraint.

### 5.4. Test 4: Two-dimensional passive transport with four components

We aim here at testing the ability of our numerical scheme to accurately transport material interfaces with a complex topology. The computational domain $D$ is the square $D = [0, 60] \times [0, 60]$. We suppose that this domain is occupied by four perfect gases numbered $k = 1, \dots, 4$. At $t = 0$, the velocity and the pressure are uniform in the domain, their value are respectively $(u_1, u_2) = (\sqrt{2}$ m/s, $\sqrt{3}$ m/s), $p = 1$ Pa. We consider the subsets $C$, $S$, $P$ of $D$ defined by

$$P = \{(x_1, x_2) \in D; (x_1, x_2) \in [27.5, 32.5]^2\},$$

$$C = \{(x_1, x_2) \in D; (x_1 - 30)^2 + (x_2 - 30)^2 \le 15^2\},$$

$$S = \{(x_1, x_2) \in D; 37.5 > x_2 > 22.5, \sqrt{3}(x_1 - 30) + 15 < x_2 < \sqrt{3}(x_1 - 30) + 45, \\ -\sqrt{3}(x_1 - 30) + 15 < x_2 < -\sqrt{3}(x_1 - 30) + 45\}.$$
<table><caption>Table 9: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Initial values for pressure and density in the domain.</caption><thead><tr><th>location</th><th>k</th><th>ρ (kg.m<sup>-3</sup>)</th><th>p (Pa)</th><th>γ<sub>k</sub></th><th>π<sub>k</sub></th></tr></thead><tbody><tr><td>0 ≤ x &lt; 0.75</td><td>1</td><td>1000</td><td>1 10<sup>9</sup></td><td>4.4</td><td>6 10<sup>8</sup></td></tr><tr><td>0.75 ≤ x ≤ 0.95</td><td>2</td><td>50</td><td>1 10<sup>5</sup></td><td>2.4</td><td>×</td></tr><tr><td>0.95 &lt; x ≤ 1</td><td>3</td><td>1</td><td>1 10<sup>5</sup></td><td>1.4</td><td>×</td></tr></tbody></table>

respectively $0.75 \le x \le 0.95$ and $0.95 < x < 1$ are occupied by two different perfect gases. The EOSs parameters as well as initial states for this test are provided in table 9. The computational domain is discretized over a 2000-cell mesh. As the acoustic solver is only first order accurate the mesh is chosen sufficiently fine to capture the waves dynamic especially in the smaller regions of the domain occupied by the perfect gases. We perform simulations with both anti-diffusive and upwind schemes using $C_{CFL} = 0.8$ until $t_{end} = 270$ µs. Transparent boundary conditions are used at each boundary of the domain.

The wave pattern for this test is the same as the one depicted in figure 8. Indeed, the initial pressure ratio between the two first material is $10^4$ which generates a strong shock that propagates towards the right at velocity $D_1 = 819.92$ m.s$^{-1}$, and a rarefaction moving to the left. At time $t_{shock} = 26.07$ µs this shock reaches the second interface (as depicted in figure 8) and generates a new shock moving two the right at the speed $D_2 = 1271$ m.s$^{-1}$ and a rarefaction wave traveling leftwards between both interfaces. The profiles of the pressure, velocity and density computed with both schemes are depicted in figure 14 along with the exact solution. Both numerical solutions are in good agreement with the exact solution and illustrate the ability of the proposed scheme to deal with high pressure ratio. For this test anti-diffusive scheme provide approximations for the color functions and mass fractions that remain very sharp (see figure 15 and 16 respectively) and that verify the maximum principle and the discrete unit constraint.

### 5.4. Test 4: Two-dimensional passive transport with four components

We aim here at testing the ability of our numerical scheme to accurately transport material interfaces with a complex topology. The computational domain $D$ is the square $D = [0, 60] \times [0, 60]$. We suppose that this domain is occupied by four perfect gases numbered $k = 1, \dots, 4$. At $t = 0$, the velocity and the pressure are uniform in the domain, their value are respectively $(u_1, u_2) = (\sqrt{2}$ m/s, $\sqrt{3}$ m/s), $p = 1$ Pa. We consider the subsets $C$, $S$, $P$ of $D$ defined by

$$P = \{(x_1, x_2) \in D; (x_1, x_2) \in [27.5, 32.5]^2\},$$

$$C = \{(x_1, x_2) \in D; (x_1 - 30)^2 + (x_2 - 30)^2 \le 15^2\},$$

$$S = \{(x_1, x_2) \in D; 37.5 > x_2 > 22.5, \sqrt{3}(x_1 - 30) + 15 < x_2 < \sqrt{3}(x_1 - 30) + 45, \\ -\sqrt{3}(x_1 - 30) + 15 < x_2 < -\sqrt{3}(x_1 - 30) + 45\}.$$
<figure><img src="image_15.png" /><figcaption>Figure 14: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Profile of the density, pressure (semi-log scale) and velocity in the domain at instant $t_{end} = 1.2 \ 10^{-4}$ s.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Profile of the density, pressure (semi-log scale) and velocity in the domain at instant $t_{end} = 1.2 \ 10^{-4}$ s.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Profile of the color functions $Z_k$, $k = 1, 2, 3$ and of $Z_1 + Z_2 + Z_3$ at instant $t_{end} = 1.2 \ 10^{-4}$ s.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 15: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Profile of the color functions $Z_k$, $k = 1, 2, 3$ and of $Z_1 + Z_2 + Z_3$ at instant $t_{end} = 1.2 \ 10^{-4}$ s.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Profile of the mass fractions $Y_k$, $k = 1, 2, 3$ and of $Y_1+Y_2+Y_3$ at instant $t_{end} = 1.2 \ 10^{-4}$ s.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 16: test 3, one-dimensional three-material Riemann problems juxtaposition involving high pressure ratios. Profile of the mass fractions $Y_k$, $k = 1, 2, 3$ and of $Y_1+Y_2+Y_3$ at instant $t_{end} = 1.2 \ 10^{-4}$ s.</figcaption></figure>
<table><caption>Table 10: test 4, two-dimensional passive transport with four components. Initial location and EOS parameters of the components.</caption><thead><tr><th>material index</th><th>γ<sub>k</sub></th><th>initial location</th><th>initial density (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>k = 1</td><td>2.2</td><td>D \ C</td><td>0.01</td></tr><tr><td>k = 2</td><td>1.6</td><td>C \ S</td><td>0.1</td></tr><tr><td>k = 3</td><td>1.4</td><td>S \ P</td><td>1.0</td></tr><tr><td>k = 4</td><td>1.2</td><td>P</td><td>10.0</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>Figure 17: test 4, two-dimensional passive transport with four components. Mapping of $\sum_{k=1}^{4} k\mathcal{Z}_k$ (left) and the density (right) at $t = 0$.</figcaption></figure>

The regions C, S, P are respectively a disk, a star shaped and a square portion of D. The EOS parameters of these fluids are given in table 10 along with their initial location within the domain and density value. The domain D is discretized over a 200 × 200-cell regular mesh and we impose periodic boundary conditions. The resulting initial condition for the material interface and the density is depicted in figure 17. Figure 18 displays a mapping of $\sum_{k=1}^{4} k\mathcal{Z}_k$ at $t = 42.5$ s obtained with both the upwind scheme and the anti-diffusive scheme. The final instant is reached after 2626 time steps for the upwind scheme and 2627 time steps for the anti-diffusive scheme. We plot in figures 18 and 19 the density profiles and the interfaces by representing $\sum_{k=1}^{4} k\mathcal{Z}_k$ at the final time $t = 42.5$ s. As in the one dimensional case, the interfaces remain quite sharp and their original topology is well preserved by the numerical approximation computed by the anti-diffusive scheme. With the upwind scheme, the interface shapes are completely spread due to numerical diffusion. Nevertheless, both schemes preserve velocity and pressure profiles as in the one-dimensional test. Indeed, at the final time the relative errors $||\bar{u} - u||_{L^1}$ and $||p - \bar{p}||_{L^1}$ for the references $\bar{u} = \sqrt{5}$ and $\bar{p} = 1$ remain bounded between $3.9 \times 10^{-16}$ and $4.9 \times 10^{-14}$ (see table 11 for more details). Finally, figure 20 displays the graphs $t \mapsto \int_D \sum_{k=1}^{4} \mathcal{Z}_k(x, t)dx$ and $t \mapsto \int_D \sum_{k=1}^{4} \mathcal{Y}_k(x, t)dx$ for the anti-diffusive scheme. Again, we observe that the scheme preserve the unity constraints (5) for both color functions and mass fractions at each time step.
<table><caption>Table 10: test 4, two-dimensional passive transport with four components. Initial location and EOS parameters of the components.</caption><thead><tr><th>material index</th><th>γ<sub>k</sub></th><th>initial location</th><th>initial density (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>k = 1</td><td>2.2</td><td>D \ C</td><td>0.01</td></tr><tr><td>k = 2</td><td>1.6</td><td>C \ S</td><td>0.1</td></tr><tr><td>k = 3</td><td>1.4</td><td>S \ P</td><td>1.0</td></tr><tr><td>k = 4</td><td>1.2</td><td>P</td><td>10.0</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Figure 17: test 4, two-dimensional passive transport with four components. Mapping of $\sum_{k=1}^{4} kZ_k$ (left) and the density (right) at $t = 0$.</figcaption></figure>

The regions C, S, P are respectively a disk, a star shaped and a square portion of D. The EOS parameters of these fluids are given in table 10 along with their initial location within the domain and density value. The domain D is discretized over a 200 × 200-cell regular mesh and we impose periodic boundary conditions. The resulting initial condition for the material interface and the density is depicted in figure 17. Figure 18 displays a mapping of $\sum_{k=1}^{4} kZ_k$ at $t = 42.5$ s obtained with both the upwind scheme and the anti-diffusive scheme. The final instant is reached after 2626 time steps for the upwind scheme and 2627 time steps for the anti-diffusive scheme. We plot in figures 18 and 19 the density profiles and the interfaces by representing $\sum_{k=1}^{4} kZ_k$ at the final time $t = 42.5$ s. As in the one dimensional case, the interfaces remain quite sharp and their original topology is well preserved by the numerical approximation computed by the anti-diffusive scheme. With the upwind scheme, the interface shapes are completely spread due to numerical diffusion. Nevertheless, both schemes preserve velocity and pressure profiles as in the one-dimensional test. Indeed, at the final time the relative errors $||\bar{u} - u||_{L^1}$ and $||p - \bar{p}||_{L^1}$ for the references $\bar{u} = \sqrt{5}$ and $\bar{p} = 1$ remain bounded between $3.9 \times 10^{-16}$ and $4.9 \times 10^{-14}$ (see table 11 for more details). Finally, figure 20 displays the graphs $t \mapsto \int_D \sum_{k=1}^{4} Z_k(x, t)dx$ and $t \mapsto \int_D \sum_{k=1}^{4} Y_k(x, t)dx$ for the anti-diffusive scheme. Again, we observe that the scheme preserve the unity constraints (5) for both color functions and mass fractions at each time step.
<figure><img src="image_19.png" /><figcaption>Figure 18: test 4, two-dimensional passive transport with four components. Mapping of $\sum_{k=1}^{4} kZ_k$ at $t = 42.5$ s obtained with both upwind scheme and the anti-diffusive scheme.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: test 4, two-dimensional passive transport with four components. Mapping of the density at $t = 42.5$ s obtained with both upwind scheme and the anti-diffusive scheme.</figcaption></figure>

<table><caption>Table 11: test 4, two-dimensional passive transport with four components. Relative error of $p$ and $u$ in $L^1$ norm to the constant profile $\bar{p} = 1$ and $\bar{u} = \sqrt{5}$ at final time $t = 42.5$ s.</caption><thead><tr><th>Scheme</th><th>||<i>p</i> - <i>p</i>||<sub><i>L</i><sup>1</sup></sub></th><th>||<i>u</i> - <i>u</i>||<sub><i>L</i><sup>1</sup></sub></th></tr></thead><tbody><tr><td>Anti-diffusive</td><td>4.88 × 10<sup>-14</sup></td><td>3.97 × 10<sup>-16</sup></td></tr><tr><td>Upwind</td><td>3.99 × 10<sup>-14</sup></td><td>1.19 × 10<sup>-15</sup></td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Figure 18: test 4, two-dimensional passive transport with four components. Mapping of $\sum_{k=1}^{4} kZ_k$ at $t = 42.5$ s obtained with both upwind scheme and the anti-diffusive scheme.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 19: test 4, two-dimensional passive transport with four components. Mapping of the density at $t = 42.5$ s obtained with both upwind scheme and the anti-diffusive scheme.</figcaption></figure>

<table><caption>Table 11: test 4, two-dimensional passive transport with four components. Relative error of $p$ and $u$ in $L^1$ norm to the constant profile $\bar{p} = 1$ and $\bar{u} = \sqrt{5}$ at final time $t = 42.5$ s.</caption><thead><tr><th>Scheme</th><th>||<i>p</i> - <i>p</i>||<sub><i>L</i><sup>1</sup></sub></th><th>||<i>u</i> - <i>u</i>||<sub><i>L</i><sup>1</sup></sub></th></tr></thead><tbody><tr><td>Anti-diffusive</td><td>4.88 &times; 10<sup>-14</sup></td><td>3.97 &times; 10<sup>-16</sup></td></tr><tr><td>Upwind</td><td>3.99 &times; 10<sup>-14</sup></td><td>1.19 &times; 10<sup>-15</sup></td></tr></tbody></table>
<figure><img src="image_21.png" /><figcaption>Figure 20: test 4, two-dimensional passive transport with four components. Graph of $t \mapsto \int_D \sum_{k=1}^4 Z_k(x, t) dx$ (left) and $t \mapsto \int_D \sum_{k=1}^4 Y_k(x, t) dx$ (right) for the approximation computed with the anti-diffusive scheme.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21: test 5: triple point problem. Initial position of the interfaces in the computational domain.</figcaption></figure>

### 5.4.1. Test 5: Triple point problem

We propose here a two-dimensional problem that involves the interaction of three Riemann problems across three initial material discontinuities and a triple point. Similar tests have been studied in the literature, see e.g. [17, 27, 40, 58]. The computation domain is $[0, 7 \text{ m}] \times [0, 3 \text{ m}]$ and it is occupied by three perfect gases located as depicted in figure 21. The medium, is initially at rest and all initial states and parameters of each fluid are given in table 12.

At the initial instant, the pressure in the fluid $k = 1$, depicted in blue in figure 21, is greater than in the rest of the domain. This generates a set of waves, including two shocks travelling towards the right end of the domain: one of these waves travels within the fluid 2 (green in figure 21), the other within the fluid 3 (red in figure 21). The jump between the densities and the material

<table><caption>Table 12: test 5: triple point problem. Initial data and EOS parameters.</caption><thead><tr><th>k</th><th>ρ (kg.m<sup>-3</sup>)</th><th>p (Pa)</th><th>u (m/s, m/s)</th><th>γ<sub>k</sub></th></tr></thead><tbody><tr><td>k = 1</td><td>1.0</td><td>1.0</td><td>(0, 0)</td><td>1.6</td></tr><tr><td>k = 2</td><td>0.125</td><td>0.1</td><td>(0, 0)</td><td>1.5</td></tr><tr><td>k = 3</td><td>1.0</td><td>0.1</td><td>(0, 0)</td><td>1.4</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>Figure 20: test 4, two-dimensional passive transport with four components. Graph of $t \mapsto \int_D \sum_{k=1}^4 Z_k(x, t) dx$ (left) and $t \mapsto \int_D \sum_{k=1}^4 Y_k(x, t) dx$ (right) for the approximation computed with the anti-diffusive scheme.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 21: test 5: triple point problem. Initial position of the interfaces in the computational domain.</figcaption></figure>

### 5.4.1. Test 5: Triple point problem

We propose here a two-dimensional problem that involves the interaction of three Riemann problems across three initial material discontinuities and a triple point. Similar tests have been studied in the literature, see e.g. [17, 27, 40, 58]. The computation domain is $[0, 7 \text{ m}] \times [0, 3 \text{ m}]$ and it is occupied by three perfect gases located as depicted in figure 21. The medium, is initially at rest and all initial states and parameters of each fluid are given in table 12.

At the initial instant, the pressure in the fluid $k = 1$, depicted in blue in figure 21, is greater than in the rest of the domain. This generates a set of waves, including two shocks travelling towards the right end of the domain: one of these waves travels within the fluid 2 (green in figure 21), the other within the fluid 3 (red in figure 21). The jump between the densities and the material

<table><caption>Table 12: test 5: triple point problem. Initial data and EOS parameters.</caption><thead><tr><th>k</th><th>ρ (kg.m<sup>-3</sup>)</th><th>p (Pa)</th><th>u (m/s, m/s)</th><th>γ<sub>k</sub></th></tr></thead><tbody><tr><td>k = 1</td><td>1.0</td><td>1.0</td><td>(0, 0)</td><td>1.6</td></tr><tr><td>k = 2</td><td>0.125</td><td>0.1</td><td>(0, 0)</td><td>1.5</td></tr><tr><td>k = 3</td><td>1.0</td><td>0.1</td><td>(0, 0)</td><td>1.4</td></tr></tbody></table>
properties across the material interface that separates fluid 2 and 3 creates an instability. For this simulation we use a 700 × 300-cell regular mesh and impose wall boundary conditions. Figures 22, 23 and 24 show respectively the color function, the density and pressure profiles at t = 5 s obtained with both the upwind and the anti-diffusive scheme. The numerical solution obtained is in good agreement with those found in the literature (see for example [27, 28]). Here again the interfaces computed by the anti-diffusive scheme remain very sharp. Finally, we observe as expected no significant difference between the pressure profiles obtained with both schemes as in [27].

<figure><img src="image_23.png" /><figcaption>Figure 22: test 5: triple point problem. Profiles of color functions for both upwind (top) and anti-diffusive (bottom) schemes at t = 5 s.</figcaption></figure>

### 5.5. Test 6: Two-dimensional shock-bubble interaction with three materials

We now consider a test inspired by the experimental results of [32]. This problem has been studied in several publications with various approaches (see e.g. [7, 27, 41]). We consider here a slightly modified set up that involves three gases: air, R22 (Chlorodifluoromethane) and Helium, that are represented here by perfect gases. The computational domain is a rectangular region P whose dimensions are $L_1 \times L_2$. Let $(X_1, X_2)$ be a point of the domain, the regions

$$ \mathcal{E} = \{(x_1, x_2) \in P \mid (x_1 - X_1)^2 + (x_2 - X_2)^2 < r_2\}, $$

$$ \mathcal{F} = \{(x_1, x_2) \in P \mid r_2 < (x_1 - X_1)^2 + (x_2 - X_2)^2 < r_1\} $$

are respectively a disc and a ring of center $(X_1, X_2)$. The region E is filled with Helium and F contains R22. The gas in the rest of the domain is air. At
properties across the material interface that separates fluid 2 and 3 creates an instability. For this simulation we use a 700 × 300-cell regular mesh and impose wall boundary conditions. Figures 22, 23 and 24 show respectively the color function, the density and pressure profiles at t = 5 s obtained with both the upwind and the anti-diffusive scheme. The numerical solution obtained is in good agreement with those found in the literature (see for example [27, 28]). Here again the interfaces computed by the anti-diffusive scheme remain very sharp. Finally, we observe as expected no significant difference between the pressure profiles obtained with both schemes as in [27].

<figure><img src="image_22.png" /><figcaption>Figure 22: test 5: triple point problem. Profiles of color functions for both upwind (top) and anti-diffusive (bottom) schemes at t = 5 s.</figcaption></figure>

### 5.5. Test 6: Two-dimensional shock-bubble interaction with three materials

We now consider a test inspired by the experimental results of [32]. This problem has been studied in several publications with various approaches (see e.g. [7, 27, 41]). We consider here a slightly modified set up that involves three gases: air, R22 (Chlorodifluoromethane) and Helium, that are represented here by perfect gases. The computational domain is a rectangular region P whose dimensions are $L_1 \times L_2$. Let $(X_1, X_2)$ be a point of the domain, the regions

$$ \mathcal{E} = \{(x_1, x_2) \in P \mid (x_1 - X_1)^2 + (x_2 - X_2)^2 < r_2\}, \\ \mathcal{F} = \{(x_1, x_2) \in P \mid r_2 < (x_1 - X_1)^2 + (x_2 - X_2)^2 < r_1\} $$

are respectively a disc and a ring of center $(X_1, X_2)$. The region E is filled with Helium and F contains R22. The gas in the rest of the domain is air. At
<figure><img src="image_24.png" /><figcaption>Figure 23: test 5: triple point problem. Profiles of density for both upwind (top) and anti-diffusive (bottom) schemes at t = 5 s.</figcaption></figure>

the beginning of the computation, a Mach 1.22 shock with initial position $x_s$ is travels through air towards the left end of the domain. For the present test, we use the following values

$$L_1 = 0.445 \text{ m}, \quad L_2 = 0.089 \text{ m}, \quad r_1 = 0.025 \text{ m}, \quad r_2 = 0.015 \text{ m}, \\ X_{\text{shock}} = 0.275 \text{ m}, \quad X_1 = 0.225 \text{ m}, \quad X_2 = 0.0445 \text{ m}.$$

The initial state of the fluids and the EOS parameters we chose for the fluids are given in table 13. We use a regular mesh composed by 1250 × 250 cells for

<table><caption>Table 13: test 6, two-dimensional shock-bubble interaction with three materials. EOS parameters and initial data.</caption><thead><tr><th>material</th><th>ρ (kg.m<sup>-3</sup>)</th><th>p (Pa)</th><th>u (m.s<sup>-1</sup>, m.s<sup>-1</sup>)</th><th>γ<sub>k</sub></th></tr></thead><tbody><tr><td>Helium</td><td>0.138</td><td>1.01325 10<sup>5</sup></td><td>(0, 0)</td><td>1.6</td></tr><tr><td>R22</td><td>3.863</td><td>1.01325 10<sup>5</sup></td><td>(0, 0)</td><td>1.249</td></tr><tr><td>Air (pre-shock)</td><td>1.686</td><td>1.59 10<sup>5</sup></td><td>(0, 0)</td><td>1.4</td></tr><tr><td>Air (post-shock)</td><td>1.225</td><td>1.01325 10<sup>5</sup></td><td>(-113.5, 0)</td><td>1.4</td></tr></tbody></table>

discretizing $P$, and the domain is closed using wall boundaries. The position of the material interfaces obtained with both the anti-diffusive scheme and the upwind scheme are displayed in figure 26 at several instants. The shock hits the air/R22 interface at $t \simeq 60 \mu s$. After the impact the bubble and the gas ring are
<figure><img src="image_23.png" /><figcaption>Figure 23: test 5: triple point problem. Profiles of density for both upwind (top) and anti-diffusive (bottom) schemes at t = 5 s.</figcaption></figure>

the beginning of the computation, a Mach 1.22 shock with initial position $x_s$ is travels through air towards the left end of the domain. For the present test, we use the following values

$$L_1 = 0.445 \text{ m}, \quad L_2 = 0.089 \text{ m}, \quad r_1 = 0.025 \text{ m}, \quad r_2 = 0.015 \text{ m}, \\ X_{\text{shock}} = 0.275 \text{ m}, \quad X_1 = 0.225 \text{ m}, \quad X_2 = 0.0445 \text{ m}.$$

The initial state of the fluids and the EOS parameters we chose for the fluids are given in table 13. We use a regular mesh composed by 1250 × 250 cells for

<table><caption>Table 13: test 6, two-dimensional shock-bubble interaction with three materials. EOS parameters and initial data.</caption><thead><tr><th>material</th><th>ρ (kg.m<sup>-3</sup>)</th><th>p (Pa)</th><th>u (m.s<sup>-1</sup>, m.s<sup>-1</sup>)</th><th>γ<sub>k</sub></th></tr></thead><tbody><tr><td>Helium</td><td>0.138</td><td>1.01325 10<sup>5</sup></td><td>(0, 0)</td><td>1.6</td></tr><tr><td>R22</td><td>3.863</td><td>1.01325 10<sup>5</sup></td><td>(0, 0)</td><td>1.249</td></tr><tr><td>Air (pre-shock)</td><td>1.686</td><td>1.59 10<sup>5</sup></td><td>(0, 0)</td><td>1.4</td></tr><tr><td>Air (post-shock)</td><td>1.225</td><td>1.01325 10<sup>5</sup></td><td>(-113.5, 0)</td><td>1.4</td></tr></tbody></table>

discretizing $P$, and the domain is closed using wall boundaries. The position of the material interfaces obtained with both the anti-diffusive scheme and the upwind scheme are displayed in figure 26 at several instants. The shock hits the air/R22 interface at $t \simeq 60 \ \mu\text{s}$. After the impact the bubble and the gas ring are
<figure><img src="image_25.png" /><figcaption>Figure 24: test 5: triple point problem. Profiles of pressure for both upwind (top) and anti-diffusive (bottom) schemes at t = 5 s.</figcaption></figure>

set in motion towards the left end of the domain. Both R22 and Helium bulk are compressed by the shock. Their shape is deformed by the motion so that two vortices appear on the right side of the gas bulk. The upwind scheme and the anti-diffusive scheme provide similar position of the material fronts, although the material interfaces obtained thanks to the anti-diffusive scheme are much sharper.

### 5.6. Test 7: Two-dimensional Kelvin-Helmoltz instability

In this last part, we want to test the capability of the proposed method to capture complex modification of the interface geometry between an arbitrary number of materials. To this end, we are interested in computing a multi-material Kelvin-Helmoltz instability inspired by [27, 62].

We consider a simple two-dimensional shear flow in a square periodic domain of dimension 1m × 1m. Here, we suppose that the domain is initially occupied by three perfect gases numbered k = 1, ..., 3 as depicted in figure 27. At t = 0, we assume the pressure constant and the lighter gas (in blue) is initially moving from the right to the left while the heavier ones (green and red) are moving in the opposite direction. The initial state of the fluid is given in the table 14. The interface is then perturbed by an imposed vertical velocity $\phi$ field that reads

$$ \phi(x_1, x_2) = \omega_0 \sin(4\pi x) \left[ e^{\frac{-(x_2-0.25)^2}{2\sigma^2}} + e^{\frac{-(x_2-0.75)^2}{2\sigma^2}} \right] $$
<figure><img src="image_24.png" /><figcaption>Figure 24: test 5: triple point problem. Profiles of pressure for both upwind (top) and anti-diffusive (bottom) schemes at t = 5 s.</figcaption></figure>

set in motion towards the left end of the domain. Both R22 and Helium bulk are compressed by the shock. Their shape is deformed by the motion so that two vortices appear on the right side of the gas bulk. The upwind scheme and the anti-diffusive scheme provide similar position of the material fronts, although the material interfaces obtained thanks to the anti-diffusive scheme are much sharper.

### 5.6. Test 7: Two-dimensional Kelvin-Helmoltz instability

In this last part, we want to test the capability of the proposed method to capture complex modification of the interface geometry between an arbitrary number of materials. To this end, we are interested in computing a multi-material Kelvin-Helmoltz instability inspired by [27, 62].

We consider a simple two-dimensional shear flow in a square periodic domain of dimension 1m × 1m. Here, we suppose that the domain is initially occupied by three perfect gases numbered k = 1, ..., 3 as depicted in figure 27. At t = 0, we assume the pressure constant and the lighter gas (in blue) is initially moving from the right to the left while the heavier ones (green and red) are moving in the opposite direction. The initial state of the fluid is given in the table 14. The interface is then perturbed by an imposed vertical velocity $\phi$ field that reads

$$ \phi(x_1, x_2) = \omega_0 \sin(4\pi x) \left[ e^{\frac{-(x_2-0.25)^2}{2\sigma^2}} + e^{\frac{-(x_2-0.75)^2}{2\sigma^2}} \right] $$
<figure><img src="image_26.png" /><figcaption>Figure 25: test 6, two-dimensional shock-bubble interaction with three materials. Mapping of $\sum_{k=1}^{3} kZ_k$ for both upwind (right) and anti-diffusive (left) schemes at different times.</figcaption></figure>

where $\omega_0 = 0.1$ and $\sigma = 0.05/\sqrt{2}$. Such a choice guarantees that only a single mode will dominate the linear growth of the instability that will develop at the interfaces. The domain is discretized by means of a 1000 x 1000-cell regular mesh.

Figure 28 represents a mapping of the computed color functions at different times between $t = 0$s and $t = 2$s for both anti-diffusive and upwind schemes. As expected in the references, the Kelvin-Helmoltz instability is illustrated by the formation for both schemes of four symmetric vertices. Once again, all the fine structures are sharply captured with the anti-diffusive scheme while it becomes very difficult to distinguish the interfaces with the upwind scheme due to numerical diffusion.

## 6. Conclusion

We propose in the present paper a method for simulating compressible flows involving an arbitrary number $m$ of materials separated by sharp interfaces.

Our approach relies on two elements. First, we introduced a $m$-component flow model that extends the two-component five-equation model with isobaric
<figure><img src="image_25.png" /><figcaption>Figure 25: test 6, two-dimensional shock-bubble interaction with three materials. Mapping of $\sum_{k=1}^{3} kZ_k$ for both upwind (right) and anti-diffusive (left) schemes at different times.</figcaption></figure>

where $\omega_0 = 0.1$ and $\sigma = 0.05/\sqrt{2}$. Such a choice guarantees that only a single mode will dominate the linear growth of the instability that will develop at the interfaces. The domain is discretized by means of a 1000 x 1000-cell regular mesh.

Figure 28 represents a mapping of the computed color functions at different times between $t = 0$s and $t = 2$s for both anti-diffusive and upwind schemes. As expected in the references, the Kelvin-Helmoltz instability is illustrated by the formation for both schemes of four symmetric vertices. Once again, all the fine structures are sharply captured with the anti-diffusive scheme while it becomes very difficult to distinguish the interfaces with the upwind scheme due to numerical diffusion.

## 6. Conclusion

We propose in the present paper a method for simulating compressible flows involving an arbitrary number $m$ of materials separated by sharp interfaces.

Our approach relies on two elements. First, we introduced a $m$-component flow model that extends the two-component five-equation model with isobaric
<figure><img src="image_27.png" /><figcaption>Figure 26: test 6, two-dimensional shock-bubble interaction with three materials. Density profiles for both upwind (right) and anti-diffusive (left) schemes at different times.</figcaption></figure>

<table><caption>Table 14: test 7, two-dimensional Kelvin-Helmoltz instability. Initial data and EOS parameters.</caption><thead><tr><th>material index</th><th>density (kg.m<sup>-3</sup>)</th><th>pressure (Pa)</th><th>velocity (m.s<sup>-1</sup>, m.s<sup>-1</sup>)</th><th>&gamma;<sub>k</sub></th></tr></thead><tbody><tr><td>k = 1</td><td>1</td><td>2.5</td><td>(-0.5, 0)</td><td>5/3</td></tr><tr><td>k = 2</td><td>2</td><td>2.5</td><td>( 0.5, 0)</td><td>1.4</td></tr><tr><td>k = 3</td><td>2</td><td>2.5</td><td>( 0.5, 0)</td><td>2.4</td></tr></tbody></table>

closure of [3, 4] to *m* components. Under simple thermodynamical assumptions this model is hyperbolic and may be expressed equivalently with a quasi-conservative form or a fully conservative form. This model is adapted to a
<figure><img src="image_26.png" /><figcaption>Figure 26: test 6, two-dimensional shock-bubble interaction with three materials. Density profiles for both upwind (right) and anti-diffusive (left) schemes at different times.</figcaption></figure>

<table><caption>Table 14: test 7, two-dimensional Kelvin-Helmoltz instability. Initial data and EOS parameters.</caption><thead><tr><th>material index</th><th>density (kg.m<sup>-3</sup>)</th><th>pressure (Pa)</th><th>velocity (m.s<sup>-1</sup>, m.s<sup>-1</sup>)</th><th>&gamma;<sub>k</sub></th></tr></thead><tbody><tr><td>k = 1</td><td>1</td><td>2.5</td><td>(-0.5, 0)</td><td>5/3</td></tr><tr><td>k = 2</td><td>2</td><td>2.5</td><td>(0.5, 0)</td><td>1.4</td></tr><tr><td>k = 3</td><td>2</td><td>2.5</td><td>(0.5, 0)</td><td>2.4</td></tr></tbody></table>

closure of [3, 4] to *m* components. Under simple thermodynamical assumptions this model is hyperbolic and may be expressed equivalently with a quasi-conservative form or a fully conservative form. This model is adapted to a
<figure><img src="image_28.png" /><figcaption>Figure 27: test 7, two-dimensional Kelvin-Helmoltz instability. Initial interface positions.</figcaption></figure>

wide range of EOSs, including tabulated laws. Second, we proposed an anti-diffusive Lagrange-Remap Finite Volume scheme combining the lines of [41] and [22]. This numerical scheme succeeds in controlling the numerical diffusion of the material fronts without any interface reconstruction process. The overall numerical scheme is conservative with respect to the mass of each component, the total mass, the momentum and the total energy. It is also endowed with stability properties under a classic CFL condition. Moreover we also proved and tested the ability of our discretization to preserve a category of solutions with constant velocity and pressure.

Our simulation strategy was tested against several one-dimensional and two-dimensional test cases involving up to five distinct materials and complex EOSs. The results show that the anti-diffusive mechanism succeeds in evolving the interfaces while preserving their sharpness. The accuracy of the overall method is enhanced at the vicinity of the contact discontinuities that drive the material fronts, while far from the interfaces the behavior of the scheme if similar to a classic upwind discretization.

The simulation approach we proposed here has been implemented in the code HERACLES [34]. Further developments will deal with enhancing the accuracy of the solver by using high order methods as in [8] and propose an extension to unstructured meshes following the lines of [22, 25].

## Appendix A. Hyperbolicity and eigenstructure analysis

This section details the hyperbolicity study of the $m$-multicomponent model as stated in proposition 3. To this goal, we start by considering some companion evolution equations of the system (8). First, combining the partial mass conservation equations (8c) we get the conservation equation for the total mass

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0. \qquad (A.1) $$
<figure><img src="image_27.png" /><figcaption>Figure 27: test 7, two-dimensional Kelvin-Helmoltz instability. Initial interface positions.</figcaption></figure>

wide range of EOSs, including tabulated laws. Second, we proposed an anti-diffusive Lagrange-Remap Finite Volume scheme combining the lines of [41] and [22]. This numerical scheme succeeds in controlling the numerical diffusion of the material fronts without any interface reconstruction process. The overall numerical scheme is conservative with respect to the mass of each component, the total mass, the momentum and the total energy. It is also endowed with stability properties under a classic CFL condition. Moreover we also proved and tested the ability of our discretization to preserve a category of solutions with constant velocity and pressure.

Our simulation strategy was tested against several one-dimensional and two-dimensional test cases involving up to five distinct materials and complex EOSs. The results show that the anti-diffusive mechanism succeeds in evolving the interfaces while preserving their sharpness. The accuracy of the overall method is enhanced at the vicinity of the contact discontinuities that drive the material fronts, while far from the interfaces the behavior of the scheme if similar to a classic upwind discretization.

The simulation approach we proposed here has been implemented in the code HERACLES [34]. Further developments will deal with enhancing the accuracy of the solver by using high order methods as in [8] and propose an extension to unstructured meshes following the lines of [22, 25].

## Appendix A. Hyperbolicity and eigenstructure analysis

This section details the hyperbolicity study of the $m$-multicomponent model as stated in proposition 3. To this goal, we start by considering some companion evolution equations of the system (8). First, combining the partial mass conservation equations (8c) we get the conservation equation for the total mass

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0. \qquad (A.1) $$
<figure><img src="image_29.png" /><figcaption>Figure 28: test 7, two-dimensional Kelvin-Helmoltz instability. Mapping of $\sum_{k=1}^{3} k\mathcal{Z}_k$ for both upwind (right) and anti-diffusive (left) schemes at different times.</figcaption></figure>

Now, introducing the material derivative $D_t$ defined for any fluid parameter $(\mathbf{x}, t) \mapsto a$ by

$$D_t a = \partial_t a + \mathbf{u} \cdot \nabla a.$$

the system (8) can be classically recast into the following form

$$\rho D_t \mathbf{u} + \nabla p = 0, \qquad (A.2a)$$

$$\rho D_t E + \nabla \cdot (p\mathbf{u}) = 0, \qquad (A.2b)$$

$$D_t(\rho_k \mathcal{Z}_k) + \rho_k \mathcal{Z}_k (\nabla \cdot \mathbf{u}) = 0, \qquad k = 1, \dots, m, \qquad (A.2c)$$

$$D_t \mathcal{Z}_k = 0, \qquad k = 1, \dots, m. \qquad (A.2d)$$

either by assuming that $(\mathbf{x}, t) \mapsto \mathbf{u}$ is a regular function or by considering that the formulation (A.2) is a formal expression. Next, the medium specific internal
<figure><img src="image_28.png" /><figcaption>Figure 28: test 7, two-dimensional Kelvin-Helmoltz instability. Mapping of $\sum_{k=1}^{3} k\mathcal{Z}_k$ for both upwind (right) and anti-diffusive (left) schemes at different times.</figcaption></figure>

Now, introducing the material derivative $D_t$ defined for any fluid parameter $(x, t) \mapsto a$ by

$$D_t a = \partial_t a + \mathbf{u} \cdot \nabla a.$$

the system (8) can be classically recast into the following form

$$\rho D_t \mathbf{u} + \nabla p = 0, \qquad (A.2a)$$

$$\rho D_t E + \nabla \cdot (p\mathbf{u}) = 0, \qquad (A.2b)$$

$$D_t(\rho_k \mathcal{Z}_k) + \rho_k \mathcal{Z}_k (\nabla \cdot \mathbf{u}) = 0, \qquad k = 1, \dots, m, \qquad (A.2c)$$

$$D_t \mathcal{Z}_k = 0, \qquad k = 1, \dots, m. \qquad (A.2d)$$

either by assuming that $(x, t) \mapsto \mathbf{u}$ is a regular function or by considering that the formulation (A.2) is a formal expression. Next, the medium specific internal
energy e verifies

$$ \rho D_t e + p(\nabla \cdot \mathbf{u}) = 0, \qquad (A.3) $$

where the pressure p is supposed to be uniquely defined by (9). Finally, we introduce the mass fraction $\mathcal{Y}_k$ of fluid $k = 1, \dots, m$ defined by

$$ \rho \mathcal{Y}_k = \rho_k \mathcal{Z}_k, \quad k = 1, \dots, m, \qquad (A.4) $$

and using (A.2c) we obtain that

$$ D_t \mathcal{Y}_k = 0, \quad k = 1, \dots, m. \qquad (A.5) $$

**Proposition 5.** The pressure p verifies

$$ D_t p = -\rho c^2 (\nabla \cdot \mathbf{u}), \qquad (A.6) $$

where $c > 0$ is the sound velocity of the m-component medium that is defined by

$$ \rho \xi c^2 = \sum_k \mathcal{Z}_k \rho_k \xi_k c_k^2. \qquad (A.7) $$

*Proof.* The specific internal energy can be expressed thanks to the phasic internal energy and the mass fractions $\mathcal{Y}_k$. Indeed,

$$ e = \sum_k \mathcal{Y}_k e_k(\rho_k, p), $$

thus

$$ \begin{aligned} \rho D_t e &= \sum_k \rho \mathcal{Y}_k \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} D_t \rho_k + \left[ \sum_k \rho \mathcal{Y}_k \left( \frac{\partial e_k}{\partial p_k} \right)_{\rho_k} \right] D_t p \\ &= \sum_k \rho_k \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} D_t (\mathcal{Z}_k \rho_k) + \left[ \sum_k \mathcal{Z}_k \rho_k \left( \frac{\partial e_k}{\partial p_k} \right)_{\rho_k} \right] D_t p \\ &= - \sum_k \mathcal{Z}_k \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} (\nabla \cdot \mathbf{u}) + \left[ \sum_k \mathcal{Z}_k \left( \frac{\partial \rho_k e_k}{\partial p_k} \right)_{\rho_k} \right] D_t p, \end{aligned} $$

which leads to

$$ \rho D_t e = -(\nabla \cdot \mathbf{u}) \sum_k \mathcal{Z}_k \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} + D_t p \sum_k \mathcal{Z}_k \xi_k. $$

Using (7) and (A.3) we obtain

$$ \xi D_t p + (\nabla \cdot \mathbf{u}) \sum_{k=1}^{m} \mathcal{Z}_k \left[ p - \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} \right] = 0. $$

One notes that (3) provides also the identity

$$ \rho_k \xi_k c_k^2 = - \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} + p_k, \quad k = 1, \dots, m, $$
energy e verifies

$$ \rho D_t e + p(\nabla \cdot \mathbf{u}) = 0, \qquad (A.3) $$

where the pressure p is supposed to be uniquely defined by (9). Finally, we introduce the mass fraction $\mathcal{Y}_k$ of fluid $k = 1, \dots, m$ defined by

$$ \rho \mathcal{Y}_k = \rho_k \mathcal{Z}_k, \quad k = 1, \dots, m, \qquad (A.4) $$

and using (A.2c) we obtain that

$$ D_t \mathcal{Y}_k = 0, \quad k = 1, \dots, m. \qquad (A.5) $$

**Proposition 5.** The pressure p verifies

$$ D_t p = -\rho c^2 (\nabla \cdot \mathbf{u}), \qquad (A.6) $$

where $c > 0$ is the sound velocity of the m-component medium that is defined by

$$ \rho \xi c^2 = \sum_k \mathcal{Z}_k \rho_k \xi_k c_k^2. \qquad (A.7) $$

*Proof.* The specific internal energy can be expressed thanks to the phasic internal energy and the mass fractions $\mathcal{Y}_k$. Indeed,

$$ e = \sum_k \mathcal{Y}_k e_k(\rho_k, p), $$

thus

$$ \begin{aligned} \rho D_t e &= \sum_k \rho \mathcal{Y}_k \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} D_t \rho_k + \left[ \sum_k \rho \mathcal{Y}_k \left( \frac{\partial e_k}{\partial p_k} \right)_{\rho_k} \right] D_t p \\ &= \sum_k \rho_k \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} D_t (\mathcal{Z}_k \rho_k) + \left[ \sum_k \mathcal{Z}_k \rho_k \left( \frac{\partial e_k}{\partial p_k} \right)_{\rho_k} \right] D_t p \\ &= - \sum_k \mathcal{Z}_k \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} (\nabla \cdot \mathbf{u}) + \left[ \sum_k \mathcal{Z}_k \left( \frac{\partial \rho_k e_k}{\partial p_k} \right)_{\rho_k} \right] D_t p, \end{aligned} $$

which leads to

$$ \rho D_t e = -(\nabla \cdot \mathbf{u}) \sum_k \mathcal{Z}_k \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} + D_t p \sum_k \mathcal{Z}_k \xi_k. $$

Using (7) and (A.3) we obtain

$$ \xi D_t p + (\nabla \cdot \mathbf{u}) \sum_{k=1}^{m} \mathcal{Z}_k \left[ p - \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} \right] = 0. $$

One notes that (3) provides also the identity

$$ \rho_k \xi_k c_k^2 = - \rho_k^2 \left( \frac{\partial e_k}{\partial \rho_k} \right)_{\rho_k} + p_k, \quad k = 1, \dots, m, $$
then

$$ \xi D_t p + (\nabla \cdot \mathbf{u}) \sum_{k=1}^{m} \mathcal{Z}_k \rho_k \xi_k c_k^2 = 0, $$

which completes the proof. □

We now turn back to the system (14). Under the assumptions of proposition 3 the matrix A satisfies the following properties

a) the matrix $\mathbf{A}(\mathbf{V})$ is diagonalizable,

b) the eigenvalues $(\lambda_k)_{k=1,...,2m+2}$ of $\mathbf{A}(\mathbf{V})$ are the real numbers

$$ \lambda_1 = u - c, \quad \lambda_2 = u + c, \quad \lambda_3 = \dots = \lambda_{2m+2} = u, $$

c) the right eigenvectors $(\mathbf{r}_k)_{k=1,...,2m+2}$ associated with the eigenvalues $(\lambda_k)_{k=1,...,2m+2}$ are

$$ \mathbf{r}_1 = (1, -c\rho, 0, 0, \dots, 0)^T, $$

$$ \mathbf{r}_2 = (1, c\rho, 0, 0, \dots, 0)^T, $$

$$ \mathbf{r}_k = (0, \dots, 0, \underbrace{1}_{k^{\text{th}} \text{ coordinate}}, 0, \dots, 0)^T, \quad k = 3, \dots, 2m+2, $$

d) the fields associated with the eigenvalues $\lambda_1$ and $\lambda_2$ are genuinely non-linear, while those associated with $\lambda_k$ are linearly degenerate for $k = 3, \dots, 2m+2$.

As a consequence, the system (8) is hyperbolic.

Finally we examine the definition of the non-conservative product $\mathbf{u} \cdot \nabla \mathcal{Z}_k$. System (8) can be expressed as a conservation law using the variables $\rho \mathbf{u}$, $\rho E$, $\rho_k \mathcal{Z}_k$, $\rho \mathcal{Z}_k$, $k = 0, \dots, m$. Under this form, when one considers a pure one-dimensional discontinuous solution of the system, a classical set of Rankine-Hugoniot conditions is available. In particular, the jump condition associated with $\rho \mathcal{Z}_k$ for $k = 0, \dots, m$ ensures that the color function variable $\mathcal{Z}_k$ and the material velocity cannot simultaneously experience a jump across the discontinuity.

## Appendix B. Expression of the m-component system in Lagrangian Coordinates

In this section we present the expression of the m-component system (8) using Lagrangian coordinates for the case of one-dimensional problems. This will be useful for the derivation of the Lagrange-Projection algorithm we aim to elaborate in the next sections for solving (8). We consider the mapping $(X, t) \in \mathbb{R} \times [0, +\infty) \mapsto \chi$ defined by the ordinary differential equation

$$ \partial\chi/\partial t = u(\chi(X, t), t), \quad \chi(X, t = 0) = X. $$
then

$$ \xi D_t p + (\nabla \cdot \mathbf{u}) \sum_{k=1}^{m} Z_k \rho_k \xi_k c_k^2 = 0, $$

which completes the proof. □

We now turn back to the system (14). Under the assumptions of proposition 3 the matrix A satisfies the following properties

a) the matrix $\mathbf{A}(\mathbf{V})$ is diagonalizable,

b) the eigenvalues $(\lambda_k)_{k=1,...,2m+2}$ of $\mathbf{A}(\mathbf{V})$ are the real numbers

$$ \lambda_1 = u - c, \quad \lambda_2 = u + c, \quad \lambda_3 = \dots = \lambda_{2m+2} = u, $$

c) the right eigenvectors $(\mathbf{r}_k)_{k=1,...,2m+2}$ associated with the eigenvalues $(\lambda_k)_{k=1,...,2m+2}$ are

$$ \mathbf{r}_1 = (1, -c\rho, 0, 0, \dots, 0)^T, $$

$$ \mathbf{r}_2 = (1, c\rho, 0, 0, \dots, 0)^T, $$

$$ \mathbf{r}_k = (0, \dots, 0, \underbrace{1}_{k^{\text{th}} \text{ coordinate}}, 0, \dots, 0)^T, \quad k = 3, \dots, 2m+2, $$

d) the fields associated with the eigenvalues $\lambda_1$ and $\lambda_2$ are genuinely non-linear, while those associated with $\lambda_k$ are linearly degenerate for $k = 3, \dots, 2m+2$.

As a consequence, the system (8) is hyperbolic.

Finally we examine the definition of the non-conservative product $\mathbf{u} \cdot \nabla Z_k$. System (8) can be expressed as a conservation law using the variables $\rho \mathbf{u}$, $\rho E$, $\rho_k Z_k$, $\rho Z_k$, $k = 0, \dots, m$. Under this form, when one considers a pure one-dimensional discontinuous solution of the system, a classical set of Rankine-Hugoniot conditions is available. In particular, the jump condition associated with $\rho Z_k$ for $k = 0, \dots, m$ ensures that the color function variable $Z_k$ and the material velocity cannot simultaneously experience a jump across the discontinuity.

## Appendix B. Expression of the $m$-component system in Lagrangian Coordinates

In this section we present the expression of the $m$-component system (8) using Lagrangian coordinates for the case of one-dimensional problems. This will be useful for the derivation of the Lagrange-Projection algorithm we aim to elaborate in the next sections for solving (8). We consider the mapping $(X, t) \in \mathbb{R} \times [0, +\infty) \mapsto \chi$ defined by the ordinary differential equation

$$ \partial\chi/\partial t = u(\chi(X, t), t), \quad \chi(X, t = 0) = X. $$
The Lagrangian system of coordinates $(X, t)$ attached to the initial instant $t = 0$ expresses that a fluid particle located at $x = \chi(X, t)$ at instant $t$ was originally at position $X$ at instant $t = 0$. This means that $X = \chi^{-1}(x, t)$. Now, let $(x, t) \mapsto q$ be a mapping that represents a medium parameter in Eulerian coordinates, we can define a Lagrangian representation $(X, t) \mapsto q^{\text{Lag}}$ of this parameter by setting $q^{\text{Lag}}(X, t) = q(\chi(X, t), t)$. If we note $\rho^{\text{Lag}}(X, t = 0) = \rho_0^{\text{Lag}}(X)$, the system (8) for one-dimensional problems is equivalent to

$$ \rho_0^{\text{Lag}} \partial_t (1/\rho^{\text{Lag}}) - \partial_X u^{\text{Lag}} = 0 \qquad (\text{B.1a}) $$

$$ \rho_0^{\text{Lag}} \partial_t u^{\text{Lag}} + \partial_X p^{\text{Lag}} = 0 \qquad (\text{B.1b}) $$

$$ \rho_0^{\text{Lag}} \partial_t E^{\text{Lag}} + \partial_X (p^{\text{Lag}} u^{\text{Lag}}) = 0 \qquad (\text{B.1c}) $$

$$ \rho_0^{\text{Lag}} \partial_t \mathcal{Y}_k^{\text{Lag}} = 0, \quad k = 1, \dots, m-1, \qquad (\text{B.1d}) $$

$$ \rho_0^{\text{Lag}} \partial_t \mathcal{Z}_k^{\text{Lag}} = 0, \quad k = 1, \dots, m-1. \qquad (\text{B.1e}) $$

## Appendix C. Preservation of iso-pressure and iso-velocity profiles

As in [41], we prove that the anti-diffusive scheme introduced in section 3 preserves a family of iso-pressure and iso-velocity profiles. Let us introduce the primitive variable $\mathbf{V}^n = (u^n, E^n, \mathcal{Y}_1^n, \dots, \mathcal{Y}_m^n)$ for one-dimensional problems.

**Proposition 6.** Let $(\rho^n \mathbf{V}^n, \tilde{\mathcal{Z}}_1^n, \dots, \tilde{\mathcal{Z}}_m^n)$ be the approximate state variable at instant $t = t^n$. Suppose that $(\rho^n \mathbf{V}^n, \tilde{\mathcal{Z}}_1^n, \dots, \tilde{\mathcal{Z}}_m^n)$ is a contact discontinuity with both uniform velocity and pressure such that

$$ p_i^n = \bar{p}, \quad u_i^n = \bar{u}, \quad (\rho_k)_i^n = \bar{\rho}_k, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

Then the approximate solution $(\rho^{n+1} \mathbf{V}^{n+1}, \tilde{\mathcal{Z}}_1^{n+1}, \dots, \tilde{\mathcal{Z}}_m^{n+1})$ computed with the anti-diffusive scheme satisfies

$$ p_i^{n+1} = \bar{p}, \quad u_i^{n+1} = \bar{u}, \quad (\rho_k)_i^{n+1} = \bar{\rho}_k, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

*Proof.* Let us first consider the following variable $(\tilde{\rho} \tilde{\mathbf{V}}^n, \tilde{\mathcal{Z}}_1^n, \dots, \tilde{\mathcal{Z}}_m^n)$ deduced from the Lagrangian step. Since the phasic densities and the pressure are constant at $t = t^n$, we have

$$ (\rho_k e_k)_i^n = \overline{\rho_k e_k}, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. \qquad (\text{C.1}) $$

In addition as the velocity is constant, the Lagrangian fluxes needed in the acoustic scheme are equal to zero i.e. $[\rho^n]_i = 0, [\rho^n u^n]_i = 0$ and $L_i = 1$. Thereby, we have $\forall i \in \mathbb{Z}$

$$ \left\{ \begin{array}{l} (\tilde{\rho}u)_i = (\rho u)_i^n, \\ (\tilde{\rho}E)_i = (\rho E)_i^n, \\ (\tilde{\rho}_k \tilde{\mathcal{Z}}_k)_i = (\rho_k \mathcal{Z}_k)_i^n, \quad 1 \le k \le m, \\ \tilde{\mathcal{Z}}_{k,i} = \mathcal{Z}_{k,i}^n, \quad 1 \le k \le m. \end{array} \right. $$
The Lagrangian system of coordinates $(X, t)$ attached to the initial instant $t = 0$ expresses that a fluid particle located at $x = \chi(X, t)$ at instant $t$ was originally at position $X$ at instant $t = 0$. This means that $X = \chi^{-1}(x, t)$. Now, let $(x, t) \mapsto q$ be a mapping that represents a medium parameter in Eulerian coordinates, we can define a Lagrangian representation $(X, t) \mapsto q^{\text{Lag}}$ of this parameter by setting $q^{\text{Lag}}(X, t) = q(\chi(X, t), t)$. If we note $\rho^{\text{Lag}}(X, t = 0) = \rho_0^{\text{Lag}}(X)$, the system (8) for one-dimensional problems is equivalent to

$$ \rho_0^{\text{Lag}} \partial_t (1/\rho^{\text{Lag}}) - \partial_X u^{\text{Lag}} = 0 \qquad (\text{B.1a}) $$

$$ \rho_0^{\text{Lag}} \partial_t u^{\text{Lag}} + \partial_X p^{\text{Lag}} = 0 \qquad (\text{B.1b}) $$

$$ \rho_0^{\text{Lag}} \partial_t E^{\text{Lag}} + \partial_X (p^{\text{Lag}} u^{\text{Lag}}) = 0 \qquad (\text{B.1c}) $$

$$ \rho_0^{\text{Lag}} \partial_t \mathcal{Y}_k^{\text{Lag}} = 0, \quad k = 1, \dots, m-1, \qquad (\text{B.1d}) $$

$$ \rho_0^{\text{Lag}} \partial_t \mathcal{Z}_k^{\text{Lag}} = 0, \quad k = 1, \dots, m-1. \qquad (\text{B.1e}) $$

## Appendix C. Preservation of iso-pressure and iso-velocity profiles

As in [41], we prove that the anti-diffusive scheme introduced in section 3 preserves a family of iso-pressure and iso-velocity profiles. Let us introduce the primitive variable $\mathbf{V}^n = (u^n, E^n, \mathcal{Y}_1^n, \dots, \mathcal{Y}_m^n)$ for one-dimensional problems.

**Proposition 6.** Let $(\rho^n \mathbf{V}^n, \tilde{\mathcal{Z}}_1^n, \dots, \tilde{\mathcal{Z}}_m^n)$ be the approximate state variable at instant $t = t^n$. Suppose that $(\rho^n \mathbf{V}^n, \tilde{\mathcal{Z}}_1^n, \dots, \tilde{\mathcal{Z}}_m^n)$ is a contact discontinuity with both uniform velocity and pressure such that

$$ p_i^n = \bar{p}, \quad u_i^n = \bar{u}, \quad (\rho_k)_i^n = \bar{\rho}_k, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

Then the approximate solution $(\rho^{n+1} \mathbf{V}^{n+1}, \tilde{\mathcal{Z}}_1^{n+1}, \dots, \tilde{\mathcal{Z}}_m^{n+1})$ computed with the anti-diffusive scheme satisfies

$$ p_i^{n+1} = \bar{p}, \quad u_i^{n+1} = \bar{u}, \quad (\rho_k)_i^{n+1} = \bar{\rho}_k, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

*Proof.* Let us first consider the following variable $(\tilde{\rho} \tilde{\mathbf{V}}^n, \tilde{\mathcal{Z}}_1^n, \dots, \tilde{\mathcal{Z}}_m^n)$ deduced from the Lagrangian step. Since the phasic densities and the pressure are constant at $t = t^n$, we have

$$ (\rho_k e_k)_i^n = \overline{\rho_k e_k}, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. \qquad (\text{C.1}) $$

In addition as the velocity is constant, the Lagrangian fluxes needed in the acoustic scheme are equal to zero i.e. $[\rho^n]_i = 0, [\rho^n u^n]_i = 0$ and $L_i = 1$. Thereby, we have $\forall i \in \mathbb{Z}$

$$ \left\{ \begin{array}{l} (\tilde{\rho}u)_i = (\rho u)_i^n, \\ (\tilde{\rho}E)_i = (\rho E)_i^n, \\ (\tilde{\rho}_k \tilde{\mathcal{Z}}_k)_i = (\rho_k \mathcal{Z}_k)_i^n, \quad 1 \le k \le m, \\ \tilde{\mathcal{Z}}_{k,i} = \mathcal{Z}_{k,i}^n, \quad 1 \le k \le m. \end{array} \right. $$
Thanks to the previous relations, we obtain $\tilde{\rho}_{k,i} = \rho_{k,i}^n = \bar{\rho}_k$, $1 \le k \le m$. Then by summation over $k = 1, \dots, m$ and using theorem 1 we have $\tilde{\rho}_i = \rho_i^n$. This leads to

$$ \tilde{u}_i = u_i^n = \bar{u} \text{ and } \tilde{E}_i = E_i^n, \quad \forall i \in \mathbb{Z}. $$

By definition of $E$, it yields $\tilde{e}_i = e_i^n$ and thus

$$ \tilde{\rho} \tilde{e}_i = (\rho e)_i^n, \quad \forall i \in \mathbb{Z}. $$

In the context of the isobaric closure we have for all $i$ in $\mathbb{Z}$

$$ \left\{ \begin{aligned} \tilde{p}_i &= p_1(\bar{p}_1, (\tilde{\rho}_1 \tilde{e}_1)_i) = \dots = p_m(\bar{\rho}_m, (\tilde{\rho}_m \tilde{e}_m)_i), \\ \tilde{\rho} e_i &= (\rho e)_i^n = \sum_{k=1}^m (\rho_k \mathcal{Z}_k)_i^n \tilde{\rho}_k \tilde{e}_k. \end{aligned} \right. \qquad (C.2) $$

However at $t = t^n$ we have $p_i^n = \bar{p}$ and $(\rho_k e_k)_i^n = \overline{\rho_k e_k}$, $1 \le k \le m$. Thus, we have

$$ \left\{ \begin{aligned} \bar{p} &= p_1(\bar{p}_1, (\tilde{\rho}_1 \tilde{e}_1)_i) = \dots = p_m(\bar{\rho}_m, (\tilde{\rho}_m \tilde{e}_m)_i), \\ (\rho e)_i^n &= \sum_{k=1}^m (\rho_k \mathcal{Z}_k)_i^n \overline{\rho_k e_k}. \end{aligned} \right. \qquad (C.3) $$

By proposition 1, we know that (C.2) and (C.3) admit a unique solution and therefore have the same solution, thus

$$ \tilde{p}_i = \bar{p}, \quad \tilde{\rho}_k \tilde{e}_k = \overline{\rho_k e_k}, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

Now let us consider the projection step, in the present case it boils down to

$$ \left\{ \begin{aligned} (\rho u)_i^{n+1} &= (\tilde{\rho} u)_i & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\rho} u \rrbracket_i, \\ (\rho E)_i^{n+1} &= (\tilde{\rho} E)_i & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\rho} E \rrbracket_i, \\ (\rho_k \mathcal{Z}_k)_i^{n+1} &= (\tilde{\rho}_k \tilde{\mathcal{Z}}_k)_i & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\rho}_k \tilde{\mathcal{Z}}_k \rrbracket_i, & 1 \le k \le m \\ \mathcal{Z}_{k,i}^{n+1} &= \tilde{\mathcal{Z}}_{k,i} & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\mathcal{Z}}_k \rrbracket_i, & 1 \le k \le m. \end{aligned} \right. \qquad (C.4) $$

Assuming without any loss of generality that $\bar{u} > 0$ and flux definitions (20), we have

$$ \tilde{\rho} \bar{u}_{i+1/2} = \bar{u} \sum_{k=0}^{m} \overline{\rho_k \tilde{\mathcal{Z}}_{k,i+1/2}} \text{ and } \tilde{\rho} E_{i+1/2} = \sum_{k=0}^{m} (\overline{\rho_k e_k} + \frac{1}{2} \overline{\rho_k} \bar{u}^2) \tilde{\mathcal{Z}}_{k,i+1/2}. $$

The update of both momentum and energy then reads

$$ (\rho u)_i^{n+1} - (\rho u)_i^n = \bar{u} \sum_{k=1}^{m} \overline{\rho_k} (\mathcal{Z}_{k,i}^{n+1} - \mathcal{Z}_{k,i}^n), \qquad (C.5) $$

and

$$ (\rho E)_i^{n+1} - (\rho E)_i^n = \bar{u} \sum_{k=1}^{m} (\overline{\rho_k e_k} + \frac{1}{2} \overline{\rho_k} \bar{u}^2) (\mathcal{Z}_{k,i}^{n+1} - \mathcal{Z}_{k,i}^n), \qquad (C.6) $$
Thanks to the previous relations, we obtain $\tilde{\rho}_{k,i} = \rho_{k,i}^n = \bar{\rho}_k$, $1 \le k \le m$. Then by summation over $k = 1, \dots, m$ and using theorem 1 we have $\tilde{\rho}_i = \rho_i^n$. This leads to

$$ \tilde{u}_i = u_i^n = \bar{u} \text{ and } \tilde{E}_i = E_i^n, \quad \forall i \in \mathbb{Z}. $$

By definition of $E$, it yields $\tilde{e}_i = e_i^n$ and thus

$$ \tilde{\rho} \tilde{e}_i = (\rho e)_i^n, \quad \forall i \in \mathbb{Z}. $$

In the context of the isobaric closure we have for all $i$ in $\mathbb{Z}$

$$ \left\{ \begin{aligned} \tilde{p}_i &= p_1(\bar{p}_1, (\tilde{\rho}_1 \tilde{e}_1)_i) = \dots = p_m(\bar{\rho}_m, (\tilde{\rho}_m \tilde{e}_m)_i), \\ \tilde{\rho} e_i &= (\rho e)_i^n = \sum_{k=1}^m (\rho_k \mathcal{Z}_k)_i^n \tilde{\rho}_k \tilde{e}_k. \end{aligned} \right. \quad (C.2) $$

However at $t = t^n$ we have $p_i^n = \bar{p}$ and $(\rho_k e_k)_i^n = \overline{\rho_k e_k}$, $1 \le k \le m$. Thus, we have

$$ \left\{ \begin{aligned} \bar{p} &= p_1(\bar{p}_1, (\tilde{\rho}_1 \tilde{e}_1)_i) = \dots = p_m(\bar{\rho}_m, (\tilde{\rho}_m \tilde{e}_m)_i), \\ (\rho e)_i^n &= \sum_{k=1}^m (\rho_k \mathcal{Z}_k)_i^n \overline{\rho_k e_k}. \end{aligned} \right. \quad (C.3) $$

By proposition 1, we know that (C.2) and (C.3) admit a unique solution and therefore have the same solution, thus

$$ \tilde{p}_i = \bar{p}, \quad \tilde{\rho}_k \tilde{e}_k = \overline{\rho_k e_k}, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

Now let us consider the projection step, in the present case it boils down to

$$ \left\{ \begin{aligned} (\rho u)_i^{n+1} &= (\tilde{\rho} u)_i & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\rho} u \rrbracket_i, \\ (\rho E)_i^{n+1} &= (\tilde{\rho} E)_i & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\rho} E \rrbracket_i, \\ (\rho_k \mathcal{Z}_k)_i^{n+1} &= (\tilde{\rho}_k \tilde{\mathcal{Z}}_k)_i & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\rho}_k \tilde{\mathcal{Z}}_k \rrbracket_i, & 1 \le k \le m \\ \mathcal{Z}_{k,i}^{n+1} &= \tilde{\mathcal{Z}}_{k,i} & - \frac{\Delta t}{\Delta x} \overline{u} \llbracket \tilde{\mathcal{Z}}_k \rrbracket_i, & 1 \le k \le m. \end{aligned} \right. \quad (C.4) $$

Assuming without any loss of generality that $\bar{u} > 0$ and flux definitions (20), we have

$$ \tilde{\rho} \bar{u}_{i+1/2} = \bar{u} \sum_{k=0}^{m} \overline{\rho_k \tilde{\mathcal{Z}}_{k,i+1/2}} \text{ and } \tilde{\rho} E_{i+1/2} = \sum_{k=0}^{m} (\overline{\rho_k e_k} + \frac{1}{2} \overline{\rho_k} \bar{u}^2) \tilde{\mathcal{Z}}_{k,i+1/2}. $$

The update of both momentum and energy then reads

$$ (\rho u)_i^{n+1} - (\rho u)_i^n = \bar{u} \sum_{k=1}^{m} \overline{\rho_k} (\mathcal{Z}_{k,i}^{n+1} - \mathcal{Z}_{k,i}^n), \quad (C.5) $$

and

$$ (\rho E)_i^{n+1} - (\rho E)_i^n = \bar{u} \sum_{k=1}^{m} (\overline{\rho_k e_k} + \frac{1}{2} \overline{\rho_k} \bar{u}^2) (\mathcal{Z}_{k,i}^{n+1} - \mathcal{Z}_{k,i}^n), \quad (C.6) $$
By combining both last relations of (C.4), we obtain

$$ \rho_{k,i}^{n+1} = \bar{\rho}_k, \quad 1 \le k \le m. $$

Using (C.5), and theorem 1 we deduce that

$$ \rho_i^n (u_i^{n+1} - \bar{u}) = 0 \Rightarrow u_i^{n+1} = \bar{u}, \quad \forall i \in \mathbb{Z} $$

similarly (C.6) leads to

$$ \sum_{k=1}^{m} ((\rho_k e_k)_i^{n+1} - \overline{\rho_k e_k}) Z_{k,i}^{n+1} = 0. $$

Using the same uniqueness arguments as for in the Lagrange step for the definition of the pressure by isobaric closure we finally obtain

$$ p_i^{n+1} = \bar{p}, \quad (\rho_k e_k)_i^{n+1} = \overline{\rho_k e_k}, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

□
By combining both last relations of (C.4), we obtain

$$ \rho_{k,i}^{n+1} = \bar{\rho}_k, \quad 1 \le k \le m. $$

Using (C.5), and theorem 1 we deduce that

$$ \rho_i^n (u_i^{n+1} - \bar{u}) = 0 \Rightarrow u_i^{n+1} = \bar{u}, \quad \forall i \in \mathbb{Z} $$

similarly (C.6) leads to

$$ \sum_{k=1}^{m} ((\rho_k e_k)_i^{n+1} - \overline{\rho_k e_k}) Z_{k,i}^{n+1} = 0. $$

Using the same uniqueness arguments as for in the Lagrange step for the definition of the pressure by isobaric closure we finally obtain

$$ p_i^{n+1} = \bar{p}, \quad (\rho_k e_k)_i^{n+1} = \overline{\rho_k e_k}, \quad \forall i \in \mathbb{Z}, \quad 1 \le k \le m. $$

□
References

[1] R. Abgrall, *How to prevent pressure oscillations in multicomponent flow calculations: a quasi-conservative approach*, J. Comput. Phys., Vol. 125, Issue 1, pp. 150-160, 1996.

[2] R. Abgrall, S. Karni, *Computations of compressible multifluids*, J. Comput. Phys., Vol. 169, Issue 2, pp. 594-623, 2001.

[3] G. Allaire, S. Clerc, S. Kokh, *A five-equation model for the numerical simulation of interfaces in two-phase flows*, C.R. Acad. Sci. Paris - Series I - Mathematics, Vol. 331, Issue 12, Part 1, pp. 1017-1022, 2000.

[4] G. Allaire, S. Clerc, S. Kokh, *A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids*, J. Comput. Phys., Vol. 181, Issue 2, pp. 577-616, 2002.

[5] M.R. Baer, J.W. Nunziato, *A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials*, Int. J. Mutliphase Flow, Vol. 12, Issue 6, pp. 861-889, 1986.

[6] P. Batten, N. Clarke, C. Lambert, D.M. Causon, *On the choice of wavespeeds for the HLLC Riemann solver*, SIAM J. Sci. Comput. 18, pp. 1553-1570, 1997.

[7] B. Braconnier, B. Nkonga, *An all-speed relaxation scheme for interface flows with surface tension*, J. Comput. Phys., Vol. 228, Issue 16, pp. 5722-5739, 2009.

[8] M. Billaud Friess, B. Boutin, F. Caetano, G. Faccanoni, S. Kokh, F. Lagoutière and L. Navoret, *A second order antidiffusive Lagrange-remap scheme for two-component flows*, ESAIM: PROCEEDINGS, Vol. 32, pp. 149-162, 2011.

[9] M. Billaud Friess, S. Kokh, *An anti-diffusive Lagrange-Remap scheme for multi-material compressible flows with an arbitrary number of component*, ESAIM: PROCEEDINGS, Vol. 35, pp. 203-209, 2012.

[10] F. Bouchut, *Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources*, Frontiers in Mathematics series, Birkhäuser, 2004.

[11] F. Boyer, C. Lapuerta, *Study of a three component Cahn-Hilliard flow model*, ESAIM: Math. Model. Num., Vol. 40, Issue 4, pp. 653-687, 2006.

[12] F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar and M. Quintard, *Cahn-Hilliard / Navier-Stokes model for the simulation of three-phase flows*, Transport Porous Med., Vol. 82, Issue 3, pp. 463-483, 2010.
References

[1] R. Abgrall, *How to prevent pressure oscillations in multicomponent flow calculations: a quasi-conservative approach*, J. Comput. Phys., Vol. 125, Issue 1, pp. 150-160, 1996.

[2] R. Abgrall, S. Karni, *Computations of compressible multifluids*, J. Comput. Phys., Vol. 169, Issue 2, pp. 594-623, 2001.

[3] G. Allaire, S. Clerc, S. Kokh, *A five-equation model for the numerical simulation of interfaces in two-phase flows*, C.R. Acad. Sci. Paris - Series I - Mathematics, Vol. 331, Issue 12, Part 1, pp. 1017-1022, 2000.

[4] G. Allaire, S. Clerc, S. Kokh, *A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids*, J. Comput. Phys., Vol. 181, Issue 2, pp. 577-616, 2002.

[5] M.R. Baer, J.W. Nunziato, *A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials*, Int. J. Mutliphase Flow, Vol. 12, Issue 6, pp. 861-889, 1986.

[6] P. Batten, N. Clarke, C. Lambert, D.M. Causon, *On the choice of wavespeeds for the HLLC Riemann solver*, SIAM J. Sci. Comput. 18, pp. 1553-1570, 1997.

[7] B. Braconnier, B. Nkonga, *An all-speed relaxation scheme for interface flows with surface tension*, J. Comput. Phys., Vol. 228, Issue 16, pp. 5722-5739, 2009.

[8] M. Billaud Friess, B. Boutin, F. Caetano, G. Faccanoni, S. Kokh, F. Lagoutière and L. Navoret, *A second order antidiffusive Lagrange-remap scheme for two-component flows*, ESAIM: PROCEEDINGS, Vol. 32, pp. 149-162, 2011.

[9] M. Billaud Friess, S. Kokh, *An anti-diffusive Lagrange-Remap scheme for multi-material compressible flows with an arbitrary number of component*, ESAIM: PROCEEDINGS, Vol. 35, pp. 203-209, 2012.

[10] F. Bouchut, *Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources*, Frontiers in Mathematics series, Birkhäuser, 2004.

[11] F. Boyer, C. Lapuerta, *Study of a three component Cahn-Hilliard flow model*, ESAIM: Math. Model. Num., Vol. 40, Issue 4, pp. 653-687, 2006.

[12] F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar and M. Quintard, *Cahn-Hilliard / Navier-Stokes model for the simulation of three-phase flows*, Transport Porous Med., Vol. 82, Issue 3, pp. 463-483, 2010.
[13] F. Boyer, S. Minjeaud, *Numerical schemes for a three component Cahn-Hilliard model*, ESAIM: Math. Model. Num., Vol. 45, Issue 4, pp. 697-738, 2011.

[14] F. Coquel, E. Godlewski, B. Perthame, A. In, P. Rascle, *Some new Godunov and relaxation methods for two-phase flow problems*, Godunov methods (Oxford, 1999), Kluwer/Plenum, New York, pp. 179-188, 2001.

[15] I.-L. Chern, J. Glimm, O. McBryan, B. Plohr, and S. Yaniv, *Front tracking for gas dynamics*, J. Comput. Phys., Vol. 62, Issue 1, pp. 83-110, 1986.

[16] S. Dellacherie, *Relaxation schemes for the multicomponent Euler system*, ESAIM: Math. Model. Num., Vol. 37, Issue 6, pp. 909-936, 2003.

[17] S. Del Pino, *Metric-based mesh adaptation for 2D Lagrangian compressible flows*, J. Comput. Phys., Vol. 230, Issue 5, pp. 1793-1821, 2011.

[18] B. Després, F. Lagoutière, *Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire*, C.R. Acad. Sci. Paris - Series I - Mathematics, Vol. 328, Issue 10, pp. 939-943, 1999.

[19] B. Després, F. Lagoutière, *Contact discontinuity capturing schemes for linear advection and compressible gas dynamics*, J. Sci Comp., Vol. 16, Issue 4, pp. 479-524, 2001.

[20] B. Després, F. Lagoutière, *Numerical resolution of a two-component compressible fluid model with interfaces*, Progr. Comput. Fluid. Dynam., Vol. 7, Issue 6, pp. 295-310, 2007.

[21] B. Després, *Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques*, Mathématiques et applications, SMAI, Springer, Vol. 68, 2010.

[22] B. Després, F. Lagoutière, E. Labourasse, I. Marmajou, *An anti-dissipative transport scheme on unstructured meshes for multicomponent flows*, Int. J. Finite Vol., Vol.7, pp. 30-65, 2010.

[23] V. Dyadechko, M. Shashkov, *Reconstruction of Multi-material Interfaces from Moment Data*, J. Comput. Phys., Vol. 227, Issue 11, pp. 5161-5184, 2008.

[24] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, *A hybrid particle level set method for improved interface capturing*, J. Comput. Phys., Vol. 183, Issue 1, pp. 83-116, 2002.

[25] V. Faucher, S. Kokh, *Extended Vofire algorithm for fast transient fluid-structure dynamics with liquid-gas flows and interfaces*, J. of Fluids and Structures, Vol. 39, pp. 102-125, 2013.
[13] F. Boyer, S. Minjeaud, *Numerical schemes for a three component Cahn-Hilliard model*, ESAIM: Math. Model. Num., Vol. 45, Issue 4, pp. 697-738, 2011.

[14] F. Coquel, E. Godlewski, B. Perthame, A. In, P. Rascle, *Some new Godunov and relaxation methods for two-phase flow problems*, Godunov methods (Oxford, 1999), Kluwer/Plenum, New York, pp. 179-188, 2001.

[15] I.-L. Chern, J. Glimm, O. McBryan, B. Plohr, and S. Yaniv, *Front tracking for gas dynamics*, J. Comput. Phys., Vol. 62, Issue 1, pp. 83-110, 1986.

[16] S. Dellacherie, *Relaxation schemes for the multicomponent Euler system*, ESAIM: Math. Model. Num., Vol. 37, Issue 6, pp. 909-936, 2003.

[17] S. Del Pino, *Metric-based mesh adaptation for 2D Lagrangian compressible flows*, J. Comput. Phys., Vol. 230, Issue 5, pp. 1793-1821, 2011.

[18] B. Després, F. Lagoutière, *Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire*, C.R. Acad. Sci. Paris - Series I - Mathematics, Vol. 328, Issue 10, pp. 939-943, 1999.

[19] B. Després, F. Lagoutière, *Contact discontinuity capturing schemes for linear advection and compressible gas dynamics*, J. Sci Comp., Vol. 16, Issue 4, pp. 479-524, 2001.

[20] B. Després, F. Lagoutière, *Numerical resolution of a two-component compressible fluid model with interfaces*, Progr. Comput. Fluid. Dynam., Vol. 7, Issue 6, pp. 295-310, 2007.

[21] B. Després, *Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques*, Mathématiques et applications, SMAI, Springer, Vol. 68, 2010.

[22] B. Després, F. Lagoutière, E. Labourasse, I. Marmajou, *An anti-dissipative transport scheme on unstructured meshes for multicomponent flows*, Int. J. Finite Vol., Vol.7, pp. 30-65, 2010.

[23] V. Dyadechko, M. Shashkov, *Reconstruction of Multi-material Interfaces from Moment Data*, J. Comput. Phys., Vol. 227, Issue 11, pp. 5161-5184, 2008.

[24] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, *A hybrid particle level set method for improved interface capturing*, J. Comput. Phys., Vol. 183, Issue 1, pp. 83-116, 2002.

[25] V. Faucher, S. Kokh, *Extended Vofire algorithm for fast transient fluid-structure dynamics with liquid-gas flows and interfaces*, J. of Fluids and Structures, Vol. 39, pp. 102-125, 2013.
[26] R. Fedkiw, T. Aslam, B. Merriman, S. Osher, *A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the Ghost Fluid Method)*, J. Comput. Phys., Vol. 152, Issue 2, pp. 457–492, 1999.

[27] S. Galéra, P.-H. Maire, J. Breil, *A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction*, J. Comput. Phys., Vol. 229, Issue 16, pp. 5755–5787, 2010.

[28] S. Galéra, J. Breil, P.-H. Maire, *A 2D unstructured multi-material Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) scheme using MOF interface reconstruction*, Comput. Fluids, Vol. 46, Issue 1, pp. 237–244, 2011.

[29] F. Daude, P. Galon, Z. Gao, E. Blaud, *Numerical experiments using a HLLC-type scheme with ALE formulation for compressible two-phase flows five-equation models with phase transition*, Computers & Fluids, Vol. 94, pp. 112–138.

[30] J. Glimm, O. McBryan, *A computational model for interfaces*, Adv. Appl. Math., Vol. 6, Issue 4, 1985.

[31] E. Godlewski, P.-A. Raviart, *Numerical approximation of hyperbolic systems of conservation laws*, Applied Mathematic Science, Springer, 1996.

[32] J.F. Haas, B. Sturtevant, *Interaction of a weak shock wave with cylindrical and spherical gas inhomogeneities*, J. Fluid Mech., Vol. 181, pp. 41–76, 1987.

[33] A. Harten, *On a class of high resolution total-variation-stable finite difference schemes*, SIAM J. Numer. Anal., Vol. 21, pp. 1–23, 1984.

[34] Code HERACLES: http://irfu.cea.fr/Projets/Site_heracles

[35] C.W. Hirt, B.D. Nichols, *Vol. of Fluid (VOF) method for the dynamics of free boundaries*, J. Comput. Phys., Vol. 39, Issue 1, pp. 201–205, 1981.

[36] S. Jaouen, F. Lagoutière, *Numerical transport of an arbitrary number of components*, Comput. Method Appl. M., Vol. 196, Issues 33-34, pp. 3127–3140, 2007.

[37] D. Juric, G. Tryggvason, *A front tracking method for dendritic solidification*, J. Comput. Phys., Vol. 123, Issue 1, pp. 127-148, 1996.

[38] S. Karni, *Multicomponent flow calculations by a consistent primitive algorithm*, J. Comput. Phys., Vol. 112, Issue 1, pp. 31–43, 1994.

[39] S. Karni, *Hybrid multifluid algorithms*, SIAM J. Sci. Comput., Vol. 17, Issue 5, pp. 1019–1039, 1996.

[40] M. Kucharik, R.V. Garimella, S.P. Schofield, M.J. Shashkov, *A comparative study of interface reconstruction methods for multi-material ALE simulations*, J. Comput. Phys., Vol. 229, Issue 7, pp. 2432–2452, 2010.
[26] R. Fedkiw, T. Aslam, B. Merriman, S. Osher, *A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the Ghost Fluid Method)*, J. Comput. Phys., Vol. 152, Issue 2, pp. 457–492, 1999.

[27] S. Galéra, P.-H. Maire, J. Breil, *A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction*, J. Comput. Phys., Vol. 229, Issue 16, pp. 5755–5787, 2010.

[28] S. Galéra, J. Breil, P.-H. Maire, *A 2D unstructured multi-material Cell-Centered Arbitrary Lagrangian-Eulerian (CCALE) scheme using MOF interface reconstruction*, Comput. Fluids, Vol. 46, Issue 1, pp. 237–244, 2011.

[29] F. Daude, P. Galon, Z. Gao, E. Blaud, *Numerical experiments using a HLLC-type scheme with ALE formulation for compressible two-phase flows five-equation models with phase transition*, Computers & Fluids, Vol. 94, pp. 112–138.

[30] J. Glimm, O. McBryan, *A computational model for interfaces*, Adv. Appl. Math., Vol. 6, Issue 4, 1985.

[31] E. Godlewski, P.-A. Raviart, *Numerical approximation of hyperbolic systems of conservation laws*, Applied Mathematic Science, Springer, 1996.

[32] J.F. Haas, B. Sturtevant, *Interaction of a weak shock wave with cylindrical and spherical gas inhomogeneities*, J. Fluid Mech., Vol. 181, pp. 41–76, 1987.

[33] A. Harten, *On a class of high resolution total-variation-stable finite difference schemes*, SIAM J. Numer. Anal., Vol. 21, pp. 1–23, 1984.

[34] Code HERACLES: http://irfu.cea.fr/Projets/Site_heracles

[35] C.W. Hirt, B.D. Nichols, *Vol. of Fluid (VOF) method for the dynamics of free boundaries*, J. Comput. Phys., Vol. 39, Issue 1, pp. 201–205, 1981.

[36] S. Jaouen, F. Lagoutière, *Numerical transport of an arbitrary number of components*, Comput. Method Appl. M., Vol. 196, Issues 33-34, pp. 3127–3140, 2007.

[37] D. Juric, G. Tryggvason, *A front tracking method for dendritic solidification*, J. Comput. Phys., Vol. 123, Issue 1, pp. 127-148, 1996.

[38] S. Karni, *Multicomponent flow calculations by a consistent primitive algorithm*, J. Comput. Phys., Vol. 112, Issue 1, pp. 31–43, 1994.

[39] S. Karni, *Hybrid multifluid algorithms*, SIAM J. Sci. Comput., Vol. 17, Issue 5, pp. 1019–1039, 1996.

[40] M. Kucharik, R.V. Garimella, S.P. Schofield, M.J. Shashkov, *A comparative study of interface reconstruction methods for multi-material ALE simulations*, J. Comput. Phys., Vol. 229, Issue 7, pp. 2432–2452, 2010.
[41] S. Kokh, F. Lagoutière, *An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model*, J. Comput. Phys., Vol. 229, Issue 8, pp. 2773-2809, 2010.

[42] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, *Modelling merging and fragmentation in multiphase flows with SURFER*, J. Comput. Phys., Vol. 113, Issue 1, pp. 134-147, 1994.

[43] F. Lagoutière, *Modélisation mathématique et résolution numérique de problèmes de fluides compressibles à plusieurs constituants*, Ph.D Thesis, Université Pierre et Marie Curie, 2000.

[44] A.-Y. Le Roux, *A numerical conception of entropy for quasi-linear equations*, Math. of Comp. 31, number 140, pp. 848-872, 1977.

[45] T.G. Liu, B.C. Khoo, K.S. Yeo, *Ghost fluid method for strong shock impact-ing on material interface*, J. Comput. Phys., Vol. 190, Issue 2, pp. 651-681, 2003.

[46] J. Massoni, R. Saurel, B. Nkonga, R. Abgrall, *Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur: Some models and Eulerian methods for interface problems between compressible fluids with heat transfer*, Int. J. Heat. Mass. Tran., Vol. 45, Issue 6, pp. 1287-1307, 2002.

[47] A. Murrone, H. Guillard, *A five equation reduced model for compressible two phase flow problems*, J. Comput. Phys., Vol. 202, Issue 2, pp. 664-698, 2005.

[48] W. Mulder, S.Osher, J.A. Sethian, *Computing interface motion in compressible gas dynamics*, J. Comput. Phys., Vol. 100, Issue 2, pp. 209-228, 1992.

[49] D. Nguyen, F. Gibou, R. Fedkiw, *A fully conservative ghost fluid method & stiff detonation waves*, In 12th Int. Detonation Symposium, San Diego, CA, 2002.

[50] S. Osher, J.A. Sethian, *Fronts propagating with curvature dependant speed: algorithms based on Hamilton-Jacobi formulations*, J. Comput. Phys., Vol. 79, Issue 1, pp. 12-49, 1988.

[51] P.L. Roe, *Approximate Riemann solvers, parameter vectors and difference schemes*, J. Comp. Phys., vol. 43, pp. 357-372, 1981.

[52] V.V. Rusanov. *Calculation of interaction of non-steady shock waves with obstacles*, J. Comp. Math. Phys. USSR, vol. 1, pp. 267-279, 1961.

[53] R. Saurel, R. Abgrall, *A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows*, J. Comput. Phys., Vol. 150, Issue 2, pp. 425-467, 1999.
[41] S. Kokh, F. Lagoutière, *An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model*, J. Comput. Phys., Vol. 229, Issue 8, pp. 2773-2809, 2010.

[42] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, *Modelling merging and fragmentation in multiphase flows with SURFER*, J. Comput. Phys., Vol. 113, Issue 1, pp. 134-147, 1994.

[43] F. Lagoutière, *Modélisation mathématique et résolution numérique de problèmes de fluides compressibles à plusieurs constituants*, Ph.D Thesis, Université Pierre et Marie Curie, 2000.

[44] A.-Y. Le Roux, *A numerical conception of entropy for quasi-linear equations*, Math. of Comp. 31, number 140, pp. 848-872, 1977.

[45] T.G. Liu, B.C. Khoo, K.S. Yeo, *Ghost fluid method for strong shock impact-ing on material interface*, J. Comput. Phys., Vol. 190, Issue 2, pp. 651-681, 2003.

[46] J. Massoni, R. Saurel, B. Nkonga, R. Abgrall, *Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur: Some models and Eulerian methods for interface problems between compressible fluids with heat transfer*, Int. J. Heat. Mass. Tran., Vol. 45, Issue 6, pp. 1287-1307, 2002.

[47] A. Murrone, H. Guillard, *A five equation reduced model for compressible two phase flow problems*, J. Comput. Phys., Vol. 202, Issue 2, pp. 664-698, 2005.

[48] W. Mulder, S.Osher, J.A. Sethian, *Computing interface motion in compressible gas dynamics*, J. Comput. Phys., Vol. 100, Issue 2, pp. 209-228, 1992.

[49] D. Nguyen, F. Gibou, R. Fedkiw, *A fully conservative ghost fluid method & stiff detonation waves*, In 12th Int. Detonation Symposium, San Diego, CA, 2002.

[50] S. Osher, J.A. Sethian, *Fronts propagating with curvature dependant speed: algorithms based on Hamilton-Jacobi formulations*, J. Comput. Phys., Vol. 79, Issue 1, pp. 12-49, 1988.

[51] P.L. Roe, *Approximate Riemann solvers, parameter vectors and difference schemes*, J. Comp. Phys., vol. 43, pp. 357-372, 1981.

[52] V.V. Rusanov. *Calculation of interaction of non-steady shock waves with obstacles*, J. Comp. Math. Phys. USSR, vol. 1, pp. 267-279, 1961.

[53] R. Saurel, R. Abgrall, *A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows*, J. Comput. Phys., Vol. 150, Issue 2, pp. 425-467, 1999.
[54] R. Saurel, R. Abgrall, *A simple method for compressible multifluid flows*, SIAM J. Sci. Comput., Vol. 21, Issue 3, pp. 1115-1145, 1999.

[55] I. Suliciu, *On modelling phase transitions by means of rate-type constitutive equations, shock wave structure*, Internat. J. Engrg. Sci. 28 pp. 829-841, 1990.

[56] R. Scardovelli, S. Zaleski, *Direct numerical simulation of free-surface and interfacial flow*, Annual Review of Fluid Mechanics, Vol. 31, pp. 567-603, 1999.

[57] J.A. Sethian, *Level set methods: evolving interfaces in geometry, fluid mechanics, computer vision and materials sciences*, Cambridge Monographs on Applied and Computational Mathematics, 1996.

[58] C.D. Sijoy, S. Chaturvedi, *Volume-of-fluid algorithm with different modified dynamic material ordering methods and their comparisons*, J. Comput. Phys., Vol. 229, Issue 10, pp. 3848-3863, 2010.

[59] K.M. Shyue, *An efficient shock-capturing algorithm for compressible multicomponent problems*, J. Comput. Phys., Vol. 142, Issue 1, pp. 208-242, 1998.

[60] K.M. Shyue, *A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with van der Waals Equation of State*, J. Comput. Phys., Vol. 156, Issue 1, pp. 43-88, 1999.

[61] K.M. Shyue, *A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with Mie-Grüneisen Equation of State*, J. Comput. Phys., Vol. 171, Issue 2, pp. 678-707, 2001.

[62] V. Springel: *E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh*, Mon. Not. R. Astron. Soc., Vol. 40, Issue 2, pp.791-851, 2010.

[63] J.J. Quirk, S. Karni, *On the dynamics of a shock-bubble interaction*, J. Fluid Mech., Vol. 328, pp. 129-163, 1996.

[64] H. Terashima, G. Tryggvason, *A front-tracking/ghost-fluid method for fluid interfaces in compressible flows*, J. Comput. Phys., Vol. 2, Issue 213, pp. 4012-4037, 2009.

[65] E.F. Toro, *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*, second edition, Springer-Verlag, 1999.

[66] S.O. Unverdi, G. Tryggvason, *A front tracking method for viscous, incompressible, multi-fluid flows*, J. Comput. Phys., Vol. 100, Issue 1, p. 25-37, 1992.
[54] R. Saurel, R. Abgrall, *A simple method for compressible multifluid flows*, SIAM J. Sci. Comput., Vol. 21, Issue 3, pp. 1115-1145, 1999.

[55] I. Suliciu, *On modelling phase transitions by means of rate-type constitutive equations, shock wave structure*, Internat. J. Engrg. Sci. 28 pp. 829-841, 1990.

[56] R. Scardovelli, S. Zaleski, *Direct numerical simulation of free-surface and interfacial flow*, Annual Review of Fluid Mechanics, Vol. 31, pp. 567-603, 1999.

[57] J.A. Sethian, *Level set methods: evolving interfaces in geometry, fluid mechanics, computer vision and materials sciences*, Cambridge Monographs on Applied and Computational Mathematics, 1996.

[58] C.D. Sijoy, S. Chaturvedi, *Volume-of-fluid algorithm with different modified dynamic material ordering methods and their comparisons*, J. Comput. Phys., Vol. 229, Issue 10, pp. 3848-3863, 2010.

[59] K.M. Shyue, *An efficient shock-capturing algorithm for compressible multicomponent problems*, J. Comput. Phys., Vol. 142, Issue 1, pp. 208-242, 1998.

[60] K.M. Shyue, *A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with van der Waals Equation of State*, J. Comput. Phys., Vol. 156, Issue 1, pp. 43-88, 1999.

[61] K.M. Shyue, *A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with Mie-Grüneisen Equation of State*, J. Comput. Phys., Vol. 171, Issue 2, pp. 678-707, 2001.

[62] V. Springel: *E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh*, Mon. Not. R. Astron. Soc., Vol. 40, Issue 2, pp.791-851, 2010.

[63] J.J. Quirk, S. Karni, *On the dynamics of a shock-bubble interaction*, J. Fluid Mech., Vol. 328, pp. 129-163, 1996.

[64] H. Terashima, G. Tryggvason, *A front-tracking/ghost-fluid method for fluid interfaces in compressible flows*, J. Comput. Phys., Vol. 2, Issue 213, pp. 4012-4037, 2009.

[65] E.F. Toro, *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*, second edition, Springer-Verlag, 1999.

[66] S.O. Unverdi, G. Tryggvason, *A front tracking method for viscous, incompressible, multi-fluid flows*, J. Comput. Phys., Vol. 100, Issue 1, p. 25-37, 1992.
[67] C.W. Wang, T.G. Liu, B.C. Khoo, *A real ghost fluid method for the simulation of multimedium compressible flow*, SIAM J. Sci. Comput. ,Vol. 28, Issue 1, pp. 278–302, 2006.

[68] Z.D. Wu, L. Sun, Z. Zong, *A mass-fraction-based interface-capturing method for multi-component flow*, Int. J. Numer. Meth. Fluids, 2013.

[69] K. Yokoi, *A variational approach to motion of triple junction of gas, liquid and solid*, Comput. Phys. Commun., Vol. 180, Issue 7, pp. 1145–1149, 2009.

[70] H.K. Zhao, T. Chan, B. Merriman, S. Osher, *A variational level set approach to multiphase motion*, J. Comput. Phys., Vol. 127, Issue 1, pp. 179–195, 1996.
[67] C.W. Wang, T.G. Liu, B.C. Khoo, *A real ghost fluid method for the simulation of multimedium compressible flow*, SIAM J. Sci. Comput. ,Vol. 28, Issue 1, pp. 278–302, 2006.

[68] Z.D. Wu, L. Sun, Z. Zong, *A mass-fraction-based interface-capturing method for multi-component flow*, Int. J. Numer. Meth. Fluids, 2013.

[69] K. Yokoi, *A variational approach to motion of triple junction of gas, liquid and solid*, Comput. Phys. Commun., Vol. 180, Issue 7, pp. 1145–1149, 2009.

[70] H.K. Zhao, T. Chan, B. Merriman, S. Osher, *A variational level set approach to multiphase motion*, J. Comput. Phys., Vol. 127, Issue 1, pp. 179–195, 1996.