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

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

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

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

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

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

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

On the other hand, another framework, that is, becoming more widespread for the simulation of multiphase flows, is the lattice Boltzmann method (LBM), as the continuum assumption is ignored, and instead, a mesoscopic kinetic operator is used to connect the molecular physics to a set of macroscopic properties obeying the Navier–Stokes (NS) equations. The LBM can be easily parallelized64,65 and adapted to complex topologies66,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 Keller69), 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 tests66,73,74 indicate that the color Gradient method provides the sharper fluid interfaces, although it is prone to numerical instabilities due to the lattice pinning, a phenomenon that also causes incorrect interface advection.73 As pointed by Qu,75 the LBM is commonly limited to incompressible flows since the model is derived from a low Mach number expansion of the Maxwellian distribution. Qu derived a compressible model valid for a wide range of Mach numbers applying a FV scheme for the solution of the discrete velocity Boltzmann equation (DVBE) as in the work of Nannelli and Succi76 and adopted an exact Riemann solver with a MUSCL scheme to capture the discontinuities. Contrarily, in the work of Joshi et al.,77 where a hybrid LBM is built for compressible multi-species flows, the intercell fluxes are evaluated through a LBM and a FV scheme is used to update the node parameters. In addition, with such blending of the LBM with the FV schemes, another limit of the former method can be overcome, that is, the difficult applicability to unstructured meshes. Applications on irregular grids are becoming, in general, a necessity in industrial applications, especially when complicated geometries are considered and present a number of advantages over structured grids, ranging from the ease of mesh smoothness requirement,78,79 to the maintenance of load balance in parallel computing.80,81 A detailed review of the application of different high-order methods on unstructured meshes is presented in Ref. 79.

Given this variety of schemes for the solution of multi-species flows and the recent introduction of hybridizations of the classical FV schemes with both the DG and LBM frameworks, this review proposes, without the ambition to cover the totality of the available computational frames, to gather the recent advancements in high-order accuracy schemes on unstructured grids and to put them in perspective for the simulation of multi-species flows with the diffuse interface method. A recent review of the main DIM models is available in Ref. 27; in this work, we focus on the implementation and adaptations required by these models within different numerical frameworks, emphasizing the challenges and advantages relative to each scheme, and comparing their efficiency and accuracy. The review is organized as follows: the DIM is presented and the variations to the complete Baer and Nunziato1 and the four-equation model of Abgrall30 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.

A cornerstone in multi-species modeling is the model introduced by Baer and Nunziato,1 from here on abbreviated with the Baer and Nunziato (BN) model, which was originally proposed for the description of reactive multiphase flow, and, in particular, the deflagration to detonation transition (DDT). In this process, the interstices in a granular reactive bed are filled with a gas phase forming an immiscible mixture. Each phase is treated as a compressible fluid in local thermodynamic equilibrium, but the mixture is allowed to be in nonequilibrium across the interface, and each phase is described through a set of inviscid conservation laws, which also accounts for the interaction between phases. As pointed by Bdzil et al.,82 this implicitly means that the self-equilibration time scale of each phase is much smaller than the equilibration time scale between the phases. Each conservation law includes a source term modeled to represent the transfer of mass, momentum, and energy between the phases, here named Ω̃,Γ̃ 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 one-dimensional system for the multi-species flow model can be expressed in vectorial form,

Ut+·F(U)+H(U)·u=S(U),
(1)

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: Ek represents the total energy with Ek=ek+12uk·uk and ek is the internal energy. Velocity and pressure of phase k will be, respectively, uk and pk.

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 Kapila2,82 reads

U=[ϕ1ϕ1ρ1ϕ2ρ2ϕ1ρ1u1ϕ2ρ2u2ϕ1ρ1E1ϕ2ρ2E2],F=[ϕ1u1ϕ1ρ1u1ϕ2ρ2u2ϕ1ρ1u12+ϕ1P1ϕ2ρ2u22+ϕ2P2ϕ1u1(ρ1E1+P1)ϕ2u2(ρ2E2+P2)],S=[Φ̃Ω̃Ω̃Γ̃Γ̃ẼẼ].
(2)

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 Ω̃,Γ̃ 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

Φ̃=ϕ1ϕ2μc(P1β1P2).
(3)

The intergranular stress β1 measures the forces between the grains and estimates of this value can be found in the work of Elban and Chiarito.83 The quantity μc is defined as the dynamic compaction viscosity, and it is an indicator of the rate at which pressure equilibrium is reached.15 A similar set of equations was derived by Saurel and Abgrall15 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 Abgrall15 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 Pg = Pl through a volume variation proportional to the compaction viscosity and to the difference between the gas phase pressure and the interfacial pressure Pi,

ϕgt=μc(PgPi).
(4)

Similar considerations hold for the velocity, whereas the relaxation process depends on the drag force Fd, which can be expressed in terms of a velocity relaxation coefficient λ and the velocity jump across the interface,

Fd=λ(ulug).
(5)

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 μc=λ+.

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,

S=[μc(PgPl)00Piϕgx+λ(ulug)Piϕgxλ(ulug)PiViϕgx+μcPi(PgPl)+λVi(ulug)PiViϕgxμcPi(PgPl)λVi(ulug)].
(6)

Questioning the acoustic properties of the original BN model1 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 Pi and Vi, 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 Abgrall15 proposed an estimate where the interfacial pressure is equal to the mixture pressure and the velocity is the one of the center of mass,

Pi=ϕkPk,Vi=ϕkρkukϕkρk.
(7)

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

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

U=[ϕ1ϕ1ρ1ϕ2ρ2ρuρE],F=[ϕ1uϕ1ρ1uϕ2ρ2uρu2+Pu(ρE+P)],H=[ϕ1ϕ2ρ1a12ρ2a22ϕ1ρ2a22+ϕ2ρ1a120000].
(8)

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

e=Y1e1(ρ1,p)+Y2e2(ρ2,p).
(9)

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ρc2=ϕ1ρ1c12+ϕ2ρ2c22,
(10)

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

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 γ 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 models10,86 provided the additional equation in terms of the specific heat ratio γ, obtained from the weighting of the mass fractions

γ=ϕ1ρ1cp1+ϕ2ρ2cp2ϕ1ρ1cv1+ϕ2ρ2cv2.
(11)

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/(γ1) instead of γ, and this four equation model reads

U=[1γ1ρρuρE],F=[1γ1uρuρu2+Pu(ρE+P)],H=[1γ1000].
(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üneisen EOS.29,31 More recently, Lei and Li87 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 γ 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's2 model resides in the source term H,

U=[ϕ1ϕ1ρ1ϕ2ρ2ρuρE],F=[ϕ1uϕ1ρ1uϕ2ρ2uρu2+Pu(ρE+P)],H=[ϕ10000].
(13)

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

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

U=[ϕ1ϕ1ρ1ϕ2ρ2ρuϕ1ρ1e1ϕ2ρ2e2],F=[ϕ1uϕ1ρ1uϕ2ρ2uρu2+Pϕ1ρ1e1uϕ2ρ2e2u],H=[ϕ1000ϕ1P1ϕ2P2],S=[μc(P1P2)000μcPi(P1P2)μcPi(P1P2)].
(14)

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

ρ(Y1e1+Y2e2+12u2)t+u(ρ(Y1e1+Y2e2+12u2)+(ϕ1p1+ϕ2p2))x=0.
(15)

Defining the acoustic impedance as Zk=ρkck, the interface pressure can be obtained from

Pi=Z2P1+Z1P2Z1+Z2,
(16)

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

c2=Y1c12+Y2c22.
(17)

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

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 model3 [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

Pg=(γg1)ρgeg,
(18)

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

pl=(γl1)ρlelγlπ,
(19)

with a good balance in accuracy and in computational requirements. The model parameters γl and π are determined through the shock wave Hugoniot data.92 The SG model accounts for molecular and repulsive effects in a simplified way, leading to a disagreement between the experimental values and those computed with Eq. (19) for the liquid specific volume. A more accurate model is the noble-able-stiffened-gas model, which accounts for covolume and repulsive effects, and was introduced by Le Métayer and Saurel in Ref. 93 and reads

pk=(γk1)ρk(ekqk)1ρkbkγkp,k.
(20)

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

p=ρ0a02f(η)+ρ0Γ0e,
(21)
f(η)=(η1)(η12γ0(η1))ηs(η1)2,
(22)

with η=ρ/ρ0,Γ0 a material constant, a0 the bulk speed of sound, and ρ0 the reference density. The parameter s represents the Hugoniot slope coefficient, defined as s=dUsdUp, given the shock wave speed Us and the particle velocity Up.

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

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

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

F(t,X)=ỹX=xX.
(23)

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

G=gij=(CL)1=(FFT)1=(FT)1F1,
(24)

which is related to the Almansi tensor εij through

G=gij=δij2εij,
(25)

with δij being the Kronecker symbol. A cobasis vector corresponding to the Lagrangian coordinates is El=Xl and satisfies the following:

FT=(E1,E2,E3),G=l=13ElEl.
(26)

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

ϕgρgt+div(ϕgρgu)=0,ϕsρst+div(ϕsρsu)=0,ρut+div(ρuu(ϕsσs+ϕgσg))=0,ρEt+div((ρE+P)u(ϕsσs+ϕgσg)u)=0,DGeDt+Ge(ux+(ux)T)=LpGe,ϕst+div(ϕsu)=ϕsϕgρsϕs2ρgϕg2ϕsρsag2+ϕsρgag2ux,Ps=Pg=P.
(27)

Here, Ge 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 Lp and the stress tensors for solid σs and gas σg are, respectively,

σs=2ρsesGG,σg=PgI.
(28)

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 Gavrilyuk8 based on the relaxation procedure, which led to the six-equation model of Saurel et al.34 A dissipation function, D̃, is introduced:

D̃=12μcϕ̇s2,ϕ̇s=DϕsDt,
(29)

and the relaxation equation reads

PsPg=D̃ϕṡ=ϕṡμc.
(30)

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

ϕ2ρ2λt+div(ϕ2ρ2uλ)=ϕ2ρ2K,
(31)

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

(ρE)t+div(ρE+p)u=KQ,
(32)

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

ρλt+div(ρ2uλ)=K.
(33)

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

Interim considerations on DI models: Some advantages and shortcomings of different DI models are discussed here, without taking into account the sharpening techniques. There is a large number of applications of the DI models on problems such as the shock–bubble interaction and spherical bubble dynamics capable of testing the properties of the models. The simplest five-equation models of Allaire et al.3 and Massoni et al.116 [Eq. (13)] are known to be very efficient37,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 materials36,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.

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

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

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

Ut+·F(U)=0,
(34)

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

U=[ρρuρvρwE],Fx=[ρuρu2+pρuvρuwu(E+p)],Fy=[ρvρuvρv2+pρvwv(E+p)],Fz=[ρwρuwρvwρw2+pw(E+p)].
(35)

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

ΩUt+ΩF(U)·ndS=0.
(36)

Assuming the domain to be a tessellation of small non-overlapping N elements Vi with volume |Vi|, where i=1,...,N, and applying the integral Euler formulation over a small control volume Vi, the semi-discrete FV scheme reads

dU¯idt|Vi|+ViF(U)·ndS=dU¯idt|Vi|+fVifF(U)·ndS=dU¯idt|Vi|+fViF¯f(U)|f|=0.
(37)

The state variable and face flux integrals have been replaced, respectively, by the volume averaged state variable U¯i and averaged normal flux F¯f,

U¯i=ViUdV|Vi|,F¯f=fF(U)·ndS|f|.
(38)

Hereafter, we will assume a domain Ω that can be any combination of hexahedral, pyramidal, or prismatic elements, with a surface area of |f| for the face f, and with the number of faces Nf. In computational fluid dynamics (CFD) codes, the averaged flux is commonly numerically approximated through a Gaussian quadrature formula.118 The accuracy of the approximation depends upon the selected number of quadrature points Nqp, denoted with xf,α, and their distribution and weights ωα depend upon the Gaussian rule. The approximate face integral is, therefore, of the form

F¯f=α=1NqpωαF(U(xf,α))·nf.
(39)

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

F(U(xf,α))F̃(UL(xf,α),UR(xf,α)).
(40)

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

Uin+1=UinΔt1|Vi|j=1Nfα=1NqpF̃(UL(xf,α),UR(xf,α))ωα|f|.
(41)

High-order space discretization: The ultimate goal in high-order FV schemes is to eliminate spurious oscillations near sharp gradients of the state variables, that is, in proximity of shocks or other sharp discontinuities, and at the same time provide high-order accuracy in smooth flows regions. Independently from the procedure to eliminate nonphysical oscillations, which in the FV frameworks comprehend a variety of methods, such as the total variation diminishing (TVD) techniques and the weighted essentially non-oscillatory schemes, the key ingredient is the reconstruction of a high-order polynomial pi(x,y,z). The way the polynomial pi 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, ui,

ui=1|Vi|Viu(x,y,z)dV=1|Vi|Vipi(x,y,z)dV.
(42)

Within the FV framework, in order to allow such reconstruction, the averages of the target cell Vi will be used as well as the averages um of the neighboring cells constituting the computational stencil S.

The central reconstruction stencil is built up adding recursively the side neighbors of an element Vi, until a desired number M of elements is reached starting from the target element for which m = 0,

S=m=0MVm.
(43)

Since in this work we want to consider unstructured arbitrary meshes, it is convenient sometimes to perform the reconstruction in a uniform reference space.49,57 Considering for example a tetrahedral element with vertices w1=(x1,y1,z1),w2=(x2,y2,z2),w3=(x3,y3,z3),w4=(x4,y4,z4), the transformation from the physical coordinate system (x, y, z) to the reference coordinate system (ξ,η,ζ) reads

(xyz)=(x1y1z1)+J·(ξηζ),
(44)

with J, the Jacobian matrix, defined as

J=[x2x1x3x1x4x1y2y1y3y1y4y1z2z1z3z1z4z1].
(45)

This transformation defines both the direct x=x(ξ) and inverse ξ=ξ(x) mappings from ξ=(ξ,η,ζ) into x=(x,y,z). Applying the inverse mapping to all the elements Vm of the reconstruction stencil S, the transformed elements Vm will constitute the transformed stencil S,

S=m=0MVm.
(46)

Over this reference system, the reconstruction polynomial is sought as expansion over local polynomial basis functions λk(ξ,η,ζ), leading to

p(ξ,η,ζ)=U0+k=1Kakλk(ξ,η,ζ),
(47)

defined U0 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

K=12(r+1)(r+2)(r+3)1.
(48)

In order to determine the degrees of freedom ak, one can impose that for each cell Vm of the stencil S, the cell average of the reconstruction polynomial to be equal to the cell average of the solution Um,

Vmp(ξ,η,ζ)dξdηdζ=U0|Vm|+k=1Kvmakλk(ξ,η,ζ)dξdηdζ=|Vm|Um.
(49)

An expression for the DOFs can be drawn in matrix form

k=1KAmkak=bm,m=1,2,,M,
(50)

where Amk represents the integral of the basis function over the cell m in the stencil and bm the right hand side vector, respectively,

Amk=Vmϕkdξdηdζ,|Vm|(UmU0)=bm.
(51)

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 Ω, and performing an integration by parts. This results in

Ωϕ(x,y,z)UtdΩ+Ωϕ(x,y,z)F(U)·ndS=Ωϕ(x,y,z)·F(U)dΩ.
(52)

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 space119 Vh={ϕhL:ϕh|ΩVk(Ω),k=0,1,2,,N}, where Vk is the space of polynomials of degree up to k. The discrete solution Uh, with expansion coefficients denoted by uh, can be seen as expansions over a finite-element basis Pjk in the aforementioned polynomial space, where d is the number of degrees of freedom

Uh(x,y,z,t)=j=1duh(t)Pjk(x,y,z).
(53)

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

ddtEϕhuhdΩ+fEfϕhF(uh)·ndS=Eϕh·F(uh)dΩ.
(54)

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

ϕh(x,y,z)=j=1dPjk(x,y,z).
(55)

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,

ddtEPikPjkuhdΩ+fEfPjkF(uh)·ndS=EPjk·F(uh)dΩ..
(56)

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

fPjkF(uh)·ndSl=1LψlF(uh(xf,l))·nPjk(xf,l)|f|,
(57)

and the volume integral with

E(Pjk)·F(uh)dΩm=1MωjF(uh(xE,m))·Pjk(xE,m)|E|,
(58)

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

ddtEPikPjkuhdΩ=m=1MωjF(uh(xE,m))·Pjk(xE,m)|E|fEl=1LψlF̃(uhL(xf,l),uhR(xf,l))·nPjk(xf,l)|f|.
(59)

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

MdUdt+R(U)=0.
(60)

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

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

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

d=(n+kk)=(n+k)!n!k!.
(61)

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

Uh=ŨB1+Ux|cΔxB2+Uy|cΔyB3+2Ux2|cΔx2B4+2Uy2|cΔy2B5+2Uxy|cΔxΔyB6,
(62)

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

B1=1,B2=xxcΔx,B3=yycΔy,B4=(xxc)22Δx21Ωeωe(xxc)2Δx2dΩ,B5=(yyc)22Δy21Ωeωe(yyc)2Δy2dΩ,B6=(xxc)(yyc)(Δx)(Δy)1Ωeωe(xxc)(yyc)ΔxΔydΩ.
(63)

The features of the FV and the DG schemes were discussed in Secs. III A and III B. The advantages of the DG formulation over the FV in higher-order applications are well known.79,118 However, it was also shown the extreme increase in computer operations when third and higher polynomials are used in three dimensions. Indeed, for the same order of accuracy, the DG methods consist of a system of equations with a higher number of unknowns compared to the FV discretization. An interesting approach, aiming to decrease the computational cost of the DG method, is the family of reconstructed DG (rDG) method, proposed by Dumbser et al. in Refs. 57 and 123. This procedure shares some similarities with the work of van Leer et al.,124,125 where a recovery method is built to compute the diffusive fluxes within the DG scheme. The recovery is to be considered as a L2 projection, which provides a polynomial solution, that is, indistinguishable in the weak sense from the underlying discontinuous solution. In the procedure proposed by Dumbser et al.,57 named PNPM schemes, the solution is initially represented by a polynomial PN of degree N. The time evolution and the intercell flux evaluation are performed with a polynomial PM, reconstructed from the underlying polynomial PN. 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 (Φl) and the hierarchical orthogonal reconstruction basis Ψl are chosen in such a way to satisfy the following properties:

Φl=Ψlfor1lLN,[Ψm,Ψm]=0m,n,withmn.
(64)

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

Ui=Ũi+UxiB2+UyiB3.
(65)

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

UiR=ŨiR+UxiRB2+UyiRB3+UxxiRB4+UyyiRB5+UxyiRB6,
(66)

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

ŨiR=Ũi,ŨxiR=Ũxi,ŨyiR=Ũyi,
(67)

the final requirement is, considering a neighboring cell j,

Uj=Ũi+UxiB2+UyiB3+UxxiRB4+UyyiRB5+UxyiRB6,Uxj=Uxi+UxxiRΔxiB2(i)+UxyiRΔyiB3(i),Uyj=Uyi+UyyiRΔyiB3(i)+UxyiRΔxiB2(i).
(68)

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

Ux=Ux|cΔx,Uy=Uy|cΔy,Uxx=2Ux2|cΔx2,Uyy=2Uy2|cΔy2,Uxy=2Uy|cΔxΔy.
(69)

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

(B4(i)B5(i)B6(i)B2(i)0B3(i)0B3(i)B2(i))(UxxiRΔxi2UyyiRΔyi2UxyiRΔxiΔyi)=(Uj(Ũi+UxiΔxiB2(i)+UyiΔyiB3(i)Δxi(UxjUxi)Δyi(UyjUyi)).
(70)

This method is indeed more efficient than the original PNPM scheme. However, since it is required that the reconstructed solution has to be equal to the underlying DG solution of the neighboring cell at the centroid, the procedure is not 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 Liu129 proposed a reconstructed DG-FV scheme in which the higher order degrees of freedom are reconstructed by minimizing the interfacial jump integration, as first introduced by Wang et al.130 for the FV framework, aiming to increase the efficiency for steady-state problems.

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

FIG. 1.

Convergence order comparison for second-order discontinuous Galerkin DG(P1), third-order green Gauss reconstructed DG P1P2(GG), least-square reconstructed DG P1P2(LSr), hybrid least-square reconstructed DG P1P2(HLSr), and third-order DG(P2), for subsonic flow past a cylinder test case. [Adapted with permission from Cheng et al., “A hybrid reconstructed discontinuous Galerkin method for compressible flows on arbitrary grids,” Comput. Fluids 139, 68–79 (2016). Copyright 2016 Elsevier Ltd.]

FIG. 1.

Convergence order comparison for second-order discontinuous Galerkin DG(P1), third-order green Gauss reconstructed DG P1P2(GG), least-square reconstructed DG P1P2(LSr), hybrid least-square reconstructed DG P1P2(HLSr), and third-order DG(P2), for subsonic flow past a cylinder test case. [Adapted with permission from Cheng et al., “A hybrid reconstructed discontinuous Galerkin method for compressible flows on arbitrary grids,” Comput. Fluids 139, 68–79 (2016). Copyright 2016 Elsevier Ltd.]

Close modal

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

TV(Un)=i=|UinUi1n|,
(71)

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

TV(Un+1)TV(Un).
(72)

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

Uij,α=Ui+Θi·k=1Kakλk(ξ,η,ζ).
(73)

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

The slope limiter Θil,α is commonly expressed as

Θil,α={min(1,UimaxUiUil,αUi)ifUil,αUi>0,min(1,UiminUiUil,αUi)ifUil,αUi<0,1ifUil,αUi=0.
(74)

The maximum Uimax and the minimum Uimin values of the state variable over the stencil are required, and the limiter Θi selects the minimum value from the slope limiter Θil,α from all the side neighbors, and quadrature points

Θi=min(Θil,α),l[1,L],α[1,Nqp].
(75)

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

pweno=m=0msωmpm(ξ,η,ζ),
(76)

given the total number of stencils ms. Substituting the expression for the single polynomial Eq. (47) in Eq. (76), one obtains the expanded expression for the WENO polynomial

pweno=u¯0+k=1K(m=0msωmak(m))λk(ξ,η,ζ),
(77)

and the non-linear weights ωm are defined as in the classical formulation in Refs. 46 and 133,

ωm=γmm=0msγm,γm=dm(ε+ISm)p.
(78)

The linear weights dm take a large value for the central stencil (102<d0<105), although in Ref. 134 it was found that higher values are better for smooth solutions and lower values works better on discontinuities, and smaller linear weights are assigned to the sectorial stencil (dm = 1). The parameters p and ε are responsible to control the non-linear weights decay for non-smooth stencils and to avoid division by zero, respectively, and typically, the values p = 4 and ε=106 are used. The smoothness indicator ISm, as the name suggests, indicates how the solution is smooth on the stencil, and is calculated with

ISm=1<|β|<rV0(Dβpm(ξ,η,ζ))2dξdηdζ,
(79)

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

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

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

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

ft+c·f=Ω.
(80)

The variable c denotes the particle velocity, and Ω is the collision operator, and as originally proposed by Bathnagar, Gross, and Krook,146 if f is close to its equilibrium state f eq, it can be approximated by a relaxation mechanism, ruled by a relaxation time τ dependent on fluid properties like the viscosity

Ω(f)=ffeqτ.
(81)

The resulting discrete velocity Boltzmann equation (DVBE) is a set of PDEs, with a number of discrete velocities ci,

fit+ci·f=fifieqτ.
(82)

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

fi(x+ciΔx,t+Δt)=fi(x,t)1τtt+Δt(fifieq)(x,t)dt.
(83)

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,

ρi=0m1fi=i=0m1fieq,ρuii=0m1cifi=i=0m1cifieq.
(84)

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

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

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

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

ft+c·f=ffeqτf,ht+c·h=hheqτh1τhf(h),
(85)

and the quantities h=c·u12(ffeq) and 1τhf=1τh1τf are defined from Guo's model.156 If one considers, the following density feq and energy heq equilibrium distributions, respectively:

feq=ρ(12πRT)D/2exp((cu)22RT),heq=ρc22(12πRT)D/2exp((cu)22RT),
(86)

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

f=ω(c,T)n=01n!a(n)H(n)(c),h=ω(c,T)n=01n!b(n)H(n)(c).
(87)

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 ω(c,T) is defined as

ω(c,T)=1(2π)D/2ec2/2.
(88)

The integration from time tn to time tn+1 leads to the so-called discrete unified gas-kinetic scheme (DUGKS).158 The feature is that like FV schemes, this can be easily applied to unstructured meshes

f(x,tn+1)=f(x,tn)ΔtVFn+12+Δt2(Ωn+1+Ωn),h(x,tn+1)=h(x,tn)ΔtVHn+12+Δt2(Ωn+1+Ωn),
(89)

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

F(n+12)=V(c·n)f(x,tn+12)dS,H(n+12)=V(c·n)h(x,tn+12)dS,
(90)

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

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

(cxk,cyk)={(0,0)for k=1,v1(cosπk2,sinπk2)for k=2,3,4,5,v2(cosπ(k+0.5)2,sinπ(k+0.5)2)for k=6,7,8,9,
(91)

with v1v20 non-zero constants, and the local equilibrium particle velocity distribution function is given by

fkeq=ρ[Ak+Bk(uxcxk+uycyk)2+Ck(uxcxk+uycyk)2],
(92)

with k[1,9]. We refer to Refs. 159 and 77 for the determination of the coefficients Ak, Bk, and Ck. 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.

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

  • Discrete conservation of mass, momentum and energy

  • High-order accuracy solution of smooth flow regions

  • Oscillation-free representation of interfaces and discontinuities.

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

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

In both Refs. 4 and 53, an Harten-Lax-van Leer contact (HLLC) approximate Riemann solver160 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:

α1t+·(α1u)=α1·u,
(93)

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

fHLLC=1+sgn(s*)2[fLa+s(q*LqL)]+1+sgn(s*)2[fRa+s(q*RqR)],
(94)

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

q*K=(sKuKsKs*)((α1ρ1)K(α2ρ2)KρKs*ρKvKEK+(s*uK)(ρKs*+pKsKuK)α1K).
(95)

The wave speed in the star region is calculated by

s*=pRpL+ρLuL(sLuL)ρRur(sruR)ρL(sLuL)ρR(sRuR),
(96)

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

s=min(0,sL),s+=max(0,sR),
(97)

with

sL=min(u¯c¯,uLcL),sR=min(u¯+c¯,uL+cL).
(98)

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

uHLLC=1+sgn(s*)2[uL+s(sLuLsLs*1)]+1+sgn(s*)2[uR+s(sRuRsRs*1)].
(99)

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

A recent investigation on the performances of high-resolution methods for shock–turbulence interaction was conducted by Johnsen et al.161 Although capable of sharp shock capturing, the well-known WENO scheme of Ref. 46 was reported to dissipate turbulent fluctuations for the higher resolvable wavenumbers and results in a degradation of spectral properties as already pointed out by Pirozzoli167 because of the upwind biased nonlinear reconstruction. In the Finite Difference context, several workarounds have been proposed, for example in form of improved nonlinear weights evaluation, by integrating the smoothness indicators Eq. (79) on a full stencil and using much smaller values for the parameter ε, as in the WENO-Z168 and in the WENO-M169 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 Lele170 and features enhanced small-scale waves resolution and spectral properties, and was adapted to flows with discontinuities by Deng and Zhang171 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-SYMBO165 and MDCD-WENO166 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-JS46 if 1020 cells per wavelength are used to compute the highest wavenumber component.

TABLE I.

Tabular with multi-species advection test164 convergence order comparison for MUSCL reconstruction, Jiang and Shu WENO (WENO-JS), bandwidth-optimized WENO165 (WENO-SYMBO), and MDCD-WENO.166 [Adapted with permission from Wang et al., “Consistent high resolution interface-capturing finite volume method for compressible multi-material flows,” Comput. Fluids 202, 104518 (2020). Copyright 2020 Elsevier Ltd.]

 MUSCLWENO-JSWENO-SYMBOMDCD-WENO
Number of cellsErrorOrderErrorOrderErrorOrderErrorOrder
7.587 × 10−02   7.012 × 10−02   6.332 × 10−02   3.688 × 10−02   
6.135 × 10−02 0.307 1.174 × 10−02 2.579 3.489 × 10−03 4.182 1.743 × 10−03 4.403 
16 2.657 × 10−02 1.207 7.089 × 10−04 4.049 4.218 × 10−04 3.048 1.614 × 10−04 3.432 
32 1.074 × 10−02 1.307 2.644 × 10−05 4.745 1.905 × 10−05 4.468 1.105 × 10−05 3.869 
64 3.658 × 10−03 1.553 8.631 × 10−07 4.937 1.121 × 10−06 4.087 6.915 × 10−07 3.998 
128 1.153 × 10−03 1.666 3.223 × 10−08 4.743 5.541 × 10−08 4.339 2.814 × 10−08 4.619 
 MUSCLWENO-JSWENO-SYMBOMDCD-WENO
Number of cellsErrorOrderErrorOrderErrorOrderErrorOrder
7.587 × 10−02   7.012 × 10−02   6.332 × 10−02   3.688 × 10−02   
6.135 × 10−02 0.307 1.174 × 10−02 2.579 3.489 × 10−03 4.182 1.743 × 10−03 4.403 
16 2.657 × 10−02 1.207 7.089 × 10−04 4.049 4.218 × 10−04 3.048 1.614 × 10−04 3.432 
32 1.074 × 10−02 1.307 2.644 × 10−05 4.745 1.905 × 10−05 4.468 1.105 × 10−05 3.869 
64 3.658 × 10−03 1.553 8.631 × 10−07 4.937 1.121 × 10−06 4.087 6.915 × 10−07 3.998 
128 1.153 × 10−03 1.666 3.223 × 10−08 4.743 5.541 × 10−08 4.339 2.814 × 10−08 4.619 

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 Kreiss174 in the context of level set methods. The advection for the interface function is therefore rewritten as

ϕt+u·ϕ=U0n·(εh|ϕ|ϕ(1ϕ)),
(100)

once defined n as the interface normal n=ϕ/|ϕ|,U as a characteristic compression velocity and εh as a grid spacing length scale. Since we are interested in sharpening the interface, small values of εh and large values of U0 are the target. Under the consideration that also the density gradients will diffuse in time, a similar correction is required for the density equation, with the only difference that the density correction has to be insensitive to the density values close to the interface, since great density ratio may be involved:

ρt+·(ρu)=H(ϕ)U0n·((εhn·ρ)(12ϕ)ρ).
(101)

The function H(ϕ) is introduced to localize the compression region and may be taken as39 

H(ϕ)=tanh[(ϕ(1ϕ)102)].
(102)

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

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

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

ϕτ=·(Dϕ),
(103)

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

D=|sLsR(ϕRϕL)sRsL|.
(104)

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:

TOLi|(·(ϕ¯)lim¯)i|Vii|((α¯)CD)i|2Vi.
(105)

This compares the gradient of the limited volume fraction (ϕ¯)lim with the gradient calculated with a central difference scheme (α¯)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 ρ across the interface is slow. In practice, to the Kapila's system [Eq. (8)], another source S(U) is added that reads

S(U)=[R̂1R̂2uR̂κR̂+(p(Γ2Γ1)+(Π2Π1)RR].
(106)

The regularization operator, in this case, reads

R(ϕ2)=L(ϕ2)U0n·(εh|ϕ|ϕ(1ϕ)),
(107)

with

L(ϕ2)={1for106<ϕ2<1106,0otherwise.
(108)

The other regularization terms are defined as

R̂1,2(ρ1,2α1,2)L(ϕ1,2)U0n·(εhn·(ρ1,2ϕ1,2))(12ϕ1,2)(ρ1,2ϕ1,2)),
(109)

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

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

FIG. 2.

Example of bubble shapes for the isolated single-bubble collapse test with pressure ratio of 10 at time t=110μs, obtained with (a) the second-order min-mod limiter, and (b) with a fifth-order WENO. [Adapted with permission from Tiwari et al., “A diffuse interface model with immiscibility preservation,” J. Comput. Phys. 252, 290–309 (2013). Copyright 2013 Elsevier Ltd.]

FIG. 2.

Example of bubble shapes for the isolated single-bubble collapse test with pressure ratio of 10 at time t=110μs, obtained with (a) the second-order min-mod limiter, and (b) with a fifth-order WENO. [Adapted with permission from Tiwari et al., “A diffuse interface model with immiscibility preservation,” J. Comput. Phys. 252, 290–309 (2013). Copyright 2013 Elsevier Ltd.]

Close modal

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 Φ is used as an indicator of the presence of an interface within a cell Ci, with an associated cell average Φ¯i, and which takes unity or zero value whether Ci is filled with a fluid a or not. For the intermediate values, the hyperbolic tangent function Φ̃i(x) mimics the jump across an interface

Φ̃i(x)=12[1+σitanh(β(xxi1/2Δxix̃i))],
(110)

given the sign function σi=sgn(ΔΦ¯i) of the jump ΔΦ¯i=Φi+1(t)Φi1(t) and a user-defined slope parameter β. Assuming conservation of volume fraction in the form

1Δxixi1/2xi+1/2Φ̃i(x)dx=Φ¯i(t),
(111)

the interface location x̃i can be determined from Eq. (110) as it is the only unknown,

x̃i=12βln[exp(β(1+σi2Φ¯i)/σi)1exp(β(1σiΦ¯i)/σi].
(112)

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

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. 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.]

Close modal

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.

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

Θ̃il,α={WimaxWi2(Wil,αWi)if(Wil,αWi)>0,WiminWi2(Wil,αWi)if(Wil,αWi)<0,1if(Wil,αWi)=0,
(113)

which allows to recover the minmod and superbee limiters from

Θil,α=max[0,min(βΘ̃il,α,1),min(Θ̃il,α,β],
(114)

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

Θil,α=max[0,min[2,2Θ̃il,α,t]],
(115)

where the term t is defined as

t=max[min(2Θ̃il,α,β),min{(2β)Θ̃il,α+2(β1),Θ̃il,α}].
(116)

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

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

BVijξ=|qij,Lξqij,Rξ|,
(117)

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

qij,Lξ=ΓijQiξ(x,y)dΓ,qij,Rξ=ΓijQijξ(x,y)dΓ.
(118)

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

TBViξ=j=1JBVijξ,
(119)

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

Ξ={Qiξ1(x,y),Qiξ2(x,y),,QiξN(x,y)}.
(120)

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

popt(ξ,η,ζ)=s=1Stλsps(ξ,η,ζ),
(121)

summing over the stencils 1<sSt with linear coefficients λs. The polynomial p1 is computed subtracting the lower-order polynomials from the optimal polynomial

p1(ξ,η,ζ)=1λ1(popt(ξ,η,ζ)s=2Stλsps(ξ,η,ζ)),
(122)

and the reconstruction polynomial is a non-linear combination of the polynomials ps as defined in Eq. (76). The central stencil λ1 is calculated by normalizing an arbitrary value assigned to the non-normalized weight λ, and subsequently, the lower-order polynomial weights are obtained from the central one:

λ1=11λ,λs=1λ1St1.
(123)

The multi-species convergence study defined in the work of Wong and Lele164 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).

FIG. 4.

Example of 2D unstructured mesh with arbitrary elements used in Ref. 178 in conjunction with CWENO reconstruction. [Reproduced with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.]

FIG. 4.

Example of 2D unstructured mesh with arbitrary elements used in Ref. 178 in conjunction with CWENO reconstruction. [Reproduced with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.]

Close modal
TABLE II.

Tabular with eL and eL2 errors achieved for the multi-species advection test of Wong and Lele.164 Convergence order and wall clock time shows the efficiency properties of the CWENO scheme over the classic WENO scheme. [Reproduced with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.]

 CWENOWENO
Order/edges per sideeLOLeL2OLTimeeLOLeL2OLTime
3rd/10 1.124 × 10−02 ⋯ 6.897 × 10−03 ⋯ 1.000 1.119 × 10−02 ⋯ 6.880 × 10−03 ⋯ 1.252 
3rd/20 1.514 × 10−03 2.89 9.535 × 10−04 2.85 8.240 1.506 × 10−03 2.89 9.492 × 10−04 2.86 8.757 
3rd/40 1.993 × 10−04 2.93 1.244 × 10−04 2.94 85.831 1.982 × 10−04 2.93 1.236 × 10−04 2.94 94.658 
3rd/80 2.487 × 10−05 3.00 1.581 × 10−05 2.98 673.159 2.473 × 10−05 3.00 1.571 × 10−05 2.98 798.829 
4th/10 1.879 × 10−03 ⋯ 5.355 × 10−04 ⋯ 1.530 1.867 × 10−03 ⋯ 5.342 × 10−04 ⋯ 1.821 
4th/20 7.487 × 10−05 4.65 2.762 × 10−05 4.28 10.864 7.427 × 10−05 4.65 2.739 × 10−05 4.29 14.081 
4th/40 4.962 × 10−06 3.92 1.857 × 10−06 3.89 115.068 4.912 × 10−06 3.92 1.839 × 10−06 3.90 144.007 
4th/80 3.482 × 10−07 3.83 1.200 × 10−07 3.95 928.223 3.451 × 10−07 3.83 1.188 × 10−07 3.95 1128.293 
5th/10 9.115 × 10−04 ⋯ 4.777 × 10−04 ⋯ 1.650 9.101 × 10−04 ⋯ 4.775 × 10−04 ⋯ 2.601 
5th/20 2.994 × 10−05 4.93 1.736 × 10−05 4.78 13.055 2.997 × 10−05 4.92 1.740 × 10−05 4.78 19.929 
5th/40 1.084 × 10−06 4.79 5.816 × 10−07 4.90 139.432 1.088 × 10−06 4.78 5.833 × 10−07 4.90 211.149 
5th/80 3.327 × 10−08 5.03 1.883 × 10−08 4.95 1119.209 3.335 × 10−08 5.03 1.889 × 10−08 4.95 1738.832 
6th/10 1.426 × 10−04 ⋯ 4.309 × 10−05 ⋯ 2.227 1.414 × 10−04 ⋯ 4.321 × 10−05 ⋯ 3.704 
6th/20 2.447 × 10−06 5.86 7.587 × 10−07 5.83 16.729 2.467 × 10−06 5.84 7.637 × 10−07 5.82 30.602 
6th/40 3.537 × 10−08 6.11 1.207 × 10−08 5.97 175.908 3.536 × 10−08 6.12 1.212 × 10−08 5.98 304.109 
6th/80 6.493 × 10−10 5.77 2.046 × 10−10 5.88 1493.933 6.572 × 10−10 5.75 2.055 × 10−10 5.88 2417.527 
 CWENOWENO
Order/edges per sideeLOLeL2OLTimeeLOLeL2OLTime
3rd/10 1.124 × 10−02 ⋯ 6.897 × 10−03 ⋯ 1.000 1.119 × 10−02 ⋯ 6.880 × 10−03 ⋯ 1.252 
3rd/20 1.514 × 10−03 2.89 9.535 × 10−04 2.85 8.240 1.506 × 10−03 2.89 9.492 × 10−04 2.86 8.757 
3rd/40 1.993 × 10−04 2.93 1.244 × 10−04 2.94 85.831 1.982 × 10−04 2.93 1.236 × 10−04 2.94 94.658 
3rd/80 2.487 × 10−05 3.00 1.581 × 10−05 2.98 673.159 2.473 × 10−05 3.00 1.571 × 10−05 2.98 798.829 
4th/10 1.879 × 10−03 ⋯ 5.355 × 10−04 ⋯ 1.530 1.867 × 10−03 ⋯ 5.342 × 10−04 ⋯ 1.821 
4th/20 7.487 × 10−05 4.65 2.762 × 10−05 4.28 10.864 7.427 × 10−05 4.65 2.739 × 10−05 4.29 14.081 
4th/40 4.962 × 10−06 3.92 1.857 × 10−06 3.89 115.068 4.912 × 10−06 3.92 1.839 × 10−06 3.90 144.007 
4th/80 3.482 × 10−07 3.83 1.200 × 10−07 3.95 928.223 3.451 × 10−07 3.83 1.188 × 10−07 3.95 1128.293 
5th/10 9.115 × 10−04 ⋯ 4.777 × 10−04 ⋯ 1.650 9.101 × 10−04 ⋯ 4.775 × 10−04 ⋯ 2.601 
5th/20 2.994 × 10−05 4.93 1.736 × 10−05 4.78 13.055 2.997 × 10−05 4.92 1.740 × 10−05 4.78 19.929 
5th/40 1.084 × 10−06 4.79 5.816 × 10−07 4.90 139.432 1.088 × 10−06 4.78 5.833 × 10−07 4.90 211.149 
5th/80 3.327 × 10−08 5.03 1.883 × 10−08 4.95 1119.209 3.335 × 10−08 5.03 1.889 × 10−08 4.95 1738.832 
6th/10 1.426 × 10−04 ⋯ 4.309 × 10−05 ⋯ 2.227 1.414 × 10−04 ⋯ 4.321 × 10−05 ⋯ 3.704 
6th/20 2.447 × 10−06 5.86 7.587 × 10−07 5.83 16.729 2.467 × 10−06 5.84 7.637 × 10−07 5.82 30.602 
6th/40 3.537 × 10−08 6.11 1.207 × 10−08 5.97 175.908 3.536 × 10−08 6.12 1.212 × 10−08 5.98 304.109 
6th/80 6.493 × 10−10 5.77 2.046 × 10−10 5.88 1493.933 6.572 × 10−10 5.75 2.055 × 10−10 5.88 2417.527 

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

Ut+Gx+HUx=0,
(124)

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

ΩϕiUtdx=ΩdϕidxGdxΩϕiHUxdx[ϕiĜ12(ϕiL+ϕiR)HUx]xk1/2xk+1/2.
(125)

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

[HUx]xk+1/2=01H(ψ(τ;UL,UR))ψτ(τ;UL,VR)dτ.
(126)

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

Another extension of the DG method to handle spatial derivatives higher than the first-order and non-conservative products is the local DG (LDG), introduced in Ref. 55, used in Ref. 184 for solving the Hamilton–Jacobi equation and extended to the DI models by Gryngarten and Menon.185 The aim of the LDG method is to rewrite the equations as a first-order system by introducing new variables, which approximate the solution derivatives and is motivated by the work of Bassi and Rebay186,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

{Ut+·F(U)=H(U)·u,q=f(U),
(127)

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αm and qαp, and is related to the source term through

q=(qαmqαp)=(α1α1).
(128)

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

qj=j=0pci,jPjk,
(129)

and multiplying Eq. (127) by a test function ϕ and integrating over the domain Ωi, one obtains the following two equations:

ΩjϕUjtdVΩjϕ·F(Uj)dV+ΩjϕF̂(U,U+)dS=ΩjϕH(Uj)dV,ΩjϕqjdV+Ωjϕf(Uj)dVΩjϕf̂(U,U+)dS=0.
(130)

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αm and the right side is used for qαp. Then, the source Sα can be approximated by the combination

Ŝα=12[(uuM)qαp+(u+uM)qαm],
(131)

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

In the work of Saleem et al.,188 the DG method is used for solving the five-equation model of Kreeft and Koren.189 This is a velocity and pressure equilibrium model with a source term, which expresses the exchange of energy between the two fluids in terms of mechanical L̇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

L̇=L̇M+L̇T=pu·ϕ1+(ϕ1ϕ1ρ1ρ)u·p+pϕ1(1ϕ1)τ2τ1τ·u,
(132)

where τ1 and τ2 are the isentropic compressibilities of the two phases

τ1=1ρ1(ρ1p)s1,τ2=1ρ2(ρ2p)s2.
(133)

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

τ=ϕ1τ1+ϕ2τ2,ρ=ϕ1ρ1+ϕ2ρ2.
(134)

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

U=[ϕ1ρ1Eϕ1ρ1ϕ2ρ2ρuρE],F=[ϕ1u(ρ1E+P)ϕ1ρ1uϕ2ρ2uρu2+Pu(ρE+P)],H=[L̇0000].
(135)

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

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

A recent attempt to extend the DG discretization for the solution of a γ-based model of the form proposed in Ref. 29 can be found in Ref. 193. This new formulation differs from the path-conservative formulation of Ref. 183 in that it reduces to the quasi-conservative FV discretization of Ref. 15 if a P0 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.

The reconstructed DG method demonstrated to be an accurate and efficient method for the solution of hyperbolic PDEs, especially under a careful selection of the polynomial basis and reconstruction procedure. Recently, Pandare et al.63 successfully solved a one-dimensional compressible multi-material problem with a high-order reconstructed DG method. The multi-species system of equation used is similar to Kapila's2 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

Sk=1τ(pkk(pkϕkρkak2)kϕkρkak2)ϕkρkak2,
(136)

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

τ=cτmaxk(hak).
(137)

Small cτ 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τ=ρmin/ρmax.

Following Luo et al.,120 the DG formulation with a Taylor basis polynomial of degree PN is employed for the solution discretization at the beginning of a time step and used to reconstruct a polynomial of degree PM before flux evaluation. Therefore, a third-order rDG(P1P2) solution can be recovered from a first order underlying solution P1, 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 P1 moments and to a WENO limiter with TVB minmod troubled cell detector for the P2 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

x(ϕkρkEk)|i=ρkEki¯ϕkx|i,
(138)

having the cell-average value of the specific total energy of the material kρkEk¯. 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).

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.]

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.]

Close modal

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 χ=γγ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),

U=[ρρuEρMχρ],F=[ρuρu2+Pu(ρE+P)ρuMχρu].
(139)

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 Maheshwari155 presented a similar approach but corrected an inconsistency of the Wang model and substituted the KT model used by Joshi et al. with the multi-entropy-level LBM of Yan and Zhang.198 The inconsistency found in Ref. 197 lies in the derivation of Eq. (139). Following the expansion of the equation for the mixture specific heat:

(1γ1){(ργ)t+(ργu)x}+(ργ){(11γ)t+u(11γ)x}=0,
(140)

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 (ργ) and the rest of the variables are first evolved in time before solving for 1γ1. Therefore, Eq. (139) is modified to be consistent with the predictor–corrector steps

U=[ρρuEρMγρ],F=[ρuρu2+Pu(ρE+P)ρuMγρu].
(141)

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

Fi1/2R=Fi1/2R++Fi1/2R,Fi1/2L=Fi1/2L++Fi1/2L,
(142)

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:

Fi1/2R+=F(i1/2),2R+F(i1/2),4R,Fi1/2R=F(i1/2),1R+F(i1/2),3R.
(143)

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

Fi1/2=Fi1/2R++Fi1/2L.
(144)

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

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

Fast transient fluid–structure dynamics was also studied by Faucher et al.201,202 using the VOFIRE (finite-volume with reconstruction) algorithm. The VOFIRE is an anti-diffusive scheme introduced by Després et al.,203 where the dissipation is controlled with a constrained downwind scheme. In illustrating the basic idea behind the VOFIRE scheme, let us define the grid coordinates for the cell center xj=(j+1/2)Δx and cell interface xj+1/2=(xj+1+xj)/2, with a uniform space interval Δx. Assuming that at a certain time instant tn an interface lies in a cell j, then in order to build the fluxes for the time update, only the discrete values hj+1/2 at the cell boundary xj+1/2 are necessary. Then, the constraints for the downwind scheme are the condition for flux consistency,

mj+1/2=min(hh,hj+1)hj+12max(hj,hj+1)=Mj+12,
(145)

and the condition for stability through advection,

mj12=min(hj1,hj))hj*max(hj1,hj)=Mj12,
(146)

with hj*=hj+uΔtΔx(hj1/2hj+1/2). The intervals for which hj+1/2 meets both stability and consistency criteria are

mj+12hj+1/2Mj+12,Mj12+ΔxuΔt(mj12hj)hj+12mj12+ΔxuΔt(mj12hj),
(147)

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

hj,k=hjn+rNpj,rμj,k,r(hknhjn),μj,k,r[0,1],
(148)

including all the incoming faces r. The upwind and the downwind values for the outgoing cell face k are referred to as hjn and hkn, respectively. Therefore, hj,kn takes the downwind value when the coefficients μ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 Nikiforakis110 and discussed in Sec. II can be coupled to a elastoplastic system of equations as shown in Ref. 9 and also to a third system of equations such the resistive MHD model.98 From a numerical standpoint, the coupling between the systems uses a diffuse interface approach, but the location of the material interfaces is retrieved with a level-set method. However, the evolution at the interface is accounted with a Riemann ghost fluid method following the procedure introduced by Sambasivan and Udaykumar.204,205 In this procedure, a band of ghost points is defined in the area adjacent to the interface and the input states at both sides for the mixed-material Riemann solver are taken at a distance 1.5Δx from the interface coordinate in order to avoid potential errors carried by the interface. The Riemann solver subsequently uses these characteristic states to return the well-known left and right star regions, similarly to Eq. (95), which in practice corresponds to the regions surrounding the material interface. The ghost cells are then updated with these values.

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

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 Sturtevant206 and numerically by Quirk and Karni207 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–Meshkov207,208 instability mechanism. The vorticity direction depends upon the Atwood number, defined as

At=ρgρairρg+ρair,
(149)

where the density of Helium and R22 has to be considered in place of ρg. For the case of Helium gas bubble, assuming a shock incident from the right, the vorticity is clockwise because the Helium is lighter than air (negative Atwood number), and the vortex roll-up move the bubble structure downstream to the surrounding fluid. The opposite mechanism happens for the R22 bubble as the higher density leads to positive Atwood numbers, and the counterclockwise vortices move the bubble structure upstream.209 In both the experiments and the computations, the Mach number for the shock wave is took as Ms=1.22, and the bubble is in initial equilibrium with the surrounding air. In the initial phases of the shock–bubble interaction, complex wave structures involving refraction, reflection, and diffraction develop, providing a challenging context for compressible solvers and shock capturing schemes, whereas higher order schemes would be beneficial for a correct evaluation of the instabilities.

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

FIG. 6.

Computational domain for the Helium bubble–shock wave test case as adopted in Ref. 200. [Reproduced with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.]

FIG. 6.

Computational domain for the Helium bubble–shock wave test case as adopted in Ref. 200. [Reproduced with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.]

Close modal

The initial conditions are defined as follows:

(ϕ1ρ1,ϕ2ρ2,u,v,p,ϕ1)={(0,1.204,0,0,101325,0)forPre-Shock,(0,1.658,114.49,0,159060,0)forPost-Shock,(0.158,0.061,0,0,101325,0.95)Insidebubble,
(150)

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-WENO165 (WENO-SYMBO) and MDCD-WENO of166 on uniform meshes, the CWENO on unstructured meshes of Ref. 200 and the hybrid LB-FV schemes of Ref. 155 in the variants with MUSCL limiter and WENO reconstruction (LAMBS-MUSCL and LAMBS-WENO), that are in good agreement with the experimental evidence of Ref. 206 and early computational results of Ref. 207.

FIG. 7.

Bubble interface position for different schemes: (a) LAMBS-MUSCL and LAMBS-WENO, Adapted with permission from Sashi Kumar and Maheshwari, “A lattice Boltzmann method for simulation of multi-species shock accelerated flows,” Int. J. Comput. Fluid Dyn. 32, 19–34 (2018). Copyright 2018 Taylor and Francis Online (b) WENO-JS, WENO-SYMBO, MDCD-WENO Adapted with permission from Wang et al., “Consistent high resolution interface-capturing finite volume method for compressible multi-material flows,” Comput. Fluids 202, 104518 (2020). Copyright 2020 Elsevier Ltd. and (c) 2D and 3D CWENO adapted with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

FIG. 7.

Bubble interface position for different schemes: (a) LAMBS-MUSCL and LAMBS-WENO, Adapted with permission from Sashi Kumar and Maheshwari, “A lattice Boltzmann method for simulation of multi-species shock accelerated flows,” Int. J. Comput. Fluid Dyn. 32, 19–34 (2018). Copyright 2018 Taylor and Francis Online (b) WENO-JS, WENO-SYMBO, MDCD-WENO Adapted with permission from Wang et al., “Consistent high resolution interface-capturing finite volume method for compressible multi-material flows,” Comput. Fluids 202, 104518 (2020). Copyright 2020 Elsevier Ltd. and (c) 2D and 3D CWENO adapted with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

Close modal

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.

FIG. 8.

Helium bubble shape at different times for test case IV G with different computational settings: (a) MUSCL-THINC/QQ-BVD, adapted with permission from Cheng et al., “Low-dissipation BVD schemes for single and multi-phase compressible flows on unstructured grids,” J. Comput. Phys. 428, 110088 (2020). Copyright 2020 Elsevier Ltd. (b) Fifth order CWENO on unstructured grid, adapted with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license. (c) Fifth-order WENO, adapted with permission from Coralic and Colonius, “Finite-volume WENO scheme for viscous compressible multicomponent flows,” J. Comput. Phys. 274, 95–121 (2014). Copyright 2014 Elsevier Ltd. (d) DG(P2) adapted with permission from Cheng et al., “A discontinuous Galerkin method for the simulation of compressible gas-gas and gas-water two-medium flows,” J. Comput. Phys. 403, 109059 (2020). Copyright 2020 Elsevier Ltd. (e) WENO-JS of Ref. 46. (f) WENO-SYMBO of Ref. 165 and (g) MDCD-WENO adapted with permission from Wang et al., “Consistent high resolution interface-capturing finite volume method for compressible multi-material flows,” Comput. Fluids 202, 104518 (2020). Copyright 2020 Elsevier Ltd.

FIG. 8.

Helium bubble shape at different times for test case IV G with different computational settings: (a) MUSCL-THINC/QQ-BVD, adapted with permission from Cheng et al., “Low-dissipation BVD schemes for single and multi-phase compressible flows on unstructured grids,” J. Comput. Phys. 428, 110088 (2020). Copyright 2020 Elsevier Ltd. (b) Fifth order CWENO on unstructured grid, adapted with permission from Tsoutsanis et al., “CWENO finite-volume interface capturing schemes for multicomponent flows using unstructured meshes,” J. Sci. Comput. 89, 64 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license. (c) Fifth-order WENO, adapted with permission from Coralic and Colonius, “Finite-volume WENO scheme for viscous compressible multicomponent flows,” J. Comput. Phys. 274, 95–121 (2014). Copyright 2014 Elsevier Ltd. (d) DG(P2) adapted with permission from Cheng et al., “A discontinuous Galerkin method for the simulation of compressible gas-gas and gas-water two-medium flows,” J. Comput. Phys. 403, 109059 (2020). Copyright 2020 Elsevier Ltd. (e) WENO-JS of Ref. 46. (f) WENO-SYMBO of Ref. 165 and (g) MDCD-WENO adapted with permission from Wang et al., “Consistent high resolution interface-capturing finite volume method for compressible multi-material flows,” Comput. Fluids 202, 104518 (2020). Copyright 2020 Elsevier Ltd.

Close modal
FIG. 9.

Helium bubble shape in terms of volume fraction isovalues. Comparison between DG0 (top half) and DG1 (bottom half) resolution. [Adapted with permission from E. Franquet and V. Perrier, “Runge-Kutta discontinuous Galerkin method for the approximation of Baer and Nunziato type multiphase models,” J. Comput. Phys. 231, 4096–4141 (2012). Copyright 2012 Elsevier Ltd.]

FIG. 9.

Helium bubble shape in terms of volume fraction isovalues. Comparison between DG0 (top half) and DG1 (bottom half) resolution. [Adapted with permission from E. Franquet and V. Perrier, “Runge-Kutta discontinuous Galerkin method for the approximation of Baer and Nunziato type multiphase models,” J. Comput. Phys. 231, 4096–4141 (2012). Copyright 2012 Elsevier Ltd.]

Close modal
FIG. 10.

Helium bubble shape computed with hybrid LB-FV scheme. Comparison of (a) LAMBS-MUSCL scheme and (b) LAMBS-WENO scheme, adapted with permission from Sashi Kumar and Maheshwari, “A lattice Boltzmann method for simulation of multi-species shock accelerated flows,” Int. J. Comput. Fluid Dyn. 32, 19–34 (2018). Copyright 2018 Taylor and Francis Online.

FIG. 10.

Helium bubble shape computed with hybrid LB-FV scheme. Comparison of (a) LAMBS-MUSCL scheme and (b) LAMBS-WENO scheme, adapted with permission from Sashi Kumar and Maheshwari, “A lattice Boltzmann method for simulation of multi-species shock accelerated flows,” Int. J. Comput. Fluid Dyn. 32, 19–34 (2018). Copyright 2018 Taylor and Francis Online.

Close modal

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 Mie–Grüneisen EOS is used for the materials with the following coefficients:

(ρ0,c0,s,Γ0)={(9.961,4.77,1.43,2.53),Molybdenum,(2.66,2.1,1.68,1.18),MORB  liquid.
(151)

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

(ρ,u,p)={(9.964,(0,0),0),pre-shockmolybdenum,(11.042,(0.543,0),30),post-shockmolybdenum,
(152)

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μs and t=110μs for density contours (Fig. 11) and pressure contours (Fig. 12). At time t=50μs, one can observe through the density plot the incident shock and the reflected shock in molybdenum and MORB, respectively, with the latter moving slower than the former. At time t=110μs, the transmitted shock still have not traveled the MORB capsule completely. Figures 11 and 12 compare results from different references computed using FV with and without sharpening techniques, and DG method with different polynomial order approximation.

FIG. 11.

Density contour plots for the shock in molybdenum-encapsulated MORB test case for different schemes: (a) FV with THINC reconstruction and (b) FV, adapted with permission from Shyue and Xiao, “An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach,” J. Comput. Phys. 268, 326–354 (2014). Copyright 2104 Elsevier Ltd. (c) FV with anti-diffusion terms, adapted with permission from So et al., “Anti-diffusion interface sharpening technique for two-phase compressible flow simulations,” J. Comput. Phys. 231, 4304–4323 (2012). Copyright 2012 Elsevier, Ltd. (d) and (e) DG with P1 and P2 order discretization, adapted with permission from Luo et al., “A quasi-conservative discontinuous Galerkin method for multi-component flow using the non-oscillatory kinetic flux,” J. Sci. Comput. 87, 69 (2020). Copyright 2020 SpringerLink.

FIG. 11.

Density contour plots for the shock in molybdenum-encapsulated MORB test case for different schemes: (a) FV with THINC reconstruction and (b) FV, adapted with permission from Shyue and Xiao, “An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach,” J. Comput. Phys. 268, 326–354 (2014). Copyright 2104 Elsevier Ltd. (c) FV with anti-diffusion terms, adapted with permission from So et al., “Anti-diffusion interface sharpening technique for two-phase compressible flow simulations,” J. Comput. Phys. 231, 4304–4323 (2012). Copyright 2012 Elsevier, Ltd. (d) and (e) DG with P1 and P2 order discretization, adapted with permission from Luo et al., “A quasi-conservative discontinuous Galerkin method for multi-component flow using the non-oscillatory kinetic flux,” J. Sci. Comput. 87, 69 (2020). Copyright 2020 SpringerLink.

Close modal
FIG. 12.

Pressure contour plots for the shock in molybdenum-encapsulated MORB test case for different schemes: (a) FV with THINC reconstruction and (b) FV, adapted with permission from Shyue and Xiao, “An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach,” J. Comput. Phys. 268, 326–354 (2014). Copyright 2014 Elsevier, Ltd. (c) FV with anti-diffusion terms, adapted with permission from So et al., “Anti-diffusion interface sharpening technique for two-phase compressible flow simulations,” J. Comput. Phys. 231, 4304–4323 (2012). Copyright 2012 Elsevier, Ltd. (d) and (e) DG with P1 and P2 order discretization, adapted with permission from Luo et al., “A quasi-conservative discontinuous Galerkin method for multi-component flow using the non-oscillatory kinetic flux,” J. Sci. Comput. 87, 69 (2020). Copyright 2020 SpringerLink.

FIG. 12.

Pressure contour plots for the shock in molybdenum-encapsulated MORB test case for different schemes: (a) FV with THINC reconstruction and (b) FV, adapted with permission from Shyue and Xiao, “An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach,” J. Comput. Phys. 268, 326–354 (2014). Copyright 2014 Elsevier, Ltd. (c) FV with anti-diffusion terms, adapted with permission from So et al., “Anti-diffusion interface sharpening technique for two-phase compressible flow simulations,” J. Comput. Phys. 231, 4304–4323 (2012). Copyright 2012 Elsevier, Ltd. (d) and (e) DG with P1 and P2 order discretization, adapted with permission from Luo et al., “A quasi-conservative discontinuous Galerkin method for multi-component flow using the non-oscillatory kinetic flux,” J. Sci. Comput. 87, 69 (2020). Copyright 2020 SpringerLink.

Close modal

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 P2 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.

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

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

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

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

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

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

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

The authors have no conflicts to disclose.

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

1.
M. R.
Baer
and
J. W.
Nunziato
, “
A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials
,”
Int. J. Multiphase Flow
12
,
861
889
(
1986
).
2.
A. K.
Kapila
,
R.
Menikoff
,
J. B.
Bdzil
,
S. F.
Son
, and
D. S.
Stewart
, “
Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations
,”
Phys. Fluids
13
,
3002
3024
(
2001
).
3.
G.
Allaire
,
S.
Clerc
, and
S.
Kokh
, “
A five-equation model for the simulation of interfaces between compressible fluids
,”
J. Comput. Phys.
181
,
577
616
(
2002
).
4.
E.
Johnsen
and
T.
Colonius
, “
Implementation of WENO schemes in compressible multicomponent flow problems
,”
J. Comput. Phys.
219
,
715
732
(
2006
).
5.
X.
Deng
and
P.
Boivin
, “
Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet
,”
Appl. Math. Modell.
84
,
405
424
(
2020
).
6.
L.
Jofre
and
J.
Urzay
, “
Transcritical diffuse-interface hydrodynamics of propellants in high-pressure combustors of chemical propulsion systems
,”
Prog. Energy Combust. Sci.
82
,
100877
(
2021
).
7.
J. T.
Oden
,
A.
Hawkins
, and
S.
Prudhomme
, “
General diffuse-interface theories and an approach to predictive tumor growth modeling
,”
Math. Models Methods Appl. Sci.
20
,
477
517
(
2010
).
8.
N.
Favrie
and
S. L.
Gavrilyuk
, “
Diffuse interface model for compressible fluid—Compressible elastic-plastic solid interaction
,”
J. Comput. Phys.
231
,
2695
2723
(
2012
).
9.
L.
Michael
and
N.
Nikiforakis
, “
A multi-physics methodology for the simulation of reactive flow and elastoplastic structural response
,”
J. Comput. Phys.
367
,
1
27
(
2018
).
10.
S.
Karni
, “
Multicomponent flow calculations by a consistent primitive algorithm
,”
J. Comput. Phys.
112
,
31
(
1994
).
11.
S. S.
Jain
,
A.
Mani
, and
P.
Moin
, “
A conservative diffuse-interface method for compressible two-phase flows
,”
J. Comput. Phys.
418
,
109606
(
2020
).
12.
B. D.
Nichols
and
C. W.
Hirt
, “
Improved free surface boundary conditions for numerical incompressible-flow calculations
,”
J. Comput. Phys.
8
,
434
448
(
1971
).
13.
J.
Glimm
,
J. W.
Grove
,
X. L.
Li
,
K. M.
Shyue
,
Y.
Zeng
, and
Q.
Zhang
, “
Three-dimensional front tracking
,”
SIAM J. Sci. Comput.
19
,
703
727
(
1998
).
14.
R.
Abgrall
and
S.
Karni
, “
Computations of compressible multifluids
,”
J. Comput. Phys.
169
,
594
623
(
2001
).
15.
R.
Saurel
and
R.
Abgrall
, “
A multiphase Godunov method for compressible multifluid and multiphase flows
,”
J. Comput. Phys.
150
,
425
467
(
1999
).
16.
J. P.
Boris
and
E. S.
Oran
, “
The numerical simulation of compressible reactive flows
,”
AlAA Paper No. 87-1323
,
1987
.
17.
R.
Saurel
and
R.
Abgrall
, “
Simple method for compressible multifluid flows
,”
SIAM J. Sci. Comput.
21
,
1115
1145
(
1999
).
18.
C.
Hirt
and
B.
Nichols
, “
Volume of fluid (VOF) method for the dynamics of free boundaries
,”
J. Comput. Phys.
39
,
201
(
1981
).
19.
R. P.
Fedkiw
,
T.
Aslam
,
B.
Merriman
, and
S.
Osher
, “
A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)
,”
J. Comput. Phys.
152
,
457
492
(
1999
).
20.
E. J.
Caramana
and
M. J.
Shashkov
, “
Elimination of artificial grid distortion and hourglass-type motions by means of Lagrangian subzonal masses and pressures
,”
J. Comput. Phys.
142
,
521
561
(
1998
).
21.
S. O.
Unverdi
and
G.
Tryggvason
, “
A front-tracking method for viscous, incompressible, multi-fluid flows
,”
J. Comput. Phys.
100
,
25
37
(
1992
).
22.
S.
Osher
and
J. A.
Sethian
, “
Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations
,”
J. Comput. Phys.
79
,
12
49
(
1988
).
23.
M.
Sussman
, “
A level set approach for computing solutions to incompressible two-phase flow
,”
J. Comput. Phys.
114
,
146
159
(
1994
).
24.
R. P.
Fedkiw
,
T.
Aslam
, and
S.
Xu
, “
The ghost fluid method for deflagration and detonation discontinuities
,”
J. Comput. Phys.
154
,
393
427
(
1999
).
25.
P.
Bigdelou
,
C.
Liu
,
P.
Tarey
, and
P.
Ramaprabhu
, “
An efficient ghost fluid method to remove overheating from material interfaces in compressible multi-medium flows
,”
Comput. Fluids
233
,
105250
(
2021
).
26.
P.
Bigdelou
, “
A numerical study of interfacial instabilities in shocked material with surface tension
,”
Electronic theses and dissertations
(The University of North Carolina at Charlotte, ProQuest Dissertation Publishing
,
2020
), No.
27959531
.
27.
R.
Saurel
and
C.
Pantano
, “
Diffuse-interface capturing methods for compressible two-phase flows
,” Annu. Rev. Fluid Mech.
50
,
105
130
(
2018
).
28.
R.
Scardovelli
and
S.
Zaleski
, “
Direct numerical simulation of free-surface and interfacial flow
,”
Annu. Rev. Fluid Mech.
31
,
567
603
(
1999
).
29.
K. M.
Shyue
, “
An efficient shock-capturing algorithm for compressible multicomponent problems
,”
J. Comput. Phys.
142
,
208
242
(
1998
).
30.
R.
Abgrall
, “
How to prevent pressure oscillations in multicomponent flow calculations: A quasi conservative approach
,”
J. Comput. Phys.
125
,
150
160
(
1996
).
31.
K. M.
Shyue
, “
A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state
,”
J. Comput. Phys.
156
,
43
88
(
1999
).
32.
J. W.
Cahn
and
J. E.
Hilliard
, “
Free energy of a nonuniform system. I. Interfacial free energy
,”
J. Chem. Phys.
28
,
258
267
(
1958
).
33.
G.
Perigaud
and
R.
Saurel
, “
A compressible flow model with capillary effects
,”
J. Comput. Phys.
209
,
139
178
(
2005
).
34.
R.
Saurel
,
F.
Petitpas
, and
R.
Abgrall
, “
Modelling phase transition in metastable liquids: Application to cavitating and flashing flows
,”
J. Fluid Mech.
607
,
313
350
(
2008
).
35.
S. K.
Godunov
, “
An interesting class of quasi-linear systems
,”
Usp. Mat. Nauk SSSR
139
,
521
523
(
1961
).
36.
S. K.
Godunov
and
E. I.
Romenskii
, “
Nonstationary equations of nonlinear elasticity theory in Eulerian coordinates
,”
J. Appl. Mech. Tech. Phys.
13
,
868
884
(
1972
).
37.
K.
Schmidmayer
,
S. H.
Bryngelson
, and
T.
Colonius
, “
An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics
,”
J. Comput. Phys.
402
,
109080
(
2020
).
38.
A.
Tiwari
,
J. B.
Freund
, and
C.
Pantano
, “
A diffuse interface model with immiscibility preservation
,”
J. Comput. Phys.
252
,
290
309
(
2013
).
39.
R. K.
Shukla
,
C.
Pantano
, and
J. B.
Freund
, “
An interface capturing method for the simulation of multi-phase compressible flows
,”
J. Comput. Phys.
229
,
7411
7439
(
2010
).
40.
S.
Ii
,
K.
Sugiyama
,
S.
Takeuchi
,
S.
Takagi
,
Y.
Matsumoto
, and
F.
Xiao
, “
An interface capturing method with a continuous function: The THINC method with multi-dimensional reconstruction
,”
J. Comput. Phys.
231
,
2328
2358
(
2012
).
41.
K. M.
Shyue
and
F.
Xiao
, “
An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach
,”
J. Comput. Phys.
268
,
326
354
(
2014
).
42.
A.
Chiapolino
,
R.
Saurel
, and
B.
Nkonga
, “
Sharpening diffuse interfaces with compressible fluids on unstructured meshes
,”
J. Comput. Phys.
340
,
389
417
(
2017
).
43.
B.
Van Leer
, “
Towards the ultimate conservative difference scheme. V. A second-order Sequel to Godunov's method
,”
J. Comput. Phys.
32
,
101
136
(
1979
).
44.
T.
Barth
and
D.
Jespersen
, “
The design and application of upwind schemes on unstructured meshes
,”
AIAA Paper
No. 89-0366,
1989
.
45.
A.
Kurganov
and
E.
Tadmor
, “
New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations
,”
J. Comput. Phys.
160
,
241
282
(
2000
).
46.
G.
Jiang
and
C.
Shu
, “
Efficient implementation of weighted ENO schemes
,”
J. Comput. Phys.
126
,
202
228
(
1996
).
47.
C.-W.
Shu
, “
Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic laws
,”
ICASE Report No. 97-65 (Brown University,
1997
), pp.
1
78
.
48.
X.
Liu
,
S.
Osher
, and
T.
Chan
, “
Weighted essentially non-oscillatory schemes
,”
J. Comput. Phys.
115
,
200
212
(
1994
).
49.
P.
Tsoutsanis
,
V. A.
Titarev
, and
D.
Drikakis
, “
WENO schemes on arbitrary mixed-element unstructured meshes in three space dimensions
,”
J. Comput. Phys.
230
,
1585
1601
(
2011
).
50.
M.
Dumbser
, “
A simple two-phase method for the simulation of complex free surface flows
,”
Comput. Methods Appl. Mech. Eng.
200
,
1204
1219
(
2011
).
51.
P.
Tsoutsanis
, “
Stencil selection algorithms for WENO schemes on unstructured meshes
,”
J. Comput. Phys.: X
4
,
100037
(
2019
).
52.
V. A.
Titarev
and
E. F.
Toro
, “
ADER: Arbitrary high order Godunov approach
,”
J. Sci. Comput.
17
,
609
618
(
2002
).
53.
V.
Coralic
and
T.
Colonius
, “
Finite-volume WENO scheme for viscous compressible multicomponent flows
,”
J. Comput. Phys.
274
,
95
121
(
2014
).
54.
B.
Cockburn
and
C.-W.
Shu
, “
TVB Runge-Kutta projection discontinuous Galerkin finite element methods. II. PDF
,”
Math. Comput.
52
,
411
435
(
1989
).
55.
B.
Cockburn
and
C. W.
Shu
, “
The local discontinuous Galerkin method for time-dependent convection-diffusion systems
,”
SIAM J. Numer. Anal.
35
,
2440
2463
(
1998
).
56.
B.
Cockburn
and
C. W.
Shu
, “
Runge-Kutta discontinuous Galerkin methods for convection-dominated problems
,”
J. Sci. Comput.
16
,
173
261
(
2001
).
57.
M.
Dumbser
,
D. S.
Balsara
,
E. F.
Toro
, and
C. D.
Munz
, “
A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes
,”
J. Comput. Phys.
227
,
8209
8253
(
2008
).
58.
H.
Luo
,
L.
Luo
,
R.
Nourgaliev
, and
V. A.
Mousseau
, “
A reconstructed discontinuous Galerkin method for the compressible Euler equations on arbitrary grids
,”
AIAA Paper No. 2009-3788
,
2009
.
59.
F.
Zhang
,
J.
Cheng
, and
T.
Liu
, “
A reconstructed discontinuous Galerkin method for incompressible flows on arbitrary grids
,”
J. Comput. Phys.
418
,
109580
(
2020
).
60.
A.
Syrakos
,
S.
Varchanis
,
Y.
Dimakopoulos
,
A.
Goulas
, and
J.
Tsamopoulos
, “
A critical analysis of some popular methods for the discretisation of the gradient operator in finite volume methods
,”
Phys. Fluids
29
,
127103
(
2017
).
61.
H.
Luo
,
Y.
Xia
,
S.
Spiegel
,
R.
Nourgaliev
, and
Z.
Jiang
, “
A reconstructed discontinuous Galerkin method based on a hierarchical WENO reconstruction for compressible flows on tetrahedral grids
,”
J. Comput. Phys.
236
,
477
492
(
2013
).
62.
J.
Qiu
and
C. W.
Shu
, “
Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method: One-dimensional case
,”
J. Comput. Phys.
193
,
115
135
(
2004
).
63.
A. K.
Pandare
,
J.
Waltz
, and
J.
Bakosi
, “
A reconstructed discontinuous Galerkin method for multi-material hydrodynamics with sharp interfaces
,”
Int. J. Numer. Methods Fluids
92
,
874
898
(
2020
).
64.
S.
Schmieschek
,
L.
Shamardin
,
S.
Frijters
,
T.
Krüger
,
U. D.
Schiller
,
J.
Harting
, and
P. V.
Coveney
, “
LB3D: A parallel implementation of the Lattice-Boltzmann method for simulation of interacting amphiphilic fluids
,”
Comput. Phys. Commun.
217
,
149
161
(
2017
).
65.
K.
Sverdrup
,
N.
Nikiforakis
, and
A.
Almgren
, “
Highly parallelisable simulations of time-dependent viscoplastic fluid flow with structured adaptive mesh refinement
,”
Phys. Fluids
30
,
093102
(
2018
).
66.
S.
Chen
and
G. D.
Doolen
, “
Lattice Boltzmann method for fluid flows
,”
Annu. Rev. Fluid Mech.
30
,
329
(
1998
).
67.
A.
Lintermann
and
W.
Schröder
, “
Lattice–Boltzmann simulations for complex geometries on high-performance computers
,”
CEAS Aeronaut. J.
11
,
745
766
(
2020
).
68.
A. K.
Gunstensen
,
D. H.
Rothman
,
S.
Zaleski
, and
G.
Zanetti
, “
Lattice Boltzmann model of immiscible fluids
,”
Phys. Rev. A
43
,
4320
4327
(
1991
).
69.
D. H.
Rothman
and
J. M.
Keller
, “
Immiscible cellular-automaton fluids
,”
J. Stat. Phys.
52
,
1119
1127
(
1988
).
70.
X.
Shan
and
H.
Chen
, “
Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation
,”
Phys. Rev. E
49
,
2941
2948
(
1994
).
71.
M. R.
Swift
,
E.
Orlandini
,
W. R.
Osborn
, and
J. M.
Yeomans
, “
Lattice Boltzmann simulations of liquid-gas and binary fluid systems
,”
Phys. Rev. E
54
,
5041
5052
(
1996
).
72.
H. W.
Zheng
,
C.
Shu
, and
Y. T.
Chew
, “
A lattice Boltzmann model for multiphase flows with large density ratio
,”
J. Comput. Phys.
218
,
353
371
(
2006
).
73.
T.
Reis
, “
A conservative interface sharpening lattice Boltzmann model
,”
SIAM J. Sci. Comput.
40
,
B1495
B1516
(
2018
).
74.
A.
De Rosis
and
C.
Coreixas
, “
Multiphysics flow simulations using D3Q19 lattice Boltzmann methods based on central moments
,”
Phys. Fluids
32
,
117101
(
2020
).
75.
K.
Qu
,
C.
Shu
, and
Y. T.
Chew
, “
Lattice Boltzmann and finite volume simulation of inviscid compressible flows with curved boundary
,”
Adv. Appl. Math. Mech.
2
,
573
586
(
2010
).
76.
F.
Nannelli
and
S.
Succi
, “
The lattice Boltzmann equation on irregular lattices
,”
J. Stat. Phys.
68
,
401
407
(
1992
).
77.
H.
Joshi
,
A.
Agarwal
,
B.
Puranik
,
C.
Shu
, and
A.
Agrawal
, “
A hybrid FVM-LBM method of single and multi-fluid compressible flow problems
,”
Int. J. Numer. Methods Fluids
62
,
403
427
(
2010
).
78.
A.
Jameson
and
D.
Mavriplis
, “
Finite volume solution of the two-dimensional Euler equations on a regular triangular mesh
,”
AIAA J.
24
,
611
618
(
1986
).
79.
Z. J.
Wang
, “
High-order methods for the Euler and Navier-Stokes equations on unstructured grids
,”
Prog. Aerosp. Sci.
43
,
1
41
(
2007
).
80.
D. F.
Harlacher
,
H.
Klimach
,
S.
Roller
,
C.
Siebert
, and
F.
Wolf
, “
Dynamic load balancing for unstructured meshes on space-filling curves
,” in
Proceedings of the IEEE 26th International Parallel and Distributed Processing Symposium Workshops (IEEE, 2012), Paper No. IPDPSW
2012
, pp.
1661
1669
.
81.
P.
Tsoutsanis
,
A. F.
Antoniadis
, and
K. W.
Jenkins
, “
Improvement of the computational performance of a parallel unstructured WENO finite volume CFD code for implicit large eddy simulation
,”
Comput. Fluids
173
,
157
170
(
2018
).
82.
J. B.
Bdzil
,
R.
Menikoff
,
S. F.
Son
,
A. K.
Kapila
, and
D. S.
Stewart
, “
Two-phase modeling of deflagration-to-detonation transition in granular materials: A critical examination of modeling issues
,”
Phys. Fluids
11
,
378
402
(
1999
).
83.
W. L.
Elban
and
M. A.
Chiarito
, “
Quasi-static compaction study of coarse HMX explosive
,”
Powder Technol.
46
,
181
193
(
1986
).
84.
R.
Saurel
,
A.
Chinnayya
, and
Q.
Carmouze
, “
Modelling compressible dense and dilute two-phase flows
,”
Phys. Fluids
29
,
063301
(
2017
).
85.
A.
Wood
,
A Textbook of Sound
(
Macmillan
,
New York
,
1930
).
86.
B.
Larrouturou
, “
How to preserve the mass fractions positivity when computing multi-component flows
,”
J. Comput. Phys.
95
,
59
84
(
1991
).
87.
X.
Lei
and
J.
Li
, “
A non-oscillatory energy-splitting method for the computation of compressible multi-fluid flows
,”
Phys. Fluids
30
,
040906
(
2018
).
88.
D. P.
Garrick
,
M.
Owkes
, and
J. D.
Regele
, “
A finite-volume HLLC-based scheme for compressible interfacial flows with surface tension
,”
J. Comput. Phys.
339
,
46
67
(
2017
).
89.
B.
Thornber
,
M.
Groom
, and
D.
Youngs
, “
A five-equation model for the simulation of miscible and viscous compressible fluids
,”
J. Comput. Phys.
372
,
256
280
(
2018
).
90.
P.
Yi
,
S.
Yang
,
C.
Habchi
, and
R.
Lugo
, “
A multicomponent real-fluid fully compressible four-equation model for two-phase flow with phase change
,”
Phys. Fluids
31
,
026102
(
2019
).
91.
M.
Pelanti
and
K. M.
Shyue
, “
A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves
,”
J. Comput. Phys.
259
,
331
357
(
2014
).
92.
O. L.
Métayer
,
J.
Massoni
, and
R.
Saurel
, “
Elaborating equations of state of a liquid and its vapor for two-phase flow models
,”
Int. J. Therm. Sci.
43
,
265
276
(
2004
).
93.
O.
Le Métayer
and
R.
Saurel
, “
The Noble-Abel stiffened-gas equation of state
,”
Phys. Fluids
28
046102
(
2016
).
94.
S.
Ndanou
,
N.
Favrie
, and
S.
Gavrilyuk
, “
Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation
,”
J. Comput. Phys.
295
,
523
555
(
2015
).
95.
N.
Favrie
,
S. L.
Gavrilyuk
, and
R.
Saurel
, “
Solid-fluid diffuse interface model in cases of extreme deformations
,”
J. Comput. Phys.
228
,
6037
6077
(
2009
).
96.
N.
Favrie
and
S.
Gavrilyuk
, “
Mathematical and numerical model for nonlinear viscoplasticity
,”
Philos. Trans. R. Soc. A
369
,
2864
2880
(
2011
).
97.
N.
Favrie
,
S.
Gavrilyuk
,
B.
Nkonga
, and
R.
Saurel
, “
Sharpening diffuse interfaces with compressible flow solvers
,”
Open J. Fluid Dyn.
04
,
44
68
(
2014
).
98.
L.
Michael
,
S. T.
Millmore
, and
N.
Nikiforakis
, “
A multi-physics methodology for four-states of matter
,”
Commun. Appl. Math. Comput.
2
,
487
514
(
2020
).
99.
S.
Kulkarni
and
A.
Tabarraei
, “
An analytical study of wave propagation in a peridynamic bar with nonuniform discretization
,”
Eng. Fract. Mech.
190
,
347
366
(
2018
).
100.
J. A.
Trangenstein
, “
A second-order Godunov algorithm for two-dimensional solid mechanics
,”
Comput. Mech.
13
,
343
359
(
1994
).
101.
G. H.
Miller
and
P.
Colella
, “
A high-order Eulerian Godunov method for elastic-plastic flow in solids
,”
J. Comput. Phys.
167
,
131
176
(
2001
).
102.
M.
Nikodemou
,
L.
Michael
, and
N.
Nikiforakis
, “
A unified multi-phase and multi-material formulation for combustion modeling
,”
Phys. Fluids
33
,
106113
(
2021
).
103.
X.
Wang
,
S. S.
Kulkarni
, and
A.
Tabarraei
, “
Concurrent coupling of peridynamics and classical elasticity for elastodynamic problems
,”
Comput. Methods Appl. Mech. Eng.
344
,
251
275
(
2019
).
104.
G. H.
Miller
and
P.
Colella
, “
A conservative three-dimensional Eulerian method for coupled solid-fluid shock capturing
,”
J. Comput. Phys.
183
,
26
82
(
2002
).
105.
S. L.
Gavrilyuk
,
N.
Favrie
, and
R.
Saurel
, “
Modelling wave dynamics of compressible elastic materials
,”
J. Comput. Phys.
227
,
2941
2969
(
2008
).
106.
L.
Li
,
Q.
Chen
, and
B.
Tian
, “
Transport diffuse interface model for simulation of solid-fluid
,”
Appl. Math. Mech.
40
,
321
330
(
2019
).
107.
P. T.
Barton
, “
An interface-capturing Godunov method for the simulation of compressible solid-fluid problems
,”
J. Comput. Phys.
390
,
25
50
(
2019
).
108.
N. S.
Ottosen
and
M.
Ristinmaa
,
The Mechanics of Constitutive Modeling
(
Elsevier Science
,
2005
).
109.
T.
Wallis
,
P. T.
Barton
, and
N.
Nikiforakis
, “
A flux-enriched Godunov method for multi-material problems with interface slide and void opening
,”
J. Comput. Phys.
442
,
110499
(
2021
).
110.
L.
Michael
and
N.
Nikiforakis
, “
A hybrid formulation for the numerical simulation of condensed phase explosives
,”
J. Comput. Phys.
316
,
193
217
(
2016
).
111.
J. W.
Banks
,
D. W.
Schwendeman
,
A. K.
Kapila
, and
W. D.
Henshaw
, “
A high-resolution Godunov method for compressible multi-material flow on overlapping grids
,”
J. Comput. Phys.
223
,
262
297
(
2007
).
112.
P. T.
Barton
,
D.
Drikakis
,
E.
Romenski
, and
V. A.
Titarev
, “
Exact and approximate solutions of Riemann problems in non-linear elasticity
,”
J. Comput. Phys.
228
,
7046
7068
(
2009
).
113.
T.
Wallis
,
P. T.
Barton
, and
N.
Nikiforakis
, “
A diffuse interface model of reactive-fluids and solid-dynamics
,”
Comput. Struct.
254
,
106578
(
2021
).
114.
I.
Peshkov
and
E.
Romenski
, “
A hyperbolic model for viscous Newtonian flows
,”
Continuum Mech. Thermodyn.
28
,
85
104
(
2016
).
115.
S. K.
Godunov
and
E. I.
Romenskii
,
Elements of Continuum Mechanics and Conservation Laws
(
Kluwer Academic–Plenum Publisher
,
2003
).
116.
J.
Massoni
,
R.
Saurel
,
B.
Nkonga
, and
R.
Abgrall
, “
Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur
,”
Int. J. Heat Mass Transfer
45
,
1287
1307
(
2002
).
117.
R.
Saurel
,
S. L.
Martelot
,
R.
Tosello
, and
E.
Lapébie
, “
Symmetric model of compressible granular mixtures with permeable interfaces
,”
Phys. Fluids
26
,
123304
(
2014
).
118.
J. A.
Ekaterinaris
, “
High-order accurate, low numerical diffusion methods for aerodynamics
,”
Prog. Aerosp. Sci.
41
,
192
300
(
2005
).
119.
D. A.
Di Pietro
and
A.
Ern
,
Mathematical Aspects of Discontinuous Galerkin Methods
(
Springer
,
Heidelberg, Dordrecht, London, New York
,
2011
).
120.
H.
Luo
,
J. D.
Baum
, and
R.
Löhner
, “
A discontinuous Galerkin method based on a Taylor basis for the compressible flows on arbitrary grids
,”
J. Comput. Phys.
227
,
8875
8893
(
2008
).
121.
H.
Luo
,
L.
Luo
, and
K.
Xu
, “
A BGK-based discontinuous Galerkin method for the Navier-Stokes equations on arbitrary grids
,”
Comput. Fluid Dyn. Rev.
1
,
103
122
(
2010
).
122.
H.
Luo
,
Y.
Xia
,
S.
Li
,
R.
Nourgaliev
, and
C.
Cai
, “
A Hermite WENO reconstruction-based discontinuous Galerkin method for the Euler equations on tetrahedral grids
,”
J. Comput. Phys.
231
,
5489
5503
(
2012
).
123.
M.
Dumbser
, “
Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier-Stokes equations
,”
Comput. Fluids
39
,
60
76
(
2010
).
124.
B.
Van Leer
and
S.
Nomura
, “
Discontinuous Galerkin for diffusion
,”
AIAA Paper No. 2005-5108
,
2005
.
125.
B.
Van Leer
,
M.
Lo
, and
M.
Van Raalte
, “
A discontinuous Galerkin method for diffusion based on recovery
,”
AIAA Paper No. 2007-4083
,
2007
, Vol.
1
, pp.
763
774
.
126.
H.
Luo
,
L.
Luo
,
R.
Nourgaliev
,
V. A.
Mousseau
, and
N.
Dinh
, “
A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids
,”
J. Comput. Phys.
229
,
6961
6978
(
2010
).
127.
J.
Cheng
,
T.
Liu
, and
H.
Luo
, “
A hybrid reconstructed discontinuous Galerkin method for compressible flows on unstructured grids
,”
AIAA Paper No. 2016-1330
,
2016
.
128.
L.
Zhang
,
L.
Wei
,
H.
Lixin
,
D.
Xiaogang
, and
Z.
Hanxin
, “
A class of hybrid DG/FV methods for conservation laws. I. Basic formulation and one-dimensional systems
,”
J. Comput. Phys.
231
,
1081
1103
(
2012
).
129.
J.
Cheng
and
T.
Liu
, “
A variational reconstructed discontinuous Galerkin method for the steady-state compressible flows on unstructured grids
,”
J. Comput. Phys.
380
,
65
87
(
2019
).
130.
Q.
Wang
,
Y. X.
Ren
,
J.
Pan
, and
W.
Li
, “
Compact high order finite volume method on unstructured grids. III. Variational reconstruction
,”
J. Comput. Phys.
337
,
1
26
(
2017
).
131.
A.
Harten
, “
ENO schemes with subcell resolution
,”
J. Comput. Phys.
83
,
148
184
(
1989
).
132.
P.
Tsoutsanis
, “
Extended bounds limiter for high-order finite-volume schemes on unstructured meshes
,”
J. Comput. Phys.
362
,
69
94
(
2018
).
133.
J.
Shi
,
C.
Hu
, and
C. W.
Shu
, “
A technique of treating negative weights in WENO schemes
,”
J. Comput. Phys.
175
,
108
127
(
2002
).
134.
M.
Dumbser
and
M.
Käser
, “
Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems
,”
J. Comput. Phys.
221
,
693
723
(
2007
).
135.
D.
Levy
,
G.
Puppo
, and
G.
Russo
, “
Compact central WENO schemes for multidimensional conservation laws
,”
SIAM J. Sci. Comput.
22
,
656
672
(
2000
).
136.
M.
Dumbser
,
I.
Peshkov
,
E.
Romenski
, and
O.
Zanotti
, “
High order ADER schemes for a unified first order hyperbolic formulation of continuum mechanics: Viscous heat-conducting fluids and elastic solids
,”
J. Comput. Phys.
314
,
824
862
(
2016
).
137.
J.
Qiu
and
C.
Shu
, “
Runge-Kutta discontinuous Galerkin method using WENO limiters
,”
SIAM J. Sci. Comput.
26
,
907
929
(
2005
).
138.
D. S.
Balsara
,
C.
Altmann
,
C. D.
Munz
, and
M.
Dumbser
, “
A sub-cell based indicator for troubled zones in RKDG schemes and a novel class of hybrid RKDG+HWENO schemes
,”
J. Comput. Phys.
226
,
586
620
(
2007
).
139.
T.
Inamuro
,
M.
Yoshino
, and
F.
Ogino
, “
Accuracy of the lattice Boltzmann method for small Knudsen number with finite Reynolds number
,”
Phys. Fluids
9
,
3535
3542
(
1997
).
140.
C. K.
Aidun
and
J. R.
Clausen
, “
Lattice-Boltzmann method for complex flows
,”
Annu. Rev. Fluid Mech.
42
,
439
472
(
2010
).
141.
A.
Scagliarini
,
L.
Biferale
,
M.
Sbragaglia
,
K.
Sugiyama
, and
F.
Toschi
, “
Lattice Boltzmann methods for thermal flows: Continuum limit and applications to compressible Rayleigh-Taylor systems
,”
Phys. Fluids
22
,
055101
055121
(
2010
).
142.
N. I.
Prasianakis
and
I. V.
Karlin
, “
Lattice Boltzmann method for thermal flow simulation on standard lattices
,”
Phys. Rev. E
76
,
016702
(
2007
).
143.
P.
Lallemand
and
L. S.
Luo
, “
Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions
,”
Phys. Rev. E
68
,
036706
(
2003
).
144.
X.
Shan
,
X. F.
Yuan
, and
H.
Chen
, “
Kinetic theory representation of hydrodynamics: A way beyond the Navier-Stokes equation
,”
J. Fluid Mech.
550
,
413
441
(
2006
).
145.
X.
He
and
L. S.
Luo
, “
Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation
,”
Phys. Rev. E
56
,
6811
6820
(
1997
).
146.
P.
Bathnagar
,
E.
Gross
, and
M.
Krook
, “
A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems
,”
Phys. Rev.
94
,
511
(
1954
).
147.
C.
Coreixas
and
J.
Latt
, “
Compressible lattice Boltzmann methods with adaptive velocity stencils: An interpolation-free formulation
,”
Phys. Fluids
32
,
116102
116123
(
2020
).
148.
F. J.
Alexander
,
S.
Chen
, and
J. D.
Sterling
, “
Lattice Boltzmann thermohydrodynamics
,”
Phys. Rev. E
47
,
R2249(R)
(
1993
).
149.
N.
Frapolli
,
S. S.
Chikatamarla
, and
I. V.
Karlin
, “
Entropic lattice Boltzmann model for compressible flows
,”
Phys. Rev. E
92
,
061301
(
2015
).
150.
X.
He
,
S.
Chen
, and
G. D.
Doolen
, “
A Novel thermal model for the lattice Boltzmann method in incompressible limit
,”
J. Comput. Phys.
146
,
282
300
(
1998
).
151.
Q.
Li
,
Y. L.
He
,
Y.
Wang
, and
W. Q.
Tao
, “
Coupled double-distribution-function lattice Boltzmann method for the compressible Navier-Stokes equations
,”
Phys. Rev. E
76
,
056705
(
2007
).