The diffuse interface models, part of the family of the front capturing methods, provide an efficient and robust framework for the simulation of multi-species flows. They allow the integration of additional physical phenomena of increasing complexity while ensuring discrete conservation of mass, momentum, and energy. The main drawback brought by the adoption of these models consists of the interface smearing, increasing with the simulation time, therefore, requiring a counteraction through the introduction of sharpening terms and a careful selection of the discretization level. In recent years, the diffuse interface models have been solved using several numerical frameworks including finite volume, discontinuous Galerkin, and hybrid lattice Boltzmann method, in conjunction with shock and contact wave capturing schemes. The present review aims to present the recent advancements of high-order accuracy schemes with the capability of solving discontinuities without the introduction of numerical instabilities and to put them in perspective for the solution of multi-species flows with the diffuse interface method.

## I. INTRODUCTION

Recent research in the area of multiphase and multi-species flows has been successfully applied to a wide range of problems in engineering and sciences, for example, detonation of high energy materials,^{1,2} shock–bubble interaction,^{3,4} propellant mixing in high-pressure combustors,^{5,6} tumor growth modeling,^{7} and fluid−solid interaction.^{8,9}

The main difficulty, especially in the context of compressible multi-species flows, is to maintain the thermodynamic consistency at material interfaces as an addition to the thermodynamic properties of the separated phases, while avoiding numerical instabilities across discontinuities.^{3,10,11} Typically, a material interface separating two different phases is on the order of few nanometers, and, in an engineering context, can be considered as a sharp discontinuity. However, the way the interface is treated and how its position is determined plays a key part in multi-species flow modeling, and differentiates the numerous methods introduced over the last three decades. Following a physical rigor principle, one can consider two families of methods, one that rigorously does not allow interface mixing, namely, the sharp interface methods (SIM),^{12,13} and the other one that considers the interface as a diffused zone, namely, the diffuse interface methods (DIM).^{14,15} Another possible classification is based on how the interface position is determined, whether through a series of interpolations following a mesh adaptation (front tracking methods^{13} and Lagrangian methods^{16}) or resolving the material interface on a fixed mesh through a reconstruction process or with the aid of a Riemann solver (front capturing methods^{2,17} and Eulerian methods^{18,19}), pure Lagrangian methods are not feasible for the large deformations often encountered in hydrodynamic flows,^{20} due to the geometric constraints affecting the mesh deformation. To circumvent these limitations, Hirt and Nichols^{18} proposed an Eulerian approach, namely, volume of fluid (VOF), where the mixture of fluid on each computational cell is represented by a volume fraction advected in time, and the interface boundaries are located through the derivatives of the volume fraction function. In the method of Unverdi and Tryggvason,^{21} developing for the computation of incompressible and immiscible flows, a moving frame overlaps a fixed mesh tracking the position of the interface through a series of marking points. The interpolation of these marking points provides the marking function used to retrieve the interface location. Although this front tracking method allows a sharp interface definition, it is still challenging to be adapted to complicated topological changes and large number of interfaces.

On the other hand, the level set method (LSM) developed in Refs. 22 and 23 is considered easier to be implemented and uses a distance function from the interface, where the zero level marks the interface location, and the positive values correspond to one of the fluids and the negative values to the other one. In the original version, this method was affected by spurious oscillations that were corrected successively in the work of Fedkiw *et al.*^{24} In this technique, named ghost fluid method, a ghost cell storing the properties of the fluid **b** is defined in correspondence of each computational cell storing the properties of the fluid **a**. After each time step, the distance function, calculated with a level set equation, determines which set of properties has to be used to compute the solution, whether the one stored in the ghost or the computational cell.^{25,26} This set of methods has the advantage to allow a sharp definition of the interface for fluids with different equation of states (EOS). However, general drawbacks they share are the difficulty to enforce the mass conservation at the discontinuities discretely and the complexity in the extension to more dimensions.^{27,28}

A method that gained considerable popularity in recent years is the DIM. This method allows an artificial mixing of the fluids in the interface region and thus requires a thermodynamically consistent mixing model to regularize the state evolution of the interface region without the introduction of spurious oscillations.^{3,29} This requirement is relaxed when the two fluid properties are similar but becomes more stringent for flows involving different phases. This problem was addressed by Karni^{10} and Abgrall^{30} who first proposed an Eulerian formulation of the DIM valid for perfect gases, while in Refs. 15 and 29 that was extended to fluids governed by stiffened EOS and in Ref. 31 for Van der Waals gases. The extension to general EOS is attributed to the work of Allaire *et al.*^{3} Two subclasses of DIM have emerged in literature, the first of which, also named phase-field method, considers a visco-capillary structure (Cahn–Hilliard^{32}) for the mixture cells and requires a fine mesh resolution to resolve the interface, being more suitable for the study of smaller scale systems. The second subclass is to some extents an evolution of the VOF method, where the classical system of conservation laws is extended with additional phase mass conservation and volume fraction advection equations. Methods based on the reduced seven-equation Baer and Nunziato model,^{1} introduced first for the simulation of detonation-to-deflagration transition of highly reactive materials, have found a wide range of applications in the context of compressible multi-species flows together with the four-equation model of Abgrall^{30} and the five-equation models of Allaire *et al.*^{3} and Kapila *et al.*^{2}

The versatility of the diffuse interface models favored the integration of additional physical phenomena such as viscosity and capillarity (Perigaud and Saurel,^{33} Saurel *et al.*^{34}) and elasto-plastic solid–fluid interaction coupling a visco-plastic model to the Euler equations as in Favrie and Gavrilyuk,^{8} recalling the work of Godunov^{35} and Godunov and Romenskii.^{36} The advantage of this setting is that the same system of equations is applied on the whole domain avoiding any need for interface tracking or reconstruction, and the entropy, mass, and momentum are discretely conserved.^{11,37} Moreover, this is one of the few methods able to recognize the formation or collapse of the interfaces.^{2,38} The inherited drawback of the DIMs is the time-increasing interface smearing, which limits the simulation time horizon before incurring in excessive diffusion. Different strategies have recently been developed to limit the interface diffusion, for example, introducing diffusion-sharpening terms as in the work of Shukla *et al.,*^{39} Tiwari *et al.,*^{38} and Jain *et al.,*^{11} or using the tangent of hyperbola for interface capturing (THINC) approach, developed first for the VOF method by Ii *et al.*^{40} and adapted to DIM by Shyue and Xiao.^{41} Contemporary, coupling of high-order reconstruction scheme, typical of single-fluid finite-volume (FV) codes, to the interface-sharpening technique became common practice,^{37,38} and more recently, applications to unstructured meshes also appeared.^{42} Classical reconstruction schemes such as monotone upstream-centered scheme for conservation law (MUSCL)^{43–45} and weighted essentially non-oscillatory (WENO)^{46,47} reached a remarkable maturity within the single-fluid finite-volume framework, allowing high-order accuracy in smooth regions on both structured^{46,48} and unstructured grids with arbitrary elements.^{49} Great effort was devoted to enhance the robustness and efficiency of the stencil selection procedure^{50,51} in order to decrease the computational load, while maintaining a high-order of accuracy.

A fifth-order WENO reconstruction for two-dimensional compressible, multiphase flows was already adopted in the work of Johnsen and Colonius^{4} and Tiwari *et al.,*^{38} and a three-dimensional extension based on the implementation of Titarev and Toro^{52} was adopted in the work of Coralic and Colonius.^{53} In the class of high-order methods, the finite-element discontinuous Galerkin (DG) formulation, introduced by Cockburn and Shu,^{54–56} is rapidly gaining popularity for its more straightforward formulation of high accuracy schemes based on the higher-order polynomial used for the representation of the piecewise solution. The initial drawback of presenting more degrees of freedom (DOFs) per cell resulting in more unknowns and higher computational times compared to the FV framework for the same grid seems to be mitigated by the introduction of the reconstructed discontinuous Galerkin (rDG) schemes by Dumbser *et al.*^{57} Here, a polynomial of degree N is used as a piecewise representation of the actual solution, while the fluxes and the time updated solution are evaluated with a polynomial of degree $M\u2265N$, reconstructed from the underlying polynomial of degree N. In this scheme, the unknown reconstructed polynomial is first projected to the solution space, that is, the DG solution in the contributing cells, making it indistinguishable from the initial solution in the weak sense, and the reconstruction is performed via least-square procedure operating on the overdetermined system obtained by extending the reconstruction stencil. This scheme, called $PNPM$, is a generalization of both the FV and DG frameworks, as in the limit of *N* = *M* the scheme is equivalent to the pure DG formulation, while for *N* = 0, this falls in a classical FV scheme. In the work of Luo *et al.,*^{58} the costly projection step that requires numerical integration is avoided by requiring the reconstructed polynomial and its first derivatives to be equal to the solution and its first derivatives of all the elements adjacent to the target cell, while in the work of Zhang *et al.,*^{59} the higher-order derivatives are obtained through a Green–Gauss procedure^{60} applied to the first-order derivatives, in order to improve the computational efficiency. Luo *et al.*^{61} extended the $PNPM$ scheme to unstructured tetrahedral grids employing an Hermite WENO^{62} reconstruction obtaining oscillation-free solutions. For $M\u2265N$, it is demonstrated that the *rDG* schemes achieve accurate solutions with a minor impact on the computational time. The application of this FV-DG hybridization to the compressible multi-species hydrodynamics was proposed by Pandare^{63} who used the approach presented in Ref. 58 in conjunction with a WENO limiter, obtaining a third-order solution.

On the other hand, another framework, that is, becoming more widespread for the simulation of multiphase flows, is the lattice Boltzmann method (LBM), as the continuum assumption is ignored, and instead, a mesoscopic kinetic operator is used to connect the molecular physics to a set of macroscopic properties obeying the Navier–Stokes (NS) equations. The LBM can be easily parallelized^{64,65} and adapted to complex topologies^{66,67} and presents various strategies for the computation of multiphase flows. The main multiphase LB models present in literature are the Color Gradient Model (Gunstensen *et al.*^{68} and Rothmann and Keller^{69}), the Shan-Chen model, also known as pseudo molecular interactions model,^{66,70} and the free energy model (Swift *et al.*^{71} and Zheng *et al.*^{72}) numerical tests^{66,73,74} indicate that the color Gradient method provides the sharper fluid interfaces, although it is prone to numerical instabilities due to the lattice pinning, a phenomenon that also causes incorrect interface advection.^{73} As pointed by Qu,^{75} the LBM is commonly limited to incompressible flows since the model is derived from a low Mach number expansion of the Maxwellian distribution. Qu derived a compressible model valid for a wide range of Mach numbers applying a FV scheme for the solution of the discrete velocity Boltzmann equation (DVBE) as in the work of Nannelli and Succi^{76} and adopted an exact Riemann solver with a MUSCL scheme to capture the discontinuities. Contrarily, in the work of Joshi *et al.,*^{77} where a hybrid LBM is built for compressible multi-species flows, the intercell fluxes are evaluated through a LBM and a FV scheme is used to update the node parameters. In addition, with such blending of the LBM with the FV schemes, another limit of the former method can be overcome, that is, the difficult applicability to unstructured meshes. Applications on irregular grids are becoming, in general, a necessity in industrial applications, especially when complicated geometries are considered and present a number of advantages over structured grids, ranging from the ease of mesh smoothness requirement,^{78,79} to the maintenance of load balance in parallel computing.^{80,81} A detailed review of the application of different high-order methods on unstructured meshes is presented in Ref. 79.

Given this variety of schemes for the solution of multi-species flows and the recent introduction of hybridizations of the classical FV schemes with both the DG and LBM frameworks, this review proposes, without the ambition to cover the totality of the available computational frames, to gather the recent advancements in high-order accuracy schemes on unstructured grids and to put them in perspective for the simulation of multi-species flows with the diffuse interface method. A recent review of the main DIM models is available in Ref. 27; in this work, we focus on the implementation and adaptations required by these models within different numerical frameworks, emphasizing the challenges and advantages relative to each scheme, and comparing their efficiency and accuracy. The review is organized as follows: the DIM is presented and the variations to the complete Baer and Nunziato^{1} and the four-equation model of Abgrall^{30} are listed. The classic FVM (finite volume method), DGM (discontinuous Galerkin method), and LBM are briefly discussed and used as a baseline for the introduction of the hybrid schemes. Finally, the available to-date applications of DI models in conjunction with sharpening techniques on different numerical frameworks and hybrid schemes are presented.

## II. DIFFUSE INTERFACE MODELS

### A. Seven-equation models

A cornerstone in multi-species modeling is the model introduced by Baer and Nunziato,^{1} from here on abbreviated with the Baer and Nunziato (BN) model, which was originally proposed for the description of reactive multiphase flow, and, in particular, the deflagration to detonation transition (DDT). In this process, the interstices in a granular reactive bed are filled with a gas phase forming an immiscible mixture. Each phase is treated as a compressible fluid in local thermodynamic equilibrium, but the mixture is allowed to be in nonequilibrium across the interface, and each phase is described through a set of inviscid conservation laws, which also accounts for the interaction between phases. As pointed by Bdzil *et al.,*^{82} this implicitly means that the self-equilibration time scale of each phase is much smaller than the equilibration time scale between the phases. Each conservation law includes a source term modeled to represent the transfer of mass, momentum, and energy between the phases, here named $\Omega \u0303,\u2009\Gamma \u0303$ and $E\u0303$, respectively.

An additional equation is necessary to evaluate the volume fraction of the phases, given the saturation condition $\varphi s+\varphi g=1$ of the solid volume fraction $\varphi s$ and gas volume fraction $\varphi g$. The one-dimensional system for the multi-species flow model can be expressed in vectorial form,

where *U* is the vector of conserved variables, *u* is the velocity field, *F* is the flux vector, and *H* and *S* are the vectors of the non-conservative quantities.

In presenting the models available in literature, the classic notation will be used: *E _{k}* represents the total energy with $Ek=ek+12uk\xb7uk$ and

*e*is the internal energy. Velocity and pressure of phase

_{k}*k*will be, respectively,

*u*and

_{k}*p*.

_{k}With this notation, a two-component flow model in the BN formulation contains six equations, one for the conservation of mass, momentum, and energy of each phase in addition to a volume fraction evolution equation, leading to the well-known seven-equation model, that in the modified version of Bdzil and Kapila^{2,82} reads

For the sake of simplicity and compactness, all the non-conservative terms have been included in the *S* vector leaving in this case *H* as a zero vector. The nature of the source terms is rather complicated and depends on the physical process under examination. For instance, in the work of Baer and Nunziato, the terms $\Omega \u0303,\u2009\Gamma \u0303$ and $E\u0303$ are the balances related to the burning and drag process of the solid phase indicating therefore the exchanges of mass, momentum, and energy-heat transfer, respectively, during flame spread. The summation of the source terms between solid and gaseous phases is null; hence, the total mass, momentum, and energy are conserved. It must be noted that in the original formulation of the BN model, the subscripts 1 and 2 are to be considered relative to the solid and gas phases, respectively. Moreover, an additional term, $\Phi \u0303$, is introduced in the volume fraction equation, that in the original work of Baer and Nunziato is referred to as the compaction equation, and includes the resistance effects of the granular configuration to the induced forces

The intergranular stress *β*_{1} measures the forces between the grains and estimates of this value can be found in the work of Elban and Chiarito.^{83} The quantity *μ _{c}* is defined as the dynamic compaction viscosity, and it is an indicator of the rate at which pressure equilibrium is reached.

^{15}A similar set of equations was derived by Saurel and Abgrall

^{15}applying an averaging process on the single-phase Navier–Stokes equations in order to build a model valid for immiscible fluids with resolved interfaces as well as two-phase mixtures. Like the BN model, this system of equations is unconditionally hyperbolic; however, in order to solve interfaces between immiscible fluids, velocity and pressure equilibrium are required across the interface. If this condition is not enforced, the characteristic directions corresponding to the volume fractions, which abruptly change from 0 to 1 across the mixing zone, would emanate degenerate waves in correspondence of the interface. Therefore, the boundary problem is ill-adapted unless a proper closure is applied at the interface. Moreover, a model of this form containing non-conservative terms does not admit Rankine–Hugoniot relations. Saurel and Abgrall

^{15}addressed these problems by introducing an interface equilibrium condition and a relaxation procedure after wave propagation that tends to the equilibrium condition.

After wave propagation in a gas phase dispersed in a liquid phase, the pressure state is in initial non-equilibrium and tends to reach the equilibrium *P _{g}* =

*P*through a volume variation proportional to the compaction viscosity and to the difference between the gas phase pressure and the interfacial pressure

_{l}*P*,

_{i}Similar considerations hold for the velocity, whereas the relaxation process depends on the drag force *F _{d}*, which can be expressed in terms of a velocity relaxation coefficient

*λ*and the velocity jump across the interface,

It is reasonable to assume that the pressure equilibrium can be reached almost instantaneously in certain applications, meaning that the compaction viscosity tends to infinity. In that case, the solution has to be obtained through an hyperbolic step resolved without source terms,^{15} followed by pressure and velocity relaxation with coefficients $\mu c=\lambda \u2192+\u221e$.

Applying the conditions above to the system of Eq. (1), the vector of source terms in Eq. (2) can be expanded in the form introduced in Ref. 15,

Questioning the acoustic properties of the original BN model^{1} in situation of dilute limit, Saurel *et al.*^{84} proposed a model in which the sound speed correctly vanishes in non-continuous phases. Regarding the choice of the interfacial pressure and velocity *P _{i}* and

*V*, in the original BN model the two variables were chosen such that the first to be equal to the pressure of the gas phase, and the latter to be equal to the velocity of the less compressible phase. However, this choice has a considerable impact on the wave pattern and Saurel and Abgrall

_{i}^{15}proposed an estimate where the interfacial pressure is equal to the mixture pressure and the velocity is the one of the center of mass,

### B. Five-equation models

It is well known that the seven-equation model is the most complete and general model although it exhibits a complex wave pattern as a consequence of the source terms nature. The numerical behavior of the system and the solution is thus heavily conditioned by the relaxation procedure employed to deal with the different length and time scales for velocity and pressure equilibration. It is estimated, however, that the time scales for velocity and pressure equilibration is small compared to the flow characteristic time,^{2} leading to stiff source terms. Kapila *et al.*^{2} proposed an asymptotic reduction of the seven-equation model under stiff mechanical relaxation in order to obtain a model with a single velocity and a single pressure. The result is a five-equation model, that is, a zero-order approximation of the BN model, with a non-conservative volume fraction equation, two mass balances, and one momentum and energy balance equation.

In this context, the vectors of system of Eq. (1) are

In this case, the only source term was incorporated in the *H* vector, and thus, the *S* vector is zero. Here, *E*, *u,* and *p* are the mixture energy, the equilibrium velocity, and pressure, respectively. The mixture internal energy is defined through the mass fraction $Yk=(\varphi \rho )k\rho $ and the mixture density $\rho =(\varphi \rho )1+(\varphi \rho )2$,

Although this system is considerably simpler than the complete BN model, it still presents a series of challenges mostly related to the presence of the non-conservative source in the volume fraction equation and to the way the sound speed is calculated. The former condition causes the model to slowly converge, or not converge at all, to the exact solution. The latter condition refers to the fact that the equilibrium sound speed is calculated with Wood's formula,^{85}

which has a non-monotonic behavior and can be deleterious in the study of wave transmission as the mixture equilibrium speed of sound can eventually become smaller than the ones in pure fluids.^{34,38}

### C. Four-equation and five-equation averaging models

Another five-equation model was derived by Allaire *et al.,*^{3} generalizing the four-equation model proposed by Abgrall,^{30} initially limited to perfect gas EOS, to arbitrary EOS. This kind of models, named $\gamma \u2212$ models, is built upon mixture rules, which usually provides an additional equation to the classic Euler system, with the role of describing the dynamics of the fluid composition. A successful mixture model is generally built in order to meet the following criteria:^{3}

in the limit of a single fluid, the model should correctly degenerate to a single fluid model,

the model should be hyperbolic,

total mass and total energy of the fluids should be conserved,

in presence of interfaces, the model should maintain numerical stability.

Early models^{10,86} provided the additional equation in terms of the specific heat ratio *γ*, obtained from the weighting of the mass fractions

However, this formulation leads to oscillations at interfaces,^{10,30} that is, in contrast with the last of the criteria previously listed. The solution, proposed by Abgrall in Ref. 30, is to average with the quantity $1/(\gamma \u22121)$ instead of *γ*, and this four equation model reads

This model was originally built for perfect gas EOS and later extended by Shyue to the stiffened-gas (SG) EOS, Van der Waals EOS, and Mie–Grüneisen EOS.^{29,31} More recently, Lei and Li^{87} presented a non-oscillatory scheme based on an energy-splitting Godunov scheme without energy exchange between materials, and using isobaric hypothesis. It was observed that the number of transport equation in the $\gamma \u2212$ models depends upon the number of parameters in the EOS, increasing rapidly in complexity for realistic material EOS. To overcome this, Allaire *et al.*^{3} proposed to use a single volume fraction transport equation in place of the algebraic closure, in addition to mass conservation equation for each phase. The resulting model is a five-equation formulation similar in philosophy to the multiphase models in mechanical equilibrium, where the main difference compared to Kapila's^{2} model resides in the source term *H*,

The system is closed with an EOS, and in Ref. 3, isobaric and isothermal closures were initially proposed. In Ref. 3, it is demonstrated that under the isobaric closure, this model can capture fluid interfaces without introducing spurious oscillations, calling this closure a consistency condition. In addition, compared to the four-equation model where only the partial thermodynamic properties of the mixture could be recovered, the introduction of separate mixture mass balances enforces phase mass conservation and allows the calculation of the properties of the separated phases, making this model a very efficient and versatile setting for the simulation of multi-species flows. Both these kind of four-equation and five-equation models have been evolved over the years in order to include additional physical phenomena. For example, the effects of capillarity were included by Garrick *et al.*^{88} by adding capillary force terms to a five-equation model. The effects of viscosity and diffusion were incorporated to a five-equation model by Thornber *et al.,*^{89} while Yi *et al.*^{90} developed a four-equation model for two-phase flow with phase change including a phase-equilibrium solver.

### D. Six-equation models

Previously, we have mentioned some of the challenges involving the Kapila's five-equation model,^{2} regarding the numerical behavior. This model is able to treat the dynamical appearance or collapse of interfaces. However, due to the presence of the non-conservative term in the volume fraction advection [Eq. (8)], as well as the sound speed following Wood's formula,^{85} it is non-trivial to solve the system numerically and especially to enforce volume positivity of the phases. In its original form, this model is not easy to adapt to account for additional physical phenomena and to be discretized on unstructured meshes. To circumvent these difficulties, Saurel *et al.*^{34} proposed a reduced BN model, where pressure non-equilibrium is maintained. This procedure is similar to the first reduction that can be found in Ref. 2 in the limit of stiff velocity relaxation [pressure equilibrium alone leads to a mixed elliptic description not suitable for evolutionary partial differential equations (PDEs)] and results in a six-equation model that in the limit of stiff pressure relaxation reduces to the Kapila's model.

Using again the notation introduced in Eq. (1), the vectors for the six-equation model are

A mixture energy equation is added to the system in order to correct the internal energy equations in the presence of shocks

Defining the acoustic impedance as $Zk=\rho kck$, the interface pressure can be obtained from

and the mixture sound of speed can be computed as a monotonic function

It is immediately noticeable that the volume fraction equation is now a transport equation with a relaxation term, and the speed of sound is now monotonic, resolving the issues introduced with the Kapila's model. The new model, with the introduction of Eq. (15), results in an overdetermined system, and as pointed in Ref. 34, this has to be regarded more as a different procedure to solve the five-equation model rather than a new physical model. Details on the solving algorithms for this model can be found in Refs. 27, 34, and 91.

### E. Choice of the closure laws

The discussion of closure laws has been so far postponed. The choice of an appropriate EOS is clearly related to the nature of the phases and their expected thermodynamic evolution during the analysis. For instance, in Ref. 1 where the subjects of the study were granular reactive materials in gaseous phases, the authors selected for the granular material a thermoelastic formulation of the Helmholtz free energy, and for the gas phase a Jones–Wilkins–Lee (JWL) EOS, that can describe the state of detonated material as well as unreacted explosives fairly accurately. On the other hand, one of the reasons for the rapid gain in popularity of the five-equation models [Eq. (13)] over the four-equation models [Eq. (12)] is that the latter type was initially limited to specific EOS such as the perfect-gas EOS. In this frame, one of the most important features carried by the Allaire's model^{3} [Eq. (13)] is the generalization to real fluids EOS, and the demonstration of the good numerical behavior of the isobaric closure.

When inert gas–fluid systems are considered, it is nowadays a common choice to have the ideal-gas law for the gas phase

and the stiffened-gas (SG) law for the liquid phase

with a good balance in accuracy and in computational requirements. The model parameters *γ _{l}* and $\pi \u221e$ are determined through the shock wave Hugoniot data.

^{92}The SG model accounts for molecular and repulsive effects in a simplified way, leading to a disagreement between the experimental values and those computed with Eq. (19) for the liquid specific volume. A more accurate model is the noble-able-stiffened-gas model, which accounts for covolume and repulsive effects, and was introduced by Le Métayer and Saurel in Ref. 93 and reads

The constants $p\u221e,k$ and *b _{k}* are characteristics of the fluid, and

*q*is a reference phase energy. The values for dodecane, water, and oxygen are available in the previous reference. Other more general EOS are available, like the Mie–Grüneisen for condensed materials upon which other EOSs are built explicitly for compressible fluid and compressible elastic–plastic solid interaction. The Mie–Gruneisen EOS is often used to model metallic materials, as well as detonation products, and is of the form

_{k}with $\eta =\rho /\rho 0,\u2009\Gamma 0$ a material constant, *a*_{0} the bulk speed of sound, and *ρ*_{0} the reference density. The parameter *s* represents the Hugoniot slope coefficient, defined as $s=dUsdUp$, given the shock wave speed *U _{s}* and the particle velocity

*U*.

_{p}### F. Fluid–solid interface models

In recent years, the diffuse interface model was successfully applied to a series of problems involving materials undergoing impact loading,^{8,94} extreme deformations,^{95} and, in general, for studying the dynamic of fracture and fragmentation,^{94,96–98} that can also be enhanced by peridynamics models.^{99} As one can imagine, due to the inherited increasing-with-time diffusion of the fluid–solid interface of the DIM, this type of setting is more suitable for rapid-dynamics problems. The idea is to extend the models discussed above for mixtures of fluids, to include a hyperelastic description of a pure solid phase. An Eulerian formulation for elastic–plastic flows in conservative form, based on the inverse of the deformation gradient, was proposed by Trangenstein,^{100} although some difficulties to verify the curl gauge constraint, $\u2207\xd7gT=0$, of the inverse of the deformation tensor *g* were experienced. A solution for this problem was presented in Ref. 101 casting the governing equations for *g* in non-conservative form. A different approach uses the non-linear elasticity formulation introduced in Godunov and Romenskii,^{36} where the limits of the visco-plastic model are able to describe solids as well as viscous and inviscid flows through a plastic relaxation term included in the equation for the deformation tensor,^{102} or in alternative with the concurrent coupling of peridynamics and classical elasticity for elastodynamic problems.^{103}

The derivation of a solid–fluid diffuse interface model is performed in Ref. 95, and the procedure is executed in two steps: first, the Hamilton principle is employed to derive a hyperbolic system of motion equations, and thereafter, the elasticity effects are incorporated through a thermodynamically consistent, displacement-based, formulation of elastic materials. Such displacement-based description was explored in Refs. 36, 104, and 105 and is advantageous with respect to the hypoelastic description widely used in commercial codes, as it allows a divergence form description, favoring convergence properties and involves only material derivatives.^{105} Plasticity effects were included subsequently in the work of Favrie and Gavrilyuk^{8} through the splitting of the internal specific energy in a hydrodynamic part dependent upon density and entropy, and a shear part function of the invariants of the Finger tensor. Thus, the pressure is related only to the hydrodynamic part of the energy and the deviatoric part of the stress tensor to the elastic energy. The model results are compatible with the von Mises yield criteria.

Here, we will discuss the DI model for the coupling of compressible fluid–compressible elastic–plastic solid interaction introduced in Ref. 8. A brief introduction to the notation follows: let the Lagrangian coordinates be $X=(Xl),l=1,2,3$ and the Eulerian coordinates be $x=(xe),e=1,2,3$. If the deformation or diffeomorphism $y\u0303$, from a reference configuration Π to the actual configuration $\Pi \u0303$ is written as $x=y\u0303(t,X)$, the deformation gradient can be expressed as

Then, the Finger tensor, that is, the inverse of the left Cachy–Green deformation tensor is

which is related to the Almansi tensor *ε _{ij}* through

with *δ _{ij}* being the Kronecker symbol. A cobasis vector corresponding to the Lagrangian coordinates is $El=\u2207Xl$ and satisfies the following:

By applying the Hamilton principle, a pressure equilibrium elastic solid–gas diffuse interface model with plasticity effects can be written

Here, *G ^{e}* has to be considered as the effective elastic Finger tensor, a modification necessary when plastic dissipative effects are included through Maxwellian terms in order to satisfy the second law of thermodynamics, and the velocity field is defined as $u(t,x)$. The plastic effects are accounted for through the term

*L*and the stress tensors for solid

^{p}*σ*and gas

_{s}*σ*are, respectively,

_{g}The above system poses two challenges. The first is related to the equilibrium formulation, and looking at Eq. (27) one can notice the similarity with the non-conservative term representing the compaction effects in the volume fraction advection equation encountered in the Kapila's model Eq. (8). This, unitedly with the sound speed which originally follows the Wood's law Eq. (10), poses potential inconsistencies similar to the ones already discussed for fluid–fluid systems. In addition, because of the nature of DI models for which pure phases do not exist, in order to avoid the degeneracy of the phase-density equation, an infinitesimal part of one phase must be artificially added to the other phase. This may not affect fluid–fluid systems as much as fluid–solid systems since the wave structure for these two cases is quite different. For the latter case, the presence of a small quantity of solid phase in the gaseous phase can introduce unphysical shear waves. A solution was suggested in the work of Favrie and Gavrilyuk^{8} based on the relaxation procedure, which led to the six-equation model of Saurel *et al.*^{34} A dissipation function, $D\u0303$, is introduced:

and the relaxation equation reads

The resulting model, with the relaxation procedure, remains hyperbolic as far as a negligible quantity of phase mixing is present. Another solution for both issues, that may be considered a too crude approximation, but seems to be effective and simple, was proposed by Li *et al.*^{106} The phase-density equations are reformulated in such a way to be unrelated to the volume fractions, which is possible by neglecting the compaction effects, that is, the Kapila's non-conservative terms. Thus, the model is able to deal with pure phases but is not indicated for the analysis of cavitating processes. A different approach was presented by Barton,^{107} where instead of using a plastic relaxation process, the model presents a mixture equation-of-state, which includes the contributions of cold compression or dilatation, temperature deviation, and shear strain. With this formulation, the fluid phase is represented as a special case of the above taking the shear modulus as zero and resulting in no shear energy contribution. The plasticity effects are included with a thermodynamically compatible source term obtained following the method of convex potentials.^{108} This approach was extended by Wallis *et al.*^{109} in order to account for interface slide and void opening, by applying modified cell interfaces boundary conditions before the standard flux calculation. The interaction between solids and fluids, and potentially a combination of the four states of matter, can be accounted for in other different ways. Michael and Nikiforakis presented in a series of papers^{9,98,110} a numerical framework for the simulation of the interaction between inert, reactive multiphase fluids and elastoplastic solids. In these publications, in order to formulate a description of condensed phase explosives, they developed a hybrid formulation of the reduced BN model and the augmented Euler equations, which considers a phase 1 made of air that interfaces with a phase 2 consisting of a homogeneous mixture of an explosive reactant *α* with products of reaction *β*. The model considers velocity and pressure equilibrium, between phase 1 and phase 2, and also between the mixture components *α* and *β*, and temperature equilibrium only between the mixture components. The resulting five-equation models are similar to Allaire's model^{3} [Eq. (13)] considered in three dimensions, with an additional transport equation

where K is a source term expressing the rate of conversion of reactants to products, and *λ* is the mass fraction of the reactants. All materials are modeled by a Mie–Grüneisen EOS and the particularity of the model is that it reduces to the five-equation model of Allaire when the phase 2 is an inert constituent (*λ* = 1 or *λ* = 0 and *K* = 0) and to a two-phase reacting five-equation model [Eq. (8)] when the combustion products are not modeled. In the latter case, the energy equation is modified with an additional source term

where *Q* includes the heat of detonation, and together with the source term of Eq. (31) represents the chemical energy release. In the case that the inert material is not present, the model should further reduce to a mixture model similar to Eq. (12) with mixture rules defined as in Ref. 111. The energy transport is taken as in Eq. (32), and instead of a transport equation for the volume fraction, a transport for the mass fraction *λ* is considered:

The elastic behavior of the solid phase is modeled by the theory developed in Ref. 112, while the inelastic follows the model introduced in Ref. 101. This hybridization between multiphase and mixture models is of interest for physical phenomena ranging on different time scales and of different dynamics. For example, in a solid-explosive-fluid interaction the multiphase type model suits the granular detonation, while the gaseous combustion is better accounted for with a mixture model. The interaction with other states of matter is treated in Refs. 9 and 113, where the previous hybrid models are coupled with the elastoplastic models of Refs. 112 and 107 for the representation of the solid phase and the two systems interact at the material boundaries via mixed-material Riemann solvers (see Sec. IV F). Magnetic effects and plasma modeling are also accounted for by adding a set of resistive magnetohydrodynamic (MHD) equations.^{98} In the work of Nikodemou *et al.,*^{102} with the attempt to allow structural response and heat conduction for all the materials in situations of detonation wave loading, the Godunov–Peshkov–Romenskii (GPR) model^{114,115} is integrated to the framework developed in Ref. 9. Being the GPR mode applied to each material of the two phases, a single formulation can be applied on the whole domain, able to model simultaneously condensed-phase explosives and material response.

*Interim considerations on DI models*: Some advantages and shortcomings of different DI models are discussed here, without taking into account the sharpening techniques. There is a large number of applications of the DI models on problems such as the shock–bubble interaction and spherical bubble dynamics capable of testing the properties of the models. The simplest five-equation models of Allaire *et al.*^{3} and Massoni *et al.*^{116} [Eq. (13)] are known to be very efficient^{37,53} and have been successfully applied with different discretization schemes and flux solvers for the shock–bubble interaction problem, with a degree of accuracy comparable with solutions obtained with non-equilibrium models. Unfortunately, the absence of a regularization source such as the one in Eq. (8), responsible for the compressibility of the mixture regions, causes this kind of models to be unsuitable for problems involving bubble collapse and dynamic interface appearance. The effects of Kapila's source term of Eq. (8) on the bubble shape during contraction are explored for instance in the work of Tiwari *et al.*^{38} and Schmidmayer *et al.*^{37} The former attributed to this source term a role in increasing the thermodynamic consistency of the model and the latter suggested how this term acts in favor of a surface sharpening, by rebalancing the pressures on the water side and on the gas side of the mixture. On the other hand, the non-conservative nature of this source term makes the model numerically unstable under the presence of strong shocks. The six-equation model of Saurel *et al.*^{34} [Eq. (14)], in this sense, was demonstrated to be an effective alternative, replacing the non-conservative term with a relaxation product and introducing an additional Energy equation. The potential of the Kapila's model [Eq. (8)] and the six-equation model [Eq. (14)] in dealing with interface dynamics inspired the work of Favrie *et al.*^{8,95} and Ndanou *et al.*^{94} that through the introduction of visco-plastic models for Maxwellian materials^{36,115} successfully simulated solid fragmentation and crack propagation. Nonetheless, additional physical phenomenon such as surface tension, phase transition, and reactive phenomenon can be easily added. The seven-equation model, although the most complete, involves a complex wave pattern, which can consist of up to seven different wave speeds,^{117} and thus requires an intricate Riemann solver in addition to the non-conservative source terms, which poses similar considerations to the ones made for the Kapila's model. Because of this, most of the effort in the literature seems to be aimed at enhancing the reduced models rather than solving the complete set of equations.

## III. NUMERICAL FRAMEWORKS

### A. Finite volume formulation

Before discussing the implementation of the DIM, the principal numerical frameworks will be briefly explained, pointing the discussion to high-order reconstruction methods on unstructured grids. This subject has seen great advancements in the past years, and it seems natural that the high order FV and DG frameworks will be commonly used in conjunction with multi-species models. The LBM will also be discussed, focusing on the compressible formulations, that is, more pertinent to the present work.

The FV framework can be considered as the most popular numerical framework for resolving conservative laws, and nowadays, the second-order FV schemes are the standard for commercial codes.^{79} The scheme relies on the integral form of the Navier–Stokes or Euler equations, reinforcing the flux conservation even for complex grid topologies and in the presence of discontinuities.

We start the discussion of the FV schemes focusing on the discretization of the convective terms. Therefore, it is convenient to recast Eq. (1) in the following simplified form:

and we limit here the discussion to the Euler system of equations. The vector of the conserved variables *U* and the components of the flux tensor *F* therefore reads

The integral form of the governing equations, over a computational domain Ω with boundaries $\u2202\Omega $ and outward normal unit $n\u2192$, is

Assuming the domain to be a tessellation of small non-overlapping *N* elements *V _{i}* with volume $|Vi|$, where $i=1,...,N$, and applying the integral Euler formulation over a small control volume

*V*, the semi-discrete FV scheme reads

_{i}The state variable and face flux integrals have been replaced, respectively, by the volume averaged state variable $U\xafi$ and averaged normal flux $F\xaff$,

Hereafter, we will assume a domain Ω that can be any combination of hexahedral, pyramidal, or prismatic elements, with a surface area of $|f|$ for the face *f*, and with the number of faces *N _{f}*. In computational fluid dynamics (CFD) codes, the averaged flux is commonly numerically approximated through a Gaussian quadrature formula.

^{118}The accuracy of the approximation depends upon the selected number of quadrature points

*N*, denoted with $xf,\alpha $, and their distribution and weights $\omega \alpha $ depend upon the Gaussian rule. The approximate face integral is, therefore, of the form

_{qp}The conservative variables are represented through piecewise functions, which are discontinuous at cell interfaces. The upwind-Godunov methods resolve the discontinuity posing a Riemann problem taking as inputs the left and right extrapolated boundary values, named $UL(xf,\alpha )$ and $UR(xf,\alpha )$. The flux function, in terms of the Riemann solver $F\u0303$ therefore, reads

Substituting Eq. (40) in Eq. (39) and then, again, the results in Eq. (37), the high-order explicit finite-volume formulation can be written as

*High-order space discretization:* The ultimate goal in high-order FV schemes is to eliminate spurious oscillations near sharp gradients of the state variables, that is, in proximity of shocks or other sharp discontinuities, and at the same time provide high-order accuracy in smooth flows regions. Independently from the procedure to eliminate nonphysical oscillations, which in the FV frameworks comprehend a variety of methods, such as the total variation diminishing (TVD) techniques and the weighted essentially non-oscillatory schemes, the key ingredient is the reconstruction of a high-order polynomial $pi(x,y,z)$. The way the polynomial *p _{i}* is built is one of the major differences between the FV and the DG, whereas the limiting techniques are mostly shared between these two frameworks with minor modifications. The procedure for building a high-order polynomial in the FV setting for unstructured meshes will now be briefly discussed.

As mentioned, the ultimate goal is to build a high-order polynomial for each element *i* that on the same element equals the average of the state quantity, *u _{i}*,

Within the FV framework, in order to allow such reconstruction, the averages of the target cell *V _{i}* will be used as well as the averages

*u*of the neighboring cells constituting the computational stencil $S$.

_{m}The central reconstruction stencil is built up adding recursively the side neighbors of an element *V _{i}*, until a desired number

*M*of elements is reached starting from the target element for which

*m*= 0,

Since in this work we want to consider unstructured arbitrary meshes, it is convenient sometimes to perform the reconstruction in a uniform reference space.^{49,57} Considering for example a tetrahedral element with vertices $w1=(x1,y1,z1),\u2009w2=(x2,y2,z2),\u2009w3=(x3,y3,z3),w4=(x4,y4,z4)$, the transformation from the physical coordinate system (*x*, *y*, *z*) to the reference coordinate system $(\xi ,\eta ,\zeta )$ reads

with *J*, the Jacobian matrix, defined as

This transformation defines both the direct $x=x(\xi )$ and inverse $\xi =\xi (x)$ mappings from $\xi =(\xi ,\eta ,\zeta )$ into $x=(x,y,z)$. Applying the inverse mapping to all the elements *V _{m}* of the reconstruction stencil $S$, the transformed elements $V\u2032m$ will constitute the transformed stencil $S\u2032$,

Over this reference system, the reconstruction polynomial is sought as expansion over local polynomial basis functions $\lambda k(\xi ,\eta ,\zeta )$, leading to

defined *U*_{0} as the vector of cell averaged conserved variables at the target cell. The upper summation bound *K* corresponds to the number of degree of freedom and is related to the degree of the polynomial *r* by the expression

In order to determine the degrees of freedom *a _{k}*, one can impose that for each cell $V\u2032m$ of the stencil $S\u2032$, the cell average of the reconstruction polynomial to be equal to the cell average of the solution

*U*,

_{m}An expression for the DOFs can be drawn in matrix form

where *A _{mk}* represents the integral of the basis function over the cell

*m*in the stencil and

*b*the right hand side vector, respectively,

_{m}### B. Discontinuous Galerkin formulation

In DG methods, higher accuracy orders are reached by means of a high-order polynomial representation of the local element solution. Similar to FV schemes, as the solution is discontinuous across the domain elements, the treatment of the numerical flux may be dealt with by introducing an upwinding technique, and through the solution of the interfacial Riemann problem.

The DG method is based on a weak formulation of the hyperbolic system of conservation laws. This form is obtained by multiplying the form Eq. (34) by a smooth test function $\varphi (x,y,z)$, integrating over the domain Ω, and performing an integration by parts. This results in

The solution is again discretely approximated by a collection of piecewise solution, but this time, on each element those are defined as a combination of *n* local polynomial basis functions $P(x,y,z)$. The discrete solution lies in a finite-element space of discontinuous functions, that is, a Sobolev space^{119} $Vh={\varphi h\u2208L\u221e:\varphi h|\Omega \u2208Vk(\Omega ),k=0,1,2,\u2026,N}$, where *V ^{k}* is the space of polynomials of degree up to k. The discrete solution

*U*, with expansion coefficients denoted by

_{h}*u*, can be seen as expansions over a finite-element basis $Pjk$ in the aforementioned polynomial space, where

_{h}*d*is the number of degrees of freedom

It is convenient at this point to express the integrals of Eq. (52) in terms of an element *E* of the domain, that in this review will be allowed to be of an arbitrary shape. The relation between the element shape and the basis functions will be clarified later. The boundary integral will be included in the summation over the faces $|f|$ of the element *E*, and the semi-discrete form will, therefore, read

In a Galerkin discretization, the test function is taken from the same polynomial space of the solution expansion

Therefore, Eq. (54) can be rewritten, and since it must be satisfied for any element *E* and any function $\varphi h$, this will result in Ref. 118,

The volume integral on the right-hand side of Eq. (56) and the surface integral on the left hand side can be calculated by an appropriate Gaussian quadrature rule. The number of quadrature points depends on the degree of the shape function used to represent the solution. Usually, three, six, and twelve points are used for linear, quadratic, and cubic basis function, respectively, for volume integrals, and two, three, and four points for the boundary integrals. The boundary integral will be approximated with

and the volume integral with

where *L* and *M* are total number of surface and volume quadrature points, respectively. The intercell fluxes in Eq. (57) are determined in the same manner as in Eq. (40) with $F\u0303$ being the resulting approximate numerical flux. Inserting this relation in Eq. (56), together with the approximated volume integral of Eq. (58), the semi-discrete system can be further extended

The term on the left-hand side represents the stiffness matrix and is frequent in literature to find Eq. (59) in a very concise form as system of ordinary differential equations

In general, *U* is the vector of the degree of freedom, *M* denotes the mass matrix, and *R*(*U*) the residual vector.

*Space discretization:* Unlike FV schemes, DG schemes enable higher-order spatial discretizations by increasing the order of the polynomial representing the solution. Thus, there is no need of an enlarged stencil like in the FV case, as the stencil remains local. This leads to great advantages in applications concerning complex, unstructured meshes, as well as scalability and parallel computing efficiency.

As already discussed, the solution is expanded over a finite basis function $Pjk$ as in Eq. (53), and the number of degrees of freedom (DOFs) *d* depends upon the degree of the polynomial *k* and the space dimension *n*. For example, for triangles and tetrahedra, the number of DOFs can be calculated through

It can be easily proved with the above that, in three dimensions, the number of polynomials increases dramatically as the order of accuracy is increased (for a second-order reconstruction, ten polynomials are required, and for a third order, the requirement is 20) with a remarkable impact on the computational efficiency. Traditionally, the basis can be represented through Lagrange finite element or node-based basis. With such a choice, the DOFs to be solved are the node variables, at the vertices for a linear reconstruction, for example. Different choices are, however, possible, and in a modal formulation, the unknowns to be solved are the polynomial expansion coefficients. Following the latter, in order to introduce the reconstructed DG schemes about which we are going to talk soon, a popular choice is the cell-centered Taylor expansion.^{120–122} In this series of papers, the solution is expressed through a quadratic expansion, here limited for compactness to two dimensions (extension to three dimensions is straightforward) and is of the form

where the basis functions are normalized in order to improve the matrix conditioning. Thus, they read as

### C. Hybrid FV-DG formulations

The features of the FV and the DG schemes were discussed in Secs. III A and III B. The advantages of the DG formulation over the FV in higher-order applications are well known.^{79,118} However, it was also shown the extreme increase in computer operations when third and higher polynomials are used in three dimensions. Indeed, for the same order of accuracy, the DG methods consist of a system of equations with a higher number of unknowns compared to the FV discretization. An interesting approach, aiming to decrease the computational cost of the DG method, is the family of reconstructed DG (rDG) method, proposed by Dumbser *et al.* in Refs. 57 and 123. This procedure shares some similarities with the work of van Leer *et al.,*^{124,125} where a recovery method is built to compute the diffusive fluxes within the DG scheme. The recovery is to be considered as a *L*_{2} projection, which provides a polynomial solution, that is, indistinguishable in the weak sense from the underlying discontinuous solution. In the procedure proposed by Dumbser *et al.*,^{57} named $PNPM$ schemes, the solution is initially represented by a polynomial *P _{N}* of degree

*N*. The time evolution and the intercell flux evaluation are performed with a polynomial

*P*, reconstructed from the underlying polynomial

_{M}*P*. The reconstruction is performed by means of a least-square procedure, operating on a overdetermined system built up by setting a weak equality between the reconstructed solution on the target cell and the DG solution in the neighboring cells. The basis function of degree N in the reference space $(\Phi l)$ and the hierarchical orthogonal reconstruction basis $\Psi l$ are chosen in such a way to satisfy the following properties:

_{N}The first equality means that up to degree *N*, the reconstruction polynomial basis and the test functions basis are set to be equal, while for reconstruction degrees higher than *N*, the two bases are set to be orthogonal. The interesting property of such schemes is that they are a unified formulation for both DG and FV schemes. The FV scheme is obtained in the limit of *N* = 0, as the solution is represented by a piecewise constant polynomial, and the classic DG scheme is obtained in the limit *M* = *N*, meaning that the reconstruction is an identity operator. Thus, to perform a reconstruction, the condition has to be necessarily *M* > *N*. The original $PNPM$ uses Lagrangian type polynomials, and thus, the polynomial depends upon the shape of the arbitrary element. In a series of works,^{58,126} Luo *et al.* proposed to use a Taylor basis, that is, in the form of Eqs. (62) and (63), with a series of benefits in terms of computational costs and robustness. They claimed in fact that the recovery-based procedure of Dumbser *et al.*,^{57} although capable of *k*th-order accuracy, is costly due to the integration operation of the *L*_{2} projection, and may also degrade to lower order solutions near boundaries if enough neighboring cells are not available. As a solution, they proposed to set the reconstructed solutions and its first derivatives to be equal to the underlying DG solutions and their respective first derivatives on the target cell and its face neighbors. The least-square procedure for the third-order, so-called reconstructed DG $(P1P2)$, or $RDG(P1P2)$, for a two-dimensional problem will be briefly introduced (for the full presentation, we refer to the original works^{58,126}). Defining the derivative operation through the subscript notation, a linear polynomial solution *U _{i}* for a cell

*i*can be expressed in a compact form

Similar to Eq. (62), referring to the elements with superscript *R* as the reconstructed DOFs, the reconstructed quadratic polynomial reads

and imposing the equality of the solutions and their first-order derivatives

the final requirement is, considering a neighboring cell *j*,

The basis functions are evaluated at the cell center, and the derivatives terms have been expressed in a compacted form as follows:

The relations in Eq. (68) can be expressed in matrix form, isolating the second-order derivatives

This method is indeed more efficient than the original $PNPM$ scheme. However, since it is required that the reconstructed solution has to be equal to the underlying DG solution of the neighboring cell at the centroid, the procedure is not *k*th-order accurate as the solution does not follow the expected order of convergence. A possible solution is to set the neighboring cell average as the approximate value of the reconstructed solution, instead of the centroid value. In this sense, the scheme uses one equation from the recovery, that is, for reconstructed solution, and the two first derivatives of the reconstructed scheme. This workaround was named hybrid least squares procedure and can be found in the work of Cheng *et al.*^{127} Another hybrid FV-DG formulation was proposed by Zhang *et al.,*^{128} where the Taylor expansion is again used for the basis function, but the reconstruction for the higher order moments is performed by means of a Green–Gauss procedure, as usually done for FV schemes. Therefore, the scheme is computationally less expensive but also less accurate.

Further improvements can be found in literature. For example, Cheng and Liu^{129} proposed a reconstructed DG-FV scheme in which the higher order degrees of freedom are reconstructed by minimizing the interfacial jump integration, as first introduced by Wang *et al.*^{130} for the FV framework, aiming to increase the efficiency for steady-state problems.

The performances of some reconstructed DG schemes were compared with traditional DG schemes in Ref. 127 and illustrated in an adapted form in Fig. 1 for the subsonic flow past a cylinder test case. From a qualitative point of view, one can observe the substantial improvement in terms of computational time in place of an accuracy comparable with the one obtained with the traditional DG method.

### D. Limiting techniques

The benefit of higher-order methods for smooth regions is well known and documented. The necessity of orders higher than the second to capture the rich structures that may be evolving, in particular, flows, however, conflicts with the requirement of oscillation-free solutions near steep gradients. To circumvent Godunov's theorem, for which a linear, monotone, scheme can be at most first-order accurate, a number of different schemes have been proposed in literature. Within this review, we are going to limit the discussion to the family of total-variation-diminishing (TVD) and monotone upstream scheme for conservation laws (MUSCL),^{43,44} and to the family of weighted essentially non-oscillatory (WENO) schemes,^{46,47} as those are the families of limiting techniques that are so far used in conjunction with diffuse interface models. The MUSCL scheme and its extensions avoid the appearance of numerical oscillations by shifting to lower order approximations if solutions bounds have been violated, while the WENO scheme provides a solution taking the data from the smoothest regions of the computational stencil. The former is robust and works on small stencils, although it is suffering from accuracy degradation and is parameter dependent, while the latter preserves the accuracy and is less parameter sensitive, but works on large stencils, that may adversely affect the computational efficiency and robustness. Based on the work of Harten,^{131} the total variation is an indicator of the oscillations in a solution and is defined as

therefore, in order for a method to be TVD, the following condition must be satisfied:

To expose the MUSCL framework, we are going to recall Eq. (47) and replacing the polynomial with the extrapolated reconstructed solution $Uij,\alpha $ of element *i* and face *j* at quadrature point $xf,\alpha $, and multiplying all the polynomials for the quadrature points by a limiting function Θ_{i},

The limiter is designed in such a way to prevent spurious oscillations for every cell *i*, and a popular one, used for unstructured meshes was developed by Barth and Jespersen,^{44} which is limited to second-order accuracy, although other variants for higher-order accuracies have been lately developed.^{132} Although the Barth and Jespersen limiter is non-differentiable, restricting the convergence properties of the resulting scheme, this limiter is one of the building blocks for other limiters developed for single and multi-species flows.

The slope limiter $\Theta il,\alpha $ is commonly expressed as

The maximum $Uimax$ and the minimum $Uimin$ values of the state variable over the stencil are required, and the limiter Θ_{i} selects the minimum value from the slope limiter $\Theta il,\alpha $ from all the side neighbors, and quadrature points

On the other hand, for a WENO reconstruction on a three-dimensional unstructured grid composed of mixed elements, a combination of central and sectorial stencils is required. The extension of WENO schemes to arbitrary elements, such hexahedrons, prisms, tetrahedrons, and pyramids, is discussed in the work of Tsoutsanis *et al.*^{49} The central stencil, built as in Eq. (43) for a physical coordinate system and in Eq. (46) for a reference coordinate system, a number of sectorial stencil is formed by adding the neighboring cells with the center lying inside a given sector, and the reader is referred to the work of Tsoutsanis^{51} for a detailed analysis of several stencil selection algorithms. Once the stencils are formed, the WENO polynomial *p _{weno}* is defined as a non-linear combination of the single reconstruction polynomial per each stencil $pm(\xi ,\eta ,\zeta )$,

given the total number of stencils *m _{s}*. Substituting the expression for the single polynomial Eq. (47) in Eq. (76), one obtains the expanded expression for the WENO polynomial

The linear weights *d _{m}* take a large value for the central stencil ($102<d0<105$), although in Ref. 134 it was found that higher values are better for smooth solutions and lower values works better on discontinuities, and smaller linear weights are assigned to the sectorial stencil (

*d*= 1). The parameters

_{m}*p*and

*ε*are responsible to control the non-linear weights decay for non-smooth stencils and to avoid division by zero, respectively, and typically, the values

*p*= 4 and $\epsilon =10\u22126$ are used. The smoothness indicator

*IS*, as the name suggests, indicates how the solution is smooth on the stencil, and is calculated with

_{m}where *D* is the derivative operator, *β* is a multiple-dimension index, and *r* is the order of the polynomial. As pointed in Ref. 134, Eq. (79) just depends on the reconstruction basis and can only be computed once, at the beginning of the calculations for each element *V _{i}*.

In order to decrease the computational footprint of high-order WENO schemes that involves large stencils, different strategies have been developed. For instance, one of the ideas is to decrease the polynomial order of the directional stencils only, as done by the so-called compact WENO (CWENO) schemes introduced by Levy *et al.*^{135} for 1D and 2D Cartesian meshes and extended by Dumbser *et al.*^{136} to unstructured simplex 3D meshes. The CWENO scheme brings robustness because the reduced stencil dimension increases the chance to find at least one smooth solution, but also requires a balancing of the parameters to blend together polynomials of different orders. When it comes to the DG method, the coupling with WENO schemes seems even more convenient because lower order, but sufficiently accurate polynomials, are evolved maintaining stencil compactness of the DG method, with the possibility to reconstruct higher-order polynomials. This can be seen as a first step toward the reconstructed hybrid DG-FV schemes. Non-linear limiting for DG schemes was introduced by Qiu and Shu,^{137} and they split the problem in two parts: the first consists of finding the cells containing the discontinuity with a TVB limiter, and in the second part, the higher moments are provided by the reconstruction, either with a classic WENO reconstruction, or with a compact Hermite WENO (HWENO) reconstruction operating on Hermite polynomials.^{62} Balsara *et al.*^{138} replaced the total variation bounded (TVB) limiter, detecting the troubled cells with a reformulated monotonicity preserving (MP) limiter, thus preserving meaningful cell substructures and providing a problem-independent model. The HWENO reconstruction was then proposed for the rDG framework by Luo *et al.,*^{122} in an attempt to eliminate the oscillations during reconstruction phase. However, in order to avoid also the oscillations carried by the underlying DG discretization, a Hierarchical Hermite reconstruction procedure is presented in Ref. 61.

### E. Lattice Boltzmann method

Between the widely known numerical frameworks, the LBM has seen an increasing number of applications thanks to the efforts by the research community in recent years. The LBM is recognized as a very efficient candidate for the simulation of incompressible and multi-species flows. In particular, thanks to the kinetic formulation of the model, the pressure is calculated from the equation of state and one is released from the inconvenience of solving the Poisson equation with clear benefits in terms of efficiency. In addition, the model is able to accurately treat flow fields in a wide range of Knudsen numbers^{139,140} and the linearity of the convective term in the phase space allows easier implementations.^{67} However, the extension to high Mach number flows is not straightforward^{141,142} and the LBM in this sense still have not reached the complete maturity. A review of the lattice Boltzmann framework for incompressible single-phase and multiphase flows can be found in Ref. 66, and the inconsistencies for thermal-acoustic properties of the discretized Boltzmann approach model are discussed in Refs. 142–144. They show that a large number of discrete speeds is necessary to recover a consistent continuum description,^{141} or alternatively, this can be achieved by adding counterterms to compensate for the lattice constraints.

The lattice Boltzmann equation can be obtained in multiple ways, either starting from the discrete kinetic equation for a particle distribution, or from the continuum Boltzmann equation.^{145} Following the latter, one can describe the time and space evolution of the velocity distribution *f* through

The variable *c* denotes the particle velocity, and Ω is the collision operator, and as originally proposed by Bathnagar, Gross, and Krook,^{146} if *f* is close to its equilibrium state *f ^{eq}*, it can be approximated by a relaxation mechanism, ruled by a relaxation time

*τ*dependent on fluid properties like the viscosity

The resulting discrete velocity Boltzmann equation (DVBE) is a set of PDEs, with a number of discrete velocities *c _{i}*,

The LBM solves this set of equations through a collide and stream process, and the explicit system reads

The target is to recover the physical properties in terms of moments of the distribution function. The zeroth-order mass and the first-order momentum read, respectively,

As pointed by Coreixas and Latt,^{147} this kind of discretization is based on static velocity stencils in time and space, and thus is fixed with respect to a population density. Moreover, the equilibrium distribution functions are usually based on the expansion over the Maxwellian status with reference state ($0,T0$). A consequence of this is the tendency to trigger oscillations if strong disequilibrium or great temperature fluctuations characterize the macroscopic quantities. For incompressible flows, this is not of great interest, and in fact, standard LBMs, based for example on 19 or 27 discrete velocities, are a much more efficient alternative to the NS equations. However, it is known that this poses a great limitation to the range of heat capacity ratio and Prandtl numbers addressable by the framework.^{143}

As anticipated previously, a straightforward solution is to enlarge the number of discrete velocities, as proposed in the multi-speed method initiated by Alexander *et al.*^{148} This method itself is isothermal and cannot be properly used as a compressible solver but within this fashion, some off-lattice methods in conjunction with an entropic formulation of the collision operator were built, resulting in a model capable of ranging in the weak compressibility.^{149} The major drawback of these methods is the increase in computational requirements in order to perform the interpolations of the non-integer velocity, which after the stream-and-collide process are not associated with any grid node, canceling the benefits of adopting such framework. This inconvenience was addressed in Ref. 147. As an alternative, He *et al.*^{150} proposed a double-distribution function, which was further simplified and elaborated in Ref. 151. In this context, an additional internal distribution function was derived, obtaining a model with a distribution for the flow field and one for the temperature, coupled by the equation of state. The method can be used for arbitrary Prandtl numbers, but the complicated recovery of the macroscopic variables limits the popularity of this method. More recently, Chen *et al.*^{152} presented an arbitrary *γ* and Prandtl model using multiple-relaxation-time schemes, where the transport process can be controlled separately, by different collision parameters associated with the transport variables.

### F. Hybrid lattice Boltzmann-finite volume schemes

Contemporary to Sec. III E, an interesting class of hybrid lattice Boltzmann–finite volume (LBFV) schemes emerged and the hybridization between the two frameworks, in the context of compressible flows, followed mainly two different paths. One type of hybridization uses a FV scheme for the temporal and spatial discretization of the lattice Boltzmann equations,^{75,153,154} and the other type uses the FV scheme for storing the node parameter data and the LBM for evolving the intercell face fluxes.^{77,155} It has to be said that the philosophies that led to these are quite different. In the first, the main goal is to reduce the dependency of the discrete velocities from the lattices, and thus to overcome the difficulties of the LBM for the compressible flow range, while the second aims to provide a robust solver for multi-species flows.

For instance, the multidimensional LBM of Feng *et al.*^{153} uses a double distribution function of the type of He *et al.*,^{150}

and the quantities $h\u2032=c\xb7u\u221212(f\u2212feq)$ and $1\tau hf=1\tau h\u22121\tau f$ are defined from Guo's model.^{156} If one considers, the following density *f ^{eq}* and energy

*h*equilibrium distributions, respectively:

^{eq}with *D* the spatial dimension, *R* the gas constant, *c* the molecule velocity, and *ρ* and *u* the fluid density and velocity, respectively, then the thermal compressible Navier–Stokes can be retrieved by a Chapman–Enskog expansion.^{156,157} As in the work of Shan and Chen,^{70} the distributions are expanded through Grad's moment in terms of Hermitian polynomials

The expansion coefficients $a(n)$ and $b(n)$ can be calculated after a projection of the distribution on the Hermitian basis. $H(n)(c)$ is n-rank tensor and the weighting function $\omega (c,T)$ is defined as

The integration from time *t _{n}* to time $tn+1$ leads to the so-called discrete unified gas-kinetic scheme (DUGKS).

^{158}The feature is that like FV schemes, this can be easily applied to unstructured meshes

with the terms *F* and *H* resembling the fluxes form of the FV schemes, expressed as

and the reconstructed interface values of gas density and energy distributions can be calculated with one of the techniques illustrated in Secs. III A–III F.

On the other hand, Joshi *et al.*^{77} used a LBM step just for the calculation of the intercell fluxes. In this work, the LBM considered is the one introduced by Kataoka and Tsutahara (KT)^{159} and that in two dimensions, it uses nine discrete velocities D2Q9, defined as

with $v1\u2260v2\u22600$ non-zero constants, and the local equilibrium particle velocity distribution function is given by

with $k\u2208[1,9]$. We refer to Refs. 159 and 77 for the determination of the coefficients *A _{k}*,

*B*and

_{k,}*C*. The variables are stored in a finite volume fashion and the time evolution is performed with a classic FV time marching scheme. The initial time step velocity distribution is obtained from the equilibrium distribution, and once the fluxes are evaluated through a decoupled procedure for the

_{k}*x*and

*y*directions using an explicit scheme of the type of Eq. (83), the macroscopic variables are retrieved through Eq. (84). Although the procedure in this form is not suitable for very high Mach number flows due to the adoption of the KT model, the benefit of such hybrid procedure consisted in the application of the LBM on unstructured meshes, and it was also claimed in Ref. 77 that the flux evaluation resulted in more accurate results compared with the classical FVM-Godunov and Flux Vector Splitting techniques, although only a qualitative study was presented.

## IV. HIGH-ORDER IMPLEMENTATION OF DIM

Synthesizing the requirements for an interface-capturing scheme to be successful in simulating compressible multi-species flows, one can say that the scheme has to satisfy the following requirements:

Discrete conservation of mass, momentum and energy

High-order accuracy solution of smooth flow regions

Oscillation-free representation of interfaces and discontinuities.

Given these requirements, it is straightforward to bring the FV framework with high-order reconstruction, and coupled to one of the techniques presented in Sec. III D as a good candidate for such purpose. Indeed, the discrete conservation is easily enforced casting the governing equations in a conservative form, typical of FV schemes. On the other hand, the advection equation for the material interface must be written in non-conservative form in order to avoid the appearance of oscillations near interfaces.^{30} One of the first applications of interface capturing schemes with high-order reconstruction on a Cartesian mesh is by Johnsen and Colonius,^{4} who used the modified version of Eq. (12) with volume fraction advection, introduced by Shyue,^{29} with up to a fifth order WENO reconstruction.

In order to avoid the appearance of oscillations, the reconstruction is performed on primitive variables projected in the characteristic space, and projected back to the physical space after reconstruction. At the beginning of each time step, the average primitive variables are obtained directly from second-order average conservative variables, deteriorating the solution with a second-order spatial error. This was discussed by Coralic and Colonius,^{53} and they proposed to use a two-point Gaussian quadrature to recover fourth-order estimate of average primitive variables that are then used for the reconstruction. This allowed for a recover of a formal fifth-order convergence property. In addition, they adopted the five-equation model of Allaire, given by Eq. (13) and included viscous effects.

In both Refs. 4 and 53, an Harten-Lax-van Leer contact (HLLC) approximate Riemann solver^{160} with an adaptation to solve advection equations was used, and in Ref. 161 the transport equation for the volume fraction equation took the following, equivalent alternative form:

where the velocity is defined as $u={u,v}$. The modified HLLC solver^{160} for a two-dimensional problem reads

where the starred values are defined as (K = L or K = R)

The wave speed in the star region is calculated by

and defining $u\xaf$ and $c\xaf$ as the average between left and right states of the Riemann problem for velocity and sound speed, the wave speeds can be calculated by

with

The velocity in the numerical source terms is taken in the same form as the advective flux

The HLLC solver is preferred to the exact Riemann solver because it is computationally more efficient, and since it preserves positivity, it is also preferred to the less dissipative Roe scheme. This last feature is of particular importance in multi-species flows where both pressure and density can reach very low values. However, in the case of strong shocks, the HLLC is known to suffer from the carbuncle phenomenon as spurious oscillations appear in the solution.^{162} Deng *et al.*^{163} proposed an alternative to the HLLC solver, by restoring the missing contact wave in the HLL solver through the introduction of a jump like in the THINC reconstruction (see Sec. IV A) and minimizing the dissipation with a boundary variation diminishing algorithm. Since the wave structure of the solver is not modified, this method retains the robustness of the original HLL solver in the case of strong shocks.

A recent investigation on the performances of high-resolution methods for shock–turbulence interaction was conducted by Johnsen *et al.*^{161} Although capable of sharp shock capturing, the well-known WENO scheme of Ref. 46 was reported to dissipate turbulent fluctuations for the higher resolvable wavenumbers and results in a degradation of spectral properties as already pointed out by Pirozzoli^{167} because of the upwind biased nonlinear reconstruction. In the Finite Difference context, several workarounds have been proposed, for example in form of improved nonlinear weights evaluation, by integrating the smoothness indicators Eq. (79) on a full stencil and using much smaller values for the parameter *ε*, as in the WENO-Z^{168} and in the WENO-M^{169} or, alternatively, by adding another stencil to the upwind-biased setting of Ref. 46 and minimizing the truncation error as done by Martin *et al.*^{165} for the bandwidth-optimized WENO scheme. A more localized dissipation was achieved later using the compact finite difference schemes. The compact FD scheme was introduced by Lele^{170} and features enhanced small-scale waves resolution and spectral properties, and was adapted to flows with discontinuities by Deng and Zhang^{171} combining it with the WENO limiting technique, leading to the so-called weighted compact nonlinear schemes (WCNS). An application of this scheme within diffuse interface models is documented in Ref. 164 and uses a hybrid interpolation technique, where the WENO-Z is used in regions with discontinuities, and a central interpolation is used in smooth regions to keep the dissipation as low as possible.

Similarly, Sun *et al.*^{172} proposed a semi-discrete FD scheme known as minimized dispersion and controllable dissipation (MDCD) scheme, where the linear weights of the linear reconstruction are controlled independently by free parameters for dispersion and dissipation, the first chosen to optimize dispersion properties and the second chosen arbitrary to provide the dissipation necessary for a specific test case. This linear reconstruction can be used inside a WENO reconstruction to provide shock capturing capability, and in Ref. 166, the resulting hybrid scheme is used for regions containing interface discontinuities and shock waves, and is coupled with a pure linear reconstruction to retrieve the most accurate solution in smooth regions following the procedure discussed by Ren *et al.*^{173}

The performances of the MDCD scheme for multi-species were compared with the MUSCL scheme and with some of the recent variations of WENO schemes in Ref. 166. From Table I, adapted from the previous reference, one can observe the convergence order of different schemes including MUSCL, WENO, and the WENO-SYMBO^{165} and MDCD-WENO^{166} variants for the multi-species advection test defined in Ref. 164. In Ref. 166, it is claimed that the MDCD-WENO is more efficient that the classical WENO-JS^{46} if $10\u221220$ cells per wavelength are used to compute the highest wavenumber component.

. | MUSCL . | WENO-JS . | WENO-SYMBO . | MDCD-WENO . | ||||
---|---|---|---|---|---|---|---|---|

Number of cells . | Error . | Order . | Error . | Order . | Error . | Order . | Error . | Order . |

4 | 7.587 × 10^{−02} | 7.012 × 10^{−02} | 6.332 × 10^{−02} | 3.688 × 10^{−02} | ||||

8 | 6.135 × 10^{−02} | 0.307 | 1.174 × 10^{−02} | 2.579 | 3.489 × 10^{−03} | 4.182 | 1.743 × 10^{−03} | 4.403 |

16 | 2.657 × 10^{−02} | 1.207 | 7.089 × 10^{−04} | 4.049 | 4.218 × 10^{−04} | 3.048 | 1.614 × 10^{−04} | 3.432 |

32 | 1.074 × 10^{−02} | 1.307 | 2.644 × 10^{−05} | 4.745 | 1.905 × 10^{−05} | 4.468 | 1.105 × 10^{−05} | 3.869 |

64 | 3.658 × 10^{−03} | 1.553 | 8.631 × 10^{−07} | 4.937 | 1.121 × 10^{−06} | 4.087 | 6.915 × 10^{−07} | 3.998 |

128 | 1.153 × 10^{−03} | 1.666 | 3.223 × 10^{−08} | 4.743 | 5.541 × 10^{−08} | 4.339 | 2.814 × 10^{−08} | 4.619 |

. | MUSCL . | WENO-JS . | WENO-SYMBO . | MDCD-WENO . | ||||
---|---|---|---|---|---|---|---|---|

Number of cells . | Error . | Order . | Error . | Order . | Error . | Order . | Error . | Order . |

4 | 7.587 × 10^{−02} | 7.012 × 10^{−02} | 6.332 × 10^{−02} | 3.688 × 10^{−02} | ||||

8 | 6.135 × 10^{−02} | 0.307 | 1.174 × 10^{−02} | 2.579 | 3.489 × 10^{−03} | 4.182 | 1.743 × 10^{−03} | 4.403 |

16 | 2.657 × 10^{−02} | 1.207 | 7.089 × 10^{−04} | 4.049 | 4.218 × 10^{−04} | 3.048 | 1.614 × 10^{−04} | 3.432 |

32 | 1.074 × 10^{−02} | 1.307 | 2.644 × 10^{−05} | 4.745 | 1.905 × 10^{−05} | 4.468 | 1.105 × 10^{−05} | 3.869 |

64 | 3.658 × 10^{−03} | 1.553 | 8.631 × 10^{−07} | 4.937 | 1.121 × 10^{−06} | 4.087 | 6.915 × 10^{−07} | 3.998 |

128 | 1.153 × 10^{−03} | 1.666 | 3.223 × 10^{−08} | 4.743 | 5.541 × 10^{−08} | 4.339 | 2.814 × 10^{−08} | 4.619 |

### A. Sharpening techniques

The strategies developed for capturing shock waves have been transposed with a degree of success for the problem of interface capturing. However, contact discontinuities pose additional challenges with respect to shock waves. In the case of shock waves, the convergence of characteristics acts as a counteraction to the numerical diffusion, differently from contact surface, whose parallel characteristics do not aid in such sense. Therefore, it is unavoidable that the function representing the interface condition artificially smears, unless some corrections are applied to the governing equations. Shukla *et al.*^{39} were one of the first groups to propose the addition of terms forcing advection toward the interface. In this work, in order to reduce the interface thickness, a combination of volume fraction and density corrections are proposed, similar in concept to the immiscibility regularization terms introduced by Olsson and Kreiss^{174} in the context of level set methods. The advection for the interface function is therefore rewritten as

once defined *n* as the interface normal $n=\u2207\varphi /|\u2207\varphi |,\u2009U$ as a characteristic compression velocity and *ε _{h}* as a grid spacing length scale. Since we are interested in sharpening the interface, small values of

*ε*and large values of $U0$ are the target. Under the consideration that also the density gradients will diffuse in time, a similar correction is required for the density equation, with the only difference that the density correction has to be insensitive to the density values close to the interface, since great density ratio may be involved:

_{h}The function $H(\varphi )$ is introduced to localize the compression region and may be taken as^{39}

These regularization terms are then coupled to the five-equation model [Eq. (13)], and the solver is presented in conjunction with a WENO limiter and a HLLC flux solver although the selection of the discretization is arbitrary. Unfortunately, the source terms turn out to be inconsistent with respect to the thermodynamic models for the mixture regions and can eventually lead to the accumulation of error even far from the interface.

In the same year, Kokh and Lagoutière^{175} proposed an anti-diffusive numerical scheme based on a Lagrange-remap algorithm. In this procedure, a five-equation model is first discretized with a Lagrangian coordinate system, and evolved in time through a Lagrangian step. During the Lagrangian step, the Lagrangian coordinate system deforms with the flow field, and in order to recover the solution in terms of Eulerian variables, the solution is remapped via limited downwind fluxes. Although conceptually simple and effective in limiting the interface diffusion, the implementation of such scheme can be intricate and not straightforward on already existing codes.

Another anti-diffusion interface sharpening approach was proposed by So *et al.,*^{176} inspired in this case by the techniques used in the VOF context. The sharpening effects are introduced as a post-processing operation, without the necessity of an interface reconstruction typical of VOF methods, through the solution of an anti-diffusion equation

where D is the anti-diffusion coefficient and *τ* is a pseudo time. The estimation of D is based on the underlying diffusivity of the numerical scheme, and as the HLL solver has been considered in this work, the resulting anti-diffusion coefficient is calculated by

The anti-diffusion equation is solved each time step with an explicit Euler scheme providing a *sharpened* value of $\varphi $. The remaining flow variables are updated according to the flux obtained for the new value of $\varphi $, and thus, the correction is applied ultimately on all flow variables, resulting in a conservative procedure. The interface profile can be corrected multiple times each step by solving recursively the diffusion equation. In this case, a stopping criteria is necessary, and based on an interface tolerance argument, the iteration stops once exceeded the following threshold:

This compares the gradient of the limited volume fraction $(\u2207\varphi \xaf)lim$ with the gradient calculated with a central difference scheme $(\u2207\alpha \xaf)CD$, and the tolerance has to be set sufficiently high, to not compromise numerical stability (TOL = 1 is used in Ref. 176).

Tiwari *et al.*^{38} introduced a new sharpening method, derived from the asymptotic reduction of a BN model corrected with a regularization operator. The result is a modified five-equation Kapila model, with a regularization operator $R(\varphi 2)$ similar to the one introduced in Ref. 39 and discussed in Eq. (100), but made more robust through the consideration that the variation of the mixture density *ρ* across the interface is slow. In practice, to the Kapila's system [Eq. (8)], another source S(U) is added that reads

The regularization operator, in this case, reads

with

The other regularization terms are defined as

and $R\u0302=R\u03021+R\u03022$. In Ref. 38, the effects of these sharpening terms are compared against the ones introduced in Refs. 39 and 176, and against the plain five-equation models of Allaire *et al.*^{3} [Eq. (13)] and Kapila *et al.*^{2} [Eq. (8)] for the isolated spherical bubble collapse problem. They measured that all the sharpening techniques consistently improves the bubble surface during collapse, although with the model introduced in Ref. 176, since there is not any parameter dependent on the grid scale, the thickness is not constant and evolves depending on the flow dynamics. It was also estimated that the sharpening techniques in Refs. 38 and 39 helps to prevent a diffusion by a factor of two respect to the baseline five-equation models, with which, in a three-dimensional case, an eight times finer mesh is required in order to have the same interface thickness.

Another interesting outcome of Ref. 38 is the demonstration of the superiority of high-order schemes (fifth-order WENO in this case) in representing the bubble radius on regular Cartesian meshes, as the directionality introduced by the mesh has a greater impact on lower order methods (second-order MUSCL with minmod limiter). A comparison of the resulting bubble surface is illustrated in Fig. 2.

An alternative kind of sharp interface recovery was proposed by Shyue and Xiao in Ref. 41. Here, a variant of the so-called tangent of hyperbola for interface-capturing (THINC) approach replaces the high-order reconstruction in proximity of the mixing zone. A function $\Phi $ is used as an indicator of the presence of an interface within a cell *C _{i}*, with an associated cell average $\Phi \xafi$, and which takes unity or zero value whether

*C*is filled with a fluid

_{i}**a**or not. For the intermediate values, the hyperbolic tangent function $\Phi \u0303i(x)$ mimics the jump across an interface

given the sign function $\sigma i=sgn(\Delta \Phi \xafi)$ of the jump $\Delta \Phi \xafi=\Phi i+1(t)\u2212\Phi i\u22121(t)$ and a user-defined slope parameter *β*. Assuming conservation of volume fraction in the form

the interface location $x\u0303i$ can be determined from Eq. (110) as it is the only unknown,

This allows the determination of the function $\Phi i$ at any distance from the interface, and coupled to the five-equation model of Allaire and a wave propagation algorithm for the time integration, it was demonstrated the effectiveness of the interface sharpening at a minimal cost in terms of implementation. The only additional step is placed after a standard MUSCL/WENO reconstruction, and the procedure can be potentially extended to arbitrary time marching schemes. However, Schmidmayer *et al.*^{37} tested the THINC scheme for the spherical bubble collapse and experienced some inconsistencies during error evaluation. Although the sharpening effect is still present, they quantified that when using different time integration schemes from wave propagation algorithm, a 44% larger error is introduced in the solution for low-pressure-ratio cases. In addition, this most likely is the reason for the less spherical shape of the bubble for the same low-pressure ratio cases. As a general rule, they found the THINC method to be better behaved when coupled to a lower order MUSCL scheme (see Fig. 3 where the irregular bubble shape at final times for lower pressure ratios is more pronounced for higher order schemes with THINC reconstruction).

Higher-order reconstructions are capable of obtaining more accurate solutions for the bubbles in initial pressure equilibrium, but in Ref. 37, it is suggested to smear the bubble interface over a few cells on the radial directions when using fifth- or higher-order reconstruction to avoid numerical instabilities. This has a major consequence when simulating bubble in initial pressure disequilibrium, and is of increasing importance as the initial pressure ratio is greater because the smeared mixing zone affects the initial wave dynamics and precludes accurate solutions. A sharpening technique is therefore necessary, although as before, for the case of the THINC approach a coupling with higher than third-order WENO reconstructions is generally avoided, to not incur into the opposite problem of sharpening the interface too much. This poses the necessity for a well-balanced combination of high-order discretization and sharpening methods.

### B. Applications on unstructured meshes

In the last few years, the diffuse interface models have been extended to unstructured meshes. This comes from the necessity of applying such models on complex and realistic engineering problems, which usually translates in geometries that are described more efficiently though unstructured grids, especially when considering three-dimensional cases. Chiapolino *et al.* introduced in Ref. 42 a novel flux limiter in the MUSCL framework for unstructured grids coupled to the pressure non-equilibrium six-equation model of Saurel *et al.,*^{34} and is presented as a new class of sharpening techniques. We recall here the MUSCL limiting procedure illustrated in Sec. III D, rewriting the slope limiter as in Ref. 42,

which allows to recover the minmod and superbee limiters from

by setting *β* = 1 and *β* = 2, respectively. It has to be noted that the reconstruction is here performed on the primitive variables $Wil,\alpha $. The newly introduced limiter, named Overbee,^{42} reads

where the term **t** is defined as

For instance, when *β* = 1 this corresponds to the Superbee limiter, while for *β* = 2, the boundary is extended to the upper first-order TVD region. Therefore, this new limiter combines the second-order TVD region of the Superbee limiter and the upper boundary of the first-order TVD region. The first-order region is of interest because it provides a compressive flux limiter for the interface discontinuity. This reconstruction is performed only during the hyperbolic step, which in Ref. 42 consist of a MUSCL-Hancock scheme, and only in proximity of the interface to not compress continuous waves in smooth regions. Taking *β* = 2, they demonstrated the capability of the method to keep the resolution of the interface discontinuity within 3 ± 1 mesh points.

A different spatial reconstruction method for unstructured grids is presented in Ref. 177 and is called MUSCL-THINC/QQ-BVD. Basically, the algorithm proposes two or more candidate reconstruction functions, for example, one linear polynomial built with a MUSCL scheme coupled to the multi-dimensional limiting process (MLP) slope limiter, and one non-polynomial reconstruction obtained with the THINC approach with quadratic surface representation and Gaussian quadrature (THINC/QQ). A boundary variation diminishing (BVD) algorithm selects the most suitable reconstruction between the two based on numerical dissipation argument. Assuming the candidate reconstruction function to be denoted with $Qi\xi n$, the first step is to calculate the boundary variation (BV) at a cell edge Γ_{ij} through

using the integrated values of the reconstruction function on the cell Ω_{i} and its neighbor Ω_{ij},

Then, the TBV of a reconstruction candidate $Qi\xi (x,y)$ is obtained summing over every edge of the cell

and the final reconstruction function is selected as the one with smallest TBV in the set of reconstruction candidates Ξ, which can be composed of an arbitrary number of reconstruction candidates *N*,

The strength of the model resides in the arbitrariness of the type and number of reconstructions that can be used to build the candidates set Ξ. For instance, in Ref. 177, a set of two and three reconstructions are compared. The first set includes one candidate built with MUSCL reconstruction and one with the THINC approach. The second sets adds to the previous another candidate, built with a smaller slope parameter *β*, which should improve the resolution of vortical structures, and compared to the previous set is proven to provide better prediction of weak discontinuities. In addition, although composed of only second-order reconstruction candidates, the results showed that such schemes using the BVD methods produced less error than higher-order schemes like the third-order WENO. On Cartesian meshes, Li *et al.*^{179} proposed an improved BVD reconstruction method employing a fifth-order modified WENO (MWENO) scheme, that when at least two adjacent target substencils are smooth, is capable of restoring the highest possible order of interpolation, improving the spectral properties of the scheme. Tsoutsanis *et al.*^{178} presented the compact WENO for multi-species flows on unstructured grids. As anticipated in Sec. III D, the CWENO uses directional stencils of reduced order with the purpose to make the scheme more efficient and robust. An optimal polynomial that uses the central stencil is defined as

summing over the stencils $1<s\u2264St$ with linear coefficients *λ _{s}*. The polynomial

*p*

_{1}is computed subtracting the lower-order polynomials from the optimal polynomial

and the reconstruction polynomial is a non-linear combination of the polynomials *p _{s}* as defined in Eq. (76). The central stencil

*λ*

_{1}is calculated by normalizing an arbitrary value assigned to the non-normalized weight $\lambda \u2032$, and subsequently, the lower-order polynomial weights are obtained from the central one:

The multi-species convergence study defined in the work of Wong and Lele^{164} is performed to prove the improvements of efficiency of the CWENO scheme over the classic WENO scheme on mixed element unstructured meshes (see Fig. 4), and an approximately 80% saving in computational time is found for the higher accuracy order (see Table II).

. | CWENO . | WENO . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Order/edges per side . | $eL\u221e$ . | $OL\u221e$ . | $eL2$ . | $OL\u221e$ . | Time
. | $eL\u221e$ . | $OL\u221e$ . | $eL2$ . | $OL\u221e$ . | Time
. |

$3rd/10$ | 1.124 × 10^{−02} | ⋯ | 6.897 × 10^{−03} | ⋯ | 1.000 | 1.119 × 10^{−02} | ⋯ | 6.880 × 10^{−03} | ⋯ | 1.252 |

$3rd/20$ | 1.514 × 10^{−03} | 2.89 | 9.535 × 10^{−04} | 2.85 | 8.240 | 1.506 × 10^{−03} | 2.89 | 9.492 × 10^{−04} | 2.86 | 8.757 |

$3rd/40$ | 1.993 × 10^{−04} | 2.93 | 1.244 × 10^{−04} | 2.94 | 85.831 | 1.982 × 10^{−04} | 2.93 | 1.236 × 10^{−04} | 2.94 | 94.658 |

$3rd/80$ | 2.487 × 10^{−05} | 3.00 | 1.581 × 10^{−05} | 2.98 | 673.159 | 2.473 × 10^{−05} | 3.00 | 1.571 × 10^{−05} | 2.98 | 798.829 |

$4th/10$ | 1.879 × 10^{−03} | ⋯ | 5.355 × 10^{−04} | ⋯ | 1.530 | 1.867 × 10^{−03} | ⋯ | 5.342 × 10^{−04} | ⋯ | 1.821 |

$4th/20$ | 7.487 × 10^{−05} | 4.65 | 2.762 × 10^{−05} | 4.28 | 10.864 | 7.427 × 10^{−05} | 4.65 | 2.739 × 10^{−05} | 4.29 | 14.081 |

$4th/40$ | 4.962 × 10^{−06} | 3.92 | 1.857 × 10^{−06} | 3.89 | 115.068 | 4.912 × 10^{−06} | 3.92 | 1.839 × 10^{−06} | 3.90 | 144.007 |

$4th/80$ | 3.482 × 10^{−07} | 3.83 | 1.200 × 10^{−07} | 3.95 | 928.223 | 3.451 × 10^{−07} | 3.83 | 1.188 × 10^{−07} | 3.95 | 1128.293 |

$5th/10$ | 9.115 × 10^{−04} | ⋯ | 4.777 × 10^{−04} | ⋯ | 1.650 | 9.101 × 10^{−04} | ⋯ | 4.775 × 10^{−04} | ⋯ | 2.601 |

$5th/20$ | 2.994 × 10^{−05} | 4.93 | 1.736 × 10^{−05} | 4.78 | 13.055 | 2.997 × 10^{−05} | 4.92 | 1.740 × 10^{−05} | 4.78 | 19.929 |

$5th/40$ | 1.084 × 10^{−06} | 4.79 | 5.816 × 10^{−07} | 4.90 | 139.432 | 1.088 × 10^{−06} | 4.78 | 5.833 × 10^{−07} | 4.90 | 211.149 |

$5th/80$ | 3.327 × 10^{−08} | 5.03 | 1.883 × 10^{−08} | 4.95 | 1119.209 | 3.335 × 10^{−08} | 5.03 | 1.889 × 10^{−08} | 4.95 | 1738.832 |

$6th/10$ | 1.426 × 10^{−04} | ⋯ | 4.309 × 10^{−05} | ⋯ | 2.227 | 1.414 × 10^{−04} | ⋯ | 4.321 × 10^{−05} | ⋯ | 3.704 |

$6th/20$ | 2.447 × 10^{−06} | 5.86 | 7.587 × 10^{−07} | 5.83 | 16.729 | 2.467 × 10^{−06} | 5.84 | 7.637 × 10^{−07} | 5.82 | 30.602 |

$6th/40$ | 3.537 × 10^{−08} | 6.11 | 1.207 × 10^{−08} | 5.97 | 175.908 | 3.536 × 10^{−08} | 6.12 | 1.212 × 10^{−08} | 5.98 | 304.109 |

$6th/80$ | 6.493 × 10^{−10} | 5.77 | 2.046 × 10^{−10} | 5.88 | 1493.933 | 6.572 × 10^{−10} | 5.75 | 2.055 × 10^{−10} | 5.88 | 2417.527 |

. | CWENO . | WENO . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Order/edges per side . | $eL\u221e$ . | $OL\u221e$ . | $eL2$ . | $OL\u221e$ . | Time
. | $eL\u221e$ . | $OL\u221e$ . | $eL2$ . | $OL\u221e$ . | Time
. |

$3rd/10$ | 1.124 × 10^{−02} | ⋯ | 6.897 × 10^{−03} | ⋯ | 1.000 | 1.119 × 10^{−02} | ⋯ | 6.880 × 10^{−03} | ⋯ | 1.252 |

$3rd/20$ | 1.514 × 10^{−03} | 2.89 | 9.535 × 10^{−04} | 2.85 | 8.240 | 1.506 × 10^{−03} | 2.89 | 9.492 × 10^{−04} | 2.86 | 8.757 |

$3rd/40$ | 1.993 × 10^{−04} | 2.93 | 1.244 × 10^{−04} | 2.94 | 85.831 | 1.982 × 10^{−04} | 2.93 | 1.236 × 10^{−04} | 2.94 | 94.658 |

$3rd/80$ | 2.487 × 10^{−05} | 3.00 | 1.581 × 10^{−05} | 2.98 | 673.159 | 2.473 × 10^{−05} | 3.00 | 1.571 × 10^{−05} | 2.98 | 798.829 |

$4th/10$ | 1.879 × 10^{−03} | ⋯ | 5.355 × 10^{−04} | ⋯ | 1.530 | 1.867 × 10^{−03} | ⋯ | 5.342 × 10^{−04} | ⋯ | 1.821 |

$4th/20$ | 7.487 × 10^{−05} | 4.65 | 2.762 × 10^{−05} | 4.28 | 10.864 | 7.427 × 10^{−05} | 4.65 | 2.739 × 10^{−05} | 4.29 | 14.081 |

$4th/40$ | 4.962 × 10^{−06} | 3.92 | 1.857 × 10^{−06} | 3.89 | 115.068 | 4.912 × 10^{−06} | 3.92 | 1.839 × 10^{−06} | 3.90 | 144.007 |

$4th/80$ | 3.482 × 10^{−07} | 3.83 | 1.200 × 10^{−07} | 3.95 | 928.223 | 3.451 × 10^{−07} | 3.83 | 1.188 × 10^{−07} | 3.95 | 1128.293 |

$5th/10$ | 9.115 × 10^{−04} | ⋯ | 4.777 × 10^{−04} | ⋯ | 1.650 | 9.101 × 10^{−04} | ⋯ | 4.775 × 10^{−04} | ⋯ | 2.601 |

$5th/20$ | 2.994 × 10^{−05} | 4.93 | 1.736 × 10^{−05} | 4.78 | 13.055 | 2.997 × 10^{−05} | 4.92 | 1.740 × 10^{−05} | 4.78 | 19.929 |

$5th/40$ | 1.084 × 10^{−06} | 4.79 | 5.816 × 10^{−07} | 4.90 | 139.432 | 1.088 × 10^{−06} | 4.78 | 5.833 × 10^{−07} | 4.90 | 211.149 |

$5th/80$ | 3.327 × 10^{−08} | 5.03 | 1.883 × 10^{−08} | 4.95 | 1119.209 | 3.335 × 10^{−08} | 5.03 | 1.889 × 10^{−08} | 4.95 | 1738.832 |

$6th/10$ | 1.426 × 10^{−04} | ⋯ | 4.309 × 10^{−05} | ⋯ | 2.227 | 1.414 × 10^{−04} | ⋯ | 4.321 × 10^{−05} | ⋯ | 3.704 |

$6th/20$ | 2.447 × 10^{−06} | 5.86 | 7.587 × 10^{−07} | 5.83 | 16.729 | 2.467 × 10^{−06} | 5.84 | 7.637 × 10^{−07} | 5.82 | 30.602 |

$6th/40$ | 3.537 × 10^{−08} | 6.11 | 1.207 × 10^{−08} | 5.97 | 175.908 | 3.536 × 10^{−08} | 6.12 | 1.212 × 10^{−08} | 5.98 | 304.109 |

$6th/80$ | 6.493 × 10^{−10} | 5.77 | 2.046 × 10^{−10} | 5.88 | 1493.933 | 6.572 × 10^{−10} | 5.75 | 2.055 × 10^{−10} | 5.88 | 2417.527 |

### C. DIM within discontinuous Galerkin framework

One of the first implementations of the DI models within the DG framework is presented in Ref. 180. The derivation of a solving algorithm is non-trivial as it is well known that for systems of equations containing non-conservative products like the ones characterizing the DI models, weak solutions in the classical form does not exist and Rankine–Hugoniot shock conditions cannot be defined. This problem was addressed for the DG formulation by Rhebergen *et al.*^{181} who used the theory developed by Dal Maso *et al.*^{182} (Dal Maso–LeFloch–Murat, DLM theory) where the left and right state across a discontinuity are path connected in phase space and developed a weak formulation and a numerical flux for space–time DG containing non-conservative products. In the work of Franquet and Perrier,^{180} this setting is used for solving multi-species problems with the BN model. However, the limiting procedure only consists of the truncation of higher order terms and lacks a proper treatment of potential oscillations. Given a non-conservative equation in the form

with $G=G(U)$ and $H=H(U)$, its equivalent weak form can be written as

The non-conservative product, following the DLM theory, can be expressed using the path function *ψ*:

de Frahan *et al.*^{183} followed this methodology and adopted a hierarchical reconstruction to limit the solution at discontinuities. The idea of the hierarchical reconstruction is to compute the polynomial coefficients within a cell using a MUSCL or WENO reconstruction hierarchically from the highest degree to the lowest, with the advantage to preserve high-order accuracy. The limiting step is then applied to the reconstructing coefficients of the total energy and to selected properties in the EOS to enforce conservation.

Another extension of the DG method to handle spatial derivatives higher than the first-order and non-conservative products is the local DG (LDG), introduced in Ref. 55, used in Ref. 184 for solving the Hamilton–Jacobi equation and extended to the DI models by Gryngarten and Menon.^{185} The aim of the LDG method is to rewrite the equations as a first-order system by introducing new variables, which approximate the solution derivatives and is motivated by the work of Bassi and Rebay^{186,187} for the discretization of the diffusive terms in the NS equation. In Ref. 185, the LDG is used to discretize Allaire's five-equation model [Eq. (13)]. Taking the volume fraction equation in the form of Eq. (93), the LDG system reads

and **q** is the auxiliary tensor used to reduce the order of the derivatives and *f*(*U*) is an arbitrary function of the state vector. The auxiliary tensor can be decomposed in $q\alpha m$ and $q\alpha p$, and is related to the source term through

Assuming that **q** can be approximated like the solution vector in Eq. (53) with the same basis function $Pjk$, and new coefficients $ci,j$,

and multiplying Eq. (127) by a test function $\varphi $ and integrating over the domain Ω_{i}, one obtains the following two equations:

Here, $F\u0302$ is the numerical flux which in Ref. 185 is approximated via an HLLC solver. When computing the flux $f\u0302$, the left side to the cell face quantities is used for $q\alpha m$ and the right side is used for $q\alpha p$. Then, the source $S\alpha $ can be approximated by the combination

with $uM$ the maximum absolute velocity within the element. The limiting procedure adopted in this work uses a troubled-cell detector with moment limiter modified for unstructured meshes with hanging nodes. The troubled cells are detected based on a TVD argument and are then *corrected* through a minmod limiter, which acts on the slopes extracted from Legendre polynomials. If after this step, a cell still has an invalid solution, this is further limited with a more diffusive procedure. In Ref. 185, it was tested that the whole solver is more accurate when reconstructing primitive variables rather than conservative variables for multiphase flows.

In the work of Saleem *et al.,*^{188} the DG method is used for solving the five-equation model of Kreeft and Koren.^{189} This is a velocity and pressure equilibrium model with a source term, which expresses the exchange of energy between the two fluids in terms of mechanical $L\u0307M$ and thermodynamical $L\u0307T$ rate of work. This formulation shares the same mass fractions, momentum, and energy conservation laws with the Kapila's model [Eq. (8)], but the volume fraction advection is substituted with an energy equation for one of the two fluids, which includes the new source term. The benefit of such substitution resides in a more straightforward integral formulation readily suitable for FV schemes that relaxes some of the numerical difficulties encountered in solving the Kapila's system of equations. In Ref. 189, the total rate of energy work is defined as

where *τ*_{1} and *τ*_{2} are the isentropic compressibilities of the two phases

*τ* and *ρ* are the bulk isentropic compressibility and the bulk density, respectively,

With the definitions above, this five-equation model in one-dimensional form is composed of the following vectors:

This system is solved in Ref. 189 in the FV context with a Osher-type Riemann solver, and the energy-exchange terms are integrated in the solution space instead of in the physical space. In Ref. 188, the DG scheme with minmod troubled cell detector and WENO limiter are used, while the flux are solved with a Lax–Friedrichs scheme.

A different approach to solve the five-equation model [Eq. (13)] in the DG framework was presented in Ref. 190, where the non-oscillatory kinetic (NOK) scheme is used to treat the Mie–Grüneisen EOS closure. The NOK schemes were introduced in Refs. 191 and 192 to solve the extended Euler equations for multi-species flows without introducing oscillations near material interfaces, by a flux splitting technique able to circumvent the instabilities near discontinuities brought by the classic gas-kinetic schemes. In Ref. 190, this flux is coupled to a troubled cell detector-nonlinear limiter similar to the one adopted in Refs. 185 and 188, and a maximum-principle-satisfying limiter built explicitly to maintain the bounds of the volume fraction between 0 and 1, applied to every Legendre–Gauss–Lobatto quadrature point.

A recent attempt to extend the DG discretization for the solution of a *γ*-based model of the form proposed in Ref. 29 can be found in Ref. 193. This new formulation differs from the path-conservative formulation of Ref. 183 in that it reduces to the quasi-conservative FV discretization of Ref. 15 if a *P*_{0} polynomial is used to represent the solution. This is obtained by a further integration by parts applied to the local discretized DG solution polynomial and carries the property of preserving uniform pressure and velocity fields at isolated material interfaces. The numerical procedure uses a flux solver based on the HLL flux and *a posteriori* correction, that is, recomputes DOFs with lower order schemes, on invalid cells, based on the procedure introduced for DG schemes in Refs. 194 and 195.

### D. DIM within hybrid DG-FV formulation

The reconstructed DG method demonstrated to be an accurate and efficient method for the solution of hyperbolic PDEs, especially under a careful selection of the polynomial basis and reconstruction procedure. Recently, Pandare *et al.*^{63} successfully solved a one-dimensional compressible multi-material problem with a high-order reconstructed *DG* method. The multi-species system of equation used is similar to Kapila's^{2} model, but with a different pressure-relaxation term in the form reported in Ref. 196, which is a more realistic modeling of situations with two material with net different bulk modulus. The pressure relaxation closure is thus defined as

with a finite pressure-equilibration time scale *τ,* which depends on a space-scale parameter *h* and can be adjusted through a constant scale parameter $c\tau $,

Small $c\tau $ values are generally avoided to not incur stiff pressure relaxation, and a problem independent choice for it is the local ratio between minimum and maximum density values of the materials $c\tau =\rho min/\rho max$.

Following Luo *et al.,*^{120} the DG formulation with a Taylor basis polynomial of degree *P _{N}* is employed for the solution discretization at the beginning of a time step and used to reconstruct a polynomial of degree

*P*before flux evaluation. Therefore, a third-order $rDG(P1P2)$ solution can be recovered from a first order underlying solution

_{M}*P*

_{1}, without increasing the total number of DOFs, and using the constrained least-square procedure exposed in Sec. III C.

In Ref. 63, the reconstruction is coupled to a Superbee limiter for the *P*_{1} moments and to a WENO limiter with TVB minmod troubled cell detector for the *P*_{2} moments. Furthermore, a consistent treatment of the higher-order DOFs for material densities and energies respect to the volume fractions is required in the mixture region to ensure properties discrete conservation. Similar to the interface treatment in the THINC approach, the slopes of bulk velocity, partial densities, and partial pressures are thus taken as zero. The Taylor basis is a convenient choice because in this case, it is sufficient to set the high-order DOFs as zero, allowing to calculate the derivatives of the material energy simply with

having the cell-average value of the specific total energy of the material *k*$\rho kEk\xaf$. This setting has shown interesting outcomes in terms of formal accuracy and computational efficiency. An example is the advection of a volume fraction Gaussian test in one-dimension for which the results FV, DG, and rDG schemes are compared (Fig. 5).

### E. DIM within hybrid FV-LB formulation

A few applications of diffuse interface models in the context of hybrid finite volume-lattice Boltzmann framework can also be found in the literature. Most of these applications rely on a conservative form of the mixture governing equations developed by Wang *et al.*^{197} This model reformulates the mixture ratio of specific heats derived from the simple imposition of mixture total energy conservation at any time, in terms of total enthalpy. This so-called thermodynamically consistent and fully conservative (TCFC) model assumes that all species move at the bulk velocity of the mixture and the complete system of equation, in addition to the NS equations, includes one additional equation for the specific heat of the mixture $\chi =\gamma \gamma \u22121$ and one for the mixture molecular mass *M*. In one-dimension, the system of Eq. (1) is defined by the only two vectors *U* and *F*(*U*),

Joshi *et al.*^{77} solved this system with the hybrid FV-LB procedure discussed in Sec. III F. This hybrid scheme overcomes the limits of the LBM, allowing the adoption of non-uniform grids and the results for the one-dimensional case shows the effectiveness of the LBM flux solver in capturing contact interfaces compared to the classic Riemann solvers. However, it would be interesting to evaluate the stability of the scheme for high Mach flows and for more stringent applications.

More recently, Sashi Kumar and Maheshwari^{155} presented a similar approach but corrected an inconsistency of the Wang model and substituted the KT model used by Joshi *et al.* with the multi-entropy-level LBM of Yan and Zhang.^{198} The inconsistency found in Ref. 197 lies in the derivation of Eq. (139). Following the expansion of the equation for the mixture specific heat:

one notices that this is the combination of the gamma-based model and the corrected model presented in Ref. 199. It is required that both terms in curly brackets approach zero to satisfy conservation of mass. However, this requirement may not be generally achieved, causing oscillations. In Ref. 155, a predictor–corrector stepping is adopted to enforce the conditions imposed by Eq. (140), where within the predictor step the quantity $(\rho \gamma )$ and the rest of the variables are first evolved in time before solving for $1\gamma \u22121$. Therefore, Eq. (139) is modified to be consistent with the predictor–corrector steps

The interface flux calculation is similar in philosophy to the flux vector splitting and to the previous procedure adopted in Ref. 77. The flux splitting is simply expressed as

where the $+,\u2212$ superscript refers to all the contributions in the stencil coming from the two sides of the stencil respect to the cell interface. Considering a four-point stencil, with nodes 1 and 3 placed at the right to the interface, and nodes 2 and 4 at the left to the interface, the contributions of each of them are split with the following convention:

The left fluxes are defined by just exchanging the *R* with the *L* superscript. The flux at the interface is then calculated by

Each of the contribution $F(i\u22121/2),iL,R$ is calculated through a LBM step, and the values are reconstructed through a high-order method. In Ref. 155, two reconstruction schemes are finally presented, namely, lattice-Boltzmann multi-species flux-based solver (LAMBS) with MUSCL (LAMBS-MUSCL) and WENO (LAMBS-WENO) reconstructions, and tested for more stringent compressible-multispecies flows, demonstrating non-oscillatory behavior and good accuracy obtained from an easy and straightforward implementation.

## F. Extension to fluid–structure interaction

Fluid–solid interaction with diffuse interface models was studied by Favrie and Gavrilyuk in a series of papers.^{8,96,97,105} As discussed in Sec. II, this augmented DI models include a hyperelastic description of the solid phase, which was previously studied in Refs. 35, 36, 101, and 104. The augmented system is generally in the form exposed in Eq. (27), and the presence of a compaction term like the one present in the Kapila's^{2} model [Eq. (8)] allows to account for material fragmentation, as studied by Ndanou *et al.*^{94} To solve numerically the visco-plastic model, a hyperbolic step is followed by a stress relaxation and a pressure relaxation. The hyperbolic stepping is performed through a Godunov type scheme decoupled in a conservative and a non-conservative part. Although the system contains seven waves, a three-wave HLLC solver is used as it is sufficient to capture longitudinal and transverse waves. However, the solution of the non-conservative part requires jump relations, which in the classical form induces entropy jumps in the presence of shocks that are not necessarily exact. Some extra steps are thus required to correct the hyperbolic step. First, a plastic relaxation correction is performed, during which the deviatoric part of the stress tensor will be forced to tend the yield surface if the Von Mises yield criteria is not initially met. Consequently, such the materials will have different pressures, and a pressure relaxation step is performed to determine a relaxed pressure, which satisfies the phase energy equation. This is then used to compute the mixture pressure from the mixture EOS and used to reset the internal energies. The source term from Kapila's model is here responsible for the mimicking of the mechanism of crack propagation. Indeed, in Ref. 94 the crack formation is explained as a strong cavitating process in solids, due to the coalescence of the gas bubbles growing from the tensile waves generated by the impact.

Fast transient fluid–structure dynamics was also studied by Faucher *et al.*^{201,202} using the VOFIRE (finite-volume with reconstruction) algorithm. The VOFIRE is an anti-diffusive scheme introduced by Després *et al.,*^{203} where the dissipation is controlled with a constrained downwind scheme. In illustrating the basic idea behind the VOFIRE scheme, let us define the grid coordinates for the cell center $xj=(j+1/2)\Delta x$ and cell interface $xj+1/2=(xj+1+xj)/2$, with a uniform space interval $\Delta x$. Assuming that at a certain time instant *t ^{n}* an interface lies in a cell

*j*, then in order to build the fluxes for the time update, only the discrete values $hj+1/2$ at the cell boundary $xj+1/2$ are necessary. Then, the constraints for the downwind scheme are the condition for flux consistency,

and the condition for stability through advection,

with $hj*=hj+u\Delta t\Delta x(hj\u22121/2\u2212hj+1/2)$. The intervals for which $hj+1/2$ meets both stability and consistency criteria are

and anti-dissipation is achieved by taking $hj+1/2$ closest to the downwind value. The flux values for the outgoing faces *k*, namely, $hj,kn$, are retrieved through a transverse reconstruction step and through a longitudinal reconstruction. The latter is the one that accounts for most of the anti-dissipation, and consists of finding the appropriate coefficients $\mu j,k,r$ for the following pseudo-one-dimensional problem:

including all the incoming faces *r*. The upwind and the downwind values for the outgoing cell face *k* are referred to as $hjn$ and $hkn$, respectively. Therefore, $hj,kn$ takes the downwind value when the coefficients $\mu j,k,r$ equal unity, and the upwind value when they equals zero. When applied to mixture problems, like in Ref. 201, the field considered in the reconstruction algorithm is naturally the volume fraction.

The hybrid multiphase-mixture model introduced by Michael and Nikiforakis^{110} and discussed in Sec. II can be coupled to a elastoplastic system of equations as shown in Ref. 9 and also to a third system of equations such the resistive MHD model.^{98} From a numerical standpoint, the coupling between the systems uses a diffuse interface approach, but the location of the material interfaces is retrieved with a level-set method. However, the evolution at the interface is accounted with a Riemann ghost fluid method following the procedure introduced by Sambasivan and Udaykumar.^{204,205} In this procedure, a band of ghost points is defined in the area adjacent to the interface and the input states at both sides for the mixed-material Riemann solver are taken at a distance $1.5\Delta x$ from the interface coordinate in order to avoid potential errors carried by the interface. The Riemann solver subsequently uses these characteristic states to return the well-known left and right star regions, similarly to Eq. (95), which in practice corresponds to the regions surrounding the material interface. The ghost cells are then updated with these values.

This kind of solver allows for communication between materials described by different systems of equations as long as they can be cast in hyperbolic form. In the works of Michael and Nikiforakis, the solver is validated for several sandwich-plate explosives tests, as well as ignition of a combustible tank by plasma arc.

### G. Helium bubble–shock wave interaction test case

A popular benchmark for multi-species solvers is the simulation of the deformation of a gas bubble following the interaction with a planar shock wave. This was investigated experimentally by Haas and Sturtevant^{206} and numerically by Quirk and Karni^{207} for studying two different bubble configurations. In one setting, the bubble is filled with Helium, and in the other with a refrigerant R22. The traveling of the shock wave across the bubble starts an impulsive mixing process and the difference in density between the gas bubble with the surrounding air will trigger vorticity production resembling the Richtmyer–Meshkov^{207,208} instability mechanism. The vorticity direction depends upon the Atwood number, defined as

where the density of Helium and R22 has to be considered in place of *ρ _{g}*. For the case of Helium gas bubble, assuming a shock incident from the right, the vorticity is clockwise because the Helium is lighter than air (negative Atwood number), and the vortex roll-up move the bubble structure downstream to the surrounding fluid. The opposite mechanism happens for the R22 bubble as the higher density leads to positive Atwood numbers, and the counterclockwise vortices move the bubble structure upstream.

^{209}In both the experiments and the computations, the Mach number for the shock wave is took as $Ms=1.22$, and the bubble is in initial equilibrium with the surrounding air. In the initial phases of the shock–bubble interaction, complex wave structures involving refraction, reflection, and diffraction develop, providing a challenging context for compressible solvers and shock capturing schemes, whereas higher order schemes would be beneficial for a correct evaluation of the instabilities.

We are now going to report some of the results for the Helium bubble case that can be found in literature on different numerical frameworks and with sharpening techniques. The computational domain traditionally consists of the one represented in Fig. 6, where the helium bubble of diameter $Db=5\u2009cm$ also contains air in a 28% mass concentration.

The initial conditions are defined as follows:

with specific heats equal to 1.4 and 1.66 for air and helium, respectively. The bubble dynamics is often described in terms of the position of the upstream, downstream, and jet interfaces, and this is used as a validation argument for the numerical schemes. In Fig. 7, the interface position retrieved with different schemes is reported, that is, the WENO-JS,^{46} bandwidth-WENO^{165} (WENO-SYMBO) and MDCD-WENO of^{166} on uniform meshes, the CWENO on unstructured meshes of Ref. 200 and the hybrid LB-FV schemes of Ref. 155 in the variants with MUSCL limiter and WENO reconstruction (LAMBS-MUSCL and LAMBS-WENO), that are in good agreement with the experimental evidence of Ref. 206 and early computational results of Ref. 207.

The qualitative analysis in terms of bubble shape through density and volume fraction is illustrated in Fig. 8 and gathers the results from previous references in addition to the results obtained with the fifth WENO scheme of Ref. 53, the MUSCL-THINC/QQ-BVD scheme of Ref. 177, and the second polynomial order DG of Ref. 193. The benefits from the adoption of high-order schemes for the reduced interface smearing and good definition of flow structures can be clearly observed, as these schemes correctly reproduce the jet and vortex rings formation. The impact of different discretization levels in the DG framework (Fig. 9) also results in a sharper interface resolution with the adoption of higher polynomial orders, and the same applies for the hybrid lattice Boltzmann scheme, Fig. 10, where the LAMBS-WENO scheme provides better interface definition and richer structures over the LAMBS-MUSCL scheme.

### H. Shock in molybdenum–encapsulated MORB interaction test case

The shock wave propagation in molybdenum through a region of encapsulated mid-ocean ridge basalt (MORB) was first introduced by Miller and Puckett in Ref. 210. For the setup of this test case, a unit square region is considered and the MORB is contained in a rectangle of dimension $[0.4,0.7]\xd7[0,0.5]$. Inflow boundary condition is imposed at the left side and outflow conditions at the right and top sides. Reflecting boundary is applied at the bottom side and a Mie–Grüneisen EOS is used for the materials with the following coefficients:

The initial conditions, with zero reference state for pressure and energy, are defined as

and generates a Mach 1.163 right moving shock traveling toward the encapsulated MORB liquid. The computations of the problem are shown at times $t=50\u2009\mu s$ and $t=110\u2009\mu s$ for density contours (Fig. 11) and pressure contours (Fig. 12). At time $t=50\u2009\mu s$, one can observe through the density plot the incident shock and the reflected shock in molybdenum and MORB, respectively, with the latter moving slower than the former. At time $t=110\u2009\mu s$, the transmitted shock still have not traveled the MORB capsule completely. Figures 11 and 12 compare results from different references computed using FV with and without sharpening techniques, and DG method with different polynomial order approximation.

The FV method, coupled to a second-order MUSCL reconstruction and without sharpening techniques [Fig. 11(b)^{,}^{41}], clearly shows the molybdenum–MORB interface diffusion at the later time, which is substantially improved by the adoption of the THINC reconstruction [Fig. 11(a)^{,}^{41}] or by the introduction of the anti-diffusion terms of So *et al.*^{176} [Fig. 11(c)]. However, some overshoots can also be observed in the pressure plot corresponding to the THINC reconstruction performed on a grid resolution of $N=200\xd7200$. On the other hand, the contours here illustrated from Ref. 176 with anti-diffusion terms are computed on a finer grid ($N=400\xd7400$), and therefore, the interface thickness is improved also by the higher grid resolution.

With the same grid size, Luo *et al.*^{190} presented this test case computed within the DG framework. The results clearly show the benefits in terms of interface thickness by adopting a higher order polynomial approximation for the solution representation. The interface thickness obtained with *P*_{2} order is comparable with the one obtained with THINC reconstruction, although similarly to the latter, some instabilities seem to affect both density and pressure contour plots at the later time.

A good indicator of the amount of interface diffusion is given by the definition of the high speed jet at the later time near the left upper tip of the MORB block. With the second-order FV discretization, this feature is indistinguishable from the interface structure, while it is clearly visible when using sharpening or higher-order discretization schemes.

## V. CONCLUSIONS

The diffuse interface models are widely used for the computation of multi-species flows. The versatility of these models is demonstrated by the variety in complexity and physical features that can be accounted for by augmenting or reducing the governing system of equations. On the other hand, the well-known drawback of the interface smearing required a considerable effort in recent years to limit the numerical error, especially for longer time scale simulations. The sharpening and anti-diffusion techniques proved to have the potential to remedy this issue, but the efficacy and the computational efficiency for three-dimensional problems is still an open question.

The recent advancements in high-order accuracy schemes can aid these techniques and are similarly effective, when the discretization level is carefully calibrated.

A wealth of schemes that allows to target high-order solution is now available although, as illustrated, the diffuse interface method have reached a fairly good level of maturity only within the FV framework. Applications on unstructured and three-dimensional meshes are nowadays well established in conjunction with MUSCL and WENO reconstructions, although only in the FV framework. Unfortunately, in the context of high-order methods, the FV schemes exhibit the limit of relying on wide stencils with the high computational requirements in three-dimensions, and the increased communication time in parallel computations. In this sense, the compact-WENO scheme, which operates on reduced stencils, seems to alleviate consistently the computational costs.

The DG formulations remove the necessity of large stencils even if the computational costs are still substantial since more degrees of freedom are associated with each cell. In addition, the treatment of the non-conservative terms is not well established and thermodynamically consistent procedures are available only for selected EOS. Therefore, the application of the DIM within the DG frameworks seems to be, so far, mostly limited to the reduced five-equation models and for the majority to structured meshes. The reconstructed DG method aims to reduce the number of DOFs to be evolved at each time step, and the results obtained in literature for classic single fluid flows are significant in terms of accuracy and computational efficiency. Applications with DIM are available, to the best of our knowledge only for one-dimensional problems. However, extensions to two-dimensional and three-dimensional problems will probably appear in the near future.

An interesting alternative appeared in the last decade in the form of a hybrid scheme, which combines the ease of implementation and the natural predisposition for multiphase and multi-species flows of the lattice-Boltzmann with the classic features of the FV schemes, alleviating the restrictions of the former with respect to the compressible flow limit and irregular meshes. Again, applications with DIM are at an early stage and cases with strong shock waves and highly irregular meshes need to be explored.

The majority of the works related to high-order applications of DIM lack convergence studies, as well as quantitative analysis to assess the performances of the schemes in terms of interface diffusion in multiple dimensions. In this sense, a library of common test problems should be established in order to provide a solid benchmark for the comparison of the existing and future numerical schemes.

## ACKNOWLEDGMENTS

The authors acknowledge the support by EPSRC (Grant No. 2497012), of Innovate UK (Grant No. 263261), and Airbus UK.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.