A seven-equation diffused interface method for resolved multiphase flows

Achyut Panchal Spencer H. Bryngelson Suresh Menon School of Aerospace Engineering,
Georgia Institute of Technology,
Atlanta, GA, 30332. School of Computational Science and Engineering,
Georgia Institute of Technology,
Atlanta, GA, 30332.

Abstract

The seven-equation model is a compressible multiphase formulation that allows for phasic velocity and pressure disequilibrium. These equations are solved using a diffused interface method that models resolved multiphase flows. Novel extensions are proposed for including the effects of surface tension, viscosity, multi-species, and reactions. The allowed non-equilibrium of pressure in the seven-equation model provides numerical stability in strong shocks and allows for arbitrary and independent equations of states. A discrete equations method (DEM) models the fluxes. We show that even though stiff pressure- and velocity-relaxation solvers have been used, they are not needed for the DEM because the non-conservative fluxes are accurately modeled. An interface compression scheme controls the numerical diffusion of the interface, and its effects on the solution are discussed. Test cases are used to validate the computational method and demonstrate its applicability. They include multiphase shock tubes, shock propagation through a material interface, a surface-tension-driven oscillating droplet, an accelerating droplet in a viscous medium, and shock–detonation interacting with a deforming droplet. Simulation results are compared against exact solutions and experiments when possible.

1 Introduction

Multiphase flows are characterized by more than one phase, each with its own material properties. Such flows are encountered in many situations, including liquid fuel-powered devices [1, 2], solid explosives [3], pharmaceutical sprays [4], and more. The design and engineering of these applications require computational modeling. From this perspective, there are two sets of methods: resolved approaches where the multiphase entities (MPE) are fully-resolved on the computational grid [5, 6], and dispersed approaches where the MPEs are unresolved on the grid and instead treated as point particles [7, 8, 9, 10].

This work focuses on resolved multiphase flow modeling, which can be further segmented into sharp- or diffused-interface methods (SIM or DIM). The SIM naturally represents materials as infinitely sharp [11, 12]. Specialized techniques are used to advect the interface, including the level-set method (LS) [11], the volume of fluid approach (VOF) [13], and the coupled LS–VOF method (CLSVOF) [14]. Example methods for handling jumps at the interface are continuum surface forces (CSF) [15] and ghost fluid methods (GFM) [16]. SIM are particularly attractive for their ease in specifying jump conditions and efficient interface tracking. However, they have historically been limited by low-speed flows, low density–viscosity ratios [17], and mass-conservation errors [11]. We note that some of these drawbacks are changing via bespoke methods [18, 19, 20].

Diffused interface methods (DIM) are hallmarked by the Baer–Nunziato (BN) “seven-equation” formulation of a compressible and mixture conservative flow [21, 6]. The equations follow from an averaging procedure to the phasic governing flow equations. The derived partial differential equation (PDE) solves for phasic volume fractions and phase-averaged quantities. This formulation evolves the conserved quantities on the grid, usually via numerical methods that allow for sharp gradients [22, 23]. As a result, the interface naturally diffuses across multiple grid cells due to numerical diffusion [6]. Besides their conservation properties, DIMs also allow for the creation or destruction of interfaces, common during supercavitation and flashing [24, 25, 26, 27].

The first numerical method for a single-component inviscid seven-equation model was developed by Saurel and Abgrall [28]. They used a two-step approach, where the two-phase equations were first solved using a hyperbolic solver, followed by a relaxation step to handle phase disequilibrium. After this initial effort, reduced forms of the Baer–Nunziato equations were developed to relieve the otherwise stiff pressure- and velocity-relaxation procedures. These methods are broadly used today. For example, Kapila et al. [29] developed a reduced five-equation model with single velocity and pressure, and its use has been widespread [30, 31, 32, 33, 34]. This is partly due to the simplicity of including additional physical phenomena like surface tension [35], mass-transfer [27], viscous effects [36], flashing [37], multiple components [38], and cavitation [39].

The drawback to the five-equation approach is that a single-pressure requires establishing a mixture EOS with a non-monotonic sound speed [40]. This results in numerical instability near strong shocks, difficulty maintaining volume fraction-positivity, and inaccurate wave transmission through the interface [41]. The five-equation model by Allaire et al. [30] is relatively robust, but it cannot model the creation and destruction of the interface due to the absence of a volume fraction source term. A six-equation model was developed to resolve this, which used a single velocity but phase-independent pressures [41]. Extensions of this work include phase-change effects [42, 43, 44, 45], viscous effects [46], and arbitrary numbers of components [38, 47].

Herein, we use the more extensive seven-equation model. The seven-equation model maintains phase-independent pressures so that it can represent with strong shocks and large density ratios: $\rho^{(2)}/\rho^{(1)},p^{(2)}/p^{(1)}\gg 1$ . There is no need to formulate a mixture EOS, so one can use phase-independent, arbitrary EOSs. It also retains phase-independent velocities, which allows representation of dispersed unresolved MPEs [8]. Considering this, the various features of SIM and DIM variants are summarized in Table 1.

After the initial development of the seven-equation model by Saurel and Abgrall [28], there have been relatively limited extensions and application use-cases [48, 49, 50, 51, 42]. This appears to be due to the complexity of formulating an efficient numerical method. In this work, we take a step towards addressing this problem. We start with the baseline numerical approach of Abgrall and Saurel [48] and provide novel extensions to develop a robust computational method for the seven-equation model that can handle strong shocks, arbitrary EOS, viscous effects, surface tension, multiple species, and reactions [52].

Novel contributions of this work:

1.

Various numerical methods have been developed for solving the BN PDE [28, 48, 49, 50, 51]. Still, the discrete equations method (DEM) [48] accurately models non-conservative terms at the interface. For 1D Riemann problems, DEM relieves the need to use a stiff pressure–velocity relaxation solver for pure fluid interfaces [48]. This effect is further studied in this work for more complex problems.
2.

Surface tension is required for modeling droplet breakup and deformation. Past works included it in the five-equation model [35] and a path-conservative framework for the seven-equation model [51] via source terms. We include it within the DEM approach via appropriate modifications to the Riemann solver.
3.

Viscous effects are necessary for applications such as drag computation of a deforming droplet or droplets interacting with turbulence. These have been previously considered in the DEM framework [36] for the five-equation model. This work extends it to the seven-equation model.
4.

The developed computational framework can deal with arbitrary and phase-independent EOS due to the pressure non-equilibrium, and this is demonstrated by using calorically perfect gas (CPG), thermally perfect gas (TPG), stiffened gas (SG), and Mie–Grüneisen (MG) EOS.
5.

Being able to handle multi-species mixtures for reactive flows, e.g., droplets in a detonation wave [53]. Each phase in this model includes an arbitrary number of species and reactions.
6.

The DIM interface can numerically diffuse across multiple cells. Considering this, an interface compression scheme, initially developed for the five-equation model [54] is modified to be used with the seven-equation model.

Table 1: Various numerical methods for modeling resolved MPEs and their features. P: Possible; NP: Not possible; PNA: Possible, not attempted. Bold text refers to cases considered in this work.

	SIM	5-eq DIM	6-eq DIM	7-eq DIM
$\rho^{(2)}/\rho^{(1)}$ $\gg 1,\quad$ $p^{(2)}/p^{(1)}$ $\gg 1$	Limited [18, 19, 20, 55]	Limited robustness [29]	P [41]	P [28, 48]
Arbitrary EOS	P [16]	NP	PNA	PNA
Viscous	P [55]	P [36]	P [38]	PNA
Surface tension	P [56]	P [35]	P [57]	P [51]
Flashing / trans-critical	NP	P [37, 39]	P [42, 58, 43]	P [42]
Evaporation	P [59, 18]	P [27]	P [42, 58, 43]	P [42]
Dispersed	NP	NP	NP	P [28]
Multi-species / reaction	P [18]	PNA	PNA	PNA
Interface Compression	Not needed	Used [54, 60]	Used [61]	PNA

Herein, Section 2 describes the seven-equation mathematical model, Section 3 formulates the numerical method, and Section 4 shows verification and validation cases.

2 Mathematical Formulation

2.1 Fine-grained description

We start from a fine-grained description of a multiphase flow before BN averaging. An indicator function for phase $\kappa=1,2$ is defined as $I^{\kappa}(x_{\textrm{i}},t)$ , which is $1$ if the location $x_{\textrm{i}}$ is in phase $\kappa$ at time $t$ , and $0$ otherwise. An evolution equation for $I^{\kappa}$ is

\frac{\partial I^{\kappa}}{\partial t}+u^{I}_{\textrm{i}}\frac{\partial I^{\kappa}}{\partial x_{\textrm{i}}}=0\quad(\textrm{i}=1,2,3),

(2.1)

where $u^{I}_{\textrm{i}}$ is the interface velocity. The single-phase Navier–Stokes equations accompany these for each phase as

\frac{\partial U}{\partial t}+\frac{\partial F_{\textrm{i},inv}}{\partial x_{\textrm{i}}}+\frac{\partial F_{\textrm{i},visc}}{\partial x_{\textrm{i}}}=\dot{S}_{\kappa}\quad({\textrm{i}}=1,2,3),

(2.2)

where,

U=\begin{pmatrix}1,\\ \rho,\\ \rho u_{\textrm{j}},\\ \rho E,\\ \rho Y_{\textrm{m}}\end{pmatrix},\,F_{{i},inv}=\begin{pmatrix}0,\\ \rho u_{\textrm{i}},\\ \rho u_{\textrm{i}}u_{\textrm{j}},\\ \rho Eu_{\textrm{i}},\\ \rho Y_{\textrm{m}}u_{\textrm{i}}\end{pmatrix}+\begin{pmatrix}0,\\ 0,\\ p\delta_{\textrm{ij}},\\ pu_{\textrm{i}},\\ 0\end{pmatrix},\,F_{{i},visc}=\begin{pmatrix}0,\\ 0,\\ \tau_{\textrm{i}\textrm{j}},\\ \tau_{\textrm{i}\textrm{j}}u_{\textrm{j}}+q_{\textrm{i}},\\ \rho Y_{\textrm{m}}v_{\textrm{im}}\end{pmatrix}.

(2.3)

where $U$ , $F_{\textrm{i},inv}$ , and $F_{\textrm{i},visc}$ are solution variables, inviscid fluxes and viscous fluxes, respectively. The total flux is $F_{\textrm{i}}=F_{\textrm{i},inv}+F_{\textrm{i},visc}$ , and $\rho$ , $u_{\textrm{j}},\ (\textrm{j}=1,2,3)$ , $p$ , and $Y_{\textrm{m}},\ (\textrm{m}=1,\dots,N_{s,\kappa})$ represent density, velocity, pressure, and mass-fractions for $\textrm{m}^{\textrm{th}}$ species, respectively. Note that $U(0)=1$ is redundant but corresponds to a volume fraction evolution equation after the averaging operation. We will show this later. The number of species that are tracked, $N^{\kappa}_{s}$ , can be different for different phases.

The total energy $E$ is related to the internal energy $e=h_{0,mix}+C_{v,mix}T$ as $E=e+1/2\ u_{\textrm{i}}u_{\textrm{i}}$ , where $C_{v,mix}$ is the mixture constant volume heat capacity and $h_{0,mix}$ is the mixture enthalpy of formation per mass. Pressure $p$ is computed as $P^{\kappa}(\rho,e,Y_{\textrm{m}})$ via an EOS. (EOS) The mixture formation enthalpy per mass is computed as $h_{0,mix}=\sum_{m=1}^{N_{s}^{\kappa}}h_{0,\textrm{m}}^{\kappa}Y_{\textrm{m}}$ with $h_{0,\textrm{m}}^{\kappa}$ being the formation enthalpy per mass of ${\textrm{m}}^{\textrm{th}}$ species. The mixture heat capacities are computed as

\displaystyle C_{v,mix}=\sum_{\textrm{m}=1}^{N_{s}^{\kappa}}C_{v,\textrm{m}}^{\kappa}Y_{\textrm{m}}\quad\text{and}\quad C_{p,mix}=\sum_{\textrm{m}=1}^{N_{s}^{\kappa}}C_{p,\textrm{m}}^{\kappa}Y_{\textrm{m}}.

(2.4)

The mixture specific enthalpy is $h_{mix}=h_{0,mix}+C_{p,mix}T$ . The species-specific heat capacities are $C_{p,\textrm{m}}^{\kappa}$ and $C_{v,\textrm{m}}^{\kappa}$ , and they can be curve-fits of the local temperature if needed. Both phases have independent EOS that determine $P^{\kappa}$ and thermodynamic properties $h_{0,\textrm{m}}^{\kappa}$ , $C_{p,\textrm{m}}^{\kappa}$ , $C_{v,\textrm{m}}^{\kappa}$ .

The viscous stress term $\tau_{\textrm{i}\textrm{j}}$ is computed from the shear stress tensor

\displaystyle\mathcal{T}^{\kappa}(u_{\textrm{i}})=\mu^{\kappa}\ {\partial u_{\textrm{i}}}/{\partial x_{j}}-2/3\mu^{\kappa}\ \delta_{\textrm{i}\textrm{j}}{\partial u_{\textrm{k}}}/{\partial x_{\textrm{k}}}.

(2.5)

The heat flux is modeled as

\displaystyle q_{\textrm{i}}=\lambda^{\kappa}\ {\partial T}/{\partial x_{\textrm{i}}}+\sum_{\textrm{m}=1}^{N_{s}^{\kappa}}h_{\textrm{m}}^{\kappa}Y_{\textrm{m}}v_{\textrm{im}}

(2.6)

where $v_{\textrm{im}}$ are species diffusive velocities, and $\mu^{\kappa}$ and $\lambda^{\kappa}$ are viscosity and thermal conductivity. Similar to the thermodynamic properties, the transport properties $\mathcal{T}^{\kappa}$ , $\mu^{\kappa}$ , $\lambda^{\kappa}$ , and species diffusivities are independent of the phases. The same form of the shear stress tensor $\mathcal{T}^{\kappa}$ is used for both the phases (corresponding to a Newtonian fluid). However, a different form can also be used, for example, to model a solid material.

The source terms $\dot{S}$ are $(0,0,g_{\textrm{i}},g_{\textrm{i}}u_{\textrm{i}},\dot{\omega}_{\textrm{m}})^{\text{T}}$ , where $g_{\textrm{i}}$ represents gravity and $\dot{\omega}_{\textrm{m}}^{\kappa}$ are the chemical reaction rates. We do not consider gravity hereon but include it for completeness. In the presence of chemical reactions, we also have the relations

\displaystyle\sum_{\textrm{{m}}}\dot{\omega}_{\textrm{m}}^{\kappa}=0\quad\text{and}\quad\sum_{\textrm{{m}}}Y_{\textrm{m}}v_{\textrm{mi}}=0.

(2.7)

The kinetics that models the reaction rates $\dot{\omega}_{\textrm{m}}^{\kappa}$ are also different and independent between the phases.

These equations are complemented by the jump conditions across the phases as

\left[(F_{\textrm{i}}-u^{I}_{\textrm{i}}U)\right]=\mathcal{J}^{U}_{\textrm{i}}\quad(\textrm{i}=1,2,3),

(2.8)

where

\mathcal{J}^{U}_{\textrm{i}}=\begin{pmatrix}0,\\ 0,\\ \sigma\beta n_{\textrm{i}}n_{\textrm{j}},\\ \sigma\beta u^{I}_{\textrm{i}},\\ 0\end{pmatrix},

(2.9)

and $[f]$ is the jump across the interface for any arbitrary field variable $f$ , e.g., $f^{(2)}-f^{(1)}$ . Here, since surface tension is considered, $n_{i}$ is the interface surface normal pointing from the liquid ( $\kappa=2$ ) to the gas ( $\kappa=1$ ) phase, and $\sigma$ and $\beta=\partial n_{\textrm{i}}/\partial x_{\textrm{i}}$ are the surface tension and the interface curvature, respectively. These jump conditions provide boundary conditions at the phase interfaces.

In the absence of any mass transfer across the phase interface and the presence of a thermal equilibrium due to heat transfer, these jump conditions transition into

[u_{\textrm{i}}]=0,\quad[p]=\sigma\beta-\left[\tau_{\textrm{ij}}n_{\textrm{i}}n_{\textrm{j}}\right],\quad[T]=0.

(2.10)

One would also maintain a Gibbs free energy equilibrium across the interface if mass transfer across the phase interface had been considered [42, 58, 44]. In the absence of mass-transfer the interface velocity $u^{I}_{\textrm{i}}$ is $u_{\textrm{i}}$ because $[u_{\textrm{i}}]=0$ .

2.2 Averaged equations

To derive the coarse-grained Baer–Nunziato formulation from the fine-grained description, an averaging operator $\langle\cdot\rangle$ is applied. The coarse-grained quantities are defined as $\mkern 1.5mu\overline{\mkern-1.5mu\alpha\vphantom{b}\mkern-1.5mu}\mkern 1.5mu^{\kappa}=\langle I^{\kappa}\rangle$ (volume fraction) and $\mkern 1.5mu\overline{\mkern-1.5muf\vphantom{b}\mkern-1.5mu}\mkern 1.5mu^{\kappa}\mkern 1.5mu\overline{\mkern-1.5mu\alpha\vphantom{b}\mkern-1.5mu}\mkern 1.5mu^{\kappa}=\langle I^{\kappa}f\rangle$ (phase-averaged quantities). Here the operator $\overline{\cdot}$ is spatial filtering as a result of the phase-averaging. However, since both the flow features and the interface are considered resolved in this work, $\mkern 1.5mu\overline{\mkern-1.5muf\vphantom{b}\mkern-1.5mu}\mkern 1.5mu^{\kappa}$ is represented as $f^{\kappa}$ from here onward. The generalized seven-equation formulation is derived from Eq. (2.1) and Eq. (2.2), by applying $\langle\cdot\rangle$ as

\begin{split}\underbrace{\frac{\partial\alpha^{\kappa}{U}^{\kappa}}{\partial t}}_{\textrm{I}}+\underbrace{\frac{\partial{\alpha}^{\kappa}{F}^{\kappa}_{\textrm{i},inv}}{\partial x_{\textrm{i}}}}_{\textrm{II}}+\underbrace{\frac{\partial{\alpha}^{\kappa}{F}^{\kappa}_{\textrm{i},visc}}{\partial x_{\textrm{i}}}}_{\textrm{III}}=&\underbrace{\left\langle\left(F^{\kappa}_{\textrm{i},inv}-u_{\textrm{i}}^{\kappa}U^{\kappa}\right)\frac{\partial I^{\kappa}}{\partial x_{\textrm{i}}}\right\rangle}_{\textrm{IV}}\\ &+\underbrace{\left\langle F^{\kappa}_{\textrm{i},visc}\frac{\partial I^{\kappa}}{\partial x_{\textrm{i}}}\right\rangle}_{\textrm{V}}+\underbrace{\left\langle U^{\kappa}(u^{\kappa}_{\textrm{i}}-u^{I}_{\textrm{i}})\frac{\partial I^{\kappa}}{\partial x_{\textrm{i}}}\right\rangle}_{\textrm{VI}}+\underbrace{\alpha^{\kappa}\dot{S}^{\kappa}}_{\textrm{VII}}.\end{split}

(2.11)

Here, term I is the time-derivative of flow-field variables, terms II and III are conservative flues, terms IV and V are non-conservative terms that are active only at the multiphase interfaces, and they are responsible for interphase exchanges. Term VI is mass transfer across the phase boundary, and they are neglected here but will be a part of future work. Term VII are source terms.

Regardless of whether the size of the MPE is larger or smaller than the filter-size of the averaging operator $\langle\cdot\rangle$ (implicitly same as the finite-volume grid), the same Eq. (2.11) can be used for modeling either dispersed or resolved multiphase flows. Numerical modeling of terms I, II, III, and VII is straightforward and they remain the same for either resolved or dispersed flows. IV and V are unclosed, and their modeling requires a specific focus. In the absence of mass transfer, these are typically modeled for the seven-equation model as [28, 51, 36, 42]

\begin{gathered}\left(IV+V\right)^{(1)}=I^{(1)}+R^{(1)}_{p}+R^{(1)}_{u}+R^{(1)}_{T},\\ I^{(1)}=\begin{pmatrix}-u_{\textrm{i}}^{I}\frac{\partial\alpha^{(1)}}{\partial x_{\textrm{j}}},\\ 0,\\ ((p^{I}-\sigma\beta\delta^{(1)})\delta_{\textrm{ij}}+\tau_{\textrm{ij}}^{I})\frac{\partial\alpha^{(1)}}{\partial x_{\textrm{j}}},\\ \left((p^{I}-\sigma\beta\delta^{(1)})u_{\textrm{i}}^{I}+\tau_{\textrm{ij}}^{I}u^{I}_{\textrm{i}}+q_{\textrm{i}}^{I}\right)\frac{\partial\alpha^{(1)}}{\partial x_{\textrm{i}}},\\ 0\end{pmatrix},\\ R^{(1)}_{p}=\begin{pmatrix}\theta^{p}\Delta p,\\ 0\\ 0,\\ (p^{I}-\sigma\beta\delta^{(1)})\theta^{p}\Delta p,\\ 0\end{pmatrix},R^{(1)}_{u}=\begin{pmatrix}0,\\ 0,\\ \theta^{u}\Delta u_{i},\\ u_{\textrm{i}}^{I}\theta^{u}\Delta u_{\textrm{i}}\\ 0\end{pmatrix},R^{(1)}_{T}=\begin{pmatrix}\frac{1}{\kappa^{I}}\theta^{T}\Delta T,\\ 0,\\ 0,\\ \theta^{T}\Delta T,\\ 0\end{pmatrix}\end{gathered}

(2.12)

where $\Delta p=p^{(1)}-p^{(2)}+\sigma\beta$ , $\Delta u_{i}=u^{(1)}_{i}-u^{(2)}_{i}$ , and $\Delta T=T^{(1)}-T^{(2)}$ are differences of pressure, velocity, and temperature between the phases. Furthermore, $p^{I}$ , $u_{{i}}^{I}$ (modeling discussed in [28]), $\tau_{{ij}}^{I}$ , $q_{{i}}^{I}$ (modeling discussed in [36]), and $\kappa^{I}$ (modeling discussed in [42]) are interface-averaged terms within the filter $\langle\cdot\rangle$ , and since they are unknown, they typically have to be modeled as listed in the references listed. For instance, $p^{I}$ , $u^{I}$ have been modeled in the past as $p^{(1)}$ and $u^{(2)}$ (with $\kappa=1$ being the gas phase and $\kappa=2$ as the liquid phase) for gas-liquid flows. In this work, however, the discrete equations method (DEM) is used [48], and its use circumvents the need for this modeling by solving Riemann problems at each discrete interface instead. A similar form as in Eq. (2.12) is also used for the phase $(2)$ , and it is not repeated here for brevity. As noted before, assuming that the phase $(1)$ is a gas-phase the phase $(2)$ is a liquid-phase, the surface tension force $\sigma\beta$ acts from the gas-phase towards the liquid-phase, and as a result, in Eq. (2.12), $\delta^{(1)}=1$ for the gas-phase, whereas $\delta^{(2)}=0$ for the liquid-phase [51].

The terms containing $\theta^{p}$ , $\theta^{u}$ , and $\theta^{T}$ are referred to as the relaxation terms for pressure, velocity, and temperature, respectively. The pressure relaxation results in volume exchange between the phases, the velocity relaxation results in interphase momentum exchanges, e.g., drag, and the temperature relaxation results in heat exchange between the phases. For resolved interfaces, these are typically modeled as $\theta^{p},\theta^{u},\theta^{T},\theta^{f}\rightarrow\infty$ , whereas, for dispersed phases, they are modeled using finite-time processes such as drag and heat transfer [28]. Enforcing this condition analytically would have resulted in a reduced five-/six-equation model [29, 41]. Past works [28, 42] had to use a stiff relaxation solver with the seven-equation model to enforce $\theta^{p},\theta^{u}\rightarrow\infty$ . However, the use of DEM circumvented this as it was shown for 1D canonical cases in [48]. This is further confirmed for more complex test cases in this work (Section 4 for more details).

2.3 Material properties

The governing equations are accompanied by an equation of state ( $P^{\kappa}$ ) and material properties for each phase. For demonstration, stiffened gas (SG, includes parameters $p_{0}^{\kappa}$ and $\gamma^{\kappa}$ ) [28], calorically perfect gas (CPG, has parameter $\gamma^{\kappa}$ ), thermally perfect gas (TPG, requires property curve-fits for individual species), and Mie–Grüneisen (MG, requires properties for individual species) [62] EOS are used in this work. These and the corresponding parameters are well-described in the available literature, and their forms are not repeated here for brevity. The corresponding parameters for various materials are described separately for each test as required in Section 4.

A Sutherland transport model is used for the gas-phase viscosity ( $\mu^{(1)}$ ). The liquid- viscosity ( $\mu^{(2)}$ ) and surface tension $\sigma$ are assumed to be constants, although it is possible to use curve-fits for those as well. The thermal conductivity $\lambda^{\kappa}$ is computed from the viscosity $\mu^{\kappa}$ using a constant Prandtl number ( $Pr=0.739$ ) as $\lambda^{\kappa}=C_{p,mix}^{\kappa}\mu^{\kappa}/Pr$ . The species diffusion is modeled using a unity Lewis number ( $Le$ ) assumption, with $Le^{\kappa}_{\textrm{m}}=\lambda^{\kappa}/\left(\rho^{\kappa}D^{\kappa}_{\textrm{m}}C_{p,\textrm{m}}^{\kappa}\right)=1$ . For the reactive simulations shown in this work, following previous work on hydrogen-air detonations [63], a 7-step, 6-species chemical mechanism [64] is used to model finite-rate kinetics.

3 Numerical Method

A three-dimensional finite-volume computational framework is used to solve Eq. (2.11). The equations are implemented in LESLIE, a parallel multi-physics solver, which has been used extensively in the past for multiphase reacting flow simulations [8, 65, 66]. The numerical approach here primarily follows published literature [48, 28, 51, 54, 36], and every detail is not repeated for brevity. However, the following novel contributions are made, and details about them are provided in the following sections.

1.

Implementation of the DEM [48] for multi-block structured curvilinear grids.
2.

Modification in HLLC [67] Riemann solver within the DEM to include surface tension.
3.

Extension of viscous effects within the DEM [36] for multiple dimensions.
4.

Improved robustness of the stiff relaxation solver [28] for dealing with arbitrary EOS.
5.

Application of an interface compression scheme [54] for the seven-equation model.

3.1 Baseline approach

A multi-block structured grid is used, wherein $\mathcal{C}_{ijk}$ is a computational cell with $i$ , $j$ , and $k$ as the indices in the three curvilinear directions. The corresponding cell volume is denoted as $V_{ijk}$ , and $S_{l,ijk}$ is used to indicate the cell faces. Note that these $i,\ j,\ k$ are different from the i, j, k that we used in Section 2 for the Einstein notation. Let $(\xi,\eta,\zeta)$ are the computational directions of increasing $i$ , $j$ and $k$ , and $\partial\xi/\partial x_{i}$ , $\partial\eta/\partial x_{i}$ , $\partial\zeta/\partial x_{i}$ are the corresponding grid metrics. Quantities computed at the face-centers are denoted with sub-scripts $i+1/2$ , $j+1/2$ and $k+1/2$ . For example, $S_{i+\frac{1}{2}jk}$ and ${\mathcal{N}}_{\textrm{i},i+\frac{1}{2}jk}$ are the face-area and the face-normal, respectively, between the cells $\mathcal{C}_{i+1jk}$ and $\mathcal{C}_{ijk}$ .

A volume-integrated form of Eq. (2.11) is numerically solved on this grid. The cell-centered quantity to be updated at each time-step is $\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk}$ , and it is also referred to as $Q_{ijk}$ . The corresponding time-update is denoted as $d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk}$ or $dQ_{ijk}$ . The increment $dQ$ (or $d\left(\alpha_{\kappa}U_{\kappa}\right)$ ) has contributions due to various terms in Eq. (2.11) and they are denoted using subscripts as $dQ_{\textrm{II}},dQ_{\textrm{III}}$ , etc.

The computation is advanced in time using an explicit two-stage MacCormack approach. The time-step $\Delta t$ is dynamically computed considering convective, acoustic, and diffusive propagation speeds of both the phases as

\Delta t^{\kappa}_{ijk}=CFL\cdot{V}_{ijk}\bigg{/}\left(\frac{1}{6}\sum_{l=1}^{6}{\mid v^{\kappa}_{l,ijk}\mid S_{{l,ijk}}+c^{\kappa}\overline{S}_{ijk}+\frac{2\gamma^{\kappa}\nu^{\kappa}}{Pr}^{\kappa}\frac{{\overline{S}_{ijk}}^{2}}{V_{ijk}}}\right),\quad\Delta t=\underset{ijk}{\textrm{min}}\ \textrm{min}\left(\Delta t^{(1)}_{ijk},\Delta t^{(2)}_{ijk}\right)

(3.1)

where, $S_{l,ijk}$ are face area and $v^{\kappa}_{l,ijk}=\sum_{\textrm{i}=1}^{3}u^{\kappa}_{\textrm{i},ijk}n_{\textrm{i},ijk}$ is the normal velocity on the $l^{\textrm{th}}$ face of the cell $\mathcal{C}_{ijk}$ . The averaged face area for cell $\mathcal{C}_{ijk}$ is given by $\overline{S}_{ijk}=\frac{1}{6}\sum_{l=1}^{6}S_{l,ijk}$ , and $\gamma^{\kappa}$ , $\nu^{\kappa}=\mu^{\kappa}/\rho^{\kappa}$ , and $Pr^{\kappa}$ are the heat capacity, the kinematic viscosity and the Prandtl number, respectively. The $CFL$ is chosen as $0.5$ for stability.

3.2 Inviscid fluxes

In the finite-volume scheme, the inviscid increments due to the terms $II$ and $IV$ in Eq. (2.11) are provided as

d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}=d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}^{i}+d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}^{j}+d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II/IV}^{k},

(3.2)

where the superscripts $i,j,k$ denote contributions from each direction. The fluxes in each direction are then computed using the DEM. Their computation is illustrated here only for the $i$ -direction, but the same forms are used in the other directions as well.

3.2.1 Discrete equations method (DEM)

Instead of directly solving for the conservative and the non-conservative terms, DEM divides the problem into three Riemann problems at each cell-face. A schematic for this is shown in Fig. 1(a), where the three Riemann problems at each face are shown corresponding to three possibilities of the two-phase boundaries, i.e., (1)–(1) (R1), (1)–(2) (R2) and (2)–(2) (R3). Note that piece-wise constant fields are shown in each cell here for illustration. However, we go up to second-order accuracy with piece-wise linear fields for the numerical method implemented. On each face, left- and right-interpolations are handled using a first/second-order MUSCL approach [23], and the interpolated values are denoted with subscripts $i+1/2,l$ , $i+1/2,r$ , etc.

Refer to caption — (a) Riemann sub-problems

3.2.2 Riemann solution

Solutions to the Riemann sub-problems (1)–(1), (1)–(2), and (2)–(2) are now needed. An HLLC approximate Riemann solver [67] has been used traditionally [28, 48]. However, modifications to it are required, including the surface tension. Note that for a (1)–(2) Riemann problem, a discontinuity in pressure would be maintained across the contact surface such as $p^{(2)}-p^{(1)}=p^{\sigma}$ , whereas, the solution should stay unaffected for (1)–(1) or (2)–(2) problems. A typical 1D Riemann problem is depicted in Fig. 2, which could be either of the three sub-problems R1, R2, or R3. The state vector $\mathcal{U}=(1,\rho,\rho u_{\textrm{i}},\allowbreak\rho E,\rho Y_{\textrm{m}})^{T}$ is set equal to either $U^{(1)}$ or $U^{(2)}$ (on both the left and the right side of the face) depending on whether it belong to the phase (1) or (2). In the case where the left and the right state vectors $\mathcal{U}^{L}$ or $\mathcal{U}^{R}$ belong to different phases, i.e., (1)–(2) Riemann problem, the contact wave separates the phases as the problem evolves.

Figure 2: A sample Riemann problem is shown here, where

\mathcal{U}^{L}

and

\mathcal{U}^{R}

are the initial left- and right-states, and

\mathcal{U}^{L*}

\mathcal{U}^{R*}

are the generated intermediate states. The waves on the extreme right/left could be a shock or an expansion wave (solid line), and the wave in the middle is a contact wave (dashed line), whose speed is denoted using

c^{*}

A sharp jump is assumed across both the left and the right shock/expansion waves in the HLLC approximate Riemann solver [67], and these waves travel with the speeds $c^{L}$ and $c^{R}$ , respectively. These are the acoustic propagation speeds corresponding to the state vectors $\mathcal{U}^{L}$ and $\mathcal{U}^{R}$ , respectively. Considering this, we have the following relations across these jumps.

\mathcal{F}^{K*}-\mathcal{F}^{K}=c^{K}(\mathcal{U}^{K*}-\mathcal{U}^{K}),\ K\in\{L,R\},\quad\mathcal{F}^{R*}-\mathcal{F}^{L*}+\mathcal{F}^{\sigma}=c^{*}(\mathcal{U}^{R*}-\mathcal{U}^{L*}).

(3.3)

where $\mathcal{F}=\left(0,\rho v_{n},\rho u_{\textrm{i}}v_{n}+pn_{\textrm{i}},\rho Ev_{n}+pv_{n},\rho Y_{\textrm{m}}v_{n}\right)$ is the flux vector, and that is defined for each state in the Riemann problem. As noted before, $v_{n}=\sum_{\textrm{i}=1}^{3}u_{\textrm{i}}n_{\textrm{i}}$ is the velocity normal to the corresponding cell-face. The term $\mathcal{F}^{\sigma}=\left(0,0,p^{\sigma}n_{\textrm{i}},p^{\sigma}c^{*},0\right)$ is the flux due to the surface tension at the interface boundary (at the contact wave in this case).

The pressures and the normal velocities across the interface (contact wave) are related as $p^{L*}=p^{R*}+p^{\sigma},\quad v^{L*}_{n}=v^{R*}_{n}$ , where, $p^{\sigma}=n^{LR}\sigma\beta$ is the pressure jump across the interface due to the surface tension. The normal vector $n^{LR}$ is $1$ if the left phase is the liquid phase and the right phase is the gas phase, i.e., (2)–(1), $-1$ if the left phase is the gas phase and the right phase is the liquid-phase, i.e., (1)–(2), and $0$ if the phases on both the side of the interface are the same, i.e., (1)–(1) or (2)–(2). Combining these equations and following the derivation in Toro [67] with these additional modifications due to the surface tension force, the speed of the contact wave and the intermediate states are given as

c^{*}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)=\frac{p^{R}-p^{L}+p^{\sigma}+\rho^{L}v^{L}_{n}(c^{L}-v^{L}_{n})-\rho^{R}v^{R}_{n}(c^{R}-v^{R}_{n})}{\rho^{L}(c^{L}-v^{L}_{n})-\rho^{R}(c^{R}-v^{R}_{n})},\quad\mathcal{F}^{*}=\rho^{K}\left(\frac{c^{K}-v^{K}_{n}}{c^{K}-c^{*}}\right)

(3.4)

\mathcal{U}^{K*}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)=\begin{pmatrix}1,\\ \mathcal{F}^{*},\\ \mathcal{F}^{*}\left(u_{\textrm{i}}^{K}+(c^{*}-v^{K}_{n})n_{\textrm{i}}\right),\\ \mathcal{F}^{*}E^{K}+\mathcal{F}^{*}\left(c^{*}-v^{K}_{n}\right)\left(c^{*}+\frac{p^{K}}{\rho^{K}(c^{K}-v^{K}_{n})}\right),\\ \mathcal{F}^{*}Y_{\textrm{m}}^{K}\end{pmatrix}.

(3.5)

Once the intermediate states are known, the conservative Riemann flux at the cell face is computed as

\begin{split}F^{cons}_{inv}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)=\begin{cases}\mathcal{F}(\mathcal{U}^{L})\quad\textrm{if }0\leq c^{L}\\ \mathcal{F}(\mathcal{U}^{L*})\quad\textrm{if }c^{L}\leq 0\leq c^{*}\\ \mathcal{F}(\mathcal{U}^{R*})\quad\textrm{if }c^{*}\leq 0\leq c^{R}\\ \mathcal{F}(\mathcal{U}^{R})\quad\textrm{if }c^{R}\leq 0\\ \end{cases}\end{split}

(3.6)

Note that the conventional single-phase Riemann solution is recovered by letting $p^{\sigma}=0$ .

3.2.3 Inviscid conservative flux (Terms $II$ )

Based on the Riemann solutions that are constructed for the sub-problems (1)–(1), (1)–(2), and (2)–(2), we can now compute the inviscid conservative fluxes. To illustrate this process with DEM, a schematic is shown in Fig. 1 (b). For this example, during the Riemann problem evolution, the interface that was initially at $i+1/2$ has now moved within the cell on the right based on the speed of the contact wave $c^{*}$ . $F^{cons}_{inv}$ are the conservative fluxes for the Riemann problems (1)–(1), (1)–(2), and (2)–(2) that act at the $i+1/2$ cell-face.

Considering the solutions from the three Riemann problems, the averaged i-directional conservative fluxes (term II in Eq. (2.11)) at face $i+1/2$ for phase (1) are given using the DEM as

\begin{split}\left\langle I^{(1)}\left(F^{(1)}_{inv}\right)\right\rangle_{i+\frac{1}{2}}=&\quad\textrm{min}\left(\alpha^{(1)}_{i+\frac{1}{2},l},\alpha^{(1)}_{i+\frac{1}{2},r}\right)F^{cons}_{inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(1)}\right)\\ &+\textrm{max}\left(\alpha^{(1)}_{i+\frac{1}{2},l}-\alpha^{(1)}_{i+\frac{1}{2},r},0\right)\left(\beta_{i+\frac{1}{2}}^{(1,2)}\right)^{+}F^{cons}_{inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(2)}\right)\\ &+\textrm{max}\left(\alpha^{(2)}_{i+\frac{1}{2},l}-\alpha^{(2)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F^{cons}_{inv}\left(U_{i+\frac{1}{2},l}^{(2)},U_{i+\frac{1}{2},r}^{(1)}\right),\end{split}

(3.7)

where $F^{cons}_{inv}$ are the fluxes based on the Riemann solutions as discussed earlier, and their coefficients are the probabilities of encountering a gas-gas, a liquid-liquid, a gas-liquid, or a liquid-gas interface on a fine-grained discrete level. Further details and justifications are provided in the original DEM work [48], and they are not repeated here for brevity. Coefficients $\beta_{i+{1}/{2}}^{(p,q)}$ are computed as $\textrm{sign}\left(c^{*}\left(U_{i+1/2,l}^{p},U_{i+1/2,r}^{q}\right)\right)$ , and $\left(\beta\right)^{+}=\mathrm{max}\left(\beta,0\right)$ .

After computing the fluxes, the increment $d\left(\alpha^{\kappa}U^{\kappa}\right)$ due to the $i-$ directional conservative inviscid flux (term II in Eq. (2.11)) is computed as

\begin{split}Vd\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,II}^{i}=&S_{i+\frac{1}{2}jk}\left\langle I^{\kappa}\left(F^{\kappa}_{inv}\right)\right\rangle_{i+\frac{1}{2}jk}-S_{i+\frac{1}{2}jk}\left\langle I^{\kappa}\left(F^{\kappa}_{inv}\right)\right\rangle_{i-\frac{1}{2}jk}.\end{split}

(3.8)

3.2.4 Inviscid non-conservative flux (Terms $IV$ )

The non-conservative fluxes are referred to as the ones that act at a (1)-(2) interface (contact) but within a finite-volume cell once the multiphase interface has moved into it, e.g., the flux $F^{lag}_{L/R}$ acting on the interface that has moved into the cell $i+1$ in Fig. 1 (b). The subscripts $L/R$ correspond to whether the flux is acting on the material on the left or right side of the interface. Following [48], second-order $i-$ directional update for phase $(1)$ due to the inviscid non-conservative flux (term IV in Eq. (2.11)) is given as

\begin{split}Vd\left(\alpha^{(1)}U^{(1)}\right)_{ijk,IV}^{i}&=\textrm{max}\left(\alpha^{(1)}_{i+\frac{1}{2},l}-\alpha^{(1)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F^{lag}_{L,inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(2)}\right)S_{i+\frac{1}{2}}\\ &-\textrm{max}\left(\alpha^{(2)}_{i+\frac{1}{2},l}-\alpha^{(2)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F^{lag}_{R,inv}\left(U_{i+\frac{1}{2},l}^{(2)},U_{i+\frac{1}{2},r}^{(1)}\right)S_{i+\frac{1}{2}}\\ &+\textrm{max}\left(\alpha^{(1)}_{i-\frac{1}{2},l}-\alpha^{(1)}_{i-\frac{1}{2},r},0\right)\left(\beta_{i+\frac{1}{2}}^{(1,2)}\right)^{+}F^{lag}_{L,,inv}\left(U_{i+\frac{1}{2},l}^{(1)},U_{i+\frac{1}{2},r}^{(2)}\right)S_{i-\frac{1}{2}}\\ &-\textrm{max}\left(\alpha^{(2)}_{i-\frac{1}{2},l}-\alpha^{(2)}_{i-\frac{1}{2},r},0\right)\left(\beta_{i-\frac{1}{2}}^{(2,1)}\right)^{+}F^{lag}_{R,inv}\left(U_{i-\frac{1}{2},l}^{(2)},U_{i-\frac{1}{2},r}^{(1)}\right)S_{i-\frac{1}{2}}\\ &+\textrm{max}\left(\alpha^{(1)}_{i-\frac{1}{2},r}-\alpha^{(1)}_{i+\frac{1}{2},l},0\right)F^{lag}_{L}\left(U_{i}^{(1)},U_{i}^{(2)}\right)\frac{1}{2}\left(S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}\right)\\ &-\textrm{max}\left(\alpha^{(2)}_{i-\frac{1}{2},r}-\alpha^{(2)}_{i+\frac{1}{2},l},0\right)F^{lag}_{R}\left(U_{i}^{(2)},U_{i}^{(1)}\right)\frac{1}{2}\left(S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}\right)\\ &+\langle N_{int}\rangle\left(F_{R}^{lag}\left(U_{i}^{(2)},U_{i}^{(1)}\right)-F_{L}^{lag}\left(U_{i}^{(1)},U_{i}^{(2)}\right)\right)\frac{1}{2}\left(S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}\right).\end{split}

(3.9)

where $F^{lag}_{L/R,inv}(U^{L},U^{R})$ are fluxes between the phases at the interface, and they are computed as part of the previously defined (1)-(2) Riemann solution as

\begin{split}F^{lag}_{L,inv}(U^{L},U^{R})=\mathcal{F}(\mathcal{U}^{L*})-c^{*}\left(\mathcal{U}^{L},\mathcal{U}^{L}\right)\mathcal{U}^{L*},\\ F^{lag}_{R,inv}(U^{L},U^{R})=\mathcal{F}(\mathcal{U}^{R*})-c^{*}\left(\mathcal{U}^{L},\mathcal{U}^{R}\right)\mathcal{U}^{R*}\end{split}

(3.10)

Note that based on the approximate Riemann solution derived earlier, we have $F^{lag}_{L}-F^{lag}_{R}=F^{\sigma}$ for flows with surface tension.

The first two terms in Eq. (3.9) correspond to the non-conservative fluxes that act within the cell $i$ as a result of a multiphase interface that would have entered through the cell-face $i+1/2$ . The following two terms correspond to an interface that would have entered through the cell-face $i-1/2$ . The subsequent two terms are for second-order accuracy when piecewise linear interpolations are used within a cell. In the last term, $\langle N_{int}\rangle$ is the average number of internal interfaces within the cell, and this term was related to the relaxation processes within a cell, which would lead to local equilibrium. For resolved interfaces, this could be modeled via the stiff relaxation solvers as described later in Section 3.5. The reader can find more details about the DEM elsewhere [48], and they are not repeated here for brevity.

3.3 Viscous fluxes

Terms III and V in Eq. (2.11) are the viscous conservative and non-conservative updates. Abgrall and Rodio [36] extended the original DEM framework [48] for viscous fluxes. However, they only did this for 1D flows, and here it is also used for 2D/3D. The viscous updates $d\left(\alpha^{\kappa}U^{\kappa}\right)_{ijk,III}$ and $d\left(\alpha_{\kappa}U_{\kappa}\right)_{ijk,V}$ follow the same formula as Eq. (3.2) for their directional contributions from $i/j/k$ . The flux computations are described here only for the $i$ -direction, but the other directions use a similar form.

3.3.1 Single-phase viscous fluxes

A Riemann problem is inherently defined for 1D inviscid flows, and therefore, it is not possible to completely replicate the inviscid computations shown earlier for the viscous terms. The computation of single-phase viscous fluxes $F^{\kappa}_{i,visc}$ is discussed first in this section, and their compilation into the two-phase conservative and non-conservative terms is described in the following sections.

Centered differencing and centered averaging are used for computing the terms $\partial_{x_{\textrm{i}}}f\,\vline_{\ i+1/2jk}$ , $f_{i+1/2jk}$ , etc. The physical derivatives $\partial_{x_{\textrm{i}}}f$ , e.g., $\partial_{x_{\textrm{j}}}u_{\textrm{i}}$ , $\partial_{x_{\textrm{i}}}T$ , etc., are obtained from the computational derivatives as

\frac{\partial f}{\partial x_{\textrm{i}}}=\frac{\partial f}{\partial\zeta}\frac{\partial\zeta}{\partial x_{\textrm{i}}}+\frac{\partial f}{\partial\eta}\frac{\partial\eta}{\partial x_{\textrm{i}}}+\frac{\partial f}{\partial\xi}\frac{\partial\xi}{\partial x_{\textrm{i}}},

(3.11)

where $\partial_{\zeta}f$ , $\partial_{\eta}f$ , and $\partial_{\xi}f$ are computed from the flow-variables as

\begin{split}&\quad\quad\frac{\partial f}{\partial\zeta}\vline_{\ i+\frac{1}{2}jk}=f_{i+1jk}-f_{ijk},\\ &\frac{\partial f}{\partial\eta}\vline_{\ i+\frac{1}{2}jk}=\frac{1}{2}\left(f_{i+\frac{1}{2}j+1k}-f_{i+\frac{1}{2}j-1k}\right),\\ &\frac{\partial f}{\partial\xi}\vline_{\ i+\frac{1}{2}jk}=\frac{1}{2}\left(f_{i+\frac{1}{2}jk+1}-f_{i+\frac{1}{2}jk-1}\right).\end{split}

(3.12)

In addition, any face-averaged quantities, e.g., $\alpha^{\kappa}$ , $\mu^{\kappa}$ , etc., are computed as

f_{i+\frac{1}{2}jk}=\frac{1}{2}(f_{i+1jk}+f_{ijk}).

(3.13)

With these, the $i$ -directional single-phase viscous fluxes are computed as

\begin{split}F^{\kappa}_{visc,i+\frac{1}{2}jk}&=F^{\kappa}_{\textrm{i},visc}\left(\mu^{\kappa}_{i+\frac{1}{2}jk},\dots,\frac{\partial u^{\kappa}_{\textrm{i}}}{\partial x_{i}}\vline\ _{i+\frac{1}{2}jk},\dots\right)\cdot n_{\textrm{i},i+\frac{1}{2}jk}.\end{split}

(3.14)

3.3.2 Viscous conservative flux (Terms $III$ )

For the DEM computation of the viscous conservative fluxes, similar forms as Eq. (3.7) and Eq. (3.8) are used. However, instead of using the $F^{cons}_{inv}$ from the Riemann solutions as before, the terms $F^{\kappa}_{visc}$ at the cell-faces (e.g., $i+1/2,j,k$ ) are now computed using Eq. (3.14). The DEM viscous fluxes are then given as

\begin{split}\left\langle I^{(1)}\left(F^{(1)}_{visc}\right)\right\rangle_{i+\frac{1}{2}}=&\textrm{min}\left(\alpha^{(1)}_{i+\frac{1}{2},l},\alpha^{(1)}_{i+\frac{1}{2},r}\right)F_{visc}^{(1)}\\ &+\textrm{max}\left(\alpha^{(1)}_{i+\frac{1}{2},l}-\alpha^{(1)}_{i+\frac{1}{2},r},0\right)\left(\beta_{i+\frac{1}{2}}^{(1,2)}\right)^{+}F_{visc}^{(1)}\\ &+\textrm{max}\left(\alpha^{(2)}_{i+\frac{1}{2},l}-\alpha^{(2)}_{i+\frac{1}{2},r},0\right)\left(-\beta_{i+\frac{1}{2}}^{(2,1)}\right)^{+}F_{visc}^{(1)},\end{split}

(3.15)

and a similar form is used for $\kappa=2$ as well.

3.3.3 Viscous non-conservative (Terms $V$ )

The non-conservative viscous terms follow the same form as Eq. (3.9), however, instead of $F^{lag}_{L/R,inv}$ coming from the 1D inviscid Riemann problems, $F^{lag}_{L/R,visc}$ are modeled as a function of the single-phase viscous terms $F^{\kappa}_{visc}$ computed at cell-faces (e.g., $i+1/2,j,k$ ) as

\begin{gathered}F^{lag}_{L,visc}(U^{(1)}_{i+\frac{1}{2},l},U^{(2)}_{i+\frac{1}{2},l})=F^{lag}_{R,visc}(U^{(2)}_{i+\frac{1}{2},l},U^{(1)}_{i+\frac{1}{2},l})=F^{(1)}_{visc,i+\frac{1}{2}},\\ F^{lag}_{L,visc}(U^{(1)}_{i},U^{(2)}_{i})=F^{lag}_{R,visc}(U^{(2)}_{i},U^{(1)}_{i})=\frac{S_{i+\frac{1}{2}}F^{(1)}_{visc,i+\frac{1}{2}}+S_{i-\frac{1}{2}}F^{(1)}_{visc,i-\frac{1}{2}}}{S_{i+\frac{1}{2}}+S_{i-\frac{1}{2}}}.\\ \end{gathered}

(3.16)

where $\kappa=1$ is assumed as the gas-phase, and therefore, it is assumed that $F^{I}_{visc}=F^{(1)}_{visc}$ . This is similar to the previous viscous DEM approach of Abgrall and Rodio [36].

3.4 Curvature computation

The interface curvature $\beta$ is required to compute $p^{\sigma}=\sigma\beta$ at the cell-centers $(i,j,k)$ . Once the curvature is computed at the cell centers, it is then interpolated at the cell faces using the MUSCL approach in the same way as the other field quantities for flux computation. The curvature for a field $f$ is defined as

\beta=\frac{\partial}{\partial x_{i}}\left(\frac{\partial f}{\partial x_{i}}\bigg{/}\bigg{|}\frac{\partial f}{\partial x_{i}}\bigg{|}\right).

(3.17)

A previous approach for curvature computation is followed that used a vertex-based staggering [51, 35], which is illustrated in 2D. The gradient of $f$ at a grid vertex is first computed as

\begin{split}\frac{\partial f}{\partial x_{1}}_{i+\frac{1}{2},j+\frac{1}{2}}&=\frac{1}{2}\left(\frac{f_{i+1,j+1}-f_{i,j+1}}{\Delta x}+\frac{f_{i+1,j}-f_{i,j}}{\Delta x}\right),\\ \frac{\partial f}{\partial x_{2}}_{i+\frac{1}{2},j+\frac{1}{2}}&=\frac{1}{2}\left(\frac{f_{i+1,j+1}-f_{i+1,j}}{\Delta x}+\frac{f_{i,j+1}-f_{i,j}}{\Delta x}\right).\end{split}

(3.18)

where $\Delta x$ and $\Delta y$ are the grid sizes in the $i$ and the $j$ directions. Based on this gradient, a unit vector normal to the interface is computed as

n_{\textrm{i},i+\frac{1}{2},j+\frac{1}{2}}=\begin{cases}\frac{\frac{\partial f}{\partial x_{\textrm{i}}}_{i+\frac{1}{2},j+\frac{1}{2}}}{\bigg{|}\bigg{|}\frac{\partial f}{\partial x_{\textrm{i}}}_{i+\frac{1}{2},j+\frac{1}{2}}\bigg{|}\bigg{|}},\ \ \ \ \ \text{if}\ \ \ \ \ \bigg{|}\bigg{|}\frac{\partial f}{\partial x_{\textrm{i}}}_{i+\frac{1}{2},j+\frac{1}{2}}\bigg{|}\bigg{|}>0,\\ 0,\ \ \ \ \ \ \text{otherwise},\end{cases}

(3.19)

where, $||\cdot||$ is the magnitude of the vector. The normal vector is then interpolated at the cell-faces as

n_{\textrm{i},i+\frac{1}{2},j}=\frac{1}{2}\left(n_{\textrm{i},i+\frac{1}{2},j+\frac{1}{2},}+n_{\textrm{i},i+\frac{1}{2},j-\frac{1}{2}}\right),

(3.20)

and finally, the curvature is computed at the cell-centers as

\kappa_{i,j}=\frac{n_{1,i+\frac{1}{2},j}-n_{1,i-\frac{1}{2},j,}}{\Delta x}+\frac{n_{2,i,j+\frac{1}{2},}-n_{i,j-\frac{1}{2}}}{\Delta y}.

(3.21)

The previous approach [51] used $f=\alpha^{(1)}$ in Eq. (3.17) as that would result in the interface curvature at the droplet surface. However, $\alpha^{(1)}$ varies sharply from 0 to 1 across an interface, and this can result in numerical errors. To address this, in this work, $f=\phi$ is used, which is defined as

\phi=\frac{\left(\alpha^{(2)}\right)^{m}}{\left(\alpha^{(2)}\right)^{m}+\left(1-\alpha^{(2)}\right)^{m}}\quad\textrm{for}\ m<1.

(3.22)

Shukla et al. [54] used this form of $\phi$ for their interface compression scheme, and it is used for the curvature computation in this work. The field $\phi$ is smoother as compared to $\alpha^{(2)}$ , and $m$ is a parameter that decides its smoothness. Note that $m=1$ would result in $\phi=\alpha^{(2)}$ . Finally, one would compute the curvature everywhere in the flow-field using Eq. (3.21). However, to reduce any numerical noise, it is only used in the regions where $0.99>\alpha^{(1)}>0.01$ , and it is set to zero otherwise.

To verify the curvature computation, a circular drop with diameter $d_{0}=9.6$ mm is placed at the center ( $x_{0}=0,y_{0}=0$ ) of a square domain of size $[6d_{0},6d_{0}]$ . The computational domain is uniformly divided into $240\times 240$ cells. The volume fraction for the droplet is initialized via

f=\sqrt{\left(\frac{x-x_{0}}{d_{0}/2}\right)^{2}+\left(\frac{y-y_{0}}{d_{0}/2}\right)^{2}},\quad\text{and}\quad\alpha^{(2)}=\frac{1}{2}\tanh\left(\frac{(f-1)d_{0}}{\delta 2\Delta x}\right)+\frac{1}{2}.

(3.23)

where $\delta$ is the thickness over which the numerical interface is diffused initially. The corresponding volume fraction fields for $\delta=0.5$ and $\delta=2.5$ , as well as the numerically computed curvature are shown in Fig. 3 along the center-line. The exact curvature for a droplet with diameter $d_{0}$ would be $2/d_{0}$ (shown by the solid horizontal line in Fig. 3). However, since the interface is diffused, if we were to compute the curvature using Eq. (3.17) and Eq. (3.23), then it would be $2/(d_{0}f)$ along the radius of the droplet (shown via dashed line in Fig. 3).

A comparison of the curvature computed using Eq. (3.21) and Eq. (3.22) at $m$ in Fig. 3 shows that for $\delta=2.5$ , when the interface is diffused across 8-10 cells, $m=0.1,0.5,1$ , provide a good match against $2/d_{0}f$ . The surface tension only acts at the sharp interface where the curvature is $2/d_{0}$ (horizontal line). Indeed, a numerically diffused volume fraction could result in further errors. Therefore, the curvature computation must work for sharp interfaces and should compute the curvature as $2/d_{0}$ . For $\delta=0.5$ , the interface is diffused over 2-3 cells. For this case, choosing $m=1$ (i.e., $f=\alpha^{(2)}$ , which was used in earlier works [51]) results in about $50$ % errors, but with $m=0.5$ , the errors are $<5$ %.

Choosing $m=0.1$ leads to slightly lower errors as compared to $m=0.5$ . However, the test considered here is an ideal scenario where the volume fraction is defined using Eq. (3.23), and as a result, $f$ is monotonic. In reality, there could be numerical oscillations in the $\alpha^{(2)}$ field, and those, when coupled with the high gradients of $\partial f/\partial\alpha$ near $\alpha=0,1$ , could result in a non-monotonic behavior of $f$ for a lower $m$ such as 0.1. This could lead to further errors in the computation, and considering this, $m=0.5$ is chosen as a balance between the two effects in this work.

3.5 Relaxation process

The quantities $R_{u}^{\kappa}$ , $R_{p}^{\kappa}$ and $R_{T}^{\kappa}$ in Eq. (2.12) are relaxation terms. Physically, they represent the interface jump conditions of Eq. (2.9) and Eq. (2.10) for a resolved interface. Since $\theta^{p},\theta^{u},\theta^{T}\rightarrow\infty$ for the seven-equation model, the relaxation processes are stiff, and they are applied separately from the inviscid and the viscous updates through a Strang splitting [28]. Correspondingly, the flow solution is updated as $Q^{n+1}=\mathcal{R}\mathcal{L}^{c}\mathcal{R}\mathcal{L}^{p}(Q^{n})$ , where $\mathcal{R}$ represents the stiff relaxation operator, $\mathcal{L}^{p}$ and $\mathcal{L}^{c}$ correspond to the predictor and the corrector steps of the time-update, and the superscript $n$ is used to represent the temporal iteration number. The pressure, velocity and the temperature relaxation operators are represented as $\mathcal{R}_{p}$ , $\mathcal{R}_{u}$ , $\mathcal{R}_{T}$ , respectively, and $\mathcal{R}=\mathcal{R}_{T}\mathcal{R}_{u}\mathcal{R}_{p}$ .

In terms of the numerical modeling of the interphase exchanges through the DEM, the terms with $\langle N_{int}\rangle$ in Eq. (3.9) correspond to the relaxation terms as shown by Abgrall and Saurel [48]. Here, $\langle N_{int}\rangle$ represents the number of unresolved interfaces within a finite-volume cell, and this would typically be finite for dispersed flows where the droplets are unresolved. Still, it is supposed to be negligible for the resolved interfaces such as those considered here. As a result of this, while using the DEM, there should be no need to use a separate stiff relaxation solver for resolved interfaces, and this is further confirmed in this work (see Section 4).

Interestingly, a recent work [68] showed that it is not necessary to model the non-conservative terms $I^{\kappa}$ in Eq. (2.12) if the relaxation terms $R^{\kappa}$ are modeled consistently. This suggests that both the non-conservative terms $I^{\kappa}$ and the relaxation terms $R^{\kappa}$ model the same interphase exchange effects but differently. While Schmidmayer et al. [68] were able to model the interphase exchanges only via $R^{\kappa}$ , our results here show that it is possible to model them only via the non-conservative terms $I^{\kappa}$ with the use of the DEM.

However, to provide a consistent comparison against the previous works that did use a relaxation solver [28], and because the seven-equation model would asymptote to the widely used five-equation model in the stiff pressure and velocity relaxation limit, the operators $\mathcal{R}_{p}$ and $\mathcal{R}_{u}$ are implemented in the current computational framework. The stiff temperature relaxation operator $\mathcal{R}_{T}$ is not included in this work, so the thermal exchanges through the interface are only via the resolved interphase exchange terms $I^{\kappa}$ in Eq. (3.9) via the DEM.

The stiff velocity relaxation operator $\mathcal{R}_{u}$ modifies the phase-averaged momentum and energy ( $\rho^{\kappa}\alpha^{\kappa}u_{i}^{\kappa}$ , $\rho^{\kappa}\alpha^{\kappa}E^{\kappa}$ ) while keeping the mixture momentum and energy conserved. The algorithm for this is described elsewhere [28] and it is not repeated here. The stiff pressure relaxation $\mathcal{R}_{p}$ results in local expansion/compression of the phases via a change in the volume fraction $\alpha^{\kappa}$ , and corresponding pressure work terms are included in the energy $\rho^{\kappa}\alpha^{\kappa}E^{\kappa}$ (see Eq. (3.9)). The algorithm matches that of [28], but with the following improvements:

The viscous and the surface tension forces are included in the pressure jump condition as shown in Eq. (2.10). As identified earlier [36], $\tau_{\textrm{i}\textrm{j}}^{(1)}n_{\textrm{i}}n_{\textrm{j}}=\tau_{\textrm{i}\textrm{j}}^{(2)}n_{\textrm{i}}n_{\textrm{j}}$ is maintained by the viscous DEM, and therefore, the pressure jump that would have to be applied at the interface via the relaxation process is $p^{(2)}-p^{(1)}=\sigma\beta$ . We can combine the volume fraction and the energy equations, resulting in

\frac{\partial\alpha^{(1)}\rho^{(1)}E^{(1)}}{\partial t}=-\left(p^{(1)}-\sigma\beta\right)\frac{\partial\alpha^{(1)}}{\partial t}\quad\text{and}\quad\frac{\partial\alpha^{(2)}\rho^{(2)}E^{(2)}}{\partial t}=p^{(1)}\frac{\partial\alpha^{(1)}}{\partial t}.

(3.24)

An iterative method similar to [28] is used to solve these equations and compute an updated volume fraction $\alpha^{(1)}$ and energies $E^{\kappa}$ to obtain $p^{(2)}-p^{(1)}=\sigma\beta$ . Unlike previous works without the surface tension, the energy exchange terms do not sum to zero, and the difference between the two, $\sigma\beta\theta^{p}\Delta p$ , is a source/sink due to the interfacial surface tension energy.

2.

Gradient descent is used as an iterative scheme. However, an adaptive stepping makes the algorithm robust for the various EOSs considered in this work and for a wide range of pressures. Specifically, it is used for the gradient descent that would reduce the step size if either the pressures computed via the EOS (either $\kappa=1,\,2$ ) become negative during the iterations or the updated volume fraction goes beyond $[\epsilon,1-\epsilon]$ with $\epsilon=10^{-10}$ .

3.6 Interface compression

The interface numerically diffuses across multiple cells as the simulation evolves. Several techniques have been proposed to alleviate this. Modifications to the MUSCL face-reconstruction [69] or the limiter [61] and applications of higher-order schemes [70, 47] have been useful for DIM. However, the DEM used in this work is a specialized technique, and we cannot directly apply it. Considering this, an interface compression scheme is borrowed here [54] which works independently of the other parts of the solver. Jain et al. [71] recently conducted a one-to-one comparison of various interface compression techniques and concluded that the choice would often depend on the modeling problem of interest. Their coupling with the seven-equation model can be considered in the future, but it is not addressed here. The employed iterative technique for interface compression requires solving the following PDE in a pseudo-time $\tau$ :

\frac{\partial\alpha^{(1)}}{\partial\tau}=n_{\textrm{i}}\frac{\partial}{\partial x_{\textrm{i}}}\left(\epsilon_{h}\ \vline\frac{\partial\alpha^{(1)}}{\partial x_{\textrm{i}}}\vline-\alpha^{(1)}(1-\alpha^{(1)})\right)=n_{\textrm{i}}\frac{\partial f_{comp}}{\partial x_{\textrm{i}}}.

(3.25)

The right-hand side of this equation contains compression and diffusion terms and the interface normal, which balance each other when a steady-state is reached in the pseudo-time. This would result in a consistent thickness across the diffused interface, which would be a function of $\epsilon_{h}$ . Shukla et al. [54] used this approach for a five-equation DIM, where they had to use an additional correction for the mixture density. This approach is not required for the seven-equation model as the densities for the phases are computed separately.

In curvilinear coordinates Eq. (3.25) is written as

\frac{\partial\alpha^{(1)}}{\partial\tau}=\left(n_{\textrm{i}}\frac{\partial\xi}{\partial x_{\textrm{i}}}\frac{\partial f_{comp}}{\partial\xi}+n_{\textrm{i}}\frac{\partial\eta}{\partial x_{\textrm{i}}}\frac{\partial f_{comp}}{\partial\eta}+n_{\textrm{i}}\frac{\partial\zeta}{\partial x_{\textrm{i}}}\frac{\partial f_{comp}}{\partial\zeta}\right),

(3.26)

where, $\alpha^{(1)}$ is defined at the cell-centers and updates to it in the pseudo-time are denoted as $d\alpha_{ijk}^{\tau}$ . First, $f_{comp}$ in Eq. (3.25) is computed at the cell-faces. To compute $f_{comp,i+\frac{1}{2}jk}$ at the cell-faces, a centered averaging is used for $\alpha^{(1)}$ following Eq. (3.13), and the derivatives $\partial\alpha^{(1)}/\partial x_{\textrm{i}}$ are computed using Eq. (3.11) and Eq. (3.12). Instead of computing the derivatives $\partial\alpha^{(1)}/\partial x_{\textrm{i}}$ directly with $\alpha^{(1)}$ , for numerical stability, they are computed as $\partial\phi/\partial x_{\textrm{i}}$ with the smooth field $\phi$ defined as previously in Eq. (3.22) [54].

As the next step, $\partial{f_{comp}}/\partial\xi$ , etc., are obtained at the cell-centers via a centered differencing such as

\displaystyle\frac{\partial f_{comp}}{\partial\xi}\ \vline_{\ ijk}=f_{comp,i+\frac{1}{2}jk}-f_{comp,i-\frac{1}{2}jk}.

(3.27)

Similar forms are used for the $j$ and the $k$ -directions. Last, for the updates $d\alpha_{ijk}^{\tau}$ in Eq. (3.26), the interface normal

\displaystyle n_{\textrm{i}}=\frac{\partial\phi/\partial x_{\textrm{i}}}{\lvert\partial\phi/\partial x_{\textrm{i}}\rvert}

(3.28)

is computed at the cell-centers via Eq. (3.11), and the corresponding cell-centered derivatives are computed as ${\partial\phi}/{\partial\zeta}\ \vline_{\ i+\frac{1}{2}jk}={1}/{2}\left(\phi_{i+1jk}-\phi_{i-1jk}\right)$ (shown here representatively in the $i$ -direction).

While applying the interface compression, the solution is evolved in the pseudo-time using explicit Euler until a steady-state solution is reached. The pseudo-time-step $d\tau$ is taken to be 0.2 for numerical stability [54]. The compression parameter $\epsilon_{h}$ is taken as $V_{ijk}^{1/3}/2$ , which would result in the interface diffusing across 2-3 cells.

Eq. (3.25) is not in a strictly conservative form. As a result, it was observed that when these pseudo-iterations are conducted for a very long time, the solution initially converges to a steady-state, but later it diverges. This is shown in Fig. 4 for example. For this test, a circular drop with diameter $d_{0}=9.6$ mm is placed at the center ( $x_{0}=0,y_{0}=0$ ) of a square domain of size $[2d_{0},2d_{0}]$ . The computational domain is uniformly divided into $30\times 30$ cells. The volume fraction for the droplet is initialized using Eq. (3.23) with $\delta=0.2$ and $\delta=2.5$ . The convergence is measured using an $L_{2}$ norm of $\alpha^{(1)}$ between the current and the previous iteration. For both the cases, initially, the iterations converge (in $\approx$ 10-30 steps), and a steady-state solution is reached. However, after running this for a longer time, i.e. $>10^{3}$ iterations, the solution diverges. The intermediate steady-state solution shows the interface diffused across 2–3 cells.

The solution to this problem in a standalone computation of the interface compression technique is straightforward. We could use an absolute or a relative tolerance limit to stop the iterations. However, while the compression scheme is being used as part of the DIM simulations, this is not possible. In addition to the interface compression, the flow fields are also being updated every time step, and even if the interface compression update was to be applied only once after every certain number (called $n^{comp}$ ) of simulation time-steps, the multiphase simulations can go on for an arbitrarily long number of iterations that could eventually cause this divergence. As a solution to this, after every $n^{comp}$ simulation time-steps, the interface compression scheme is iterated for at least 3 pseudo-time-steps, and the volume fraction field is updated only if the $L_{2}$ norms during this monotonically decrease. If not, the solution is reverted to its original state without any interface compression updates. The numerical results show that the interface compression did not cause any divergence by using this approach even when we ran the simulations for an arbitrarily long time (up to $10^{7}$ flow-field iterations). The interface compression was still active but acted only when necessary, and it maintained the interface thickness diffused only across 2-3 cells. $n^{comp}=100$ is used in this work. The results for only a circular drop are shown here, but the interface compression scheme has also been validated for more complex shapes, e.g., ellipse, two lobes, etc. Those results are not shown here for brevity.

3.7 Boundary Conditions

Seven different treatments are used in this work at the computational domain boundaries: slip walls, no-slip walls, subsonic inflow, subsonic outflow, supersonic inflow, supersonic outflow, periodic. For the no-slip walls, the wall-normal and the tangential velocities of both the phases are set to zero at the boundary face. Only the wall-normal component is set to zero for the slip walls, whereas a zero gradient Neumann condition is applied for the tangential component. All the other quantities, i.e., the temperature, the species mass-fractions, and the volume fraction, are set using a zero gradient Neumann condition at the walls.

The eigenvalues of the seven-equation system being solved here are $\{u^{(1)}-c^{(1)},u^{(1)},u^{(1)}+c^{(1)},u^{(2)}-c^{(2)},u^{(2)},u^{(2)}+c^{(2)},u^{I}\}$ [28]. Correspondingly, seven waves are propagating in this hyperbolic system of equations. The inflow/outflow treatment is dependent on the wave-propagation directions at the boundaries. For supersonic flows, i.e., $u^{(1)},u^{(2)}>c^{(1)}$ and $u^{(1)},u^{(2)}>c^{(2)}$ , all the waves propagate in the same direction. Therefore, a supersonic inflow, where all the waves enter the domain, is set using a Dirichlet condition. A supersonic outflow, where all the waves leave the domain, is set using a zero-gradient Neumann condition.

For subsonic inflow/outflows, where certain waves are entering the domain and certain that are leaving the domain, characteristic-based boundary conditions [72] are used. These were initially developed for multi-species single-phase flows, but the same form is used here for the two-phase flows. Since the wave-propagation speeds considered within this single-phase framework only correspond to the gas-phase wave speeds $\{u^{(1)}-c^{(1)},u^{(1)},u^{(1)}+c^{(1)}\}$ , it is implicitly assumed that the other phase properties, as well as the volume fraction, are propagating at speed $u^{(1)}$ as if they are passive scalars. This is not a valid approach in general. Still, it is retained here, considering that only a single phase (typically the gas) is present at the inflow/outflow boundaries for the configurations studied. The development of a characteristic-based boundary condition for fully compressible two-phase flows is the subject of future work.

4 Results

Various features of the seven-equation DIM framework are evaluated in this section. Table 2 lists the tests that are considered in this work and the features that they validate. For the seven-equation model, the phase-averaged flow fields for both the phases (i.e., $\rho^{\kappa}$ , $p^{\kappa}$ , etc.) are defined everywhere in the computational domain. For instance, even in a region where only the phase (1) is present, $\rho^{(2)}$ , $p^{(2)}$ etc. are still mathematically defined. The corresponding volume fraction of the phase (2) would be very small, i.e., $\alpha^{(2)}=\epsilon$ , and this would result $f^{mix}=\sum_{\kappa=1}^{2}\alpha^{\kappa}f^{\kappa}\approx f^{(1)}$ , so the phase (2) values, i.e., $\rho^{(2)}$ , $p^{(2)}$ , etc., aren’t expected to affect the mixture flow-fields with $\epsilon\rightarrow 0$ , but they still have to be defined and initialized. To avoid a division by zero and to make sure that the numerical solution does not get corrupted, the corresponding $\epsilon$ has to be larger than the machine precision (double precision used here, $10^{-15}$ ), and in this work it is chosen as $10^{-6}$ unless explicitly mentioned otherwise.

Table 2: List of verification and validation test cases for the seven-equation DIM.

Feature	Test-case and approach
Baseline solver	Periodic pulse convection - Order of scheme - Non-disturbing condition [28, 48]
$p^{(2)}/p^{(1)},\rho^{(2)}/\rho^{(1)}>>1$	Multiphase shock tube - Comparison against exact solution [48]
Arbitrary EOS & shocks	Shock transmission through - Air-Aluminum (CPG/SG) [73] - Air-HMX (CPG/MIEG) interface - Comparison against an exact solution
Surface tension	$\Delta p$ in a circular drop [51, 35] - Comparison against exact jump - Analysis of parasitic currents
Surface tension	Oscillating elliptical droplet - Compared against exact frequency [51, 35]
Viscous effects	Two-phase lid-driven cavity - Comparison against exact solution [74]
Viscous effects	Cylindrical and deforming droplet drag - Comparison against $C_{d}$ correlations [75]
Droplet deformation	Shock-droplet interaction - Comparison against experiments [76, 77, 31]
Reactive flows	Detonation-droplet interaction demonstration

4.1 Convection

A 1D periodic computational domain of length $L=1$ m is used for this test. The volume fraction $\alpha^{(1)}$ is initialized (at $t=t_{0}$ ) as a sine function varying between 0 to 1 as $\alpha^{(1)}(x)=(1/4)\sin\left(x/(2\pi L)\right)+1/2$ . It is capped between $\epsilon$ and $1-\epsilon$ to prevent any division by zero in the code. Volume fraction for the other phase is set as $\alpha^{(2)}=1-\alpha^{(1)}$ . Uniform pressure and velocity are specified as $p^{(1)}=p^{(2)}=10^{5}$ Pa and $u^{(1)}=u^{(2)}=100$ m/s. The densities are set as $\rho^{(1)}=1.0$ kg/m³ and $\rho^{(2)}=10^{3}$ kg/m³. The EOS for phase (1) is specified as CPG, with $\gamma^{(1)}=1.4$ , and for phase (2) it is specified as SG with $\gamma^{(2)}=4.4$ and $p_{0}^{(2)}=6\times 10^{8}$ . The surface tension, the viscous terms, and the interface compression are absent for this test, and a single specie is used for both the phases. Four different grids are tested with the number of cells varying as $N_{x}=20,\,50,\,100,\,200$ as well as the first- and the second-order accurate schemes. Simulations are conducted both with and without stiff relaxation.

The volume fraction pulse is allowed to convect in the periodic domain for one flow-through time ( $t=t_{f}$ ). The solution at the end is compared against the initial solution ( $t=t_{0}$ ), which serves as a ‘truth’ solution for this case. The ‘truth’ solutions are denoted with a superscript ‘ $truth$ ’. The $L_{2}$ error norms for any field-variable $f$ are computed as

||f||_{2}=\sqrt{\frac{1}{N_{x}}\sum_{i_{x}=1}^{N_{x}}\left({f(i_{x},t=t_{f})-f^{truth}(i_{x})}\right)}.

(4.1)

As noted in earlier studies [48], a “non-disturbing” condition for this case means that regardless of the volume fraction in the domain, since the simulation starts with a uniform pressure/velocity, one should maintain a uniform pressure/velocity at all times with machine accuracy. Both the pressure and the velocity maintain uniform values with $||p^{\kappa}||_{2}/p^{(1)}<10^{-15}$ and $||u^{\kappa}||_{2}/u^{(1)}<10^{-15}$ . The volume fraction $L_{2}$ error norms are plotted in Fig. 5. The volume fraction suffers numerical diffusion. However, the scheme order of accuracy is maintained as expected. The differences in the reference and the numerical order of accuracy slopes are less than 5% suggesting a reasonable match. The results shown here are without using the stiff relaxation solver. Since the uniform pressure/velocity condition is maintained between both the phases anyway, using the stiff relaxation solver did not make any difference. Those results are not shown here for brevity.

4.2 Multiphase shock tube

The objective of this test is to evaluate the method’s ability to deal with sharp jumps across a material interface [48]. A 1D computational domain of length $L=1$ m is used. A phase interface is initialized at $x/L=0.8$ . The flow variables at the initial time ( $t_{0}$ ) are set as $\alpha^{(1)}=\epsilon$ , $p^{(2)}=p^{(1)}=2\times 10^{8}$ for $x/L<0.8$ , and $\alpha^{(2)}=\epsilon$ , $p^{(2)}=p^{(1)}=10^{5}$ for $x/L\geq 0.8$ . The initial velocities are $u^{(1)}=u^{(2)}=0$ everywhere. The initial densities are set as $\rho^{(1)}=50$ and $\rho^{(2)}=1000$ everywhere. Based on this, the pressure and the velocity equilibrium are already maintained at $t_{0}$ . The phase (1) is modeled as CPG with $\gamma^{(1)}=1.4$ to mimic compressed air, whereas the phase (2) is modeled as SG with $\gamma^{(2)}=4.4$ , $p_{0}^{(2)}=6\times 10^{8}$ to mimic liquid water. Supersonic outflows are used at both ends of the tube to let the waves pass without reflection. Surface tension, viscous terms, and interface compression are absent for this test. Four different grids with the number of cells $N_{x}=200,\,500,\,1000,\,2000$ are tested with both the first- and the second-order-accurate schemes.

The simulation is conducted for 0.2 ms, and the final solution is compared against an analytical ‘truth’ solution computed via an exact Riemann solver as shown in Fig. 6. A shock propagates on the right side of the interface (in air), and a rarefaction travels towards the left (in water). The contact wave, and therefore, the interface, also moves slightly towards the right. Note that the flow variables for both phases ( $\kappa=1,\,2$ ) are evolved for the seven-equation model. However, for the resolved interface simulations that are considered here, only the mixture properties computed as $f^{mix}$ are of relevance, and only those are compared. The simulation results match the exact solution with increasing accuracy with an increase in the grid resolution and the order of the scheme (see Fig. 6). For instance, the normalized $L_{2}$ error norm of the pressure ( $||p^{mix}||_{2}/p^{(2)}_{x/L<0.8}$ ) reduces from $0.07$ to $0.03$ going from $N_{x}=200$ to $N_{x}=1000$ for the first order scheme, and almost a similar improvement in the accuracy is obtained for $N_{x}=200$ by going from the first to the second order scheme.

Results with and without the stiff relaxation (SR) solver match each other closely at all grid resolutions and scheme orders. For example, the normalized $L_{2}$ norm of the difference in $p^{mix}$ with and without the SR for $N_{x}=200$ and the second-order-accurate scheme is 0.002. This shows that the use of the DEM relieves the need to use a separate stiff relaxation solver. The DEM accurately models the interphase exchanges. The pressure and the velocity equality between the phases are obtained within the numerically diffused interface ([48] for further details) even without enforcing it. When a stiff pressure relaxation is applied the pressure and the velocity equality are enforced everywhere, and therefore, $p^{mix}=p^{1}=p^{2}$ , $u^{mix}=u^{1}=u^{2}$ . However, without it, $p^{1}$ , $p^{2}$ and $u^{1}$ , $u^{2}$ can be different away from the numerically diffused interface, but $f^{mix}\approx f^{\kappa}$ is still be maintained in the pure fluid regions as $\alpha^{\kappa}\approx 1$ .

4.3 Shock interaction with material interface

The passage of a shock through a material interface is simulated in this test. A 1D computational domain of length $L=1$ m is used. The region $x/L<0.5$ is filled with air ( $\alpha^{(2)}=\epsilon$ , CPG, $\gamma^{(1)}=1.4$ ), whereas, $x/L>0.5$ is Aluminum ( $\alpha^{(1)}=\epsilon$ , SG, $p_{0}^{(2)}=21.13\times 10^{9}$ , $\gamma^{(2)}=3.8$ ). A shock of strength Mach 2 is initialized in the air at $x/L=0.3$ to propagate towards the material interface. The pre-shock conditions are set to $p^{(1)}=p^{(2)}=10^{5}$ Pa, $u^{(1)}=u^{(2)}=0$ m/s, and $\rho^{(1)}=1.2$ kg/m³, $\rho^{(2)}=2784$ kg/m³. The post-shock conditions for $x/L<0.3$ are computed using the Rankine-Hugoniot relations for air as $p^{(1)}=4.56\times 10^{5}$ Pa, $u^{(1)}=429.0$ m/s, and $\rho^{(1)}=3.211$ kg/m³.

Only pure air is present for $x/L<0.3$ in reality. However, we also need to specify the flow fields for phase (2) for the seven-equation numerical modeling. The different initial conditions (IC) attempt to understand their effects on the final solution.

IC-1

: $\epsilon=10^{-3}$ , $p^{(2)}=p^{(1)}$ and $u^{(2)}=u^{(1)}$ for $x/L<0.3$
IC-2

: $\epsilon=10^{-6}$ , $p^{(2)}=p^{(1)}$ and $u^{(2)}=u^{(1)}$ for $x/L<0.3$
IC-3

: $\epsilon=10^{-9}$ , $p^{(2)}=p^{(1)}$ and $u^{(2)}=u^{(1)}$ for $x/L<0.3$
IC-4

: $\epsilon=10^{-3}$ , $p^{(2)}=10^{5}$ Pa and $u^{(2)}=0$ for $x/L<0.3$ .

For the first three IC, the pressure and the velocity equilibria are already attained between the phases at $t_{0}$ . However, this can be a potential problem since the pre/post-shock conditions that are specified correspond to a Mach 2 shock in the air, and they are not for the liquid, which still has a volume fraction of $\epsilon$ in this region. This would create a Riemann problem for the liquid resulting in artificial waves even before the shock hits the material interface. IC-4 aims at solving this problem by initializing phase (2) with uniform pressure/velocity even though phase (1) is initialized as a shock. This means that the SR cannot be enforced for this, but there would not be any artificial waves at the start.

The results at $t=0.3$ ms, after the shock has imparted the material interface, are analyzed and compared against an analytical solution [73] in Fig. 7. Depending on the impendence of the materials and the shock strength, the imparted shock would reflect/transmit as either a shock or a rarefaction. In this case, both the reflected and the transmitted waves are shocks. Due to the shock impact, the material interface also moves, although for this case, its velocity is only 0.095 m/s, and it is not noticeable from these plots.

In terms of the effect of the IC, due to the inconsistent pre/post-shock conditions as discussed before and a higher $\epsilon$ , as expected, IC-1 shows significant errors resulting from the artificial waves created in the phase (2) near the initial shock location ( $x/L=0.3$ ). The errors are already present at t=0.04 ms, well before the shock hits the interface, and therefore, the results at later times are also corrupted (not shown here). Both IC-2 and IC-3 show a significant improvement over IC-1. Since $\epsilon$ is smaller for them, even though the artificial waves are created in phase (2) initially, their effects on the mixture quantities, i.e., $p^{mix}$ , $u^{mix}$ seem to be minimal. IC-4 also seems to be a valid initialization approach for the seven-equation model and captures the shock transmission/reflection through the material interface, even with $\epsilon=10^{-3}$ as there are no artificial waves due to the initialization anymore. However, as noted before, this IC cannot be used with an SR solver since it is not enforced at the $t_{0}$ .

The results shown here are only for the second-order scheme. They improve accuracy with grid refinement while comparing against the reference exact solution (except for IC-1). Like the multiphase shock tube test earlier, the results with and without a stiff relaxation (SR) show a similar trend and reasonably match. Another test, where a shock initiated in the Aluminum enters the air, was also simulated. For this, the transmitted wave is a shock, whereas the reflected wave is a rarefaction due to the lower impendence of the air. These results are not shown here for brevity.

Next, to demonstrate the ability of this method to handle arbitrary EOS, the Al that was previously modeled using SG is now replaced with HMX, which is modeled using an MG EOS with the same parameters as used in an earlier work [62]. Unfortunately, an analytical solution to compare is not available here, but the obtained results are shown in Fig. 8. The results look qualitatively similar to the Al/air test, where the imparted shock resulted in a transmitted and a reflected shock at the interface. These results are obtained with the IC-2 approach, and they prove the robustness of this method to handle arbitrary EOS and large pressure and density ratios.

4.4 Surface tension effects

Several tests are conducted for validating the surface tension modeling following Nguyen and Dumbser [51]. A 2D square computational domain of width $L=1$ m is filled with air for the first test. A circular liquid droplet of radius $r_{0}=0.1549$ m is placed at its center. Slip walls are used as boundary conditions. As a result of the droplet curvature, higher pressure should be maintained within the droplet. The objective of the first test is to evaluate the ability of the DEM surface tension modifications to handle and maintain this pressure jump. Instead of computing the curvature from the volume fraction field, it is set to $\beta=1/r_{0}$ everywhere. The volume fractions inside and outside the droplet are initialized as $\alpha^{(1)}=\epsilon$ and $\alpha^{(2)}=\epsilon$ . The velocities are initialized as zero everywhere. The pressures are set as $p^{(1)}=10^{3}$ Pa and $p^{(2)}=p^{(1)}+\sigma\beta$ everywhere, and this would result in a jump in $p^{mix}$ at the droplet interface. In terms of the DEM flux computation, this is appropriate since this means that the pressure jump is only present at a phase (1)–(2) interface, where we included modifications for the surface tension. Whereas the pressure is uniform across the (1)–(1) and the (2)–(2) interfaces. The gas and the liquid phases are modeled using CPG and SG, respectively, with $\gamma^{(1)}=1.4$ , $\gamma^{(2)}=4.4$ , $p_{0}^{(2)}=6\times 10^{8}$ . The surface tension is set to $\sigma=342$ N/m. The gas-phase and the liquid-phase densities are $\rho^{(1)}=1.0$ kg/s and $\rho^{(2)}=100$ kg/s, respectively. This test is conducted both with and without the stiff relaxation solver. However, the pressure difference $p^{(2)}-p^{(1)}$ is already maintained to $p_{\sigma}=\sigma\beta$ from the start, which would mean that the relaxation solver does not have to act.

Three different grids, $200\times 200$ , $400\times 400$ , and $800\times 800$ and both first and the second-order schemes are tested. The simulations evolve until a physical time of 0.085 s, and at the end of it, the $L_{2}$ error norms for $p^{mix}$ , $u^{mix}$ and $\alpha^{(2)}$ are computed considering the initial conditions as the truth. For all the grids, for all the orders of the scheme, with or without the stiff relaxation solver, and with or without the interface compression, the normalized errors are maintained to machine accuracy ( $<10^{-15}$ ), suggesting that the DEM modifications for the surface tension can handle the surface jump condition exactly.

For the second test, the same setup is used. However, the curvature is computed from the volume fraction field as described in Section 3.4 instead of setting it as a constant. The pressure within the domain is initially set as $p^{(1)}=p^{(2)}=10^{3}$ Pa, and as a result of the surface tension force, it is expected to rise within the droplet (ideally to $\sigma\beta=2207$ Pa). The simulations are run for 0.034 s ( $\sim$ 40,000 steps for the 800 $\times$ 800 grid), and the pressure jump within the droplet at the end of it is measured by integrating over the entire computational domain $\Omega$ as

\Delta p=\int_{\Omega}p^{(2)}\alpha^{(2)}\bigg{/}\int_{\Omega}\alpha^{(2)}-\int_{\Omega}p^{(1)}\alpha^{(1)}\bigg{/}\int_{\Omega}\alpha^{(1)}.

(4.2)

This is compared against the exact pressure difference based on Young-Laplace relation as $\Delta p^{err}=|\Delta p-\sigma\beta|$ . Even though the background flow should have been at rest, parasitic currents are observed due to the numerical errors in the curvature computation. They are consistent with the previous works [51]. These are quantified using the kinetic energy as

KE^{mix}=\int_{\Omega}\left(\frac{1}{2}\alpha^{(1)}\rho^{(1)}{u^{(1)}}^{2}+\frac{1}{2}\alpha^{(2)}\rho^{(2)}{u^{(2)}}^{2}\right)\bigg{/}\int_{\Omega}\left(\alpha^{(1)}+\alpha^{(2)}\right).

(4.3)

Values of both the $\Delta p^{err}$ and the $KE^{mix}$ are reported in Table 3 for various grids and with/without the stiff relaxation (SR) and with/without the interface compression (comp.). Only the second-order scheme results are shown here. However, the results with the first-order scheme are also consistent with the conclusions here. The lowest errors are observed for the case without the SR but with the interface compression. The errors in the computed pressure jump are less than 2% for all the grids for this case, and the parasitic velocities are $<0.2$ m/s. The parasitic currents and the $p^{mix}$ within the droplet for this case are shown in Fig. 9 for reference.

Regardless of the interface compression, the errors with the SR are at least an order of magnitude higher than those without the SR. Further insights into this are provided below. The errors do not necessarily reduce with the grid, and this could be due to the contrasting effects, as we would compute the curvature more accurately for the finer grids. However, the coarser grids would help the parasitic currents to diffuse numerically. The errors with the interface compression are also significantly lower than those without it since the droplet interface becomes irregular due to the parasitic currents. Still, the interface compression helps keep it regular, resulting in a more accurate curvature computation.

Table 3: Errors

\Delta p^{err}

[Pa] and

KE^{mix}

[kg m²/s²] for the circular droplet surface tension test are shown for various grids and various simulation options.

Case	200 $\times$ 200		400 $\times$ 400		800 $\times$ 800
	$\Delta p^{err}$	$KE^{mix}$	$\Delta p^{err}$	$KE^{mix}$	$\Delta p^{err}$	$KE^{mix}$
w/o SR, w/o comp.	$1.06\times 10^{2}$	$7.99\times 10^{-6}$	$2.64\times 10^{2}$	$5.92\times 10^{-5}$	$2.09\times 10^{2}$	$4.35\times 10^{-5}$
w SR, w/o comp.	$8.84\times 10^{2}$	$1.58\times 10^{-4}$	$1.18\times 10^{3}$	$1.56\times 10^{-3}$	$1.25\times 10^{3}$	$2.18\times 10^{-3}$
w/o SR, w comp.	$3.52\times 10^{1}$	$2.26\times 10^{-6}$	$1.49\times 10^{1}$	$2.06\times 10^{-6}$	$5.86\times 10^{1}$	$2.10\times 10^{-6}$
w SR, w comp.	$1.44\times 10^{2}$	$1.68\times 10^{-4}$	$1.01\times 10^{3}$	$4.59\times 10^{-4}$	$4.66\times 10^{2}$	$1.49\times 10^{-3}$

To understand the reasons for the large errors with the SR, $p^{\kappa}$ , $\rho^{\kappa}$ , and $\alpha^{(2)}$ are plotted along the centerline for the 800 $\times$ 800 test with and without the SR in Fig. 10. As noted before, the DEM modifications in the surface tension are designed so that a pressure jump is handled across the (1)-(2) interface, whereas the (1)–(1) and the (2)–(2) interfaces should behave as before. As expected, the case with the SR maintains pressure equality away from the interface and enforces a pressure jump between $p^{(1)}$ and $p^{(2)}$ within the numerically diffused interface. However, because of this, pressure jumps are present in both $p^{(1)}$ and $p^{(2)}$ near the boundaries of the diffused interface. They result in pressure jumps at the (1)–(1) and the (2)–(2) interfaces within the DEM, possibly causing these larger errors.

On the other hand, for the case without the SR, $p^{(1)}$ stays constant, and only the $p^{(2)}$ changes within the numerically diffused interface. The liquid pressure $p^{(2)}$ reduces almost to the same value as $p^{(1)}$ right before the interface ( $x/d_{0}>0.5$ ) and then increases back due to the surface tension terms at the interface. The initial reduction in it is still in the gas phase (i.e., $\alpha^{(2)}$ is small), and therefore, it does not adversely affect the overall mixture ( $mix$ ) flow fields. The surface tension is accurately captured without the SR within the numerically diffused interface due to the DEM modifications in the (1)–(2) Riemann problem. The errors mentioned above due to the SR are also apparent in the density fields, as $\rho^{(1)}$ in the presence of the SR increases by almost $100$ % within the droplet. In contrast, the oscillations in both $\rho^{(1)}$ and $\rho^{(2)}$ are $<0.5$ % for the case without the SR.

For the third test, instead of starting from a circular droplet, an elliptical drop is used with the following shape [35, 51]

\frac{(x-x_{0})^{2}}{0.2^{2}}+\frac{(y-y_{0})^{2}}{0.12^{2}}=1.

(4.4)

The ellipsoidal droplet starts converting into a circular shape because of the initial higher/lower curvature on the horizontal/vertical ends of the droplet. The potential energy stored in the droplet’s surface tension converts into kinetic energy. The droplet oscillates from an initial ellipse into an intermediate circular shape and back to an ellipse. The kinetic energy oscillates during this process, and we can compute its frequency against an analytical value based on the Rayleigh formula [78]

\omega=(o^{3}-o)\frac{\sigma}{(\rho^{(1)}+\rho^{(2)})r_{0}^{3}}.

(4.5)

where $o$ is the oscillation mode. Corresponding to $o=2$ and $r_{0}=0.1549$ m, an oscillation frequency is determined as $t_{osc}=2\pi/\omega=0.085$ s.

Based on the conclusions from the previous test, this test is simulated without the SR. The evolution of the droplet shape at various times is shown Fig. 11, and corresponding $KE^{mix}$ is shown in Fig. 12. Similar to the previous test, interface compression is necessary for accurate curvature and surface tension computation. As expected, the initial elliptical droplet converts to a circular droplet by $t/t_{osc}=0.3$ . This is accompanied by an increase in the $KE^{mix}$ . The droplet regains an elliptical shape by $t/t_{osc}=0.5$ and $KE^{mix}$ shows reduces to zero. Ideally, this process should continue in the absence of any viscous effects. However, because of the numerical diffusion here, the peaks of $KE^{mix}$ reduce with time. The 400 $\times$ 400 grid only captures two peaks of $KE^{mix}$ , whereas the finer grids can capture more. The computed oscillation frequency matches the exact value with $<$ 15% accuracy.

4.5 Viscous effects

To verify the implementation of the viscous terms, a lid-driven cavity with a Reynolds number (Re) of 100 is simulated first. Initially, the fluid is at rest in a square domain of size $L$ . The top wall is accelerated to a velocity $u^{\infty}$ at $t=0$ , and all the other walls stay at rest. The flow within the domain starts moving and rotating due to the no-slip walls and the viscous effects. The simulations are run until a steady-state solution is achieved. Single-phase simulations [79] for this case have been established before, and the results match with an exact solution [74]. Two-phase results are obtained here with the DIM, where both the phases are set as CPG and with the same properties. The volume fraction $\alpha^{(2)}$ is initialized corresponding to a circular droplet within the domain using Eq. (3.23). Both the single-phase and the two-phase simulation results with and without the SR match the exact solution as shown in Fig. 13, verifying the implementation of the viscous terms. These results are with the second-order scheme.

Next, to validate the viscous effects in two-phase flows, the drag of a droplet in a viscous medium is computed. A 2D domain of length $L_{x}=15r_{0}$ and width $L_{y}=10r_{0}$ is used, where $r_{0}=0.5$ mm is the radius of an initially circular droplet initially placed at $(x,y)=(4r_{0},5r_{0})$ . The gas and the liquid phases are modeled using CPG and SG, respectively, with $\gamma^{(1)}=1.4$ , $\gamma^{(2)}=4.4$ , $p_{0}^{(2)}=6\times 10^{8}$ . The velocity field surrounding the droplet is initiated using a potential flow solution, using a stream function

\psi=U_{0}r\rm{sin}(\theta)[1-R_{0}^{2}/r^{2}],\quad u^{(1)}_{1}=u^{(2)}_{1}=-\partial\psi/\partial y,\quad u^{(1)}_{2}=u^{(2)}_{2}=\partial\psi/\partial x.

(4.6)

where $U_{0}$ is the free-stream velocity, and $\theta=1^{0}$ is used to trigger vortex shedding. The velocities inside the droplet are initialized as zero. The dynamic viscosity for the gas phase is set as $\mu^{(1)}=1.8\times 10^{-5}$ kg/(m-s), and for the liquid it is $\mu^{(1)}=1.14\times 10^{-3}$ kg/(m-s). The surface tension is specified as $\sigma=0.0728$ N/m. The Reynolds number ( $Re$ ) is defined as $Re=2\rho^{(1)}r_{0}U_{0}/\mu_{g}$ , and the Weber number ( $We$ ) is defined as $We=2\rho^{(1)}r_{0}U_{0}^{2}/\sigma$ .

Two conditions are simulated following a previous work using the five-equation model [80]. One with $U_{0}=5$ m/s, that corresponds to $Re=81.3$ and $We=0.1$ , and another with $U_{0}=50$ m/s, that corresponds to $Re=813$ and $We=10$ . A time-scale for normalization is defined as $t_{0}=r_{0}/U_{0}$ . Three grids are used with the number of cells across the droplet diameter $N_{x}=20,50,100$ . The simulations are conducted with and without the SR and with and without the interface compression. The interface diffuses irregularly without interface compression, similar to the previous surface tension test. The boundary and shear layers that are supposed to grow surrounding the sharp gas-liquid interface are wrongly predicted. Therefore, only the results with the interface compression are discussed further in this section.

The flow fields for the two simulated conditions are shown in Fig. 14. The droplet at $We=0.1$ maintains a circular shape, but it deforms for $We=10$ . The low $Re$ case develops a separation bubble behind the droplet, and the flow field remains relatively steady. On the other hand, the high $Re$ case demonstrates a vortex shedding process. The positive and negative velocities surrounding the droplet are due to the viscous effects. For further processing, the droplet velocity and a drag coefficient $C_{D}$ are computed as

U_{c}={\int_{\Omega}\alpha^{(2)}u^{(2)}_{1}}\bigg{/}{\int_{\Omega}\alpha^{(2)}},\quad C_{D}=\frac{\rho^{(2)}}{\rho^{(1)}}\frac{dU_{c}}{dt}\frac{\pi r_{0}}{(U_{0}-U_{c})^{2}}.

(4.7)

The computed $C_{D}$ is compared against a steady-state drag correlation for a rigid cylinder $C_{D,S}(Re)$ [75] for $We=0.1$ in Fig. 15. The droplet with $We=10$ starts deforming, and therefore, a comparison against the correlation is not possible for that. These results show that establishing the viscous flow surrounding the droplet takes approximately $5t_{0}$ , and $C_{D}$ reaches a steady value after that. The $C_{D}$ is overpredicted for the coarser grids. However, a better match against the correlation $C_{D,S}(Re)$ is achieved with grid-refinement. The $C_{D}$ for the finest grid is also $\approx$ 10% higher than the correlation $C_{D,S}(Re)$ , but this could be attributed to the relative acceleration effect [81] that is not included in the correlation $C_{D,S}(Re)$ , but are present in the current simulations since the relative velocity $U_{r}=U_{c}-U_{0}$ decreases with time. These results also show that even though the $C_{D}$ with and without the SR do not match exactly, they are very similar, and this further indicates the efficacy of the viscous DEM in accurately modeling the non-conservative interface terms such that there is no need to use the SR with the seven-equation model.

4.6 Shock-droplet interaction

This test simulates the early-stage deformation of a cylindrical water droplet when hit by a shock wave. The setup follows the experiments by Igra et al. [76] and other numerical works that used a five-equation DIM [77, 31]. A circular droplet of diameter ( $d_{0}$ ) 4.8 mm is placed in a 2D computational domain of size 23 $d_{0}$ $\times$ 12 $d_{0}$ . The background gas-phase is initialized with a shock of strength $M=1.47$ that is located at a distance $d_{0}$ ahead of the droplet center The pre-shock (background) pressure and velocity are $p^{(1)}=p^{(2)}=101325$ Pa and $u^{(1)}=u^{(2)}=0$ m/s, respectively. The corresponding post-shock conditions are subsonic, and they are set as $p^{(1)}=p^{(2)}=238558$ Pa and $u^{(1)}=u^{(2)}=u_{0}=226$ m/s. The water density is set as $\rho^{(2)}=1000$ m/s. The pre- and the post-shock gas-phase densities are $\rho^{(1)}=1.204$ kg/m³ and $\rho^{(1)}=2.18$ kg/m³, respectively. Note that these initial conditions are equivalent to the IC-2 as described earlier in Section 4.3. The liquid water is modeled with a SG EOS with $\gamma^{(2)}=6.12$ and $p_{0}^{(2)}=343.44\times 10^{6}$ Pa. The gas phase is modeled with a CPG EOS with $\gamma^{(1)}=1.4$ . The volume fraction is initialized using Eq. (3.23) with $\delta=0.5$ .

The computational domain is discretized using a uniform grid $N_{x}=17,\,42,\,83,\,170,\,255,\,340$ cells across the droplet diameter. The top and the bottom boundaries are set as slip walls. The left boundary behind the shock is set as a non-reflective subsonic inflow at the post-shock velocity $u_{0}$ . The right boundary, ahead of the shock, is set as a supersonic outflow to let the incoming shock wave pass through without any reflections. The Reynolds number ( $Re$ ) and Weber number ( $We$ ) corresponding to this condition are $Re=4860$ and $We=72320$ , respectively. Under these conditions, the droplet deformation is primarily driven by inertia during the early stages. Therefore, the surface tension and viscous effects can be neglected, as confirmed in the earlier numerical studies [77, 31]. Note that these effects would still be relevant near the final stages of the breakup and at much smaller length scales that are not resolved here. The simulations are conducted and compared with and without the surface tension and with and without the interface compression.

A non-dimensional time is defined as $\tau=t/\left(u_{0}/d_{0}\sqrt{\rho^{(1)}/\rho^{(2)}}\right)$ , where the denominator characterizes time for breakup via Rayleigh–Taylor or Kelvin–Helmholtz instabilities [82]. The simulation results at various times are shown in Fig. 16 for the case with $N_{x}=83$ . These are with the interface compression and with/without the stiff relaxation. Pressure contours and Numerical Schliren ( $|\partial\rho^{(1)}/\partial x_{\textrm{i}}|$ ) are used to identify the shock structure. The droplet is using a $\alpha^{(2)}=0.5$ contour. As the initial shock hits the droplet, there are transmitted and reflected waves. Unlike the previous 1D case, the shock structures are no longer planar because of the circular droplet. The waves that propagate inside the cylinder transmit and reflect through the material interface. The droplet remains circular during the earlier stages (up to $\tau\approx 0.17$ ). However, it starts to deform due to the vortical structures created around it later. These features are qualitatively similar to the experimental observations [76] and previous numerical results [77, 31]. The results with and without the SR match closely, reconfirming the ability of the DEM to alleviate the need to use an SR.

The experiments measured the droplet motion and the droplet deformation in four geometrical parameters. 1) Leading edge drift ( $x_{L}$ ), that is, the x-location of the droplet leading edge measured with respect to its initial location, 2) droplet centerline width ( $w$ ), which measures the width of the cylinder along the $x$ -axis, 3) droplet spanwise diameter ( $d$ ), that is the width of the droplet in the $y$ -direction, and 4) droplet area ( $A$ ). All these are computed from the simulation flow-fields using a threshold volume fraction $\alpha^{(2)}_{T}$ , considering that $\alpha^{(2)}>\alpha^{(2)}_{T}$ is the droplet shape whose geometrical parameters have to be determined.

Because of the uncertainty in the experimental measurements, it is unclear what value $\alpha^{(2)}_{T}$ should take to best match the data. Even though the gas-liquid interface is supposed to remain sharp in reality, this region may appear diffused depending on the measurement technique because of the small droplets strip surrounding the parent droplet surface. In addition to this, the interface also numerically diffuses for the simulations. The results are analyzed for $\alpha^{(2)}_{T}=0.1,\,0.5,\,0.9$ , but only those with $\alpha^{(2)}_{T}=0.9$ are shown here. Sensitivity of five-equation DIM results in the choice of $\alpha^{(2)}_{T}$ can be found elsewhere [31] and our observations are similar for the cases without the interface compression. The simulation results with the interface compression do not show a sensitivity to $\alpha^{(2)}_{T}$ since a regular interface thickness is maintained. However, those results are considered erroneous for this test, and further details are discussed later.

The simulations with the SR and the interface compression take longer to run, e.g., the simulations are 1.6 $\times$ slower with the interface compression and 5.6-times slower with the SR. Because of the cost of these simulations, the finer grids are not simulated until completion with the SR and the interface compression. As discussed next, the results are still considered sufficient to draw valuable conclusions.

The droplet leading edge ( $x_{L}/d_{0}$ ) evolution is shown in Fig. 17. All grid resolutions show a similar trend, but a grid convergence is achieved for $N_{x}\geq 170$ . Stiff relaxation or interface compression does not result in a noticeable difference in leading-edge predictions. The results match the experimental data for $\tau<0.7$ , however, the difference eventually grows to almost 100% by $\tau=1.3$ . This difference could be due to the subsonic inflow boundary conditions used on the left boundary for the simulations to maintain the post-shock velocity. The experiments were conducted within a longer shock tube where We did not need this.

The centerline droplet width ( $w$ ) decreases with time due to the droplet compression, and the spanwise droplet diameter ( $d$ ) increases as the droplet gets elongated. These are plotted in Fig. 18 and Fig. 19, respectively. Without the SR and the interface compression, the predicted droplet width $w$ matches well against the data ( $<$ 5% error) for all grid resolutions. On the other hand, the droplet height $d$ is grossly underpredicted at coarser resolutions, and a grid convergence is obtained only for $N_{x}\geq 255$ . The fine grid results capture the same trend as the experimental measurements [76]. However, they show about 30-40% overprediction. Similar errors have also been observed in earlier numerical studies using the five-equation DIM, but the reasons are unclear. The droplet area is expected to reduce with time for the experiments due to the stripping of small-sized droplets. The simulations capture this without the SR and the interface compression as shown in Fig. 20. The coarsest grid shows up to 25% underprediction. However, a grid convergence is obtained for $N_{x}\geq 170$ , and the finer grids’ results match the measurements with ¡5% error.

Regarding the stiff relaxation solver usage, similar to the leading-edge droplet predictions, the predictions for all three geometrical features ( $d$ , $w$ , $A$ ) are similar with/without the SR. The most significant differences are observed later in $d$ and $A$ predictions, which are still $<$ 10-15%. This further confirms the efficacy of the DEM in alleviating the need to use a separate SR solver, which, as demonstrated for this case, could be as much as 5.6-times as costly.

The effect of the interface compression on the prediction is widely apparent on $w$ , $d$ , and $A$ . In reality, the droplet area would reduce with time due to the stripping of small-sized droplets. One could estimate the size of these stripped droplets as $\sim 6\ \mu$ m based on $We\sim 1$ . It is not possible to resolve these small droplets with our computational grid, and in the DIM simulations considered here, they would be represented as intermediate volume fractions ( $\epsilon<\alpha^{(2)}<1-\epsilon$ ). The interface compression approach that is used here can only work with the resolved fluids, and it wrongly interprets the regions with $\epsilon<\alpha^{(2)}<1-\epsilon$ as numerically diffused regions, applying the compression there. As a result of this, the results with the interface compression are insensitive to the threshold volume fraction $\alpha^{(2)}_{T}$ , but they do not match the data. Transition to a dispersed phase Eulerian–Eulerian and Eulerian–we can consider lagrangian modeling approaches such as [8] to model these regions with $\epsilon<\alpha^{(2)}<1-\epsilon$ , however, that is not the focus of this work, and it would be considered in the future.

4.7 Detonation-droplet interaction

One could encounter a detonation-droplet interaction in practical applications such as a liquid-fueled or liquid-assisted rotating detonation engine (RDE) [83, 84] or an energetic material packed with metal particles [85]. Some recent efforts have focused on studying this phenomenon experimentally [53]. However, the simulations here serve as a qualitative demonstration without detailed measurements for quantitative validation.

A 2D computational domain of size $(1500,100)$ mm is used. A liquid droplet ( $\kappa=2$ ) of diameter $d_{0}$ mm is initialized at $(x,y)=(750,50)$ mm. A uniform grid with $N_{x}=100$ cells across the droplet diameter is used. The left and the right boundaries are set as supersonic outflows to let the internal waves pass through without reflection. The top and the bottom boundaries are set as slip walls. The entire domain is at rest initially ( $u^{(1)}=u^{(2)}=0$ ) and the pressure are set as $p^{(1)}=p^{(2)}=10^{5}$ Pa. The densities are $\rho^{(1)}=1.204$ kg/m³, and $\rho^{(2)}=1000$ kg/m³. The gas-phase ( $\kappa=1$ ) contains gaseous $\rm{H_{2}}$ and air at stoichiometry to mimic conditions in a typical rotating detonation engine [63].

In order to initiate a 2D detonation wave, $x<10$ mm is initially filled with detonation products at $T^{(1)}=3500$ K and $p^{(1)}=5\times 10^{6}$ Pa. A circular disturbance of radius 1 mm is asymmetrically added in the domain to trigger the transverse detonation waves [86].The droplet volume fraction $\alpha^{(2)}$ is initialized using Eq. (3.23) with $\delta=0.5$ . The droplet diameter $d_{0}=40$ mm is relatively large compared to real applications. However, it is chosen here considering the cost and the grid requirements for resolving a droplet.

The liquid phase is modeled using an SG EOS with $\gamma^{(2)}=4.1$ and $p^{0}=6\times 10^{8}$ , and the gas phase is modeled using a TPG EOS that uses NASA polynomial curve-fits for the species. The droplet would eventually vaporize because of the high temperature behind the detonation. However, for the current conditions, the ratio of the droplet vaporization time ( $t_{drop}$ ) to the detonation passage time ( $t_{deton}$ ) is estimated to be substantial, $t_{drop}/t_{deton}\sim 10^{6}$ , suggesting that one can simulate the early stages of this detonation-droplet interaction without using any heat/mass-transfer models. Similar to the previous test in Section 4.6, the inertia primarily drives the droplet deformation at these high-speed conditions. Therefore, the viscous and the surface tension effects are neglected for this test.

The simulation results at various time instances are shown in Fig. 21. The results shown here are with the interface compression scheme and the SR. However, simulations without them have also been conducted. As a result of the initiation, a 2D detonation wave propagates in the gas phase before reaching the droplet. The structure of this detonation wave is confirmed by a pressure spike (of 2.55 MPa), an increase in the temperature (to 3078 K) due to coupled shock, and burning at the front. The detonation speed before reaching the droplet is measured as 1968 m/s, which is extremely close to the Chapman–-Jouguet velocity of 1964 m/s computed using a Chemical Equilibrium Applications (CEA) code at the same conditions.

Similar to the shock-droplet interaction, as the detonation wave hits the droplet, it results in the transmission and reflection of non-planar waves. A circular reflected shock is observed propagating back into the detonation products. As the wave propagates after crossing the droplet, it eventually regains its reactive detonation structure. The droplet remains circular at initial times. However, it deforms and elongates, as was also the case with the shock-droplet interaction.

5 Conclusions

A seven-equation diffused interface method (DIM) framework is developed in this work. Unlike the commonly used reduced DIM models, the seven-equation model does not enforce a velocity/pressure equilibrium. Alleviating a strict pressure relaxation allows this method for arbitrarily high pressure/density ratio and arbitrary EOS of the phases, and alleviating the strict velocity equilibrium would enable it to be used in a dispersed phase regime where the MPE are unresolved. The seven-equation model has been used for simplistic 1D tests [48, 28]. However, it is extended to be used in 2D/3D with arbitrary EOS, surface tension, viscous effects, and multi-species reacting flows. Various tests, i.e., multiphase shock tube, shock-material interaction, oscillating droplet, drag in a viscous medium, shock-droplet interaction, and detonation-droplet interaction, are conducted to validate different parts of the method. Following are some novel contributions and conclusions of this work

1.

A discrete equations method (DEM) [48] is used in this work to compute the non-conservative interphase exchange terms accurately. It is demonstrated via 1D and 2D tests that this relieves the need to use a stiff relaxation solver, which is otherwise required to maintain pressure–velocity equilibrium between the phases [28] and could make the simulations as much as 5 $\times$ slower. For the cases where the initial conditions are phase-specific, e.g., a shock-material interaction problem with the shock initialized in one of the phases, the use of the pressure/velocities equilibrium at the start results in errors due to mismatching pre-/post-shock conditions between the phases.
2.

A stiff relaxation solver is still consistent with the previous works that used a seven-equation model, but certain modifications have been made. Adaptive-step gradient descent is used for pressure relaxation to make the scheme robust and useful for arbitrary EOS. The developed seven-equation framework is robust for arbitrarily high density/pressure ratios ( $\rho^{(2)}/\rho^{(1)}\sim 1-10^{3}$ , $\rho^{(2)}/\rho^{(1)}\sim 1-10^{5}$ ), and arbitrary EOS (SG, CPG, TPG, MG).
3.

To include surface tension, a modification is proposed in the inviscid fluxes of DEM and the corresponding HLLC Riemann solver. This proposed method can handle a surface-tension-related pressure jump across a material interface in an exact numerical manner. The pressure jump within a circular droplet is captured with $<2$ % error with numerically computed curvature. An oscillating droplet where the surface potential energy converts back and forth into the kinetic energy is simulated, and the oscillation period matches the exact value with a $<15$ % error. The subsequent peaks of kinetic energy reduce due to numerical diffusion. However, the results improve with grid refinement.
4.

Modifications for modeling 2D and 3D viscous effects within the DEM framework have been included. Droplet evolution in a viscous medium is simulated. The computed drag coefficient $C_{D}$ matches a steady-state drag correlation for a low Weber number case where the droplet does not deform. Vortex shedding is observed for a higher $Re$ case, but the droplet also deforms due to the higher $We$ .
5.

An interface compression scheme initially used for the five-equation model is modified for the seven-equation model. The interface compression was necessary for the surface tension tests and the viscous drag computation. Otherwise, the droplet interface gets irregular, leading to difficulties in curvature computation and boundary and shear layer predictions surrounding the droplet surface.
6.

A shock-droplet interaction case is simulated at experimental conditions, and the computed droplet deformation parameters are shown to match against available data. This would create many small droplets due to the atomization process at such high $Re$ and $We$ in reality. However, since the grid does not resolve these, they manifest as regions of intermediate volume fractions ( $\epsilon>\alpha^{(2)}>1-\epsilon$ ) with DIM. As a result, the interface compression scheme was detrimental for this case since it would treat these regions as numerically diffused regions and apply compression. We can consider modeling these regions using a dispersed-phase approach in the future.
7.

Finally, the scheme has been demonstrated to work for multi-species, reacting flows via a detonation-droplet interaction problem. This test is a qualitative demonstration without any validation/verification data.

This work does not consider heat and mass transfer at the droplet interface, which is an obvious extension of this work. Infinite [42, 43] or finite-rate [87] thermochemical relaxation solvers can be used. A key benefit of the seven-equation model is its possible extension to dispersed-phase regions where an equilibrium between the velocities is not maintained. It is shown in this work that with the use of DEM, there is no need to enforce a stiff relaxation for resolved multiphase simulations, and the next step could be to use drag and heat and mass transfer laws for interphase exchange in the regions where there are unresolved droplets. Such coupling of the seven-equation DIM with the dense-to-dilute dispersed model [8] will be considered in the future.

Acknowledgment

This work was supported by NASA Glenn Research Center [grant number NNX15AU91A].

References

Lefebvre [1998] A. H. Lefebvre, Gas turbine combustion, CRC press, 1998.
Sutton [2006] G. P. Sutton, History of liquid propellant rocket engines, AIAA, 2006.
Tarver et al. [1996] C. M. Tarver, S. K. Chidester, A. L. Nichols, Critical conditions for impact-and shock-induced hot spots in solid explosives, The Journal of Physical Chemistry 100 (1996) 5794–5799.
Vehring [2008] R. Vehring, Pharmaceutical particle engineering via spray drying, Pharmaceutical research 25 (2008) 999–1022.
Gorokhovski and Herrmann [2008] M. Gorokhovski, M. Herrmann, Modeling primary atomization, Annu. Rev. Fluid Mech. 40 (2008) 343–366.
Saurel and Pantano [2018] R. Saurel, C. Pantano, Diffuse-interface capturing methods for compressible two-phase flows, Annual Review of Fluid Mechanics 50 (2018) 105–130.
Balachandar and Eaton [2010] S. Balachandar, J. K. Eaton, Turbulent dispersed multiphase flow, Annual Review of Fluid Mechanics 42 (2010) 111–133.
Panchal and Menon [2021] A. Panchal, S. Menon, A hybrid Eulerian-Eulerian/Eulerian-Lagrangian method for dense-to-dilute dispersed phase flows, Journal of Computational Physics 439 (2021) 110339.
Bryngelson et al. [2019] S. H. Bryngelson, K. Schmidmayer, T. Colonius, A quantitative comparison of phase-averaged models for bubbly, cavitating flows, International Journal of Multiphase Flow 115 (2019) 137–143.
Bryngelson et al. [2021] S. H. Bryngelson, R. O. Fox, T. Colonius, Conditional moment methods for polydisperse cavitating flows, arXiv preprint arXiv: 2112.14172 (2021).
Sethian and Smereka [2003] J. A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annual Review of Fluid Mechanics 35 (2003) 341–372.
Sussman et al. [2007] M. Sussman, K. M. Smith, M. Y. Hussaini, M. Ohta, R. Zhi-Wei, A sharp interface method for incompressible two-phase flows, Journal of computational physics 221 (2007) 469–505.
Hirt and Nichols [1981] C. W. Hirt, B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1981) 201–225.
Sussman and Puckett [2000] M. Sussman, E. G. Puckett, A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows, Journal of Computational Physics 162 (2000) 301–337.
Brackbill et al. [1992] J. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of Computational Physics 100 (1992) 335–354.
Fedkiw et al. [1999] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method), Journal of Computational Physics 152 (1999) 457–492.
Elghobashi [2019] S. Elghobashi, Direct numerical simulation of turbulent flows laden with droplets or bubbles, Annual Review of Fluid Mechanics 51 (2019) 217–244.
Das and Udaykumar [2020] P. Das, H. Udaykumar, A sharp-interface method for the simulation of shock-induced vaporization of droplets, Journal of Computational Physics 405 (2020) 109005.
Terashima and Tryggvason [2009] H. Terashima, G. Tryggvason, A front-tracking/ghost-fluid method for fluid interfaces in compressible flows, Journal of Computational Physics 228 (2009) 4012–4037.
Barton et al. [2011] P. T. Barton, B. Obadia, D. Drikakis, A conservative level-set based method for compressible solid/fluid problems on fixed grids, Journal of Computational Physics 230 (2011) 7867–7890.
Baer and Nunziato [1986] M. Baer, J. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials, International journal of multiphase flow 12 (1986) 861–889.
Liu et al. [1994] X.-D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, Journal of computational physics 115 (1994) 200–212.
Van Leer [1979] B. Van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, Journal of Computational Physics 32 (1979) 101–136.
Petitpas et al. [2009] F. Petitpas, J. Massoni, R. Saurel, E. Lapebie, L. Munier, Diffuse interface model for high speed cavitating underwater systems, International Journal of Multiphase Flow 35 (2009) 747–759.
Schmidmayer et al. [2020] K. Schmidmayer, F. Petitpas, S. Le Martelot, É. Daniel, ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows, Computer Physics Communications 251 (2020) 107093.
Saurel et al. [2016] R. Saurel, P. Boivin, O. Le Métayer, A general formulation for cavitating, boiling and evaporating flows, Computers & Fluids 128 (2016) 53–64.
Rodio and Abgrall [2015] M. G. Rodio, R. Abgrall, An innovative phase transition modeling for reproducing cavitation through a five-equation model and theoretical generalization to six and seven-equation models, International Journal of Heat and Mass Transfer 89 (2015) 1386–1401.
Saurel and Abgrall [1999] R. Saurel, R. Abgrall, A multiphase godunov method for compressible multifluid and multiphase flows, Journal of Computational Physics 150 (1999) 425–467.
Kapila et al. [2001] A. Kapila, R. Menikoff, J. Bdzil, S. Son, D. S. Stewart, Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations, Physics of Fluids 13 (2001) 3002–3024.
Allaire et al. [2002] G. Allaire, S. Clerc, S. Kokh, A five-equation model for the simulation of interfaces between compressible fluids, Journal of Computational Physics 181 (2002) 577–616.
Meng and Colonius [2015] J. Meng, T. Colonius, Numerical simulations of the early stages of high-speed droplet breakup, Shock Waves 25 (2015) 399–414.
Liu et al. [2019] N. Liu, Z. Wang, M. Sun, R. Deiterding, H. Wang, Simulation of liquid jet primary breakup in a supersonic crossflow under adaptive mesh refinement framework, Aerospace Science and Technology 91 (2019) 456–473.
Murrone and Guillard [2005] A. Murrone, H. Guillard, A five equation reduced model for compressible two phase flow problems, Journal of Computational Physics 202 (2005) 664–698.
Petitpas et al. [2007] F. Petitpas, E. Franquet, R. Saurel, O. Le Metayer, A relaxation-projection method for compressible flows. Part II: Artificial heat exchanges for multiphase shocks, Journal of Computational Physics 225 (2007) 2214–2248.
Perigaud and Saurel [2005] G. Perigaud, R. Saurel, A compressible flow model with capillary effects, Journal of Computational Physics 209 (2005) 139–178.
Abgrall and Rodio [2014] R. Abgrall, M. G. Rodio, Discrete equation method (DEM) for the simulation of viscous, compressible, two-phase flows, Computers & Fluids 91 (2014) 164–181.
Saurel et al. [2008] R. Saurel, F. Petipas, R. Abgrall, Modelling phase transition in metastable liquids: application to cavitating and flashing flows, Journal of Fluid Mechanics 607 (2008) 313–350.
Bryngelson et al. [2021] S. H. Bryngelson, K. Schmidmayer, V. Coralic, J. C. Meng, K. Maeda, T. Colonius, MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver, Computer Physics Communications 266 (2021) 107396.
Schmidmayer et al. [2020] K. Schmidmayer, S. H. Bryngelson, T. Colonius, An assessment of multicomponent flow models and interface capturing schemes for spherical bubble dynamics, Journal of Computational Physics 402 (2020) 109080.
Wood [1930] A. B. Wood, A Textbook of Sound, G. Bell and Sons Ltd., 1930.
Saurel et al. [2009] R. Saurel, F. Petitpas, R. A. Berry, Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures, Journal of Computational Physics 228 (2009) 1678–1712.
Zein et al. [2010] A. Zein, M. Hantke, G. Warnecke, Modeling phase transition for compressible two-phase flows applied to metastable liquids, Journal of Computational Physics 229 (2010) 2964–2998.
Pelanti and Shyue [2019] M. Pelanti, K.-M. Shyue, A numerical model for multiphase liquid–vapor–gas flows with interfaces and cavitation, International journal of multiphase flow 113 (2019) 208–230.
Demou et al. [2022] A. D. Demou, N. Scapin, M. Pelanti, L. Brandt, A pressure-based diffuse interface method for low-mach multiphase flows with mass transfer, Journal of Computational Physics 448 (2022) 110730.
Rodriguez et al. [2021] M. Rodriguez, S. H. Bryngelson, S. Cao, T. Colonius, Acoustically-induced bubble growth and phase change dynamics near compliant surfaces, in: 11th International Symposium on Cavitation, 2021.
Dorschner et al. [2020] B. Dorschner, L. Biasiori-Poulanges, K. Schmidmayer, H. El-Rabii, T. Colonius, On the formation and recurrent shedding of ligaments in droplet aerobreakup, Journal of Fluid Mechanics 904 (2020).
Coralic and Colonius [2014] V. Coralic, T. Colonius, Finite-volume weno scheme for viscous compressible multicomponent flows, Journal of Computational Physics 274 (2014) 95–121.
Abgrall and Saurel [2003] R. Abgrall, R. Saurel, Discrete equations for physical and numerical compressible multiphase mixtures, Journal of Computational Physics 186 (2003) 361–396.
Chang and Liou [2007] C.-H. Chang, M.-S. Liou, A robust and accurate approach to computing compressible multiphase flow: Stratified flow model and AUSM+-up scheme, Journal of Computational Physics 225 (2007) 840–873.
Tokareva and Toro [2010] S. Tokareva, E. F. Toro, HLLC-type riemann solver for the Baer-Nunziato equations of compressible two-phase flow, Journal of Computational Physics 229 (2010) 3573–3604.
Nguyen and Dumbser [2015] N. T. Nguyen, M. Dumbser, A path-conservative finite volume scheme for compressible multi-phase flows with surface tension, Applied Mathematics and Computation 271 (2015) 959–978.
Panchal [2022] A. Panchal, Modeling moderately dense to dilute multiphase flows, Ph.D. thesis, Georgia Institute of Technology, 2022.
Dyson et al. [2022] D. Dyson, S. Vasu, A. Arakelyan, N. Berube, S. Briggs, J. Ramirez, E. M. Ninnemann, K. Thurmond, G. Kim, W. H. Green, et al., Detonation wave-induced breakup and combustion of RP-2 fuel droplets, in: AIAA SciTech 2022 Forum, 2022, pp. AIAA–2022–1453.
Shukla et al. [2010] R. K. Shukla, C. Pantano, J. B. Freund, An interface capturing method for the simulation of multi-phase compressible flows, Journal of Computational Physics 229 (2010) 7411–7439.
Dodd and Ferrante [2016] M. S. Dodd, A. Ferrante, On the interaction of Taylor length scale size droplets and isotropic turbulence, Journal of Fluid Mechanics 806 (2016) 356–412.
Fechter et al. [2017] S. Fechter, C.-D. Munz, C. Rohde, C. Zeiler, A sharp interface method for compressible liquid–vapor flow with phase transition and surface tension, Journal of Computational Physics 336 (2017) 347–374.
Schmidmayer et al. [2017] K. Schmidmayer, F. Petitpas, E. Daniel, N. Favrie, S. Gavrilyuk, A model and numerical method for compressible flows with capillary effects, Journal of Computational Physics 334 (2017) 468–496.
Pelanti and Shyue [2014] M. Pelanti, K.-M. Shyue, A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves, Journal of Computational Physics 259 (2014) 331–357.
Tanguy et al. [2007] S. Tanguy, T. Ménard, A. Berlemont, A level set method for vaporizing two-phase flows, Journal of Computational Physics 221 (2007) 837–853.
Jain et al. [2020] S. S. Jain, A. Mani, P. Moin, A conservative diffuse-interface method for compressible two-phase flows, Journal of Computational Physics 418 (2020) 109606.
Chiapolino et al. [2017] A. Chiapolino, R. Saurel, B. Nkonga, Sharpening diffuse interfaces with compressible fluids on unstructured meshes, Journal of Computational Physics 340 (2017) 389–417.
Akiki et al. [2017] M. Akiki, T. P. Gallagher, S. Menon, Mechanistic approach for simulating hot-spot formations and detonation in polymer-bonded explosives, AIAA Journal 55 (2017) 585–598.
Salvadori et al. [2022] M. Salvadori, P. Tudisco, D. Ranjan, S. Menon, Numerical investigation of mass flow rate effects on multiplicity of detonation waves within a H2/Air rotating detonation combustor, International Journal of Hydrogen Energy 47 (2022) 4155–4170.
Baurle et al. [1994] R. Baurle, G. Alexopoulos, H. Hassan, Assumed joint probability density function approach for supersonic turbulent combustion, Journal of Propulsion and Power 10 (1994) 473–484.
Patel and Menon [2008] N. Patel, S. Menon, Simulation of spray–turbulence–flame interactions in a lean direct injection combustor, Combustion and Flame 153 (2008) 228–257.
Balakrishnan et al. [2010] K. Balakrishnan, D. Nance, S. Menon, Simulation of impulse effects from explosive charges containing metal particles, Shock Waves 20 (2010) 217–239.
Toro [2009] E. F. Toro, The hll and hllc riemann solvers, in: Riemann solvers and numerical methods for fluid dynamics, Springer, 2009, pp. 315–344.
Schmidmayer et al. [2021] K. Schmidmayer, J. Cazé, F. Petitpas, E. Daniel, N. Favrie, Modelling interactions between waves and diffused interfaces, 2021. URL: https://hal.archives-ouvertes.fr/hal-03387818, working paper or preprint.
Shyue and Xiao [2014] K.-M. Shyue, F. Xiao, An eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic thinc approach, Journal of Computational Physics 268 (2014) 326–354.
Gryngarten and Menon [2013] L. D. Gryngarten, S. Menon, A generalized approach for sub-and super-critical flows using the local discontinuous galerkin method, Computer Methods in Applied Mechanics and Engineering 253 (2013) 169–185.
Jain et al. [2021] S. S. Jain, M. C. Adler, J. R. West, A. Mani, P. Moin, S. K. Lele, Assessment of diffuse-interface methods for compressible multiphase fluid flows and elastic-plastic deformation in solids, arXiv preprint arXiv:2109.09729 (2021).
Poinsot and Lele [1992] T. J. Poinsot, S. K. Lele, Boundary conditions for direct simulations of compressible viscous flows, Journal of Computational Physics 101 (1992) 104–129.
Sridharan et al. [2015] P. Sridharan, T. L. Jackson, J. Zhang, S. Balachandar, Shock interaction with one-dimensional array of particles in air, Journal of Applied Physics 117 (2015) 075902.
Ghia et al. [1982] U. Ghia, K. Ghia, C. Shin, High-Re solutions for incompressible flow using the navier-stokes equations and a multigrid method, Journal of Computational Physics 48 (1982) 387–411.
White and Corfield [2006] F. M. White, I. Corfield, Viscous fluid flow, volume 3, McGraw-Hill New York, 2006.
Igra et al. [2002] D. Igra, T. Ogawa, K. Takayama, A parametric study of water column deformation resulting from shock wave loading, Atomization and Sprays 12 (2002).
Chen [2008] H. Chen, Two-dimensional simulation of stripping breakup of a water droplet, AIAA Journal 46 (2008) 1135–1143.
Rayleigh et al. [1879] L. Rayleigh, et al., On the capillary phenomena of jets, Proc. R. Soc. London 29 (1879) 71–97.
Gallagher et al. [2017] T. Gallagher, M. Akiki, S. Menon, V. Sankaran, V. Sankaran, Development of the generalized maccormack scheme and its extension to low mach number flows, International Journal Numerical Methods in Fluids 85 (2017) 165–188.
Patel [2007] N. V. Patel, Simulation of hydrodynamic fragmentation from a fundamental and an engineering perspective, Ph.D. thesis, Georgia Institute of Technology, 2007.
Temkin and Kim [1980] S. Temkin, S. S. Kim, Droplet motion induced by weak shock waves, Journal of Fluid Mechanics 96 (1980) 133–157.
Pilch and Erdman [1987] M. Pilch, C. Erdman, Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop, International Journal of Multiphase Flow 13 (1987) 741–757.
Salvadori et al. [2022] M. Salvadori, A. Panchal, D. Ranjan, S. Menon, Numerical study of detonation propagation in H2-Air with kerosene droplets, in: AIAA SciTech 2022 Forum, 2022, pp. AIAA–2022–0394.
Kindracki [2015] J. Kindracki, Experimental research on rotating detonation in liquid fuel–gaseous air mixtures, Aerospace Science and Technology 43 (2015) 445–453.
Gogulya et al. [2004] M. Gogulya, M. Makhov, A. Y. Dolgoborodov, M. Brazhnikov, V. Arkhipov, V. Shchetinin, Mechanical sensitivity and detonation parameters of aluminized explosives, Combustion, Explosion, and Shock Waves 40 (2004) 445–457.
Kailasanath et al. [1985] K. Kailasanath, E. Oran, J. Boris, T. Young, Determination of detonation cell size and the role of transverse waves in two-dimensional detonations, Combustion and Flame 61 (1985) 199–209.
Pelanti [2021] M. Pelanti, Arbitrary-rate relaxation techniques for the numerical modeling of compressible two-phase flows with heat and mass transfer, arXiv preprint arXiv:2108.00556 (2021).