2.1. Solution difficulties

The underlying difficulties in developing an accurate set of governing equations for coupled gas and granular flow are based on the difficulties of representing important physical interactions, some of which are not defined as clearly as they are for single-phase gaseous flows. This is the topic addressed in the next section (§ 3) on model formulation. Once we have such a model, there remains the issue of how we can use and adapt existing numerical simulation tools to solve specific problems of interest, and what new issues arise from the formulation. In the case of coupled gas–granular flows, the underlying new difficulties are essentially physical, but these also have counterparts in how the solutions must be implemented numerically.

The difficulties in actually solving the equations shown below in § 3 result from a combination of non-conservative interphase coupling terms and numerical stiffness in dilute and dense granular regions. More specifically, we need to deal with: (1) non-conservative ‘nozzling’ terms (Schwendeman, Wahle & Kapila Reference Schwendeman, Wahle and Kapila2006; Karni & Hernández-Dueñas Reference Karni and Hernández-Dueñas2010; Crochet & Gonthier Reference Crochet and Gonthier2013) arising from variations of the particle volume fraction within a control volume; (2) a $pDV$ work term that affects the gas as particles enter and leave a control volume (Abgrall & Karni Reference Abgrall and Karni2010); (3) hyperbolic degeneracy when the intergranular stress in the granular phase is zero; (4) stiffness of the intergranular stress in highly packed regions; and (5) potential loss of hyperbolicity from complex-valued wave speeds.

When the pressure is uniform and the background fluid is stationary (zero velocity), the gas-phase momentum equations with a variation of particle volume fraction reduce to the ‘non-disturbing condition’ (Saurel & Abgrall Reference Saurel and Abgrall1999),

(2.3) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}=p_{g}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}.\end{eqnarray}$$

Here ${\it\alpha}_{g}$ is the gas-phase volume fraction, $p_{g}$ is the gas-phase pressure, and the term on the right-hand side is the ‘nozzling term’ due to its analogy with the $pDA$ term in the quasi-one-dimensional gas dynamic nozzle flow equations. Computing the left- and right-hand sides of equation (2.3) independently (Toro Reference Toro, A. and B.1989; Nussbaum et al. Reference Nussbaum, Helluy, Hérard and Carriére2006) can lead to situations where, for purely numerical reasons, the terms may not balance. This induces unphysical flow and pressure oscillations that can cause numerical solutions to fail (Liou et al. Reference Liou, Chang, Nguyen and Theofanous2008).

Another complication is that the compaction wave speed of the particles (analogous to the sound speed of a gas) ranges between zero, when ${\it\Theta}_{s}$ is zero or at low particle volume fractions, and approaches infinity at the packing limit. If the compaction wave speed is zero, the result is a pressureless granular gas (sometimes called ‘sticky particle’ models (LeVeque Reference LeVeque2004)), which causes numerical problems (specifically, hyperbolic degeneracy (Pelanti & LeVeque Reference Pelanti and LeVeque2006)). This means that characteristic-based methods, such as many weighted essentially non-oscillatory (WENO) schemes (Balsara & Shu Reference Balsara and Shu2000), or the standard piecewise parabolic method (PPM) (Colella & Woodward Reference Colella and Woodward1984) cannot be used to solve these equations in a straightforward manner. Godunov-based methods would require a Riemann solver for the granular phase that is valid in the pressureless state and transitions to a high-pressure state. Despite hyperbolic degeneracy, approximate Riemann solvers were previously designed for pressureless dust (Collins et al. Reference Collins, Ferguson, Chien, Kuhl, Krispin and Glaz1994; Pelanti & LeVeque Reference Pelanti and LeVeque2006), but not for kinetic-theory-based granular models.

The intergranular stress is extremely sensitive to even very small fluctuations of solid volume fraction in dense regions that approach the packing limit. This makes it difficult to develop solution approaches that are both low-dissipation (high-order) and robust. A difficulty arises because the volume fraction of particles influences the gas and granular phases differently. Particles restrict the area where fluid can flow, and so they accelerate the gas by a nozzling term. On the other hand, the volume fraction controls the compressibility of the granular phase in the same way that density does for the gas phase.

To deal with these numerical problems, the coupled system of equations may be solved using highly dissipative numerical methods. In these cases, refining the grid or using any high-order (therefore more accurate) numerical methods can trigger instabilities that can be masked by excessive numerical diffusion (Nussbaum et al. Reference Nussbaum, Helluy, Hérard and Carriére2006). Recent methods appear to have overcome this for the BN equations (Karni & Hernández-Dueñas Reference Karni and Hernández-Dueñas2010; Crochet & Gonthier Reference Crochet and Gonthier2013). Another solution method, used in interior ballistics, is to split the gas and granular conservation equations and solve a Riemann problem for each phase independently and then to compute the non-conservative nozzling and $pDV$ coupling terms separately (Toro Reference Toro, A. and B.1989), but this approach can violate the non-disturbing condition.

Now high-order, low-dissipation methods for solving the compressible Euler and Navier–Stokes equations are becoming widespread due, in part, to their excellent properties for describing shock waves and compressible turbulence (Grinstein, Margolin & Rider Reference Grinstein, Margolin and Rider2007; Thornber et al. Reference Thornber, Mosedale, Drikakis, Youngs and Williams2008; Poludnenko & Oran Reference Poludnenko and Oran2010, Reference Poludnenko and Oran2011). The nonlinear filtering effect and added dissipation near sharp gradients that is inherent in shock capturing schemes acts, in many ways, like the explicit filtering and subgrid turbulence models used in traditional large-eddy simulation. In addition, the turbulence in many high-speed shock-laden flows is non-Kolmogorov (Oran & Gamezo Reference Oran and Gamezo2007), that is, inhomogeneous, non-isotropic and driven by shock–shock interactions and Richtmyer–Meshkov and Rayleigh–Taylor instabilities. Physically accurate subgrid turbulence models usually do not exist for such conditions. Often the only choice in these situations is to use implicit large-eddy simulation, which requires low-dissipation numerical algorithms (Grinstein et al. Reference Grinstein, Margolin and Rider2007; Thornber et al. Reference Thornber, Mosedale, Drikakis, Youngs and Williams2008).

This presents an interesting quandary for multiphase models based on kinetic theory or any other multiphase model where the particles are supposed to be subgrid. In such models, the computational cell size is limited by the particle diameter. It is necessary to have enough particles in a computational cell so that the granular phase can be represented as some kind of fluid. Thus, there are grid limitations on the size of the finest mesh. This limiting computational cell size may be too coarse to capture features of the gas phase, such as the viscous sublayer of a turbulent boundary layer.

In this paper, we first present a model that solves two sets of Euler equations, one for the gas phase and one for the particles, with coupling terms that connect them. We then show how the two sets of equations can be solved using a high-order, low-dissipation numerical method. The technique uses components from existing Riemann solvers and edge interpolation schemes, and incorporates a relatively small modification to the AUSM $^{+}$ -up flux (Liou Reference Liou1996, Reference Liou2006) for the granular phase. The fluxes and computed states obtained from solving the gas- and granular-phase Riemann problems are assembled to form the convective fluxes and non-conservative terms. Tests of the method that show its ability to preserve the non-disturbing condition, compute multiphase shock-tube problems and simulate strong shocks interacting with dense layers of dust are discussed. Finally, we present comparisons from this model to experimental measurements of shock and particle-curtain interactions (Rogue et al. Reference Rogue, Rodriguez, Haas and Saurel1998; Ling et al. Reference Ling, Wagner, Beresh, Kearney and Balachandar2012; Wagner et al. Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012) and transmitted oblique granular shocks and granular contact surfaces in layers of dust (Fan et al. Reference Fan, Chen, Jiang and Li2007). In future papers, we will discuss extensions that include effects of granular viscosity and heat and mass diffusion, which are not included in this paper.

3.1. Euler equations for a single-phase, multispecies flow

The Euler equations for a multicomponent, reacting, compressible gas flow with $N_{g}$ chemical species can be written as

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{v})=0, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}Y_{j}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}Y_{j}\boldsymbol{v})=\dot{{\bf\omega}}_{j},\quad j=1,\ldots ,N_{g}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}\boldsymbol{v}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{v}\boldsymbol{v})+\boldsymbol{{\rm\nabla}}p={\it\rho}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}E}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[\boldsymbol{v}({\it\rho}E+p)]={\it\rho}\,\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{v}, & \displaystyle\end{eqnarray}$$

$Y_{j}$

${\it\rho}$

$p$

$T$

$E$

$\boldsymbol{v}$

$\boldsymbol{g}$

$j$

$\dot{{\bf\omega}}_{j}$

The ideal gas equation of state is used to relate the pressure, chemical composition, temperature and density of the gas phase,

(3.5) $$\begin{eqnarray}p={\it\rho}R_{u}T\mathop{\sum }_{j=1}^{N_{g}}\frac{Y_{j}}{M_{j}},\end{eqnarray}$$

where $M_{j}$ is the molecular weight of species $j$ and $R_{u}$ is the universal gas constant. The total energy of the gas phase, $E$ , is given by

(3.6) $$\begin{eqnarray}E=H-\frac{p}{{\it\rho}}=\mathop{\sum }_{j=1}^{N}Y_{j}\left(h_{f\,j}^{0}+\int _{T_{0}}^{T}C_{Pj}(s)\,\text{d}s\right)-\frac{p}{{\it\rho}}+\frac{\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}}{2},\end{eqnarray}$$

where $H$ is the total enthalpy, $h_{f\,j}^{0}$ is the enthalpy of formation at reference temperature $T_{0}$ and $C_{Pj}$ is the constant-pressure specific heat of species $j$ . Thermodynamic data may be taken from a number of sources. Here data are taken from Goos, Burcat & Rusnic (Reference Goos, Burcat and Rusnic2010). The sound speed, $c$ , is given by

(3.7) $$\begin{eqnarray}c^{2}={\it\gamma}\frac{p}{{\it\rho}},\end{eqnarray}$$

where ${\it\gamma}$ is the ratio of specific heats.

3.2. Euler equations for a gas in a multiphase flow

Now we modify (3.1)–(3.4) so they can be used to describe the gas phase in the gas–granular multiphase problem by accounting for the volume taken up by the particles through the gas-phase volume fraction ${\it\alpha}_{g}$ (Ishii & Hibiki Reference Ishii and Hibiki2006),

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g})=S^{M}, & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}Y_{g,j}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}Y_{g,j}\boldsymbol{v}_{g})={\it\alpha}_{g}\dot{{\bf\omega}}_{g,j}+S_{g,j}^{Y},\quad j=1,\ldots ,N_{g}, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}\boldsymbol{v}_{g})+\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}=p_{int}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}+\boldsymbol{S}^{\boldsymbol{p}}, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}E_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[{\it\alpha}_{g}\boldsymbol{v}_{g}({\it\rho}_{g}E_{g}+p_{g})]=-p_{int}\frac{\partial {\it\alpha}_{g}}{\partial t}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{v}_{g}+S_{g}^{E}, & \displaystyle\end{eqnarray}$$

$g$

${\it\alpha}_{g}$

(3.12a,b ) $$\begin{eqnarray}{\it\alpha}_{g}+{\it\alpha}_{s}=1\quad \text{and}\quad \boldsymbol{{\rm\nabla}}{\it\alpha}_{g}=-\boldsymbol{{\rm\nabla}}{\it\alpha}_{s},\end{eqnarray}$$

where ${\it\alpha}_{s}$ is the volume fraction of particles. The quantities $Y_{g,j}$ , ${\it\rho}_{g}$ , $p_{g}$ , $T_{g}$ , $E_{g}$ and $\boldsymbol{v}_{g}$ are the mass fraction of species $j$ , density, pressure, temperature, total energy and velocity vector of the gas, respectively. The pressure at the gas–solid interface around the particles, $p_{int}$ , is used in the nozzling ( $p_{int}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}$ ) and $pDV$ work ( $-p_{int}(\partial {\it\alpha}_{g}/\partial t)$ ) terms. The reaction rate due to homogeneous chemical reactions is denoted by $\dot{{\bf\omega}}_{g,j}$ . Additional interphase coupling terms describing transfer of mass, momentum and energy between the gas and the particles are denoted by $S^{M}$ , $S_{g,j}^{Y}$ , $\boldsymbol{S}^{\boldsymbol{p}}$ and $S_{g}^{E}$ . These are described in § 3.5.

3.3. Euler equations for the particulate phase in a multiphase flow

These equations are similar in form to the Euler equations for the gas phase (Gidaspow Reference Gidaspow1994). The main difference is that compressibility is introduced by changes in solid volume fraction rather than the material density of the particles, so that compaction waves (analogous to acoustic waves, rarefactions and shocks in a molecular gas) travelling through granular media directly change the particle volume fraction, ${\it\alpha}_{s}$ . The granular-phase governing equations for mass, momentum, pseudo-thermal energy (PTE) and internal energy with $N_{s}$ species within each particle are

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s})=-S^{M}, & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}Y_{s,j}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}Y_{s,j}\boldsymbol{v}_{s})={\it\alpha}_{s}\dot{{\bf\omega}}_{s,j}+S_{s,j}^{Y},\quad j=1,\ldots ,N_{s}, & \displaystyle\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s}\boldsymbol{v}_{s})+\boldsymbol{{\rm\nabla}}p_{s}+\boldsymbol{{\rm\nabla}}p_{fric}=-\boldsymbol{{\rm\nabla}}{\it\alpha}_{s}p_{g}-p_{int}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}+{\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{g}-\boldsymbol{S}^{\boldsymbol{p}}, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}E_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}E_{s}\boldsymbol{v}_{s})=-p_{s}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}_{s}+S_{s}^{PTE}, & \displaystyle\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}e_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}e_{s}\boldsymbol{v}_{s})=S_{s}^{E}, & \displaystyle\end{eqnarray}$$

$p_{s}$

$p_{fric}$

$e_{s}$

$E_{s}$

$\dot{{\it\gamma}}$

$E_{s}$

$i$

$\dot{{\it\omega}}_{s,i}$

$S_{s,j}^{Y}$

$S_{s}^{PTE}$

$S_{s}^{E}$

The PTE, $E_{s}$ , represents the energy due to random translational motion of the particles,

(3.18) $$\begin{eqnarray}E_{s}={\textstyle \frac{3}{2}}{\it\Theta}_{s}.\end{eqnarray}$$

It is a function of the granular temperature, ${\it\Theta}_{s}$ , defined as the mean square of the particle velocity fluctuations. The granular temperature is usually presented in units of $\text{m}^{2}~\text{s}^{-2}$ rather than K (Gidaspow Reference Gidaspow1994; Van Wachem et al. Reference Van Wachem, Schouten, van den Bleek, Krishna and Sinclair2001; Brilliantov & Pöschel Reference Brilliantov and Pöschel2004).

The ‘normal’ temperature of the particles, $T_{s}$ , is determined from the solid sensible internal energy, $e_{s}$ , through

(3.19) $$\begin{eqnarray}e_{s}=\mathop{\sum }_{j=1}^{N_{s}}Y_{s,j}\left(e_{f\,j}^{0}+\int _{T_{0}}^{T_{s}}C_{V,s,j}(s)\,\text{d}s\right),\end{eqnarray}$$

where $N_{s}$ is the number of species in the solid phase, $e_{f\,j}^{0}$ is an internal energy of formation and $C_{V,s,j}$ is the constant-volume specific heat of species $j$ in the particle phase. This equation assumes that the temperature distribution within each particle is uniform, and so the Biot modulus is small. More complex equations that include the effect of non-uniform temperature distributions inside the particles could be used if the Biot number were large (Koo & Kuo Reference Koo and Kuo1977; Markatos & Kirkcaldy Reference Markatos and Kirkcaldy1983; Nussbaum et al. Reference Nussbaum, Helluy, Hérard and Carriére2006). (Note that the granular temperature, ${\it\Theta}_{s}$ , is not related to ‘normal’ temperature $T_{s}$ .)

The density of the solid particles is given by

(3.20) $$\begin{eqnarray}{\it\rho}_{s}=\mathop{\sum }_{j=1}^{N_{s}}Y_{s,j}{\it\rho}_{s,j}(T_{s}),\end{eqnarray}$$

where ${\it\rho}_{s,j}$ is the density of species $j$ . In this paper, we take ${\it\rho}_{s}$ to be constant. Relaxing this assumption will be a topic of future work.

The solids pressure, $p_{s}$ , represents collisional effects of the particles and is given by an equation of state for a granular gas (Gidaspow Reference Gidaspow1994),

(3.21) $$\begin{eqnarray}p_{s}={\it\alpha}_{s}{\it\rho}_{s}{\it\Theta}_{s}[1+2(1+e){\it\alpha}_{s}g_{0}],\end{eqnarray}$$

where $e$ is the coefficient of restitution. The radial distribution function, $g_{0}$ , is defined by (Gidaspow Reference Gidaspow1994)

(3.22) $$\begin{eqnarray}\frac{1}{g_{0}}=1-\left(\frac{{\it\alpha}_{s}}{{\it\alpha}_{s,max}}\right)^{1/3},\end{eqnarray}$$

where ${\it\alpha}_{s,max}$ is the packing limit, which is an input parameter commonly set to 0.65 (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Igci et al. Reference Igci, Andrews, Sundaresan, Pannala and O’Brien2008), which is representative of the random close-packing limit of uniform spheres (Andreotti et al. Reference Andreotti, Forterre and Pouliquen2013). Other expressions for $g_{0}$ and $p_{s}$ could be used (see e.g. Kamenetsky et al. Reference Kamenetsky, Goldshtein, Shapiro and Degani2000; Van Wachem et al. Reference Van Wachem, Schouten, van den Bleek, Krishna and Sinclair2001).

When particles occupy a high volume fraction, frictional–collisional pressure must be added to $p_{s}$ to describe pressure-like forces that arise as layers of granular material slide over each other. The inclusion of frictional–collisional pressure has the effect of limiting compaction when ${\it\Theta}_{s}$ approaches zero. In highly packed conditions, the binary collision assumption in the Boltzmann equation, which is used to derive the granular-phase equation of state and transport properties, becomes invalid since particles are in contact with several other particles. In these conditions, high rates of inelastic particle collisions can rapidly reduce ${\it\Theta}_{s}$ to zero. This results in $p_{s}=0$ , which allows particles to pack to unphysically high levels. In such dense regions, frictional–collisional pressure is needed to provide the intergranular stress necessary to limit compaction. In this work we use the correlation of Johnson & Jackson (Reference Johnson and Jackson1987),

(3.23) $$\begin{eqnarray}p_{fric}=\left\{\begin{array}{@{}ll@{}}0\quad & \text{if }{\it\alpha}_{s}<{\it\alpha}_{s,crit},\\ 0.1{\it\alpha}_{s}\displaystyle \frac{({\it\alpha}_{s}-{\it\alpha}_{s,crit})^{2}}{({\it\alpha}_{s,max}-{\it\alpha}_{s})^{5}}\quad & \text{if }{\it\alpha}_{s}\geqslant {\it\alpha}_{s,crit},\end{array}\right.\end{eqnarray}$$

where ${\it\alpha}_{s,crit}$ is 0.5 unless otherwise noted and $p_{fric}$ is in units of Pa. Other expressions for frictional–collisional pressure could be used as well (Koo & Kuo Reference Koo and Kuo1977; Markatos & Kirkcaldy Reference Markatos and Kirkcaldy1983; Jenike Reference Jenike1987; Gidaspow Reference Gidaspow1994; Saurel & Abgrall Reference Saurel and Abgrall1999; Van Wachem et al. Reference Van Wachem, Schouten, van den Bleek, Krishna and Sinclair2001; Nussbaum et al. Reference Nussbaum, Helluy, Hérard and Carriére2006; Ling et al. Reference Ling, Wagner, Beresh, Kearney and Balachandar2012; Schneiderbauer, Aigner & Pirker Reference Schneiderbauer, Aigner and Pirker2012). Here the sum of the solids and frictional–collisional pressures is called the total intergranular stress,

(3.24) $$\begin{eqnarray}p_{s,tot}=p_{s}+p_{fric},\end{eqnarray}$$

to avoid ambiguity with the terminology used for interior ballistics.

(Another approach to limiting compaction of the solid phase is to limit the minimum granular temperature near the packing limit (Kamenetsky et al. Reference Kamenetsky, Goldshtein, Shapiro and Degani2000). This approach, while seemingly ad hoc, has a physical basis. The coefficient of restitution, $e$ , is not a constant, as assumed in many kinetic-theory granular multiphase models, but is a function of impact velocity. Higher impact velocities (higher ${\it\Theta}_{s}$ ) convert more PTE, $E_{s}$ , into internal energy, $e_{s}$ , through viscoelastic deformation (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004). Collisions become increasingly elastic ( $e\rightarrow 1$ ) as the granular temperature decreases. This effect limits both the rate at which ${\it\Theta}_{s}$ can decrease and its lower value. Recent granular gas models account for this to some extent (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004), but, to our knowledge, these have not been used in multiphase flow calculations.)

We do not include the mean kinetic energy in the conservation equation for PTE (equation (3.16)). This avoids small truncation errors of kinetic energy, which can lead to unphysical values of $E_{s}$ . Dissipation of $E_{s}$ due to inelastic particle collisions (see § 3.5) often reduces $E_{s}$ to the point where it is small compared to the mean kinetic energy, $\boldsymbol{v}_{s}\boldsymbol{\cdot }\boldsymbol{v}_{s}/2$ . In a numerical solution of these equations, small errors in the mean kinetic energy can cause small fluctuations of ${\it\Theta}_{s}$ , which, in dense-particle regimes, produces severe oscillations of $p_{s}$ and $p_{fric}$ . Such oscillations can degrade calculations to the point of failure. Decoupling kinetic energy from the PTE equation makes the numerical algorithm more robust. A similar approach is often used in astrophysical calculations (e.g. ASC FLASH Center (2012)). Developing a numerical algorithm that allows us to put the kinetic energy in the PTE equation while maintaining robustness is a topic of ongoing research.

The pressure term in the $E_{s}$ equation (3.16) does not include flow work from frictional–collisional pressure. This is consistent from a thermodynamics point of view. The $p_{fric}$ models used here are analogous to a barotropic equation of state for a fluid. It is possible to show, using thermodynamic arguments, that the total energy equation for a barotropic fluid simplifies to the mechanical energy equation (see appendix A). Performing a similar analysis for granular mixtures reveals that only the solids pressure, $p_{s}$ , contributes to changes in $E_{s}$ from compression and expansion (changes in volume fraction) of the particles. If the mean kinetic energy were included in the definition of $E_{s}$ , the friction pressure would enter the total PTE equation (3.18) as $\boldsymbol{v}_{s}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}p_{fric}$ , while the solids pressure would be present in the usual manner, $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{v}_{s}p_{s})$ (see (A 9)).

The Euler equations for the granular particles are compressible, similar to those for the gas phase. One key difference, mentioned above, is that this compressibility manifests itself through changes in particle volume fraction in the form of compaction waves. This granular-compaction wave speed can be derived for the particulate phase in the same way that the sound speed can be derived for a real gas (Serna & Marquina Reference Serna and Marquina2005). The compaction wave speed for a general granular gas, where friction pressure is only a function of ${\it\alpha}_{s}$ , is

(3.25) $$\begin{eqnarray}c_{s}^{2}=\frac{1}{{\it\rho}_{s}}\left[\left.\frac{\partial p_{s}^{tot}}{\partial {\it\alpha}_{s}}\right|_{{\it\Theta}_{s}}+\frac{2}{3}\frac{{\it\Theta}_{s}\left(\left.\displaystyle \frac{\partial p_{s}^{tot}}{\partial {\it\Theta}_{s}}\right|_{{\it\alpha}_{s}}\right)^{2}}{{\it\rho}_{s}{\it\alpha}_{s}^{2}}\right].\end{eqnarray}$$

Then by using the definitions for $p_{s}$ and $p_{fric}$ (see (3.21) and (3.23)), we find

(3.26) $$\begin{eqnarray}c_{s}^{2}={\it\Theta}_{s}(A+{\textstyle \frac{2}{3}}A^{2}+{\it\alpha}_{s}B)+c_{fric}^{2},\end{eqnarray}$$

where

(3.27a−c ) $$\begin{eqnarray}A=1+2(1+e){\it\alpha}_{s}g_{0},\quad B=2(1+e)(g_{0}+{\it\alpha}_{s}g_{0}^{\prime }),\quad g_{0}^{\prime }=\frac{g_{0}^{2}}{3{\it\alpha}_{s,max}}\left(\frac{{\it\alpha}_{s,max}}{{\it\alpha}_{s}}\right)^{2/3}\!.\quad\end{eqnarray}$$

The frictional contribution to the compaction wave speed (in MKS units) is

(3.28) $$\begin{eqnarray}c_{fric}^{2}=\left\{\begin{array}{@{}ll@{}}0\quad & \text{if }{\it\alpha}_{s}<{\it\alpha}_{s,crit},\\ \displaystyle \frac{1}{{\it\rho}_{s}}\frac{({\it\alpha}_{s}-{\it\alpha}_{s,crit})}{({\it\alpha}_{s,max}-{\it\alpha}_{s})^{5}}\left[{\it\alpha}_{s}\left(\frac{1}{5}+\frac{1}{2}\frac{{\it\alpha}_{s}-{\it\alpha}_{s,crit}}{{\it\alpha}_{s,max}-{\it\alpha}_{s}}\right)+\frac{{\it\alpha}_{s}-{\it\alpha}_{s,crit}}{10}\right]\quad & \text{if }{\it\alpha}_{s}\geqslant {\it\alpha}_{s,crit}.\end{array}\right.\end{eqnarray}$$

Several models are available for $p_{s}$ , $p_{fric}$ and $g_{0}$ , as mentioned above (see e.g. Kamenetsky et al. Reference Kamenetsky, Goldshtein, Shapiro and Degani2000; Van Wachem et al. Reference Van Wachem, Schouten, van den Bleek, Krishna and Sinclair2001; Nussbaum et al. Reference Nussbaum, Helluy, Hérard and Carriére2006). Typically, the coefficient of restitution, $e$ , the random close-packing limit, ${\it\alpha}_{s,max}$ , and the threshold volume fractions above which $p_{fric}$ becomes non-zero, ${\it\alpha}_{s,crit}$ , are input parameters. Exploring alternative intergranular stress models and the influence of input parameters is ongoing work.

3.4. The nozzling and $pDV$ work terms

In granular multiphase flows, it is common to choose $p_{int}=p_{g}$ (Gidaspow Reference Gidaspow1994). If we further assume that the particle density ( ${\it\rho}_{s}$ ) is not a function of pressure or intergranular stress, the $pDV$ work term reduces to

(3.29) $$\begin{eqnarray}-p_{int}\frac{\partial {\it\alpha}_{g}}{\partial t}=p_{g}\frac{\partial {\it\alpha}_{s}}{\partial t}=-p_{g}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\alpha}_{s}\boldsymbol{v}_{s}-p_{g}\frac{S^{M}}{{\it\rho}_{s}}.\end{eqnarray}$$

With these expressions, the final momentum equations for the gas and granular phase are

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}\boldsymbol{v}_{g})+\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}=p_{g}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}+\boldsymbol{S}^{\boldsymbol{p}}, & \displaystyle\end{eqnarray}$$

(3.31) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s}\boldsymbol{v}_{s})+\boldsymbol{{\rm\nabla}}p_{s}+\boldsymbol{{\rm\nabla}}p_{fric}=-{\it\alpha}_{s}\boldsymbol{{\rm\nabla}}p_{g}+{\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{g}-\boldsymbol{S}^{\boldsymbol{p}}, & \displaystyle\end{eqnarray}$$

(3.32) $$\begin{eqnarray}\frac{\partial {\it\alpha}_{g}{\it\rho}_{g}E_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[{\it\alpha}_{g}\boldsymbol{v}_{g}({\it\rho}_{g}E_{g}+p_{g})]=-p_{g}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\alpha}_{s}\boldsymbol{v}_{s}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{v}_{g}+S_{g}^{E},\end{eqnarray}$$

where the $pDV$ work associated with phase change is absorbed into $S_{g}^{E}$ .

3.5. Interphase coupling terms

Source terms that couple the governing equations for the gas and granular phases (aside from the nozzling and $pDV$ work terms) include drag, convective heat transfer, phase change and dissipation of PTE from inelastic collisions, and lifting forces. We neglect the virtual mass and Basset (history) forces (even though they can be important when computing particle trajectories in some high-speed applications (Ling, Haselbacher & Balachandar Reference Ling, Haselbacher and Balachandar2011 a ,Reference Ling, Haselbacher and Balachandar b )). Selected variable definitions and the source terms and their physical meaning are given in table 1.

Table 1.

Abbreviated list of variables and their physical interpretation.

The mass exchange term due to phase change is

(3.33) $$\begin{eqnarray}S^{M}={\dot{M}},\end{eqnarray}$$

where the total mass rate of phase change is

(3.34) $$\begin{eqnarray}{\dot{M}}=\mathop{\sum }_{i=1}^{N_{g}}S_{g,i}^{Y}=-\mathop{\sum }_{i=1}^{N_{s}}S_{s,i}^{Y},\end{eqnarray}$$

and $S_{g,i}^{Y}$ and $S_{s,i}^{Y}$ are the mass rates of production for species $i$ for the gas and granular phases due to phase change. The momentum source terms resulting from phase change, drag and lift are

(3.35) $$\begin{eqnarray}\boldsymbol{S}^{\boldsymbol{p}}=\boldsymbol{f}_{Drag}+\boldsymbol{f}_{Lift}+\boldsymbol{v}_{int}{\dot{M}}.\end{eqnarray}$$

The drag force for the particulate mixture is given by the correlation of Gidaspow (Reference Gidaspow1994),

(3.36) $$\begin{eqnarray}\boldsymbol{f}_{Drag}=K_{sg}(\boldsymbol{v}_{s}-\boldsymbol{v}_{g}),\end{eqnarray}$$

where

(3.37) $$\begin{eqnarray}K_{sg}=\left\{\begin{array}{@{}ll@{}}0.75C_{d}\displaystyle \frac{{\it\rho}_{g}{\it\alpha}_{g}{\it\alpha}_{s}|\boldsymbol{v}_{s}-\boldsymbol{v}_{g}|}{d_{s}{\it\alpha}_{g}^{2.65}}\quad & \text{if }{\it\alpha}_{g}\geqslant 0.8,\\ \displaystyle 150\frac{{\it\alpha}_{s}^{2}{\it\mu}_{g}}{{\it\alpha}_{g}d_{s}^{2}}+1.75\frac{{\it\rho}_{g}{\it\alpha}_{s}|\boldsymbol{v}_{g}-\boldsymbol{v}_{s}|}{d_{s}}\quad & \text{if }{\it\alpha}_{g}<0.8,\end{array}\right.\end{eqnarray}$$

the drag coefficient, $C_{d}$ , is

(3.38) $$\begin{eqnarray}C_{d}=\left\{\begin{array}{@{}ll@{}}\displaystyle 24({\it\alpha}_{g}Re)^{-1}[1+0.15({\it\alpha}_{g}Re)^{0.687}]\quad & \text{if }{\it\alpha}_{g}Re<1000,\\ 0.44\quad & \text{if }{\it\alpha}_{g}Re\geqslant 1000,\end{array}\right.\end{eqnarray}$$

and the Reynolds number is defined by

(3.39) $$\begin{eqnarray}Re=\displaystyle \frac{{\it\rho}_{g}|\boldsymbol{v}_{g}-\boldsymbol{v}_{s}|d_{s}}{{\it\mu}_{g}}.\end{eqnarray}$$

Other forms of the drag coefficient are available in the literature (Syamlal et al. Reference Syamlal, Rogers and O’Brien1993; Van Wachem et al. Reference Van Wachem, Schouten, van den Bleek, Krishna and Sinclair2001; Parmar, Haselbacher & Balachandar Reference Parmar, Haselbacher and Balachandar2010), and should be tested, but this is not the subject of this paper

The lift force, $\boldsymbol{f}_{Lift}$ , is modelled using the Magnus lift force defined by

(3.40) $$\begin{eqnarray}\boldsymbol{f}_{Lift}=C_{l}{\it\alpha}_{s}{\it\rho}_{g}(\boldsymbol{v}_{s}-\boldsymbol{v}_{g})\times (\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}_{g}),\end{eqnarray}$$

where the lift coefficient, $C_{l}$ , has a value of 0.5 (Drew & Lahey Reference Drew and Lahey1987). The influence of the Saffman lift force is currently neglected.

The interfacial velocity used for momentum exchange and energy during phase change are defined using the value from the phase that is losing mass,

(3.41) $$\begin{eqnarray}\boldsymbol{v}_{int}=\left\{\begin{array}{@{}ll@{}}\boldsymbol{v}_{s}\quad & \text{if }{\dot{M}}>0,\\ \boldsymbol{v}_{g}\quad & \text{if }{\dot{M}}<0.\end{array}\right.\end{eqnarray}$$

The source terms in the energy equations are

(3.42) $$\begin{eqnarray}\displaystyle & \displaystyle S_{g}^{E}=-q_{conv}+{\it\phi}_{visc}-{\it\phi}_{slip}+E_{g,int}{\dot{M}}+(\boldsymbol{f}_{Drag}+\boldsymbol{f}_{Lift})\boldsymbol{\cdot }\boldsymbol{v}_{s}, & \displaystyle\end{eqnarray}$$

(3.43) $$\begin{eqnarray}\displaystyle & \displaystyle S_{s}^{E}=q_{conv}-e_{s,int}{\dot{M}}+\dot{{\it\gamma}}, & \displaystyle\end{eqnarray}$$

(3.44) $$\begin{eqnarray}\displaystyle & \displaystyle S_{s}^{PTE}=-\dot{{\it\gamma}}-{\it\phi}_{visc}+{\it\phi}_{slip}-E_{s,int}{\dot{M}}, & \displaystyle\end{eqnarray}$$

$E_{g,int}=H_{g,int}-p_{int}/{\it\rho}_{s}$

$\dot{{\it\gamma}}$

${\it\phi}_{visc}$

${\it\phi}_{slip}$

Figure 1.

Relationship between the energy transfer mechanisms modelled in the source terms to the energy equations, where $\boldsymbol{f}_{Net}=\boldsymbol{f}_{Drag}+\boldsymbol{f}_{Lift}$ and ${\rm\Delta}\boldsymbol{v}=\boldsymbol{v}_{s}-\boldsymbol{v}_{g}$ .

The convective heat flux between the gas and particles is

(3.45) $$\begin{eqnarray}q_{conv}=6\frac{{\it\alpha}_{s}{\it\lambda}_{g}Nu}{d_{s}^{2}}(T_{g}-T_{s}),\end{eqnarray}$$

where the Nusselt number is estimated using the correlation of Gunn (Reference Gunn1978),

(3.46) $$\begin{eqnarray}Nu=(7-10{\it\alpha}_{g}+5{\it\alpha}_{g}^{2})(1+0.7Re^{0.2}Pr_{g}^{1/3})+(1.33-2.4{\it\alpha}_{g}+1.2{\it\alpha}_{g}^{2})Re^{0.7}Pr_{g}^{1/3},\end{eqnarray}$$

and $Pr_{g}$ is the gas-phase Prandtl number.

The model for viscous damping of PTE is adopted from Gidaspow (Reference Gidaspow1994),

(3.47) $$\begin{eqnarray}{\it\phi}_{visc}=3K_{sg}{\it\Theta}_{s}.\end{eqnarray}$$

The rate of internal energy production into the gas phase due to frictional heating during particle acceleration is

(3.48) $$\begin{eqnarray}q_{Drag}={\rm\Delta}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{f}_{Net},\end{eqnarray}$$

where $\boldsymbol{f}_{Net}=\boldsymbol{f}_{Drag}+\boldsymbol{f}_{Lift}$ and ${\rm\Delta}\boldsymbol{v}=\boldsymbol{v}_{s}-\boldsymbol{v}_{g}$ . Not all of this internal energy heats the gas as it flows around the particles. Some of this heating is diverted into producing random translational motion (PTE) of the particles. The production rate of PTE during particle acceleration, ${\it\phi}_{slip}$ , is given by Koch & Sangani (Reference Koch and Sangani1999) with their correction factor fixed at unity,

(3.49) $$\begin{eqnarray}{\it\phi}_{slip}=\frac{81{\it\alpha}_{s}{\it\mu}_{g}^{2}}{g_{0}d_{s}^{3}{\it\rho}_{s}\sqrt{{\rm\pi}}}\frac{|\boldsymbol{v}_{g}-\boldsymbol{v}_{s}|^{2}}{\sqrt{{\it\Theta}_{s}}}.\end{eqnarray}$$

Alternative models for ${\it\phi}_{slip}$ could also be used (e.g. Koch & Sangani Reference Koch and Sangani1999), but again these are not the subject of this paper.

The granular dissipation model adopted for this study (Lun et al. Reference Lun, Savage, Jeffrey and Chepurniy1984) is a variant of Haff’s cooling law (Haff Reference Haff1983; Brilliantov & Pöschel Reference Brilliantov and Pöschel2004),

(3.50) $$\begin{eqnarray}\dot{{\it\gamma}}=\frac{12(1-e^{2})g_{0}{\it\alpha}_{s}^{2}}{\sqrt{{\rm\pi}}}\displaystyle \frac{{\it\rho}_{s}{\it\Theta}_{s}^{3/2}}{d_{s}}.\end{eqnarray}$$

Similar to $\boldsymbol{v}_{int}$ , the interfacial energy exchange terms for phase change are defined by the phase that is losing mass,

(3.51a,b ) $$\begin{eqnarray}E_{g,int}=\left\{\begin{array}{@{}ll@{}}e_{s}+E_{s}+\boldsymbol{v}_{s}\boldsymbol{\cdot }\boldsymbol{v}_{s}/2\quad & \text{if }{\dot{M}}>0,\\ E_{g}\quad & \text{if }{\dot{M}}<0,\end{array}\right.\quad E_{s,int}=\left\{\begin{array}{@{}ll@{}}E_{s}\quad & \text{if }{\dot{M}}>0,\\ 0\quad & \text{if }{\dot{M}}<0,\end{array}\right.\end{eqnarray}$$

(3.52) $$\begin{eqnarray}e_{s,int}=\left\{\begin{array}{@{}ll@{}}e_{s}\quad & \text{if }{\dot{M}}>0,\\ E_{g}-\boldsymbol{v}_{g}\boldsymbol{\cdot }\boldsymbol{v}_{g}/2\quad & \text{if }{\dot{M}}<0.\end{array}\right.\end{eqnarray}$$

4.1. Solution of the hyperbolic terms

The gas-phase hyperbolic terms, which include gravitational acceleration to ensure hydrostatic equilibrium, are

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}Y_{g,j}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}Y_{g,j}\boldsymbol{v}_{g})=0, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}\boldsymbol{v}_{g})+\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}=p_{g}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}+\boldsymbol{f}_{Lift}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{g}{\it\rho}_{g}E_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}H_{g})=-p_{g}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\alpha}_{s}\boldsymbol{v}_{s}+\boldsymbol{f}_{Lift}\boldsymbol{\cdot }\boldsymbol{v}_{s}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{v}_{g}. & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s})=0, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}Y_{s,j}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}Y_{s,j}\boldsymbol{v}_{s})=0, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{v}_{s}\boldsymbol{v}_{s})+\boldsymbol{{\rm\nabla}}p_{s}+\boldsymbol{{\rm\nabla}}p_{fric}=-{\it\alpha}_{s}\boldsymbol{{\rm\nabla}}p_{g}-\boldsymbol{f}_{Lift}+{\it\alpha}_{s}{\it\rho}_{s}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}E_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}E_{s}\boldsymbol{v}_{s})=-p_{s}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}_{s}, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{s}{\it\rho}_{s}e_{s}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{s}{\it\rho}_{s}e_{s}\boldsymbol{v}_{s})=0. & \displaystyle\end{eqnarray}$$

The semi-discrete form of these equations, implementation details and step-by-step solution procedure are given in appendix B. The lift force is discretized using standard second-order differencing if $\min ({\it\alpha}_{s})>{\it\alpha}_{s,min}$ for all points in the stencil where ${\it\alpha}_{min}=10^{-10}$ . A fifth-order numerical method based on Houim & Kuo (Reference Houim and Kuo2011) is used for spatial discretization. Overcoming the challenges associated with solving the equations with non-conservative terms and the granular fluxes is now discussed.

4.1.1. Non-conservative nozzling and $pDV$ terms

The left-hand side of the governing equation for each phase is in a form that is solvable by many standard numerical algorithms, such as Godunov-style methods, as long as they are solved independently for each phase and also provided that the nozzling term ( $p_{g}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}$ ) in the momentum equation and the $pDV$ work term ( $p_{g}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\alpha}_{s}\boldsymbol{v}_{s}$ ) for the energy equation are treated as separate source terms. This is a straightforward approach that works well enough if granular interfaces vary smoothly in space, but can fail when there is a steep granular interface (e.g. the edge of a dust layer). The major challenge to obtaining a good numerical solution centres around problems with including the non-conservative nozzling and $pDV$ work terms and enforcing the ‘pressure non-disturbing condition’ (Saurel & Abgrall Reference Saurel and Abgrall1999). This usually cannot be maintained rigorously with a straightforward approach (as discussed in the Introduction) and unphysical pressure oscillations may develop.

An approach that has been shown to work for coupling compressible gases to compressible liquids with sharp material interfaces is the stratified-flow concept developed by Chang & Liou (Reference Chang and Liou2007). They consider control volumes for the gas and liquid phases separately. In this approach each computational cell edge is divided into gas–gas, liquid–liquid and gas–liquid sections. The hyperbolic fluxes in each section are computed independently with an appropriate Riemann solver. The flux divergence and pressure gradient terms are assembled from the solution all of these Riemann problems. The application of this concept to granular multiphase flows is illustrated in figure 2.

Figure 2.

Application of the stratified-flow concept to granular multiphase flows in a control volume used to sum the pressure forces. The computational cell index is $i$ and the right and left edges are denoted by $i\pm 1/2$ , respectively. (a) Actual control volume. (b) Control volume with particles collapsed.

The stratified-flow concept can be simplified substantially for a granular flow by using the incompressibility assumption of the particles. The volume fraction (while integrating the hyperbolic terms) is controlled by motion of the incompressible particles. This implies that the volume fraction is not a function of the gas-phase state. Thus, the product rule can be applied to the pressure gradient term ( $\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}$ ), which simplifies the gas-phase momentum equation to

(4.10) $$\begin{eqnarray}\frac{\partial {\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{v}_{g}\boldsymbol{v}_{g})=-{\it\alpha}_{g}\boldsymbol{{\rm\nabla}}p_{g}+\boldsymbol{f}_{Lift}+{\it\alpha}_{g}{\it\rho}_{g}\boldsymbol{g}.\end{eqnarray}$$

It is straightforward to satisfy the pressure non-disturbing condition with this simplification because there is only one pressure term to compute rather than having to perfectly balance two separate pressure terms.

This simplification used to produce (4.10) is only rigorously possible to make if the particles are assumed to be incompressible. If the particles are compressible, the solid material density would be a function of $T_{s}$ and the total isotropic stress acting on the particle (i.e. ${\it\rho}_{s}={\it\rho}_{s}(T_{s},p_{g},p_{s,tot})$ ). Conservation of particle mass would then dictate that the solid volume fraction is implicitly a function of ${\it\rho}_{s}$ , since ${\it\alpha}_{s}{\it\rho}_{s}$ is a conserved variable. For example, the high pressure behind a gas-phase shock would increase ${\it\rho}_{s}$ , with a resulting decrease in particle diameter, and, in turn, decrease ${\it\alpha}_{s}$ . Thus, a gas-phase shock would cause a discontinuous decrease of ${\it\alpha}_{s}$ in addition to a pressure jump. As a result, the quantity ${\it\alpha}_{g}p_{g}$ is, in general, non-differentiable and $\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}\neq {\it\alpha}_{g}\boldsymbol{{\rm\nabla}}p_{g}+p_{g}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}$ across a shock. However, if particle incompressibility is assumed, then ${\it\rho}_{s}={\it\rho}_{s}(T_{s})$ and $\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}p_{g}={\it\alpha}_{g}\boldsymbol{{\rm\nabla}}p_{g}+p_{g}\boldsymbol{{\rm\nabla}}{\it\alpha}_{g}$ , since ${\it\alpha}_{s}$ would then not be an implicit function of the gas-phase pressure and, as a result, would not change discontinuously due to compression behind a gas-phase shock.

One issue that merging the pressure gradient and nozzling terms introduces in (4.10) is that there are no portions of the gas-phase hyperbolic terms that resemble the Euler equations. As a result, standard Godunov methods cannot be used. Despite this, the control volume shown in figure 2 is still valid and the solution to the Riemann problem on the gas–gas fraction of the cell face where gas flows between the particles is physically correct. We use this observation and solve a gas-phase Riemann problem at the section marked ${\rm\Delta}{\it\alpha}_{flow}$ in figure 2. The pressure computed from this Riemann problem ( $p_{g,i+1/2}$ ) is stored at each cell edge and used to compute the non-conservative pressure term ${\it\alpha}_{g}\boldsymbol{{\rm\nabla}}p_{g}$ at each computational cell edge. This is computed in the $x$ -direction by

(4.11) $$\begin{eqnarray}{\it\alpha}_{g}\frac{\partial p_{g}}{\partial x}\approx {\it\alpha}_{g,i}\frac{p_{g,i+1/2}-p_{g,i-1/2}}{{\rm\Delta}x},\end{eqnarray}$$

where ${\rm\Delta}x$ is the width of a computational cell.

The only requirement for evaluating the gas-phase Riemann problem is that pressure is explicitly calculated and is not directly included in the momentum flux. The HLLC (Harten, Lax & van Leer Reference Harten, Lax and van Leer1983; Toro, Spruce & Speares Reference Toro, Spruce and Speares1994) and AUSM (Liou Reference Liou2006) flux schemes can be modified to match this condition. The HLLC flux is used in this work (see appendix B.2). The remaining flux derivatives are computed in a similar manner. For example, the enthalpy flux divergence is

(4.12) $$\begin{eqnarray}\frac{\partial {\it\alpha}_{g}{\it\rho}_{g}u_{g}H_{g}}{\partial x}\approx \frac{{\it\alpha}_{g,i+1/2}{\it\rho}_{g,i+1/2}u_{g,i+1/2}H_{g,i+1/2}-{\it\alpha}_{g,i-1/2}{\it\rho}_{g,i-1/2}u_{g,i-1/2}H_{g,i-1/2}}{{\rm\Delta}x},\end{eqnarray}$$

where ${\it\rho}_{g,i+1/2}$ , $u_{g,i+1/2}$ and $H_{g,i+1/2}$ come from the HLLC Riemann solver and ${\it\alpha}_{g,i+1/2}$ is computed using the flux scheme for the granular phase discussed below. More details of the gas-phase solution for the hyperbolic terms can be found in appendix B.

4.1.2. Granular fluxes

Calculation of the granular flux has its own particular set of complications. One issue is its extreme stiffness as the packing limit is approached and the friction pressure and the compaction wave speed increase towards infinity. In such regions, any minute fluctuation of solid volume fraction can produce extremely large changes in friction pressure. The solution we are using to avoid numerical instability is to smoothly degrade the order of the numerical method from fifth-order to first-order as the packing limit is approached. This helps to maintain numerical stability and monotonic solutions of the granular conservative variables. The details of this process are discussed in appendix B.4.

One of the most difficult issues is that the compaction wave speed and solid pressure become zero (1) when the granular temperature is zero or (2) at low volume fractions, resulting in hyperbolic degeneracy (as discussed previously). One approximate Riemann solver used for the particulate phase for a pressureless ( $p_{s,tot}=0$ ) dust model is from Collins et al. (Reference Collins, Ferguson, Chien, Kuhl, Krispin and Glaz1994):

(4.13) $$\begin{eqnarray}\boldsymbol{F}_{s,i+1/2}=\left\{\begin{array}{@{}ll@{}}\boldsymbol{F}_{s}(\boldsymbol{U}_{s}^{L})\quad & \text{if }u_{s}^{L}\geqslant 0\text{ and }u_{s}^{R}>0,\\ \boldsymbol{F}_{s}(\boldsymbol{U}_{s}^{R})\quad & \text{if }u_{s}^{L}<0\text{ and }u_{s}^{R}\leqslant 0,\\ \boldsymbol{F}_{s}(\boldsymbol{U}_{s}^{L})+\boldsymbol{F}_{s}(\boldsymbol{U}_{s}^{R})\quad & \text{if }u_{s}^{L}>0\text{ and }u_{s}^{R}<0,\\ 0\quad & \text{if }u_{s}^{L}<0\text{ and }u_{s}^{R}>0,\end{array}\right.\end{eqnarray}$$

where $\boldsymbol{U}_{s}^{L,R}$ are the left- and right-reconstructed conserved variable vector for the solid phase at edge $i+1/2$ and $\boldsymbol{F}_{s}$ is the flux vector of the particle phase, excluding solids pressure and the non-conservative terms. In other locations, where the granular phase is not pressureless, (4.13) is incorrect and a more usual Riemann solver will work. Ideally, however, the same Riemann solver would handle both cases without any ad hoc switches, which is the case for the AUSM $^{+}$ -up solver (Liou Reference Liou2006). In addition to reducing the order of the numerical method to first order in highly packed regions, an additional dissipative term is added to the AUSM flux for the granular phase. Numerical details of the modified AUSM flux for the granular phase are given in § B.4.

$\text{AUSM}^{+}$ -up is also modified to return $u_{s,i+1/2}$ needed for $p_{s}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}_{s}$ , the quantity ${\it\alpha}_{s,i+1/2}u_{s,i+1/2}$ for the $pDV$ work term ( $p_{g}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\alpha}_{s}\boldsymbol{v}_{s}$ ), and ${\it\alpha}_{g,i+1/2}$ needed to complete the gas-phase fluxes. Computation of the gas-phase edge volume fraction is based on where the particles have moved at $t=0^{+}$ , as shown in figure 3, and the gas-phase volume fraction at the cell edge is computed by

(4.14) $$\begin{eqnarray}{\it\alpha}_{g,i+1/2}=\left\{\begin{array}{@{}ll@{}}{\it\alpha}_{g,i+1/2}^{L}\quad & \text{if }{\dot{m}}_{s,i+1/2}>0,\\ {\it\alpha}_{g,i+1/2}^{R}\quad & \text{if }{\dot{m}}_{s,i+1/2}\leqslant 0,\end{array}\right.\end{eqnarray}$$

where ${\dot{m}}_{s,i+1/2}$ is the granular mass flux calculated from the AUSM $^{+}$ -up flux (see § B.4 for details). Similarly,

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\alpha}_{s,i+1/2}u_{s,i+1/2}=\left\{\begin{array}{@{}ll@{}}{\dot{m}}_{s,i+1/2}/{\it\rho}_{s,i+1/2}^{L}\quad & \text{if }{\dot{m}}_{s,i+1/2}>0,\\ {\dot{m}}_{s,i+1/2}/{\it\rho}_{s,i+1/2}^{R}\quad & \text{if }{\dot{m}}_{s,i+1/2}\leqslant 0,\end{array}\right. & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle u_{s,i+1/2}=\left\{\begin{array}{@{}ll@{}}{\dot{m}}_{s,i+1/2}/{\it\alpha}_{s,i+1/2}^{L}{\it\rho}_{s,i+1/2}^{L}\quad & \text{if }{\dot{m}}_{s,i+1/2}>0,\\ {\dot{m}}_{s,i+1/2}/{\it\alpha}_{s,i+1/2}^{R}{\it\rho}_{s,i+1/2}^{R}\quad & \text{if }{\dot{m}}_{s,i+1/2}\leqslant 0.\end{array}\right. & \displaystyle\end{eqnarray}$$

Then the pressure term in (4.8) is approximated as

(4.17) $$\begin{eqnarray}p_{s}\frac{\partial u_{s}}{\partial x}\approx p_{s,i}\frac{u_{s,i+1/2}-u_{s,i-1/2}}{{\rm\Delta}x},\end{eqnarray}$$

and the $pDV$ term can be approximated as

(4.18) $$\begin{eqnarray}p_{g}\frac{\partial {\it\alpha}_{s}u_{s}}{\partial x}\approx p_{g,i}\frac{{\it\alpha}_{s,i+1/2}u_{s,i+1/2}-{\it\alpha}_{s,i-1/2}u_{s,i-1/2}}{{\rm\Delta}x}.\end{eqnarray}$$

Figure 3.

Schematic of the gas-phase volume fraction at a cell face at $t=0^{+}$ .

5.1. Advection of a particle curtain

Here we examine the ability of the model and numerical algorithm to preserve the pressure non-disturbing condition near a steep granular interface. We advect a curtain of particles with a volume fraction of 40 % in a domain with uniform velocity, pressure and temperature of $100~\text{m}~\text{s}^{-1}$ , 1 atm and 300 K, respectively. To increase the complexity of the calculation, we also place a bubble of pure helium where the particles are located. The gas surrounding the He bubble is pure nitrogen. Specifically, the initial conditions are: $Y_{\text{N}_{2}}=1$ and ${\it\alpha}_{s}=0.4$ if $0.4<x<0.6$ and $Y_{\text{He}}=1$ and ${\it\alpha}_{s}=0$ otherwise. The domain was discretized with 200 grid points per metre. Inflow and outflow boundary conditions were used on the left and right boundaries, respectively. We advect the combined granular and gaseous material interface. The particle and helium curtain should be passively advected downstream without changing the pressure and temperature fields. The computed solution, shown in figure 4, indicates that temperature and pressure error are negligible and that the pressure non-disturbing condition is met.

Figure 4.

Advection of a granular and multicomponent gaseous interface with uniform temperature for a distance of 1 m. (a) Comparison of gas-phase density and ${\it\alpha}_{s}$ to the exact solution and (b) gas-phase temperature and error in pressure.

5.2. Two-phase granular shock-tube problem

In this section we test the ability of the model and solution procedure to solve shock-tube problems with low particle volume fraction. The granular-phase solids and frictional–collisional pressures are nearly zero due to the low particle volume fraction. As discussed above, this creates hyperbolically degenerate conditions that are difficult for most standard Riemann solvers to deal with. It is a situation that we must test to ensure that the computed results agree with reference calculations that used specialized Riemann solvers developed for pressureless dust. The initial conditions follow (Saito et al. Reference Saito, Marumoto and Takayama2003; Fedorov et al. Reference Fedorov, Kharlamova and Khmel’2007):

(5.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}p_{g}^{L}=10~\text{atm},\quad p_{g}^{R}=1~\text{atm},\\ T^{L}=270~\text{K},\quad T^{R}=270~\text{K},\\ Y_{g,air}^{L}=1,\quad Y_{g,air}^{R}=1,\\ {\it\alpha}_{s}^{L}=0,\quad {\it\alpha}_{s}^{R}=5.172\times 10^{-4},\\ {\it\Theta}_{s}^{L}=0,\quad {\it\Theta}_{s}^{R}=0.\end{array}\right\}\end{eqnarray}$$

The density, diameter and specific heat of the particles are ${\it\rho}_{s}=2500~\text{kg}~\text{m}^{-3}$ , $d_{s}=10~{\rm\mu}\text{m}$ and $C_{V,s}=718~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ , respectively. The domain is $0.257798$  m in length and the diaphragm is placed at $0.129$  m. The solution at $184~{\rm\mu}$ s after the diaphragm is ruptured is shown in figure 5 for two grid spacings. The physics of dust–gas shock-tube problems are well known (Miura & Glass Reference Miura and Glass1982; Saito et al. Reference Saito, Marumoto and Takayama2003); nevertheless, there are no exact solutions that can be used for verification. Instead, the computed results can be compared with reference solutions from Saito et al. (Reference Saito, Marumoto and Takayama2003) and Fedorov et al. (Reference Fedorov, Kharlamova and Khmel’2007). The results, shown in figure 5, compare well with computed results from these references, which indicates that our modified AUSM $^{+}$ -up scheme used to compute the granular flux (see § B.4) is accurate when the dust is nearly pressureless.

Figure 5.

Computed (a) gas-phase pressure and (b) bulk densities for the gas and granular phases for the two-phase granular Riemann problem at $184~{\rm\mu}\text{s}$ . The densities were scaled by the initial gas-phase density to the right of the interface. The length scale was scaled so that the domain ranges from 0 to 10 to match the output used by Fedorov et al. (Reference Fedorov, Kharlamova and Khmel’2007).

The burst of the diaphragm sends a shock wave into the low-pressure air and particle mixture, a rarefaction wave into the high-pressure side, and a gas-phase contact surface that passes through the particles. Drag forces between the post-shock gas and particles locally decelerate the gas and produce a subsequent rise in pressure. The pressure increase is due to high-speed gas from regions without particles flowing into and filling a region where gas has been decelerated by drag. The local increase in pressure produces acoustic waves that propagate towards the rarefaction and cause the hook-like structure at $x\sim 5$ . Mechanical and thermal losses of the post-shock gases due to drag also have the effect of weakening the shock, which is illustrated in the next test problem. The two-phase shock wave formed in this test problem corresponds to a Type I structure characterized by Khmel’ & Fedorov (Reference Khmel’ and Fedorov2014 a ,Reference Khmel’ and Fedorov b ), where there is a discontinuous gas-phase shock that propagates through the particles with smooth relaxation of granular properties, such as volume fraction, behind the shock.

5.3. Two-phase granular shock-tube problem with a gaseous He– $N_{2}$ interface

This problem is similar to the test presented in the previous section, but now we add a gaseous material interface. This examines whether the model and numerical method can produce the correct weak solution for a multiphase shock-tube problem that contains a gaseous material interface with varying ${\it\gamma}$ . Such a situation may exist in a dust explosion when shocks interact with thin dust flames that can separate reactants and products, which have different specific heat ratios. Here the high-pressure side contains pure helium and the low-pressure side contains a mixture of particles and pure N $_{2}$ . The initial states for this test case are

(5.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}p_{g}^{L}=10~\text{atm},\quad p_{g}^{R}=1~\text{atm},\\ T^{L}=300~\text{K},\quad T^{R}=300~\text{K},\\ Y_{g,\text{He}}^{L}=1,\quad Y_{g,\text{He}}^{R}=0,\\ Y_{g,\text{N}_{2}}^{L}=0,\quad Y_{g,\text{N}_{2}}^{R}=1,\\ {\it\alpha}_{s}^{L}=0,\quad {\it\alpha}_{s}^{R}=4.555\times 10^{-4},\\ {\it\Theta}_{s}^{L}=0,\quad {\it\Theta}_{s}^{R}=0.\end{array}\right\}\end{eqnarray}$$

The density, specific heat and restitution coefficient of the particles are $2500~\text{kg}~\text{m}^{-3}$ , $745~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ and $0.999$ , respectively. The domain is $1$  m in length and the diaphragm is placed at $0.5$  m. The volume fraction of the right state was chosen so that the bulk densities of the gas and granular phases were the same. Calculations were done for two different particle diameters, 10 and $25~{\rm\mu}\text{m}$ . The equations were integrated up to $400~{\rm\mu}\text{s}$ of physical time and 400 grid points were used in each calculation. A comparison of the computed results to the exact solution for this Riemann problem when particles are not present is shown in § B.2.

The results, shown in figure 6, are qualitatively similar to the previous test problem. The gas-phase pressure in multiphase shock tubes is actually higher if particles are present, even though the shock is weaker (smaller pressure jump and lower speed). Both of these are side effects of mechanical losses to the post-shock gas from drag as discussed earlier. The locally higher value of the solid bulk density and velocity at ${\sim}0.65$  m is caused by a sharp change in the drag force at the He– $\text{N}_{2}$ contact surface. This causes a pile-up of particles on the nitrogen side of the gas-phase contact surface where the drag force is higher. The solids pressure is very small and has a negligible influence for the low volume fractions in this problem. Spurious oscillations are not present near the granular or He– $\text{N}_{2}$ contact surfaces, which indicates that the model and solution method can accurately solve two-phase granular shock-tube problems that contain mixtures of gases with varying specific heats.

Figure 6.

Computed (a) gas-phase and solids pressure, (b) velocity, (c) temperature and (d) bulk density for multiphase shock-tube problem with a He– $\text{N}_{2}$ interface at $400~{\rm\mu}\text{s}$ . The solid temperature and velocity are not shown if the particle volume fraction lies below $10^{-8}$ .

5.4. Dense granular two-phase shock-tube problem

The solids and frictional–collisional pressures were negligible in the previous problems. This test is for a two-phase shock-tube problem with a high driving pressure of 100 atm and a dense particle volume fraction of 40 %. This situation tests whether the model can predict relaxation structures of granular shocks (Kamenetsky et al. Reference Kamenetsky, Goldshtein, Shapiro and Degani2000). This case also stresses capabilities of the solution method to compute high-pressure gases interacting with dense granular regions that approach the packing limit with steep granular interfaces. The initial states for this problem are

(5.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}p_{g}^{L}=100~\text{atm},\quad p_{g}^{R}=1~\text{atm},\\ T^{L}=300~\text{K},\quad T^{R}=300~\text{K},\\ Y_{g,Air}^{L}=1,\quad Y_{g,Air}^{R}=1,\\ {\it\alpha}_{s}^{L}=0,\quad {\it\alpha}_{s}^{R}=0.4,\\ {\it\Theta}_{s}^{L}=0,\quad {\it\Theta}_{s}^{R}=0.\end{array}\right\}\end{eqnarray}$$

The domain measures 0.06 m in length and the diaphragm is placed at 0.03 m. The resulting shock structure is advanced for 100  ${\rm\mu}\text{s}$ on a grid with 1200 computational cells. The particle diameter, density, specific heat and coefficient of restitution were $5~{\rm\mu}\text{m}$ , $1470~\text{kg}~\text{m}^{-3}$ , $987~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ and 0.999, respectively. The volume fraction in this problem approaches the packing limit and additional dissipation was added to the granular phase by setting the hyperbolic dissipation parameter ( $\mathscr{D}$ ) to 1 (see (B 25) in § B.4). The solution at $100~{\rm\mu}\text{s}$ is shown in figure 7. The influence of $\mathscr{D}$ and coefficient of restitution (e) on the granular shock structure are shown in figures 8 and 9.

Figure 7.

Computed (a) pressure, (b) velocity, (c) temperature and (d) bulk density for the dense granular shock-tube problem at $100~{\rm\mu}\text{s}$ with $e=0.999$ and $\mathscr{D}=1$ .

Figure 8.

Effect of $\mathscr{D}$ and $e$ on calculated profiles of ${\it\alpha}_{s}$ and ${\it\Theta}_{s}$ near a granular shock: (a)  $e=0.999$ , (b)  $e=0.9$ .

Figure 9.

Effect of $\mathscr{D}$ and $e$ on calculated profiles of pressure and $p_{fric}/p_{s,tot}$ near a granular shock: (a)  $e=0.999$ , (b)  $e=0.9$ .

Unlike previous problems shown in §§ 5.1–5.3, particle–particle interactions are now dominant. After the diaphragm is ruptured, a strong granular shock is transmitted into the granular phase due to the combination of Archimedes forces ( ${\it\alpha}_{s}\boldsymbol{{\rm\nabla}}p_{g}$ ) and particle drag from the high-pressure gas flowing into the dense layer of particles. Frictional heating from gas initially penetrating the particle bed increases the gas and particle temperature by ${\sim}2.5$  K near the contact surface. The computed granular shock has a relaxation structure similar to those reported by Kamenetsky et al. (Reference Kamenetsky, Goldshtein, Shapiro and Degani2000). Sources of granular energy, $E_{s}$ , decrease rapidly away from the granular shock. Without these sources, collisional dissipation of $E_{s}$ ( $\dot{{\it\gamma}}$ ) decreases ${\it\Theta}_{s}$ and $p_{s}$ . Nevertheless, mechanical equilibrium needs to be maintained and the solid volume fraction rises to compensate for the reduction in ${\it\Theta}_{s}$ by packing even further to increase the friction pressure. Eventually, the granular temperature decreases to zero, and the total intergranular stress, $p_{s,tot}$ , is solely from friction pressure, as shown in figure 9. The granular temperature decreases much faster and the resulting granular shock structure is almost entirely due to friction pressure when the coefficient of restitution is lowered to 0.9.

A weak gas-phase shock-like structure is present at the location of the granular shock. This is due to the increase in particle volume fraction, which compresses the gas and produces a sharp rise in temperature and pressure along with the propagating granular shock. This corresponds to a Type II two-phase granular shock structure characterized by Khmel’ & Fedorov (Reference Khmel’ and Fedorov2014 a ,Reference Khmel’ and Fedorov b ), a case where a propagating granular shock drives a weak propagating gas-phase discontinuity along with it. This structure is in contrast to the Type I structure mentioned above, where a gas-phase shock drives a smooth compression of particles (increase of  ${\it\alpha}_{s}$ ).

A peculiarity of granular systems is that shocks form in granular systems even though $c_{s}$ and $p_{s,tot}$ can be zero in front of the shock. Assuming that the intergranular stress and particle velocity are zero upstream of the shock and $p_{s,tot}=p_{fric}$ , it can be shown through a Rankine–Hugoniot analysis that the shock and post-shock particle velocities for a purely granular system are

(5.4a,b ) $$\begin{eqnarray}u_{shock}=\sqrt{\frac{{\it\alpha}_{s,2}}{{\it\alpha}_{s,2}-{\it\alpha}_{s,1}}\frac{p_{fric}({\it\alpha}_{s,2})}{{\it\alpha}_{s,1}{\it\rho}_{s}}},\quad u_{s,2}=\sqrt{\frac{{\it\alpha}_{s,2}-{\it\alpha}_{s,1}}{{\it\alpha}_{s,2}}\frac{p_{fric}({\it\alpha}_{s,2})}{{\it\alpha}_{s,1}{\it\rho}_{s}}},\end{eqnarray}$$

where states 1 and 2 refer to conditions upstream and downstream of the shock, respectively. Applying the post-shock volume fraction shown in figure 8(b) ( ${\it\alpha}_{s,2}=0.639$ and ${\it\alpha}_{s,1}=0.4$ ) produces a shock velocity of $186~\text{m}~\text{s}^{-1}$ (corresponding to a position of 0.0486 m after $100~{\rm\mu}\text{s}$ ) and a particle velocity of $69.7~\text{m}~\text{s}^{-1}$ . These analytical values are close enough to the computed results shown in figures 7(b) and 9(b), considering that the gas phase was not taken into account when computing the analytical post-shock condition.

Small gas-phase temperature oscillations of the order of 1 K are present near the granular shock. These do not cause any numerical instabilities and quickly damp from heat transfer between the gas and granular phases. The bump in gas-phase velocity near the solid contact surface in figure 7 is caused by a reduction in flow area when the gas flows into the particles. Oscillations of pressure, velocity and temperature near the granular interface are not present.

The effectiveness of increasing $\mathscr{D}$ to control oscillations of intergranular stress in dense granular regions is shown in figure 9(b). Severe oscillations of friction pressure are produced when $\mathscr{D}$ is zero, but increasing it to 1 or 10 effectively eliminates the oscillations at the cost of some numerical diffusion of granular interfaces.

5.5. High-pressure outgassing

The previous problem simulated high-pressure gas entering a region of dense particles. The next problem, on the other hand, simulates high-pressure gas leaving a region of dense particles. Consider high-pressure gas that has been injected into a granular layer as a shock passes over. This gas will be ejected when a subsequent expansion wave reduces the pressure at the surface. The outflow induced by the expansion wave will entrain some of the particles. The same initial conditions as the previous problem are used, with the exception that the particles are now located in the high-pressure region of the shock tube,

(5.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}p_{g}^{L}=100~\text{atm},\quad p_{g}^{R}=1~\text{atm},\\ T^{L}=300~\text{K},\quad T^{R}=300~\text{K},\\ Y_{g,Air}^{L}=1,\quad Y_{g,\mathit{Air}}^{R}=1,\\ {\it\alpha}_{s}^{L}=0.4,\quad {\it\alpha}_{s}^{R}=0.0,\\ {\it\Theta}_{s}^{L}=0,\quad {\it\Theta}_{s}^{R}=0.\end{array}\right\}\end{eqnarray}$$

The domain is 0.6 m long and the diaphragm was placed at 0.3 m. The solution at $400~{\rm\mu}\text{s}$ is shown in figure 10 for three grid resolutions.

A right-running shock wave and left-running rarefaction wave are produced, similar to the other shock-tube problems. Now, however, drag and buoyancy forces are induced by the escaping gas from the granular bed, and this ejects particles. This process reduces the particle volume fraction to a profile that is similar in form to the gaseous rarefaction wave. Intense mechanical and thermal losses to the gas phase by entraining particles severely weaken the transmitted gaseous shock. The spikes in temperature near the granular interface are physical and result from particles heating the gas at the gas-phase contact surface and a slight gap between the gas and granular contact surfaces (see figure 10 d). If the particles were not present, the temperature would be much higher in the post-shock state and much lower behind the gas-phase contact surface.

Figure 10.

Computed (a) gas-phase pressure, (b) velocity, (c) gas-phase temperature and (d) bulk density for high-pressure gaseous ejection from a dense layer of particles at a time of $400~{\rm\mu}\text{s}$ . The solid velocity is not shown if the particle volume fraction lies below $10^{-8}$ .

5.6. Interaction of a shock wave with a low-volume-fraction dust layer

This problem models the interaction of a shock wave with a loose dust layer. A schematic diagram of the geometrical set-up and initial and boundary conditions is shown in figure 11. The heights of the channel and dust layer are 6 cm and 2 cm, respectively. The left edge of the dust pile is placed at 2 cm, and the shock is initially located at 1 cm. A shock of Mach number 1.6 propagates into air at 1 atm and 288 K. The initial volume fraction of the dust layer is 0.04 %. The particle diameter, density and specific heat are $d_{s}=5~{\rm\mu}\text{m}$ , ${\it\rho}_{s}=1470~\text{kg}~\text{m}^{-3}$ and $C_{V,s}=987~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ , respectively. Figure 12 shows the computed mixture density at $900~{\rm\mu}\text{s}$ for grid spacings of 1.67, 0.83 and 0.42 mm at the finest AMR refinement level. The computed solutions, including the particle and gas-phase vortices, converge with increasing grid resolution.

Figure 11.

Schematic diagram of the initial conditions used for multidimensional test problems.

Figure 12.

Computed mixture density ( ${\it\rho}_{m}={\it\alpha}_{g}{\it\rho}_{g}+{\it\alpha}_{s}{\it\rho}_{s}$ , in units of $\text{kg}~\text{m}^{-3}$ ) for the low-volume-fraction shock wave and dust layer interaction problem at $900~{\rm\mu}\text{s}$ for three different grids: (a)  ${\rm\Delta}x_{min}=1.67$  mm, (b) 0.83 mm and (c) 0.42 mm.

Figure 12 shows that a vortex is introduced in the gas phase from a shear layer produced as particle drag locally decreases the fluid velocity in the granular layer. The gas-phase vortex entrains some of the particles. The particle vortex does not exactly coincide with the gas-phase vortex due to particle inertia. Mechanical and thermal losses from entraining particles cause the shock near the top of the layer to curve and weaken with time. Drag from high-pressure gas entering the layer from the $y$ -direction compresses the particles, which thins the dust layer immediately behind the shock. These results agree well with those of Fedorov & Fedorchenko (Reference Fedorov and Fedorchenko2010), in spite of the use of different drag and convective heat transfer models.

6.1. Shock and particle-curtain interaction

Here we test the model against nearly one-dimensional multiphase shock-tube experimental measurements in which a shock wave interacts with a curtain of particles. Two conditions are tested: a lower-volume-fraction particle curtain (Wagner et al. Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012) and a high-volume-fraction particle curtain (Rogue et al. Reference Rogue, Rodriguez, Haas and Saurel1998). A schematic diagram of the initial conditions and qualitative $x$ – $t$ diagram of the wave structures produced from a shock interacting with a particle curtain are shown in figure 13. A shock travelling at velocity $u_{sh}$ and strength $M_{sh}$ impacts a curtain of particles of length $L$ and initial volume fraction ${\it\alpha}_{s,0}$ . Drag forces resulting from gases that enter and flow through the particle bed and buoyancy forces ( ${\it\alpha}_{s}\boldsymbol{{\rm\nabla}}p_{g}$ ) from gas pushing on the granular interface accelerate the particles. Gas-phase energy that is not reflected or dissipated to the particles during the acceleration stage is transmitted to the downstream side of the particle bed in the form of a weak shock. The computational domain measured from $-100L$ to $100L$ is discretized with 4000 computational cells. Non-reflecting conditions were used on the left and right boundaries. The upstream edge of the particle curtain is placed at $x=0$ .

Figure 13.

Initial conditions and qualitative $x{-}t$ diagram for the shock wave interacting with a particle curtain.

6.1.1. Particle curtain with ${\it\alpha}_{s,0}=21\,\%$

Here we compare results of computations to experimental measurements made using the multiphase shock-tube apparatus described in Wagner et al. (Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012). Specifically, we simulate case 1 of Ling et al. (Reference Ling, Wagner, Beresh, Kearney and Balachandar2012), which presents additional data and compares the results of a Lagrangian multiphase model to these measurements. The length of the particle curtain is 2 mm. The initial volume fraction, diameter, density and specific heat of the soda lime particles is $21\,\%$ , $115~{\rm\mu}\text{m}$ , $2520~\text{kg}~\text{m}^{-3}$ and $840~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ , respectively. A Mach 1.66 shock propagates into air initially at 82.7 kPa and 296.4 K. Comparisons of the computed results and experimental data are shown in figure 14.

Figure 14.

Computed $x$ – $t$ diagram of (a) gas-phase density (in units of $\text{kg}~\text{m}^{-3}$ ) and (b) solid volume fraction for Mach 1.66 shock interacting with a particle curtain with $L=2$  mm. The particles were $115~{\rm\mu}\text{m}$ in diameter and had an initial volume fraction of 21 %. Experimental data presented in Ling et al. (Reference Ling, Wagner, Beresh, Kearney and Balachandar2012) and Wagner et al. (Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012) are represented by symbols.

The locations of the reflected shock, the transmitted shock and the upstream and downstream fronts of the particle curtain are in close agreement with experimental data. The largest discrepancy is the location of the upstream front of particles. The magnitude of this error depends on a precise definition of the particle-curtain edge, which is not stated in Ling et al. (Reference Ling, Wagner, Beresh, Kearney and Balachandar2012) or Wagner et al. (Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012). The computed position of the upstream edge of particles lags experimental results from $x/L=1.5$ to $2$ if a threshold of ${\it\alpha}_{s}=0.1\,\%$ or 0.01 % is used to define the edge of the particle curtain. On the other hand, the computed position of the upstream particle edge is much closer if a threshold value of ${\it\alpha}_{s}=1\,\%$ is used. The computed downstream edge of the particle curtain agrees well with measurements regardless of how the edge of the particle curtain is defined due to the locally high gradient of volume fraction.

6.1.2. Particle curtain with ${\it\alpha}_{s,0}=65\,\%$

The experiment of Rogue et al. (Reference Rogue, Rodriguez, Haas and Saurel1998) has a much larger volume fraction than the previous test ( ${\it\alpha}_{s,0}=65\,\%$ ). The high initial volume fraction is at a common definition of the random close-packing limit (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Igci et al. Reference Igci, Andrews, Sundaresan, Pannala and O’Brien2008). As a result, the parameters ${\it\alpha}_{s,crit}$ and ${\it\alpha}_{s,max}$ were raised to $65\,\%$ and $70\,\%$ , respectively, to keep the friction and solids pressures from becoming singular (see (3.21)–(3.23)). The length of the particle curtain is $2$  cm. The diameter, density and specific heat of the glass beads that were used to make the curtain are 1.5 mm, $2600~\text{kg}~\text{m}^{-3}$ and $840~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ , respectively. A Mach 1.3 shock propagates into air at 1 atm and 294.2 K. Comparisons of the computed results and experimental data are shown in figure 15.

Figure 15.

Computed $x$ – $t$ diagram of (a) gas-phase density (in units of $\text{kg}~\text{m}^{-3}$ ) and (b) solid volume fraction for Mach 1.3 shock interacting with a particle curtain with $L=2$  cm. The particles were 1.5 mm in diameter and had an initial volume fraction of 65 %. Experimental data of Rogue et al. (Reference Rogue, Rodriguez, Haas and Saurel1998) are represented by symbols.

The results of this test case are qualitatively similar to the previous case. One notable exception is that the transmitted shock is much weaker. This is primarily due to the much higher mechanical and thermal losses created as the shock propagates through the highly packed glass beads. The calculated results match closely with the experimental measurement of the downstream particle front. The upstream front of the particles does not agree as well. This is due, in part, to the ambiguity of defining the edge of particles in the computation, which is undefined in Rogue et al. (Reference Rogue, Rodriguez, Haas and Saurel1998). The upstream front of particles is smoother than the downstream edge, and, as a result, the modelling error is highly sensitive to how the particle edge is defined. If a threshold volume fraction of 1 % is used (see figure 15 b), the error is quite large and the front barely appears to move. The agreement between the measurements and the computed results is much better if a threshold volume fraction between 7 % and 10 % is used to define the curtain edge. This ambiguity could be avoided if the particle-curtain edges were precisely defined (based, perhaps, on a threshold volume fraction) in experimental measurements, but this is difficult due to, among other issues, uncertainties of extinction coefficients through particle clouds and details of the experimental apparatus.

Figure 16.

Comparison between computed and experimental (Rogue et al. Reference Rogue, Rodriguez, Haas and Saurel1998) pressure–time traces for a Mach 1.3 shock interacting with a dense particle curtain.

Rogue et al. (Reference Rogue, Rodriguez, Haas and Saurel1998) also measured the pressure–time ( $p$ – $t$ ) traces at physical distances of $-11$  cm ( $x/L=-5.5$ ) and 4.3 cm ( $x/L=2.15$ ) with respect to the upstream edge of the particle curtain. A comparison of the computed to the measured pressure–time traces at these locations is shown in figure 16. The calculations closely match the measurements, and also match the timing and strength of the reflected shock, as well as the subsequent pressure decay as the particle curtain (which acts like a porous moving piston) accelerates. The model results also agree with the strength and timing of the transmitted shock.

The computed $p$ – $t$ traces do not match as well when the particle cloud arrives at $x/L=2.15$ and the measured pressure rises rapidly. The sharp pressure rise is produced as the particle bed moves past the measurement location. The $p$ – $t$ trace is essentially a record of the pressure drop as a function of position in the particle curtain as it moves past the measurement location, which is highly sensitive to the volume fraction and drag, and, to some extent, to the models selected for solids and frictional–collisional pressure. (A detailed comparison of the effect of drag and frictional–collisional pressure models on the structure of the particle curtain is beyond the scope of this paper and can be undertaken as a separate study.)

6.2. Transmitted granular oblique shock and contact angles

In the experiment of Fan et al. (Reference Fan, Chen, Jiang and Li2007), a shock interacts with a dense layer of dust. The angle that a transmitted compaction wave and granular contact surface make with the horizontal when a shock wave passes over the dust layer, as shown in figure 17(a), was measured. An oblique granular shock is transmitted into the dense layer of dust by forces from high-pressure gas behind the shock pushing on the granular interface as well drag from gases penetrating into the layer. These forces compact the layer and induce an oblique granular contact surface in addition to a transmitted granular shock.

The initial conditions follow the diagram in figure 11. For this test, the layer thickness and domain heights are 2 and 10 cm, respectively. The initial position of the shock is 1 cm, and the left edge of the dust layer is at 0 m. The initial volume fraction was 47 %, and the density, diameter, coefficient of restitution and specific heat of the particles were ${\it\rho}_{s}=1100~\text{kg}~\text{m}^{-3}$ , $d_{s}=15~{\rm\mu}\text{m}$ , $e=0.9$ and $C_{V,s}=987~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ , respectively. The velocity of the shock wave, which propagated into air at 1 atm and 300 K, is $990~\text{m}~\text{s}^{-1}$ . The bottom boundary was changed to a non-reflecting condition to avoid reflections and increase the measurement distance for the transmitted compaction wave and granular contact angles. The equations were integrated for $100~{\rm\mu}\text{s}$ with ${\rm\Delta}x=347~{\rm\mu}\text{m}$ at the smallest level of AMR refinement.

The computed solid volume fraction at $100~{\rm\mu}\text{s}$ is shown in figure 17(b). The computed transmission angles of the granular contact surface and compaction wave are $1.18^{\circ }$ and $4.78^{\circ }$ , respectively, which are close to the measured values of $1^{\circ }$ – $2^{\circ }$ for the contact surface and $4^{\circ }$ for the granular shock (Fan et al. Reference Fan, Chen, Jiang and Li2007).

Figure 17.

Angles of the transmitted compaction wave and granular contact from a shock wave propagating across the dust layer: (a) definition of the angles and (b) computed solution at $100~{\rm\mu}\text{s}$ and a grid spacing of $347~{\rm\mu}\text{m}$ at the finest level of refinement.