High-order methods for diffuse-interface models in compressible multi-medium flows: A review

The diffuse interface models, part of the family of the front capturing methods, provide an efﬁcient and robust framework for the simulation of multi-species ﬂows. They allow the integration of additional physical phenomena of increasing complexity while ensuring discrete conservation of mass, momentum


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,9he 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,11Typically, 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,15nother 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 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. 73As 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,81A 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 selfequilibration 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 X; C and Ẽ , respectively.
An additional equation is necessary to evaluate the volume fraction of the phases, given the saturation condition / s þ / g ¼ 1 of the solid volume fraction / s and gas volume fraction / g .The onedimensional 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 2 u k Á u k and e k is the internal energy.Velocity and pressure of phase k will be, respectively, u k and 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 X; C and Ẽ 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, Ũ, 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 b 1 measures the forces between the grains and estimates of this value can be found in the work of Elban and Chiarito. 83The quantity l c is defined as the dynamic compaction viscosity, and it is an indicator of the rate at which pressure equilibrium is reached. 15A 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 l through a volume variation proportional to the compaction viscosity and to the difference between the gas phase pressure and the interfacial pressure P i , @/ g @t ¼ l c ðP g À 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 k 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 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 i , 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 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 fiveequation 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 Y k ¼ ð/qÞ k q and the mixture density q ¼ ð/qÞ 1 þ ð/qÞ 2 , e ¼ Y 1 e 1 ðq 1 ; pÞ þ Y 2 e 2 ðq 2 ; pÞ: 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 1 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 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 cÀ 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 c, 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=ðc À 1Þ instead of c, and this four equation model reads  : (12)   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€ uneisen EOS. 29,31More 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 cÀ 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,

:
(13 The system is closed with an EOS, and in Ref.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 phaseequilibrium 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 Z k ¼ q k c k , 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.

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 c l and p 1 are determined through the shock wave Hugoniot data. 92The 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-stiffenedgas model, which accounts for covolume and repulsive effects, and was introduced by Le M etayer and Saurel in Ref. 93 and reads The constants p 1;k and b k are characteristics of the fluid, and q k 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€ uneisen 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 with g ¼ q=q 0 ; C 0 a material constant, a 0 the bulk speed of sound, and q 0 the reference density.The parameter s represents the Hugoniot slope coefficient, defined as s ¼ dUs dUp , 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][97][98] that can also be enhanced by peridynamics models. 99As 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, r Â g T ¼ 0, of the inverse of the deformation tensor  36 where the limits of the viscoplastic 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. 103he 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, displacementbased, 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. 105Plasticity 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 ¼ ðX l Þ; l ¼ 1; 2; 3 and the Eulerian coordinates be x ¼ ðx e Þ; e ¼ 1; 2; 3.If the deformation or diffeomorphism ỹ, from a reference configuration P to the actual configuration P is written as x ¼ ỹð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 e ij through with d ij being the Kronecker symbol.A cobasis vector corresponding to the Lagrangian coordinates is E l ¼ rX l and satisfies the following: By applying the Hamilton principle, a pressure equilibrium elastic solid-gas diffuse interface model with plasticity effects can be written @/ g q g @t þ divð/ g q g uÞ ¼ 0; @/ s q s @t þ divð/ s q s uÞ ¼ 0; 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 p and the stress tensors for solid r s and gas r g are, respectively, 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, 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. 108This 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 a with products of reaction b.The model considers velocity and pressure equilibrium, between phase 1 and phase 2, and also between the mixture components a and b, 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 k is the mass fraction of the reactants.All materials are modeled by a Mie-Gr€ uneisen 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 (k ¼ 1 or k ¼ 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 k 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. 98In 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 condensedphase 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.

Physics of Fluids
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. 79he 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 X with boundaries @X and outward normal unit ñ, is ð Assuming the domain to be a tessellation of small non-overlapping N elements V i with volume jV i j, where i ¼ 1; :::; N, and applying the integral Euler formulation over a small control volume V i , the semidiscrete FV scheme reads The state variable and face flux integrals have been replaced, respectively, by the volume averaged state variable U i and averaged normal flux Hereafter, we will assume a domain X that can be any combination of hexahedral, pyramidal, or prismatic elements, with a surface area of jf j 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. 118The accuracy of the approximation depends upon the selected number of quadrature points N qp , denoted with x f ;a , and their distribution and weights x a depend upon the Gaussian rule.The approximate face integral is, therefore, of the form 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 U L ðx f ;a Þ and U R ðx f ;a Þ.The flux function, in terms of the Riemann solver F 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 p i ð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 m of the neighboring cells constituting the computational stencil S .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,

Physics of Fluids
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,57Considering for example a tetrahedral element with vertices Þ, the transformation from the physical coordinate system (x, y, z) to the reference coordinate system ðn; g; fÞ reads with J, the Jacobian matrix, defined as This transformation defines both the direct x ¼ xðnÞ and inverse n ¼ nðxÞ mappings from n ¼ ðn; g; fÞ into x ¼ ðx; y; zÞ.Applying the inverse mapping to all the elements V m of the reconstruction stencil S , the transformed elements V 0 m will constitute the transformed stencil S 0 , Over this reference system, the reconstruction polynomial is sought as expansion over local polynomial basis functions k k ðn; g; fÞ, leading to pðn; g; fÞ 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 0 m of the stencil S 0 , 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 m the right hand side vector, respectively,

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 /ðx; y; zÞ, integrating over the domain X, and performing an integration by parts.This results in ð X /ðx; y; zÞ 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 V h ¼ f/ h 2 L 1 : / h j X 2 V k ðXÞ; k ¼ 0; 1; 2; …; Ng, where V k is the space of polynomials of degree up to k.The discrete solution U h , with expansion coefficients denoted by u h , can be seen as expansions over a finite-element basis P k j in the aforementioned polynomial space, where d is the number of degrees of freedom U h ðx; y; z; tÞ ¼ X d j¼1 u h ðtÞP k j ðx; y; zÞ: 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 jf j 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 / 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 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 P k j 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.1][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,118However, 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 P N P M 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 M , reconstructed from the underlying polynomial P N .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 ðU l Þ and the hierarchical orthogonal reconstruction basis W l are chosen in such a way to satisfy the following properties:

Physics of Fluids
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 P N P M 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 kth-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 ðP 1 P 2 Þ, or RDGðP 1 P 2 Þ, 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 P N P M 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 kth-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 steadystate 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 U ij;a of element i and face j at quadrature point x f ;a , and multiplying all the polynomials for the quadrature points by a limiting function H 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. 132Although 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 H il;a is commonly expressed as The maximum U max i and the minimum U min i values of the state variable over the stencil are required, and the limiter H i selects the minimum value from the slope limiter H il;a from all the side neighbors, and quadrature points On the other hand, for a WENO reconstruction on a threedimensional 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 p m ðn; g; fÞ, x m p m ðn; g; fÞ; 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 and the non-linear weights x m are defined as in the classical formulation in Refs.46 and 133, The linear weights d m take a large value for the central stencil (10 2 < d 0 < 10 5 ), although in Ref. 134 where D is the derivative operator, b 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. 62Balsara 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 multispecies 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. 67However, 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. 145ollowing the latter, one can describe the time and space evolution of the velocity distribution f through The variable c denotes the particle velocity, and X 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 s 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; T 0 ).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. 143s 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. 149The 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 c 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,155It 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

@f
and the quantities h 0 ¼ c Á u À 1 2 ðf À f eqÞ and 1 shf ¼ 1 sh À 1 sf are defined from Guo's model. 156If one considers, the following density f eq and energy h eq equilibrium distributions, respectively: ; with D the spatial dimension, R the gas constant, c the molecule velocity, and q and u the fluid density and velocity, respectively, then the thermal compressible Navier-Stokes can be retrieved by a Chapman-Enskog expansion. 156,157As in the work of Shan and Chen, 70 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 xðc; TÞ is defined as xðc; TÞ ¼ 1 The integration from time t n to time t nþ1 leads to the so-called discrete unified gas-kinetic scheme (DUGKS). 158The 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 with v 1 6 ¼ v 2 6 ¼ 0 non-zero constants, and the local equilibrium particle velocity distribution function is given by with k 2 ½1; 9. We refer to Refs.159 and 77 for the determination of the coefficients A k , B k, and C k .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 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. 30One 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 twopoint 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 ¼ fu; vg.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 and c 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. 162Deng 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 e, 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 À 20 cells per wavelength are used to compute the highest wavenumber component.

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 ¼ r/=jr/j; U as a characteristic compression velocity and e h as a grid spacing length scale.Since we are interested in sharpening the interface, small values of e h and large values of U 0 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: The function Hð/Þ 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 Lagoutie `re 175 proposed an antidiffusive 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 postprocessing 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 s 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 /.The remaining flow variables are updated according to the flux obtained for the new value of /, 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 ðr/Þ lim with the gradient calculated with a central difference scheme ðraÞ 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ð/ 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 q 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 ¼ R1 þ R2 .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 threedimensional 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 U is used as an indicator of the presence of an interface within a cell C i , with an associated cell average U i , and which takes unity or zero value whether C i is filled with a fluid a or not.For the intermediate values, the hyperbolic tangent function Ũi ðxÞ mimics the jump across an interface

Physics of Fluids
given the sign function ÀU iÀ1 ðtÞ and a user-defined slope parameter b.Assuming conservation of volume fraction in the form the interface location xi can be determined from Eq. ( 110) as it is the only unknown, This allows the determination of the function U 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-pressureratio 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 threedimensional 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 H il;a ¼ max 0; minðb Hil;a ; 1Þ; minð Hil;a ; b h i ; by setting b ¼ 1 and b ¼ 2, respectively.It has to be noted that the reconstruction is here performed on the primitive variables W il;a .The newly introduced limiter, named Overbee, 42 reads where the term t is defined as t ¼ max minð2 Hil;a ; bÞ; min ð2 À bÞ Hil;a þ 2ðb À 1Þ; Hil;a n o h i : For instance, when b ¼ 1 this corresponds to the Superbee limiter, while for b ¼ 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 b ¼ 2, they demonstrated the capability of the method to keep the resolution of the interface discontinuity within 3 6 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 Q nn i , the first step is to calculate the boundary variation (BV) at a cell edge C ij through

Physics of Fluids
using the integrated values of the reconstruction function on the cell X i and its neighbor X ij , Then, the TBV of a reconstruction candidate Q n i ð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 N, 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 N. 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 b, 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 and the reconstruction polynomial is a non-linear combination of the polynomials p s as defined in Eq. ( 76).The central stencil k 1 is calculated by normalizing an arbitrary value assigned to the nonnormalized weight k 0 , 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).

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 w: 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. 185The 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 @U @t þ r Á FðUÞ ¼ ÀHðUÞr Á u; q ¼ rf ðUÞ; 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 am and q ap , 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 P k j , and new coefficients c i;j , and multiplying Eq. ( 127) by a test function / and integrating over the domain X i , one obtains the following two equations: Here, F is the numerical flux which in Ref. 185 is approximated via an HLLC solver.When computing the flux f , the left side to the cell face quantities is used for q am and the right side is used for q ap .Then, the source S a can be approximated by the combination with u M 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. 189This 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 M and thermodynamical _ L T 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 s 1 and s 2 are the isentropic compressibilities of the two phases s and q are the bulk isentropic compressibility and the bulk density, respectively, With the definitions above, this five-equation model in onedimensional 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 nonoscillatory kinetic (NOK) scheme is used to treat the Mie-Gr€ uneisen 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 c-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 onedimensional 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 s, which depends on a space-scale parameter h and can be adjusted through a constant scale parameter c s , Small c s 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 s ¼ q min =q 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 M before flux evaluation.Therefore, a third-order rDGðP 1 P 2 Þ solution can be recovered from a first order underlying solution 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 q k E k .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 v ¼ c cÀ1 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-entropylevel LBM of Yan and Zhang. 198The 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 ðqcÞ and the rest of the variables are first evolved in time before solving for 1 cÀ1 .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 þ; À 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 L;R ðiÀ1=2Þ;i 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,105As 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 viscoplastic 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 nonconservative part.Although the system contains seven waves, a threewave 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 es 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 x j ¼ ðj þ 1=2ÞDx and cell interface x jþ1=2 ¼ ðx jþ1 þ x j Þ=2, with a uniform space interval Dx.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 h jþ1=2 at the cell boundary x jþ1=2 are necessary.Then, the constraints for the downwind scheme are the condition for flux consistency, ) and the condition for stability through advection, The intervals for which h jþ1=2 meets both stability and consistency criteria are and anti-dissipation is achieved by taking h jþ1=2 closest to the downwind value.The flux values for the outgoing faces k, namely, h n j;k , 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 l j;k;r for the following pseudo-onedimensional problem: including all the incoming faces r.The upwind and the downwind values for the outgoing cell face k are referred to as h n j and h n k , respectively.Therefore, h n j;k takes the downwind value when the coefficients l 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. 98From 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,205In 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:5Dx 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 wellknown 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 q 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. 209In both the experiments and the computations, the Mach number for the shock wave is took as M s ¼ 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 D b ¼ 5 cm also contains air in a 28% mass concentration.
The initial conditions are defined as follows: ð/ 1 q 1 ; / 2 q 2 ; u; v; p; / 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  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 Â ½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  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 ls and t ¼ 110 ls for density contours (Fig. 11) and pressure contours (Fig. 12).At time t ¼ 50 ls, 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 ls, 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 Â 200.On the other hand, the contours here illustrated from Ref. 176 with anti-diffusion terms are computed on a finer grid (N ¼ 400 Â 400), 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  V C Author(s) 2022

Physics of Fluids
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 nonconservative 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 threedimensional 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.

FIG. 3 .
FIG. 3. Example of bubble shapes at different times for the single-bubble collapse test resulted for (a) a pressure ratio of 10 and 25 nodes on each coordinate direction, and (b) pressure ratio of 1427 and 50 nodes.[Adapted with permission from Schmidmayer et al., "An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics," J. Comput.Phys.402, 109080 (2020).Copyright 2020 Elsevier Ltd.]

FIG. 5 .
FIG. 5. Convergence order comparison for the one-dimensional advection test of a Gaussian volume fraction.[Reproduced with permission from Pandare et al., "A reconstructed discontinuous Galerkin method for multi-material hydrodynamics with sharp interfaces," Int.J. Numer.Methods Fluids 92, 8 (2020).Copyright 2020 John Wiley and Sons.]
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 m ¼ 1).The parameters p and e 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 e ¼ 10 À6 are used.The smoothness indicator IS m , as the name suggests, indicates how the solution is smooth on the stencil, and is calculated with ) Phys.Fluids 34, 021301 (2022); doi: 10.1063/5.007731434,021301-14V C Author(s) 2022