Sharp Interface Capturing in Compressible Multi-Material Flows with a Diffuse Interface Method

: Compressible multi-materialﬂows are encountered in a wide range of natural phenomena and industrial applications, such as supernova explosions in space, high speed ﬂows in jet and rocket propulsion, underwater explosions, and vapor explosions in post accidental situations in nuclear reactors. In the numerical simulations of these ﬂows, interfaces play a crucial role. A poor numerical resolution of the interfaces could make it difﬁcult to account for the physics, such as material separation, location of the shocks and contact discontinuities, and transfer of the mass, momentum and heat between different materials/phases. Owing to such importance, sharp interface capturing remains an active area of research in the ﬁeld of computational physics. To address this problem in this paper we focus on the Interface Capturing (IC) strategy, and thus we make use of a newly developed Diffuse Interface Method (DIM) called Multidimensional Limiting Process-Upper Bound (MLP-UB) . Our analysis shows that this method is easy to implement, can deal with any number of material interfaces, and produces sharp, shape-preserving interfaces, along with their accurate interaction with the shocks. Numerical experiments show good results even with the use of coarse meshes.


Introduction
Multi-phase compressible flows are present in a wide and ever expanding range of natural and industrial applications. Among them we can mention supernova explosions [1], high speed flows in jet and rocket propulsion [2], underwater explosions [3,4], and the scenario of vapor explosions in post accidental situations in nuclear reactors [5,6]. To compute these flows, two main families of numerical methods are used: Lagrangian and Eulerian. In the Lagrangian approach, the objective is to describe the phenomena that govern the state of an individual fluid element, while in the Eulerian approach, the objective is to describe the phenomena that govern the state of the fluid in a fixed volume element. In these modeling approaches, the main concerns are the treatment of the interfaces and the phenomena associated with them.
In the case of a Lagrangian approach, the interfaces are followed by moving and deforming the mesh [7,8]. For the Eulerian approach, it is possible to capture or reconstruct the interface by adding a dedicated equation. This is the approach adopted in the methods such as Volume Of Fluid (VOF) [9][10][11] and Level Set [12,13]. Front tracking methods [14,15] also exist, coupling some advantages of the Lagrangian and Eulerian representations. These methods consider the interface as a local discontinuity. However, algorithms dedicated for the determination of the position of interfaces are often expensive. Moreover, strong topology changes and formation/disappearance of the interfaces are particularly complex to deal with, especially with the increase in the number of phases present.
There is another class of Eulerian methods called algebraic methods, which belong to the class of sharp Interface Capturing (IC) approach [16] . These methods are often referred to as Diffuse Interface Methods (DIM) [17], an exhaustive up-to-date review of which can be found in [18]. In this family of methods, the diffusion of the interface is allowed. The interfaces are not subjected to any particular geometrical treatment, which makes their implementation easier. Moreover, these methods can handle the dynamic formation, disappearance and large deformation of the interfaces in flows with several phases. The propagation of the acoustic waves can also be naturally taken into account. However, one of the drawbacks of these methods is the smearing of the already diffused interfaces, due to the numerical diffusion. Significant smearing of the interfaces can then result in an extended zone of a pure numerical (non-physical) mixture of the different materials in the vicinity of the interfaces. Thus, the accurate description of any physical phenomena at the interface may be deteriorated by this wide diffusion zone.
For our own end-use applications, we have decided to focus in this paper on the IC approach to solve the multi-material compressible problem. To ensure the sharp interface capturing in the solution of a compressible multi-phase model, a second-order multidimensional MUSCL [19] reconstruction is used with a recently developed compressive, shapepreserving gradient limiter, referred to as Multidimensional Limiting Process-Upper Bound (MLP-UB) [18], which is an updated variant of the classical MLP approach demonstrated in [20]. Sharp interface capturing methods are generally Finite Volume (FV) methods, either based on the high order flux evaluations or on the slope reconstruction, but in all the cases they make use of the limiters (flux or slope). Besides the MLP-UB approach adopted here, other popular methods/schemes found in the literature are: limited downwind scheme [21], THINC scheme [22] or the method using correction algorithms [23].
The one-pressure, compressible system of the continuous equations considered in this work is focused on the phase volume fractions together with the evolution of the mean mass, momentum and energy of each phase [24][25][26]. This system is a well-known problem allowing the assessment of the good behavior of our IC approach together with the use of a compressive, shape-preserving gradient limiter for an arbitrary number of phases. This study is the sequel to the work done in [18], which was restricted to the cases modeled with a pure advection equation. Thus, the starting point model proposed in this paper is intended to evolve using the appropriate closure terms in order to generate the models of industrial relevance. These added terms could be obtained by applying a statistical averaging operator on the corresponding elementary conservation laws written for each fluid. This produces many fluctuation correlation terms that must be closed, usually by introducing the supplementary physical properties of the system. This paper is organized as follows: after this introduction, we present in Section 2, modelling and governing equations for the compressible multi-phase flows considered in this work. Next, Sections 3-5 introduce the time discretization, and interface sharpening strategy within the framework of space discretization. Following this we present the numerical results obtained for the two-material flows in Section 6, and for the multimaterial flows in Section 7. Finally, the conclusions are drawn and scope for the future work is briefly discussed.

Physical Model and Governing Equations
Considering a compressible multi-phase flow with N phases, where each phase is indexed by φ ∈ {1, ..., N}, and is modelled as a compressible medium with a thermodynamic state defined by its: 1. mass density ρ φ ; 2. specific internal energy e φ ; 3.
Further, each phase φ obeys a coupled system of governing equations, composed of the balance equations of mass, momentum and energy for each phase. We consider here the situation where the pressures are at the instantaneous equilibrium everywhere, so that we have: where p is an equilibrium pressure. We also consider that every control volume (in the terminology of the Eulerian framework) is filled by different phases, such that: with α φ ∈ [0, 1] is the volume fraction of any phase φ. Additionally, the system is completed with two equations of state (EOS), relating the four thermodynamic variables. The wellknown complete set of volume-averaged equations for any phase φ is written here in terms of the internal energies, as shown below: The above system is identified by the unknowns u φ , ρ φ , α φ , e φ , T φ , p φ . u φ is the phase velocity and we denote the components of the vector u φ = (u φ , v φ ) in 2D and u φ = (u φ , v φ , w φ ) in 3D. EOS i stands for the Equation Of State number i; i = 1, 2. The form of EOS i mentioned above is generic. The specific definitions are given in Sections 6 and 7.
Furthermore, in the present study, as evident from the above system of equations, we have separate momentum equations for each phase with independent velocities per phase. It is known that without the velocity relaxation terms, this system is conditionally hyperbolic with the loss of hyperbolicity under certain conditions [27,28]. This issue can be fixed by using a large friction source term F φ for balancing velocities between the phases.
In this work we deal with the situation where there is a full mechanical equilibrium. This implies that the phases share a common velocity, which is achieved through the use of a term F φ . The form of which is: where D is a friction constant which determines the rate at which velocities equilibrate, and φ * is a conjugate phase i.e., φ * = φ.

Time Discretization
The time segment [0, T] is split into successive time-steps [t n , t n+1 ] with ∆t ≡ ∆t n = t n+1 − t n > 0 and n ∈ N. We further identify any variable or parameter ψ at time t n as ψ n . In this study we first consider a semi-discrete, first-order, semi-implicit time discretization scheme, where discrete spatial operators are written in an explicit form, except in the pressure terms. The pressure terms are made implicit for numerical stability reasons and to make use of larger time steps, choice of which is restricted by CFL-like [29] conditions based on material velocities. This gives rise to a coupled, nonlinear system of semi-discrete equations: Here superscripts n and n + 1 denote the term at previous discrete time t n (explicit terms) and next discrete time t n+1 (implicit terms), respectively. The above system of partial differential Equations (11)-(16) is in the semi-implicit form in time. With regard to the time scheme, we have taken a two-step approach for the numerical solution of the system (11)- (16). These are:

1.
Explicit convection/advection step: This step deals with the numerical calculation of all explicit terms in the system (11)- (16). This step is also called a convection/advection step (referred to as advection step hereafter) because these explicit terms can be determined by simply solving the advection problems for u φ , α φ , α φ ρ φ and α φ ρ φ e φ , as shown below: Time discrete system for the explicit advection step where superscript c denotes the advected quantities.

Remark 1.
In Section 5 we will focus on the application of the sharp interface capturing strategy for the space discretization of divergence operators in Equations (18)-(20).

2.
Implicit step: This step deals with the numerical solution of only implicit terms in the system (11)-(16) after calculating the explicit ones by solving the explicit advection step. To rewrite the system (11)- (16) in form of only implicit quantities, we perform the following substitution: (19) in (11); (17) in (12); (20) and (18) in (13). We consequently get the system of equations containing only implicit terms at time t n+1 : At this point it is worth mentioning that the time scheme presented here is a modified version of the so-called Implicit Continuous Eulerian (ICE) scheme [30,31]. But unlike ICE, where velocity variables (u φ ) are made implicit in the advection equations of volume fractions (α φ ), partial masses (ρ φ α φ ) and partial energies (ρ φ α φ e φ ), we have used instead explicit time discretization. Indeed, in this first approach, we made the choice to solve only pressure terms in an implicit way because, as explained in [32] for a two-phase flow system, the characteristic time scale for the pressure propagation is much smaller than that of advection. Thus, choosing an implicit time discretization for pressure will avoid the choice of very small time steps for its numerical resolution.
Next, rewriting Equations (21)- (26) in the form of phase volumes (V φ = α φ V, where V denotes the sum of volumes of all the phases) and mass (m φ ) is preferred: where The first relation given by Equation (27) will be solved after solving the following ones in order to compute partial mass density for each phase (see Section 4). Notice that the operation m n+1 forbids the use of a zero phase volume. Thus, we summarize the final time discrete system for the implicit step with the new Time discrete system for the implicit step This time discrete system in this implicit step is coupled and non-linear. We generically denote it for any phase φ by: where Q time represents the system of Equations (33)- (37) in the implicit step, and discretized only with respect to time.

Space Discretization
In our study we deal with the Finite Volume method (FVM) [33] discretization of Equations (33)- (37). It is important to remark at this point that here in this section, for the sake of conciseness, we prefer not to highlight the numerical strategy for the FV discretization of the divergence operators in the advection problems' set (17)-(20), we reserve that analysis for Section 5. Further, the two-dimensional (2D) and three-dimensional (3D) versions have been implemented, but at this step, only a 2D version has been validated and thus presented in this paper. We thus consider a two-dimensional (2D) uniform Cartesian mesh [34] with staggered variables for space discretization (in the spirit of Raviart-Thomas finite elements), in order to avoid the pressure oscillations due to the checkerboard problem [33]. Thus, as usual, the normal velocity variables are stored at faces of control volumes, whereas all other scalars (volume fractions and thermodynamic variables, generically noted by Π hereafter) are stored at cell centers, as shown in Figure 1 [35].
After performing a finite volume space discretization of time discrete model (33)-(37), we get a non-linear and coupled algebraic model, which generically can be written as: where Q time space represents the system of Equations (33)- (37) in the implicit step, and discretized with respect to time and space.
Further, to find the numerical solution of the above model on 2D Cartesian meshes, we have used the Newton's method [36] from PETSc library (see Section 4), and thus for any phase φ we can express it in the general form of a solvable linear system at time t n+1 : where with superscript k refering to the k th Newton iteration, and X 0 is an initial approximation given by the values: is a block and sparse Jacobian matrix of the non-linear and coupled system (39). Again, block nature of the Jacobian matrix highlights the fact that the algebraic system (39) is coupled.

Software Architecture for Numerical Resolution of The Model
To compute the solution of the nonlinear, coupled multi-phase model, we use our in-house high-performance software platform MMICADO: Mesh based Multi-phase Interaction Capabilities for severe Accident Dedicated cOdes. MMICADO is a CEA's proprietary software designed by following the Object-Oriented Programming (OOP) principles in the modern C++. MMICADO provides a flexible framework for meshes, parallelism (using DUNE library [37]), numerical methods, input data parsers (TinyXML library [38]), solvers (mainly with PETSc library [39]), primarily for developing any new multi-phase application. An overview of such a platform is shown in Figure 2. For example, to solve system (39) we have used PETSc's Linear Search Newton (SNESNEWTONLS) [40] method (a non-linear solver) within the framework of MMICADO, where the Jacobian matrix is computed numerically using a finite difference approximation. This is followed by the use of PETSc's linear solver KSPPREONLY [41] for solving the final global linear matrix, which uses the LU factorization (specifically the PETSc's PCLU [41]) as a pre-conditioner only once. Having said that, we now summarize the time algorithm in Figure 3 showing the sequence in which the model is solved in MMICADO.

Remark 2.
To get the numerical solution of the system (39), the following CFL condition based on the material velocities should be fulfilled: where η ≤ 1 2 , ∆t is the discrete time step and h is the constant distance between cell centers in a 2D Cartesian mesh.
In the next section we will show the details of the interface sharpening approach for the advection step shown in Figure 3. To ensure the sharp interface capturing in the solution of the compressible multimaterial model, the reconstruction of the scalars Π(α φ , ρ φ , e φ ) in each finite volume cell during the advection step depicted in time algorithm (see Figure 3) is done using a secondorder multidimensional MUSCL slope reconstruction strategy [19], where the predicted slope/gradient is limited using a recently developed, compressive, shape-preserving limiter, referred to as Multidimensional Limiting Process-Upper Bound (MLP-UB) (see [18]), which is an updated version of the MLP methodology detailed in [20]. We briefly elaborate this new MLP-UB method for a generic scalar variable z(x, t) advected in a medium with velocity u: Let us use K as a generic notation for a given finite volume cell, A as an edge (or a face in 3D) of this volume K, and a vector ν A being a normal unit vector to A pointing out of the cell K (see Figure 4). We look for the semi-discrete finite volume schemes of the form: where z K is the value of z defined at the center of a cell K, and z A (resp. u A ) denotes a certain value of z (resp. u) at the midpoint x A of an edge A.
As usual for a finite volume scheme, flux is linearly reconstructed in each cell, and a slope limiter is applied to limit this reconstructed gradient. The general expression of a MUSCL [19] linear reconstruction of a discrete field z K in cell K is: where x K is the center of mass of cell K, (∇ h z) K is a predicted estimate of the gradient, and Θ K is a positive scalar intended for limiting this gradient. Next the reconstructed value z A at the midpoint x A of an edge A in cell K is given by: which, in case of a Cartesian mesh with h as a constant distance between cell centers, can be simplified as: Main steps for the calculation of gradient limiter (Θ K ): Here we summarize main steps of the MLP-UB method. Detailed analysis of these steps can be found in [18].
Step 1: Prediction of the gradient In this study, gradient is predicted with a 8-point formula (see [18]). This 8-point formula for calculating the gradient of a scalar z K in cell K of a 2D Cartesian mesh (see Figure 5) is: Step 2: Calculate limiters at cell corners (c) In order to define the scalar limiting functions Θ c for each corner of a cell K, we compute a non-limited (NL) reconstructed value with the predicted centered gradient: We denote it byẑ K (x) = I NL K z(x) for convenience in the next formulae, and we define the extremum value z extr c as: depending on whether the non-limited reconstructed value is greater or smaller than the cell value. Here,z c and z c are the local extremum values at each corner (c) determined from the different cell neighbors around a cell corner (c) (see Figure 5) by the formulae: where Cells(c) refers to the set of neighbor cells around a cell corner c. The formula for the corner limiting function is then given as (see [18]): In case of a zero denominator in the above formula, Θ c is set to zero, which corresponds to zero gradient. In order to provide the method an ability to limit compressive effects, a threshold value of β is introduced (see [18]), such that the centered gradient can be multiplied at most by β. If β = 1, no compressive effect is allowed, otherwise β ∈ (1, 2], and some compressive effect is allowed; the value β = 2 corresponds to the most compressive effect. Remark 3. β ≤ 2 bound can be interpreted as the multidimensional equivalent of Ultrabee [42] stability bounds in 1D.

Step 3: Calculate limiters at cell centers (K)
Finally, limiters at the cell centers are given by the following formula (see [18]):

Reformulation of the Advection
Step with the MLP-UB Method for the Sharp Interface Capturing In this section we describe the methodology for the application of the above discussed MLP-UB [18] method to the advection Equations (18)- (20) of scalars Π(α φ , ρ φ , e φ ). The main points are discussed in the following sub-sections.

Reconstruction of Primitive Form of Scalars Instead of Their Conservative Form
In the present study, although we work with the conservative form of advection equations for α φ , ρ φ , e φ (see Equations (19) and (20)), we perform the reconstruction of primitive scalars α φ , ρ φ , e φ with MUSCL followed by the gradient limiting with MLP-UB to get their convective fluxes. This is mainly done to serve two purposes:

1.
As shown by [27,43,44], the application of the MUSCL on the primitive form of scalars leads to a more accurate reconstruction at the interface.

2.
Reconstructing α φ , ρ φ , e φ in their primitive form provides us the freedom to use different types of gradient limiters for α φ and ρ φ , e φ . Here, we prefer to use a compressive flux limiter (MLP-UB in present case) for α φ , and a less compressive one for ρ φ and e φ individually. For instance we have used MINMOD [45] type limiter for ρ φ and e φ .
Further elaborating on point 2 given above, we choose to do so because regular wave solutions shown by ρ φ and e φ in multi-phase test problems are prone to suffer from staircase effect (demonstrated in Section 5.2.2) when reconstructed with compressive limiters. However, the evolution of α φ maintains the initial discontinuous solution, and thus to preserve this sharp discontinuity, we use a compressive limiter like MLP-UB.

Modification of MLP-UB Method to Include the Functionality of the Minmod Limiter: 1D Analysis
Here we show a 1D analysis which shows the modification in Formula (45) in order to integrate both MINMOD and MLP-UB methods into a single formula. The 1D version of (45) reads: We introduce a new multiplicative parameter θ in the above formula which results in: In our analysis, we have tested several values of β and θ, and have observed that β = 1 and θ = 0.5 corresponds to MINMOD limiter. In Figure 6 we can see that a combination of a smooth and hat function is chosen as a initial condition, and subsequently it is linearly advected, where for gradient limiting after MUSCL reconstruction in each finite volume cell, we use relation (48) with two set of values for (β, θ) ∈ {(1, 0.5), (2, 1)}. The time scheme used is an explicit second order accurate Heun's RK2 method [46]. We see that the curve produced with MINMOD limiter overlaps with the curve obtained by using relation (48) with β = 1 and θ = 0.5.
It is also interesting to see that the advection of smooth function with a compressive limiter (using β = 2 and θ = 1 in relation (48)) leads to the staircase effect (see green curve on left-hand side in Figure 6). Thus, we choose to only apply MLP-UB (β = 2 and θ = 1) on α φ , which initially is discontinuous, and it is necessary to maintain this discontinuity as it gets advected at later times. However, for ρ φ , e φ , we perform their individual reconstruction by using MINMOD limiter (here MLP-UB with β = 1 and θ = 0.5).
Further, we also extend this analysis to 2D, and subsequently modify the advection step (see Figure 3) with the application of MLP-UB on scalars Π(α φ , ρ φ , e φ ), as shown in Figure 7.

Two-Dimensional (2D) Air-Helium Shock-Bubble Interaction: Experimental Validation
We focus here on the comparison of the numerical results produced with our method for the case of an air-helium shock-bubble interaction with one of the classical experiments shown in [47]. The goal, however, is not to have a precise quantitative comparison, but a qualitative one, to show that the method provides physically relevant results with a good resolution of the interfaces, especially when the coarse meshes are used. To study the case, a suitable computational domain is designed (see Figure 8). The dimensions of the computational domain are given in Table 1.
In order to simulate a Mach 1.22 shock wave similar to the one generated in the experiment, we consider the initial conditions shown in Table 2. We also assume that both air and helium are governed by a perfect gas EOS: where  Further, the CFL number for this problem is 0.2. We have used a coarse Cartesian mesh of size 294 × 42. Homogeneous Neumann boundary conditions are used for scalars (Π), and slip boundary condition is used for velocities. It the present case we have assumed symmetry in y-direction, thus only the upper half of the computational domain is considered.
Focusing on the numerical results shown in Figure 9, we observe the interaction of a Mach 1.22 shock wave with helium bubble while the shock travels from right to left in the computational domain. The numerical results shown in Figure 9a,c,e are at the three interaction times.These approximately correspond to the same time when shadowphotographs (see Figure 9b,d,f) are taken in the experiment.
Qualitatively, the numerical results correctly depict the deformation of the lighter helium bubble by an air-shock until the end of simulation time. It is interesting to see that even for a relative coarse mesh of size 294 × 42, we get a satisfactory phase separation. This is mainly because of the good interface sharpening effect produced by using MLP-UB method [18]. The results are also free from any interface oscillations or artifacts. We also observe that the values of cell volume fraction (α) for both air and helium are in the range [0, 1] and thus bounded. Some slight differences in the topology could be observed between the numerical results and experimental graphs. This difference can be explained by the absence of a surface tension force in the current numerical model.  (Figures (a,c,e)) and experimental results (Figures (b,d,f)) [47]. The scale of (α) is shown in Figure (g).

Two-Dimensional (2D) Triple Point Case: Two Materials and Three Phases
In this section we show a 2-material and 3-phase triple point test case. This is a standard test case corresponding to a 2D Riemann problem in a rectangular computational domain, which together with the dimensions and initial conditions is displayed in Figure 10. We also consider that the two materials are governed by an ideal gas EOS given by Equation (49). We choose a coarse Cartesian mesh of size 140 × 60. The boundary conditions for this case are Homogeneous Neumann for scalars (Π), and slip for velocities.
Our analysis for this test case primarily focuses on the comparison between the numerical results produced by using two different numerical recipes for the advection of scalars Π(α φ , ρ φ , e φ ): 1.

2.
MINMOD for all Π(α φ , ρ φ , e φ ). Firstly, we can observe in Figure 11 that two shocks propagate with different speeds due to the initial density difference between the left, top-right and top-bottom domains (see Figure 11). This creates the shear along the initial phase interfaces, and results in the formation of a vorticity. Regarding the coarse mesh used to compute this vorticity, we observe in Figure 11b,f,j that the numerical solution is diffused due to the use of MINMOD limiter in the numerical advection step of scalars Π(α φ , ρ φ , e φ ). However, if we focus on Figure 11a,e,i; a significantly better capture of this vorticity can be seen. Numerical results also show an accurate capture of the shocks and multi-phase interfaces. Again, like the previous test case, we observe that the values of cell volume fraction (α) for all the phases are in the range [0, 1].

Two-Dimensional (2D) Sod-like Shock Tube: 3-Material Shock-Bubble Interaction
In this 3-material test case we have simulated a situation where a bubble of some weakly compressible heavy material is surrounded by a vapor film. We have numerically studied the interaction of a shock, generated in the Sod-like shock tube [48], with this bubble-film configuration. The 2D computational domain for this case is displayed in Figure 12, where zone 4 is a bubble of heavy material governed by a stiffened gas EOS: where zone 3 here represents a vapor film which, along with zone 1 and zone 2, is governed by a perfect gas EOS Equation (49). The initial conditions and dimensions of the problem are shown in Tables 3 and 4, respectively. The boundary conditions for this case are same as the previous ones. The size of Cartesian mesh for this case is 800 × 400.    From the results shown in Figure 13 we can see that after an interaction with the shock, vapor film detaches from the heavy material bubble, and consequently overturns on the right-hand side of it. It then mixes with the heavy bubble interface as time progresses. We can also see a complete dissociation of two small parts of the vapor film. These two breakaway structures can be seen to have both translational and circular trajectories as they travel towards the right boundary. At the end of simulation they break into two more parts due to their interaction with the reflected shock wave from the right boundary.
In this case, mesh is sufficiently fine to capture the occurrence of the Richtmyer-Meshkov [49,50] instabilities between the interface of the light vapor film and heavy material bubble. This is evident from the mixing zone created between the two interfaces at later times. These instabilities occur in the present case because of the perturbed shape of the interface, which is mainly due to the initial representation of a circular interface by a piecewise constant profile, in each finite volume cell. After a shock interacts with multi-material interfaces, this initial numerical perturbation grows in time in a similar fashion to the growth of perturbations in the physical Richtmyer-Meshkov instability.
Thus, in this complex 3-material case as well we see the robustness of our method, which ensures the sharp interface capturing for the intricate evolution of the multi-material interfaces. Furthermore, in this test case, the values of cell volume fraction (α) for all the three materials is observed to be bounded in the range [0, 1]. (e) t = 250 µs.
(k) Figure 13. 2D numerical results for the 3-material shock-bubble interaction test case using a Cartesian mesh of size 800 × 400. Values of α in zone 3 and zone 4 are shown together. Thus, its scale (see Figure (k)) is scaled up accordingly.

Conclusions and Scope for the Future Work
In the present work, we have considered the hyperbolic compressible Euler system of equations, which can serve as the basis for numerical simulations of the several industrial applications in the domain of the compressible multi-material flows. We have worked with this model in the limit of a full mechanical equilibrium, and also considered an instantaneous pressure equilibrium. Within the scope of such a mathematical model, we focused on the problem of the interface sharpening, where we chose to work with the Interface Capturing approach. We performed the analysis of interface sharpening using a recently developed Diffuse Interface Method: Multidimensional Limiting Process-Upper Bound (MLP-UB) [18]. This method has been proposed in a recent paper [18], where the authors have demonstrated its robustness for the pure advection cases. The present paper is thus considered as a natural sequel to the work presented in [18].
For performing the numerical simulations in this work we have used our in-house, high-performance, modern software: MMICADO; within its framework, we have first proposed a robust time scheme, where we first explicitly calculate the advection/convection terms. Within the framework of Finite Volume space discretization on 2D Cartesian meshes, the divergence terms in these advection/convection equations are computed by employing a second order accurate, multidimensional MUSCL [19] slope reconstruction strategy along with the MLP-UB gradient limiting method. A 1D analysis is also presented to highlight the choice made for scalar reconstruction in their primitive form (instead of their conservative form) by using the most and the least compressive form of MUSCL+MLP-UB approach.
Finally, we have provided three complex multi-material test cases, where our aim was to produce accurate results on the coarse 2D Cartesian meshes. Looking at the results we can conclude that our approach has produced:

1.
Accurate and sharp capturing of the intricate interface evolution without numerical artifacts.

2.
Accurate interaction of the shocks with the multi-material interfaces.

3.
Bounded numerical results for any number of materials.
A direct extension of the work presented in this paper is to take into account the phenomena of phase change by considering the heat and mass transfer across the multi-material interfaces. Another work which is conceptualized is the study of sliding interfaces, where we want to consider the relative motion between the multi-material interfaces. Another immediate work concerns the verification and validation of the 3D version, and theoretical demonstration of the positivity preserving property (boundedness) of the method. We plan to address these issues in the future research articles.  Data Availability Statement: The present paper contains all the necessary data needed for the thorough understanding of the work. In case a requirement for any supplementary material arises, the readers are encouraged to contact the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.