Turbulence transport in moderately dense gas–particle compressible flows

Abstract

This study employs three-dimensional particle-resolved simulations of planar shocks passing through a suspension of stationary solid particles to study wake-induced gas-phase velocity fluctuations, termed pseudo-turbulence. Strong coupling through interphase momentum and energy exchange generates unsteady wakes and shocklets in the interstitial space between particles. A Helmholtz decomposition of the velocity field shows that the majority of pseudo-turbulence is contained in the solenoidal component from particle wakes, whereas the dilatational component corresponds to the downstream edge of the particle curtain where the flow chokes. One-dimensional phase-averaged statistics of pseudo-turbulent kinetic energy (PTKE) are quantified at various stages of flow development. Reduction in PTKE is observed with increasing shock Mach number due to decreased production, consistent with single-phase compressible turbulence. The anisotropy in Reynolds stresses is found to be relatively constant through the curtain and consistent over all the conditions simulated. Analysis of the budget of PTKE shows that the majority of turbulence is produced through drag and balanced by viscous dissipation. The energy spectra of the streamwise gas-phase velocity fluctuations reveal an inertial subrange that begins at the mean interparticle spacing and decays with a power law of $-5/3$ and steepens to $-3$ at scales much smaller than the particle diameter. A two-equation model is proposed for PTKE and its dissipation. The model is implemented within a hyperbolic Eulerian-based two-fluid model and shows excellent agreement with the particle-resolved simulations.

1. Introduction

High-speed flows through particulate media occur in diverse applications, such as detonation blasts (Zhang et al. 2001), volcanic eruptions (Chojnicki, Clarke & Phillips 2006; Lube et al. 2020), coal-dust explosions (Sapko et al. 2000; Zheng et al. 2009), pulsed-detonation engines (Chang & Kailasanath 2003; Roy et al. 2004) and plume–surface interactions during interplanetary landings (Plemmons et al. 2009; Morris et al. 2011; Capecelatro 2022). In these examples, turbulence plays a crucial role in governing processes like reactant mixing and particle dispersion. However, the nature of this turbulence is distinct from both single-phase compressible turbulence and low-speed multiphase turbulence, posing a challenge to the accuracy of existing models.

Compressibility effects in turbulent flows are often characterized using the turbulent Mach number (Sagaut & Cambon 2008; Jagannathan & Donzis 2016). For values of $M_t\leqslant 0.3$ , large-scale separation exists between acoustics and turbulence. This results in a nearly incompressible flow called the quasi-isentropic regime. For higher values of $M_t$ (i.e. $0.3\lt M_t\leqslant 0.6$ ), dilatational effects are significant, leading to a nonlinear subsonic regime. The flows considered in the present study predominantly fall within this regime.

Since the 1970s, numerous studies have investigated the role of compressibility in the development of turbulent mixing layers and the generation of turbulent kinetic energy (Brown & Roshko 1974; Bradshaw 1977; Sarkar et al. 1991). Early work by Zeman (1990) and Sarkar et al. (1991) examined the effects of dilatational dissipation, $\epsilon _d$ , finding that its increase with $M_t$ leads to a reduction in turbulent kinetic energy, thereby decreasing turbulent mixing. They suggested that the suppression of growth rate is linked to increased $\epsilon _d$ due to shocklets. They developed a mathematical model to incorporate this effect into Reynolds stress closure models. However, Sarkar (1995) later showed, using direct numerical simulations of turbulent homogeneous shear flow, that the reduction of turbulent kinetic energy is primarily due to decreased turbulence production, rather than directly caused by dilatational dissipation. Subsequent studies by Vreman, Sandham & Luo (1996) and Pantano & Sarkar (2002) corroborated this finding, showing that dilatational dissipation is negligible. Instead, the reduced growth rate of turbulence is linked to diminished pressure fluctuations and, consequently, lower turbulence production resulting from a reduction in the pressure-strain term.

Kida & Orszag (1990) were among the first to analyse the kinetic energy spectrum in forced compressible turbulence, observing that its scaling is largely independent of Mach number. Donzis & Jagannathan (2013) also found that the turbulent kinetic energy spectrum in compressible isotropic turbulence follows a $-5/3$ power law in the inertial range for $0.1 \leqslant M_t \leqslant 0.6$ , consistent with the classical Kolmogorov scaling for incompressible flows (Kolmogorov, 1941b ). Further insights into compressibility scaling emerge from a Helmholtz decomposition of the velocity field $\boldsymbol{u}$ into its solenoidal component $\boldsymbol{u}_s$ and dilatational component $\boldsymbol{u}_d$ (Kida & Orszag 1990; Donzis & Jagannathan 2013; Wang et al. 2011, 2012; San & Maulik 2018). Compressibility effects are typically attributed to $\boldsymbol{u}_d$ , and both Wang et al. (2011) and Donzis & Jagannathan (2013) observed that the majority of turbulent kinetic energy resides in the solenoidal component, with $\boldsymbol{u}_d$ increasing with $M_t$ . However, all of these studies have focused on single-phase compressible turbulent flows in the absence of particles.

In multiphase flows, interphase coupling introduces additional complexity that significantly influences energy transfer and turbulence characteristics. Fluid velocity fluctuations induced by particle wakes are referred to as pseudo-turbulence (Lance & Bataille 1991; Mehrabadi et al. 2015), a term also applied to bubble-induced turbulence (BIT) in liquid flows (Risso 2018). Lance & Bataille (1991) first demonstrated that a homogeneous swarm of bubbles generates pseudo-turbulence with a spectral subrange exhibiting a $-3$ power law. They showed that at a statistically steady state, this spectral scaling results from a balance between viscous dissipation and energy production due to drag forces from rising bubbles. Similar scaling has since been observed in other bubbly flows (Mercado et al. 2010; Risso 2018; Mezui et al. 2022, 2023). Subsequent experimental studies coupling BIT with shear-induced turbulence have found that the spectra of liquid velocity fluctuations follow a $-3$ scaling at small wavenumbers, transitioning to a $-5/3$ scaling at higher wavenumbers, suggesting a single-phase signature is preserved at the smallest scales (Risso 2018). Numerical simulations of gas–particle turbulent channel flow reveal that two-way coupling between the phases results in reduction in fluid-phase turbulent kinetic energy at the scale of individual particles, while a broadband reduction over all scales is observed at moderate-to-high mass loading (Capecelatro, Desjardins & Fox 2018).

Over the past few decades, turbulence models have evolved to incorporate the effects of particles (Troshko & Hassan 2001; Fox 2014; Ma et al. 2017). A production term must be included to account for generation of turbulence through drag. A dissipation time scale is often employed based on the relative velocity between the phases ( $u_r$ ) and particle diameter ( $d_p$ ), given by $\tau = d_p / u_r$ . The use of two-equation transport models for gas–solid flows dates back to the work of Elghobashi & Abou-Arab (1983), who derived a rigorous set of equations for dilute concentrations of particles in incompressible flow using a two-fluid approach. Since then, models have been proposed for denser regimes in shear turbulence (Ma & Ahmadi 1990). Crowe, Troutt & Chung (1996) provided a review of numerical models for turbulent kinetic energy in two-phase flows. However, these models are limited to intrinsic turbulence whereby the carrier-phase turbulence would exist even in the absence of particles, as opposed to pseudo-turbulence that is entirely generated by the particle phase. Mehrabadi et al. (2015) recently developed an algebraic model for pseudo-turbulent kinetic energy (PTKE) based on particle-resolved simulation data that depends on the slip Reynolds number and particle volume fraction. A limitation of algebraic models is that PTKE can only be predicted in regions of finite volume fraction. In cases where turbulence is generated within a suspension of particles and advects downstream into the surrounding gas, transport equations for PTKE are more appropriate (Shallcross, Fox & Capecelatro 2020).

Particle-laden compressible flows challenge numerical models due to the strong coupling between shock waves, particles and turbulence over a wide range of scales. Using particle-resolved simulations of compressible homogeneous flows past random arrays of particles, Khalloufi & Capecelatro (2023) found that both $M_t$ and PTKE increase with particle volume fraction for a fixed free stream Mach number. Experimental and numerical studies of particle-laden underexpanded jets have demonstrated significant modification of shock structures due to the two-way coupling between the gas and particles even at low volume fractions where one-way coupling would be deemed appropriate for single-phase flow (Sommerfeld 1994; Patel et al. 2024). Two-dimensional particle-resolved simulations of shock–particle curtain interactions revealed PTKE magnitudes comparable to the resolved kinetic energy (Regele et al. 2014; Hosseinzadeh-Nik, Subramaniam & Regele 2018). In three-dimensional inviscid simulations, Mehta, Jackson & Balachandar (2020) reported velocity fluctuations reaching up to 50 % of the kinetic energy based on the mean flow, with increasing velocity fluctuations observed at higher shock Mach numbers, $M_s$ , and particle volume fractions, $\varPhi _p$ . It should be noted that shock-driven multiphase flows in radial configurations are prone to instabilities not captured in the planar geometries considered in the aforementioned studies (McFarland et al. 2016; Middlebrooks et al. 2018). However, because turbulence transport occurs on much shorter time scales than particle dispersal or instability growth, such effects are not relevant to the present focus on gas-phase turbulence.

Models for PTKE in compressible gas–particle flows have only recently begun to emerge. Osnes et al. (2019) proposed an algebraic model for PTKE based on particle-resolved simulations of shock–particle interactions that depends on the mean flow speed and particle volume fraction. Shallcross et al. (2020) proposed a one-equation model for PTKE containing a production term due to drag and an algebraic closure for dissipation. The dissipation model employs a time scale based on the particle diameter and local slip velocity – consistent with that used in BIT models (Ma et al. 2017) – augmented with a blending function to account for regions devoid of particles. However, the results were found to be highly sensitive to the closure applied to dissipation, limiting its applicability.

Building on these efforts, the present study aims to advance understanding of turbulence transport in compressible gas–particle flows, particularly at moderate volume fractions and Mach numbers. The paper is organized as follows. In § 2, the simulation configuration and governing equations are presented. Simulation results are provided in § 3, starting with a qualitative assessment of the flow, followed by one-dimensional phase-averaged statistics of the gas-phase velocity. The budget of PTKE is presented next, revealing key production and dissipation mechanisms. The energy spectra within the particle curtain are then presented and separate contributions from solenoidal and dilatational components highlight the sources of PTKE. In § 4, a two-equation turbulence model for PTKE and its dissipation is proposed and implemented within a hyperbolic two-fluid model. Results from § 3 are used to guide closure. An a posteriori analysis is performed and first- and second-order statistics are compared. Key findings and results are summarized in § 5.

2. Simulation set-up and methods

2.1. Flow configuration

To isolate shock–particle–turbulence interactions, we perform three-dimensional, particle-resolved simulations of a planar shock propagating through a suspension of stationary, monodisperse particles. The assumption of frozen particles is justified, as the acoustic time scale is several orders of magnitude shorter than the particle response time for the high density ratios ( $\rho _p/\rho \gt 10^3$ ) typical of gas–solid flows (Ling, Haselbacher & Balachandar 2011). The simulations are designed to emulate the multiphase shock-tube experiments of Wagner et al. (2012). Figure 1 shows a volume rendering of the gas-phase velocity magnitude within the simulation domain at a moment when the shock has advanced significantly beyond the curtain and exited the domain. The velocity increases across the particle curtain with maximum values at the downstream curtain edge where the flow chokes due to the sudden change in volume fraction.

Figure 1. The simulation domain showing particle position and a volume rendering of the gas-phase velocity magnitude after the shock has passed the curtain ( $t/\tau _L = 2$ ) with $\varPhi _p=0.2$ and $M_s=1.66$ .

Particles with diameter $D=115$ $\unicode {x03BC}$ m and density $\rho _p=2520 \ \rm {kg}\,\rm {m}^{-3}$ are randomly distributed within a curtain of thickness $L=2$ mm ( $L=17.4D$ ). A minimum of two grid points is maintained between particle surfaces. A planar shock is initially placed at a non-dimensional length of $x=5.5D$ with the flow direction parallel to the $x$ -axis. The upstream edge of the curtain is placed at $x=7D$ . Periodic boundary conditions are imposed in the two spanwise ( $y$ and $z$ ) directions. The domain size for all but one case is $[L_x\times L_y\times L_z]=[30\times 12\times 12] D$ . Here $L_y$ and $L_z$ were chosen based on a domain size independence study summarized in Appendix A. The domain is discretized with uniform grid spacing $\Delta x=D/40$ , corresponding to $[1201\times 480\times 480]$ grid points.

The pre-shock gas-phase density is $\rho _\infty =0.987 \ \rm {kg}\,\rm {m}^{-3}$ , pressure $P_\infty =82.7\ \textrm {kPa}$ , sound speed $c_\infty =343 \ \rm {m}\,\rm {s}^{-1}$ and velocity $u_\infty =0 \ \rm {m}\,\rm {s}^{-1}$ . Post-shock properties, denoted by the subscript ps, are obtained via the Rankine–Hugoniot conditions. The shock Mach number is defined as ${M_s} = u_s/c_\infty$ , where $u_s$ is the shock speed. A reference time scale based on the distance (in terms of particle curtain length) that the shock travels is defined as $\tau _L=L/u_s$ . The particle Reynolds number based on post-shock properties is defined as $Re_{ps} = \rho _{ps}u_{ps}D / \mu _{ps}$ , where $\mu _{ps}$ is the gas-phase viscosity at temperature $T_{ps}$ . The number of particles $N_p$ within the curtain is determined from the average volume fraction, $\varPhi _p$ . A summary of the cases considered in this study is given in table 1. Cases $1{-}9$ represent different combinations of $M_s$ and $\varPhi _p$ . Case $10$ exhibits a longer domain length to study turbulence transport downstream of the particle curtain.

Table 1. Parameters for the various runs used in this study.

2.2. Governing equations

The gas-phase is governed by the viscous compressible Navier–Stokes, given by

(2.1) \begin{equation} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}) = 0, \end{equation}

(2.2) \begin{equation} \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u} \otimes \boldsymbol{u} + p \mathbb{I} - \boldsymbol{\sigma })=0 \end{equation}

and

(2.3) \begin{equation} \frac {\partial \rho E}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\{ \rho E + p \} \boldsymbol{u} + \boldsymbol{q} - \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\sigma }) = 0, \end{equation}

where $\rho$ is the gas-phase density, $\boldsymbol{u}=(u,v,w)$ is the velocity and $E$ is the total energy. The viscous stress tensor is

(2.4) \begin{equation} \boldsymbol{\sigma } = \mu (\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^T) + \mu ^{\prime} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} \end{equation}

and the heat flux is

(2.5) \begin{equation} \boldsymbol{q} = - k \boldsymbol{\nabla }T \end{equation}

where $k$ is the thermal conductivity. The dynamic viscosity is modelled as a power law, $\mu = \mu _0 [(\gamma -1)T/T_0]^n$ , where $\gamma =1.4$ is the ratio of specific heats and $n=0.666$ . The second coefficient of viscosity is $\mu ^{\prime} = \mu _B - 2/3 \mu$ where the bulk viscosity $\mu _B=0.6 \mu$ is chosen as a model for air. The thermal conductivity is varied with a similar power law as viscosity to maintain a constant Prandtl number of $0.7$ . Thermodynamic relations for temperature and pressure are given by

(2.6) \begin{equation} T = \frac {\gamma p}{(\gamma -1)\rho } \quad \text{and}\quad p= (\gamma -1)\bigg(\rho E - \frac {1}{2} \rho \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{u}\bigg). \end{equation}

No-slip, adiabatic boundary conditions are enforced at the particle surfaces. Details on the numerical implementation are provided in the following section.

2.3. Numerics

The simulations are performed using the compressible multiphase flow solver jCode (Capecelatro 2023). Spatial derivatives are approximated using narrow-stencil finite-difference operators that satisfy the summation by parts (SBP) property (Strand 1994; Svärd et al. 2007). A sixth-order centred finite-difference scheme is used for the interior points, and a fourth-order, one-sided finite difference is applied at the boundaries. Kinetic energy preservation is achieved using a skew-symmetric-type splitting of the inviscid fluxes (Pirozzoli 2011), providing nonlinear stability at low Mach number. To ensure proven temporal stability, the SBP scheme is combined with the simultaneous approximation-term treatment that weakly enforces characteristic boundary conditions at the inflow and outflow (Svärd et al. 2007). High-order SBP dissipation operators (Mattsson, Svärd & Nordström 2004) are employed to dampen spurious high-wavenumber modes. Localized artificial diffusivity is used as a means of shock capturing by following the ‘LAD-D2-0’ formulation in Kawai, Shankar & Lele (2010). To limit the artificial diffusivity to regions of high compression (shocks), we employ the sensor originally proposed by Ducros et al. (1999) and later improved by Hendrickson, Kartha & Candler (2018). More details can be found in Khalloufi & Capecelatro (2023). The equations are advanced in time using a standard fourth-order Runge–Kutta scheme, with a constant Courant–Friedrichs–Lewy number of 0.8.

A ghost-point immersed boundary method originally proposed by Mohd-Yusof (1997) and later extended to compressible flows by Chaudhuri, Hadjadj & Chinnayya (2011) is employed to enforce boundary conditions at the surface of the particles. Values of the conserved variables at ghost points residing within the solid are assigned after each Runge–Kutta subiteration to enforce no-slip, adiabatic boundary conditions. The framework was validated in our previous study (Khalloufi & Capecelatro 2023), demonstrating that 40 grid points across the particle diameter is sufficient to capture drag and PTKE. An assessment of the domain size and sensitivity to random particle placement is reported in Appendix A.

3. Results

3.1. Flow visualization

Instantaneous snapshots of the flow field (see supplementary Movie S1) corresponding to Case 5 ( ${M}_s=1.66, \varPhi _p=0.2$ ) at $t/\tau _L=0.5$ , $1$ and $2$ are presented in figure 2. A two-dimensional slice in the $x{-}y$ plane shows the local gas-phase Mach number and numerical schlieren in the vicinity of the particles. The incident shock travels in the positive $x$ direction and impinges the particle curtain located at $x_0$ at time $t=0$ . Upon impact, the shock splits into a weaker transmitted shock that penetrates the curtain, as shown in figure 2 $(a)$ . At the upstream edge of the curtain, the arrival of the shock generates multiple shocklets at the surface of each particle, which coalesce into a reflected shock wave. Shock–particle interactions are seen to generate significant fluctuations in the gas-phase velocity. Contour lines of $M=1$ (shown in purple) demarcate local supersonic regions. In figure 2 $(b)$ , the shock has nearly reached the downstream curtain edge, and the local supersonic regions move downstream with the flow. Figure 2 $(c)$ shows that the flow has stabilized with both the transmitted and reflected shocks having exited the domain boundaries. The particles restrict the area of the transmitted shock, causing the gas phase to choke near the downstream edge of the curtain due to the abrupt change in volume fraction, followed by a supersonic expansion. Velocity fluctuations induced by the particles advect downstream from the curtain and dissipate, akin to grid-generated turbulence.

Figure 2. Two-dimensional planes at the centre ( $z/D=0$ ) showing schlieren (top-half of each panel) and local Mach number $M=\lVert \boldsymbol{u}\rVert /c$ (bottom-half of each panel) at (a) $t/\tau _L=0.5$ , (b) $t/\tau _L=1$ and (c) $t/\tau _L=2$ for ${M}_s=1.66$ and $\varPhi _p=0.2$ . Contour lines of $M=1$ shown in purple. Blue circles depict particle cross-sections.

3.2. Averaging operations

The flows under consideration are unsteady and statistically homogeneous in the two spanwise directions. Averaged quantities depend solely on one spatial dimension ( $x$ ) and time. Due to the presence of particles and gas-phase density variations, special attention must be given to the averaging process. To facilitate statistical phase-averaging, an indicator function is defined as

(3.1) \begin{equation} \mathcal{I}(\boldsymbol{x}) = \begin{cases} 1 & \text{if } \boldsymbol{x} \in \text{gas phase},\\ 0 & \text{if } \boldsymbol{x} \in \text{particle}. \end{cases} \end{equation}

Spatial averages are taken as integrals over $y{-}z$ slices. The integration of the indicator function yields a volume fraction $\alpha$ (or area fraction in this case) that depends solely on $x$ (time is omitted since the particles being stationary), given by

(3.2) \begin{equation} \alpha _g (x) =\langle \mathcal{I}\rangle \equiv \frac {1}{L_y L_z}\int _{L_y}\int _{L_z}\mathcal{I} \ \textrm {d}y\,\textrm {d}z, \end{equation}

where angled brackets denote a spatial average. Two other important averaging operations that will be used throughout this study are phase averages and density-weighted (Favre) averages. If $\psi (\boldsymbol{x},t)$ represents a random field variable, these averages are defined as

(3.3) \begin{equation} \begin{aligned} &\text{Spatial-average: } \langle \psi \rangle (x,t) \equiv \frac {1}{L_y L_z}\int _{L_y}\int _{L_z} \psi \ \textrm {d}y\,\textrm {d}z, \\ &\text{Phase-average: } \overline {\psi }(x,t) \equiv \frac {\langle \mathcal{I}\psi \rangle }{\langle \mathcal{I}\rangle }\equiv \frac {\langle \mathcal{I}\psi \rangle }{\alpha _g}, \\ &\text{Favre-average: } \widetilde {\psi }(x,t) \equiv \frac {\langle \mathcal{I} \rho \psi \rangle }{\langle \mathcal{I}\rho \rangle } \equiv \frac {\overline {\rho \psi }}{\overline {\rho }}. \end{aligned} \end{equation}

Spatial averages and phase averages are related via $\langle {\mathcal{I}\psi } \rangle = \alpha _g \overline {\psi }$ and similarly $\overline {\rho \psi }=\overline {\rho }\widetilde {\psi }$ . A field variable can be decomposed into its phase-average and a fluctuating quantity as $\psi = \overline {\psi } + \psi ^{\prime}$ . Similarly, the Favre decomposition is $\psi = \widetilde {\psi } + \psi ^{\prime\prime}$ .

3.3. Mean velocity, fluctuations and anisotropy

The Favre-averaged gas-phase velocity, $\widetilde {u}$ , as a function of the streamwise direction at three different time instances ( $t/\tau _L = 0.5$ , $1$ and $2$ ) is shown in figure 3. The abrupt drop in velocity observed at early times ( $t/\tau _L = 0.5$ and $1$ ) marks the location of the transmitted shock. The flow decelerates significantly as it approaches the particle curtain due to drag, with greater reduction in velocity relative to the post-shock velocity at higher volume fractions. The flow then accelerates as it traverses the curtain. At the latest time ( $t/\tau _L = 2$ ), a sharp increase in $\widetilde {u}$ at the downstream edge of the curtain is seen across all cases, indicating a region of choked flow transitioning to supersonic velocities. Similar trends in the velocity field have been reported previously (e.g. Mehta et al. 2018; Theofanous, Mitkin & Chang 2018; Osnes et al. 2019). Mehta et al. (2018) obtained an analytical solution of the Riemann problem for a duct with a sudden change in cross-sectional area as a simpler means of predicting the flow through a particle curtain. The solution was found to compare well with inviscid simulations of shock–particle interactions, though it is unable to predict the choking behaviour leading to supersonic velocities observed here.

Figure 3. Mean gas-phase velocity profiles. Darker lines indicate higher volume fractions: $\varPhi _p=0.1$ (light pink), $\varPhi _p=0.2$ (pink), $\varPhi _p=0.3$ (purple). Here (a,d,g) $t/\tau _L=0.5$ , (b,e,h) $t/\tau _L=1$ and (c,f,i) $t/\tau _L=2$ ; $(a)$ – $(c)$ $M_s=1.2$ , $(d)$ – $(f)$ $M_s=1.66$ , $(g)$ – $(i)$ $M_s=2.1$ . The grey-shaded region indicates the location of the particle curtain.

The amplitude and speed of the reflected shocks, indicated by the abrupt increase in velocity upstream of the particle curtain, increase with $\varPhi _p$ . The transmitted shock travels faster through the curtain at lower $\varPhi _p$ where the flow is less obstructed. For a given volume fraction, the magnitude of $\widetilde {u}$ decreases with increasing $M_s$ , and the flow-expansion region at the downstream edge rises sharply with increasing $M_s$ .

The root-mean-square (r.m.s.) gas-phase velocity fluctuations in the streamwise direction is defined as $u_{{rms}}^2 = \widetilde {u^{\prime\prime}u^{\prime\prime}}$ . Due to symmetry, the spanwise fluctuations are taken as $v_{{rms}}^2 = (\widetilde {v^{\prime\prime}v^{\prime\prime}\,}+\widetilde {w^{\prime\prime}w^{\prime\prime}})/2$ . Figure 4 shows these components at $t/\tau _L=1$ and $2$ . All values are normalized by the post-shock kinetic energy, $u_{ps}^2$ . Velocity fluctuations originate almost immediately within the particle curtain. The magnitude of the streamwise fluctuations is nearly twice the spanwise components. The fluctuations are higher at initial times, shortly after the shock passes over the particles. The maximum velocity fluctuations occur at the downstream edge where the flow chokes. Overall, the fluctuations decrease in magnitude with increasing $M_s$ . This reduction can be attributed to an increase in compressibility effects with higher $M_s$ . The precise dissipation mechanisms will be quantified in § 3.4, where individual terms of the PTKE budget are reported.

Figure 4. Velocity fluctuations at (a,c,e) $t/\tau _L=1$ and(b,d,f) $t/\tau _L=2$ . Here $(a,b)$ $M_s=1.2$ , $(c,d)$ $M_s=1.66$ and $(e,f)$ $M_s=2.1$ . Streamwise fluctuations $u_{{rms}}^2$ (pink/purple), spanwise fluctuations $v_{{rms}}^2$ (shades of green). Here $\varPhi _p=0.1$ (light shade), $\varPhi _p=0.1$ (intermediate shade), $\varPhi _p=0.3$ (dark shade).

It is interesting to note that the normalized fluctuations are nearly invariant with volume fraction except for the lowest shock Mach number case at early times (see figure 4 $a$ ). Previous studies by Mehta et al. (2020) observed an increase in velocity fluctuations with $\varPhi _p$ . However, we only observe significant variation due to $\varPhi _p$ at the downstream edge of the curtain.

To better quantify the level of anisotropy, we define the gas-phase anisotropy tensor as

(3.4) \begin{equation} b_{ij} = \frac {R_{ij}}{2k_g} - \frac {\delta _{ij}}{3}, \end{equation}

where $R_{ij}=\widetilde {u_i^{\prime\prime}u_j^{\prime\prime}}$ is the pseudo-turbulent Reynolds stress, $k_g=\widetilde {u_i^{\prime\prime}u_i^{\prime\prime}}/2$ (repeated indices imply summation) is the PTKE and $\delta _{ij}$ is the Dirac delta function. The streamwise component $b_{11}$ is dominant compared with the components perpendicular to the flow direction $b_{22}$ and $b_{33}$ . The cross-correlation of velocity fluctuations, $b_{12}$ , is often negligible in gas–solid flows (Mehrabadi et al. 2015). Due to symmetry in the flow, only $b_{11}$ and $b_{22}$ are reported.

All nine cases are overlaid in figure 5 at $t/\tau _L=2$ with each line style representing a distinct volume fraction and each shade of colour representing a distinct shock Mach number. The values for anisotropy obtained here are in close agreement with those reported by Osnes et al. (2019). Interestingly, the level of anisotropy remains approximately constant across the curtain, with $b_{11}\approx 0.2$ and $b_{22}\approx -0.1$ for all cases, regardless of $\varPhi _p$ and $M_s$ . Variations are confined to the edges of the curtain, where the streamwise component becomes increasingly dominant. Upstream of the curtain, the anisotropy tensor becomes ill-defined as $k_g$ vanishes in the absence of turbulence production. At the curtain boundaries, peaks in the anisotropy tensor emerge, resulting from localized turbulence generation driven by mean gradients of streamwise velocity fluctuations and particle-induced drag. As turbulence propagates downstream through the curtain, velocity fluctuations are redistributed into the spanwise directions through the pressure–dilatation term (see § 3.4). At the downstream edge, flow expansion triggers renewed turbulence production in the streamwise direction, further amplifying anisotropy. Overall, the level of anisotropy within the curtain indicates that approximately 50 % of PTKE is contained in $\widetilde {u^{\prime\prime}u^{\prime\prime}}$ , with the remaining 50 % equally partitioned between $\widetilde {v^{\prime\prime}v^{\prime\prime}\,}$ and $\widetilde {w^{\prime\prime}w^{\prime\prime}}$ .

Figure 5. Components of the Reynolds stress anisotropy tensor for Cases 1–9 at $t/\tau _L=2$ . Here $\varPhi _p=0.1$ ( ), $\varPhi _p=0.2$ ( ) and $\varPhi _p=0.3$ ( ). Parallel component $b_{11}$ (light pink, pink, purple) and perpendicular component $b_{22}$ (light green, green, dark green) for $M_s=1.2$ , $1.66$ and $2.1$ (light to dark).

The PTKE advects and decays downstream of the particle curtain. Case $10$ extends $3L$ downstream of the particle curtain to examine this behaviour in greater detail. Figure 6 shows the r.m.s. velocity components at $t/\tau _L=5$ , where the flow reaches a steady state. It can be seen that the flow remains anisotropic beyond the curtain and eventually the fluctuations completely decay. This is analogous to grid-generated turbulence (Batchelor & Townsend 1948; Mohamed & Larue 1990; Kurian & Fransson 2009). According to Batchelor & Townsend (1948), the decay of turbulence intensity downstream of a grid (or screen) with mesh width $\Delta M$ follows a power law, given by

(3.5) \begin{equation} \Bigg (\frac {u_{{rms}}}{u_0}\Bigg )^2 = A \Bigg ( \frac {x-x_L}{\Delta M} \Bigg )^{n}, \end{equation}

where $u_0$ is the velocity of the gas phase at a point of virtual origin of turbulence $x_0$ and $A$ is an empirical constant. An analogy can be drawn to our shock–particle configuration by setting the mesh width to the average interparticle spacing, $\lambda$ , which can be defined as

(3.6) \begin{equation} \lambda = D \left (\frac {\pi }{6 \varPhi _p} \right )^{1/3}. \end{equation}

Additionally, we set the point of origin of turbulence decay to the location of the downstream curtain edge, $x_L=x_0+L$ , and consider the velocity at this point $u_L=\tilde {u}(x=x_L;t=5\tau _L)$ when normalizing the turbulence intensity.

Figure 6. $(a)$ Components of r.m.s. velocities $u_{{rms}}^2$ (pink) and $v_{{rms}}^2$ (green) for Case 10 downstream of the particle curtain at $t/\tau _L=5$ . The inset shows the components in log-scale. $(b)$ Ratio of the spanwise to streamwise r.m.s. velocities as a function of streamwise distance normalized by interparticle spacing $\lambda$ .

Figure 6 $(a)$ shows the decay of streamwise and spanwise components of r.m.s. velocities as a function of the downstream distance normalized by $\lambda$ . The inset illustrates this decay in log scale, from which we conclude that the decay does indeed follow a power-law behaviour with an exponent $n=-1.7$ . This value is slightly higher than the reported values for $n$ in incompressible, single-phase grid-generated turbulence reported in the literature, which range from $-1.13$ to $-1.6$ (Mohamed & Larue 1990; Kurian & Fransson 2009). The ratio $v_{{rms}}/u_{{rms}}$ shown in figure 6 $(b)$ highlights significant anisotropy of approximately $0.7$ , while downstream it reduces to $\approx 0.5$ suggesting that the flow remains anisotropic even at later time periods.

Figure 7. Turbulent Mach number for Cases $1-9$ at $t/\tau _L=2$ . Same legend as $b_{11}$ in figure 5.

The turbulent Mach number, defined as $M_t = \sqrt {2k_g}/\overline {c}$ , is shown at $t/\tau _L=2$ for Cases $1{-}9$ in figure 7. Here $M_t$ tends to increase rapidly at the upstream edge of the curtain where turbulence is first generated, then gradually increases throughout the curtain and peaks at the downstream edge where the flow chokes. The turbulent Mach number increases monotonically with the incident shock speed. Within the curtain, $M_t$ is relatively independent of $\varPhi _p$ , but increases with increasing $\varPhi _p$ at the downstream edge. For the cases with the lowest shock Mach number ( $M_s=1.2$ ), $M_t\approx 0.2$ , which falls in the quasi-isentropic regime, as classified by Sagaut & Cambon (2008), where pressure fluctuations are not significant. These cases are distinct from the higher $M_s$ cases in that the velocity does not rapidly increase at the downstream edge of the curtain (see figure 3 $c$ ) and the mean sound speed remains relatively constant (not shown) and thus the trends in $M_t$ are qualitatively different from the two higher $M_s$ cases. For the two higher shock Mach number cases, $M_t$ varies between $0.3$ and $0.8$ , placing them in the nonlinear subsonic regime where dilatational fluctuations are expected to be important.

3.4. Budget of PTKE

The presence of particles in the flow generates local gas-phase velocity fluctuations characterized by the PTKE, defined as $k_g = (\widetilde {{u}_i^{\prime\prime} {u}_i^{\prime\prime}})/2$ . Reynolds-averaged transport equations for compressible flows have previously been derived by Sarkar et al. (1991), among others. In this study, the transport equation for PTKE is derived in a similar manner, but the presence of particles is accounted for by including the indicator function in the averaging process as defined in § 3.2. Multiplying through the Navier–Stokes equations in § 2.2 by the indicator function and averaging over the homogeneous $y$ - and $z$ -directions yields a one-dimensional, time-dependent transport equation for PTKE (a similar derivation is given by Vartdal & Osnes (2018)), given by

(3.7) \begin{equation} \begin{aligned} \frac {\partial }{\partial t}(\alpha _g \, \overline {\rho } \, k_g) + \frac {\partial }{\partial x}(\alpha _g \, \overline {\rho } \, \widetilde {u} \, k_g) &= \mathcal{P}_S + \mathcal{T} + \varPi + \alpha _g \overline {\rho } \epsilon + \mathcal{M} + \mathcal{P}_D. \end{aligned} \end{equation}

The terms on the right-hand side represent various mechanisms for producing, dissipating and transporting PTKE. Here $\mathcal{P}_{s}$ is production due to mean shear, $\mathcal{T}$ is a term akin to diffusive transport, $\varPi$ is the pressure-dilatation correlation term and $\epsilon$ is the viscous dissipation tensor. The trailing terms arising from the averaging procedure are lumped into $\mathcal{M}$ . Here $\mathcal{P}_D=\mathcal{P}_D^P+\mathcal{P}_D^V$ is production due to drag that contains contributions from pressure and viscous stresses, respectively. These terms are each defined as

(3.8) \begin{align} \mathcal{P}_S &= - \alpha _g \overline {\rho } \widetilde {u^{\prime\prime} u^{\prime\prime}} \frac {\partial \widetilde {u}}{\partial x}, \end{align}

(3.9) \begin{align} \mathcal{T}&= -\frac {1}{2} \frac {\partial }{\partial x}(\alpha _g \overline {\rho {u_i^{\prime\prime} u_i^{\prime\prime} u^{\prime\prime}}} ) - \frac {\partial }{\partial x}(\alpha _g \overline {p^{\prime}u^{\prime\prime}}) + \frac {\partial }{\partial x}(\alpha _g \overline {u_i^{\prime\prime} \sigma _{i1}^{\prime}}), \end{align}

(3.10) \begin{align} {\varPi } &= \alpha _g \ \overline { p^{\prime} \frac {\partial u_i^{\prime\prime}}{\partial x_i} }, \end{align}

(3.11) \begin{align} \alpha _g \overline {\rho } \epsilon & = - \alpha _g \ \overline { \sigma _{ik}^{\prime} \frac {\partial u_i^{\prime\prime}}{\partial x_k} }, \end{align}

(3.12) \begin{align} \mathcal{M} &= -\frac {\partial } {\partial x_i}(\alpha _g \overline {p} \overline {u_i^{\prime\prime}}) + \alpha _g \overline {p} \frac {\partial \overline {u_i^{\prime\prime}}}{\partial x_i} + \frac {\partial } {\partial x_i} (\alpha _g \overline {\sigma _{11}}\overline {u_i^{\prime\prime}}) - \alpha _g \overline {\sigma _{11}} \frac {\partial \overline {u_i^{\prime\prime}}}{\partial x_i}, \end{align}

(3.13) \begin{align} \mathcal{P}_D &= \mathcal{P}_D^P + \mathcal{P}_D^V = \overline {p^{\prime} u_i^{\prime\prime} \frac {\partial \mathcal{I}}{\partial x_i}} - \overline {\sigma _{ik}^{\prime} u_i^{\prime\prime} \frac {\partial \mathcal{I}}{\partial x_k}}. \end{align}

Figure 8 shows the budget of PTKE at $t/\tau _L=0.5$ and $2$ for different $M_s$ and $\varPhi _p$ . The terms are normalized by post-shock quantities: $\rho _{ps}u_{ps}^3/D$ . The statistics from the particle-resolved simulations are noisy due to the indicator function and to provide reliable data, a low-pass (Gaussian) filter is applied in the streamwise direction with a standard deviation $3D$ after averaging in the periodic directions. It should be noted that most coarse-grain simulations of particle-laden flows use grid spacing larger than $D$ . Also, the resulting profiles were found to be insensitive to a wide range of filter sizes. Note in figures 8(c), 8(e) and 8(h), small oscillations upstream of the curtain indicate the location of the reflected shock.

Figure 8. Budgets of PTKE at (a,c,e) $t/\tau _L=0.5$ and (b,d,f) $t/\tau _L=2$ . Here $(a,b)$ $M_s=1.2$ , $(c,d)$ $M_s=1.66$ and $(e,f)$ $M_s=2.1$ . Here $\varPhi _p=0.1$ ( ), $\varPhi _p=0.2$ ( ) and $\varPhi _p=0.3$ ( ).

The majority of PTKE is generated via drag production, which is balanced by viscous dissipation. The remaining terms are negligible except for shear production, $\mathcal{P}_s$ , and the pressure-dilatation correlation term $\varPi$ near the shock and at the edge of the curtain where the volume fraction gradient is large. Here $\mathcal{M}_s$ is omitted from the plots since it was found to be negligible. At later times after the shock has passed through the curtain, mean-shear production and the pressure-strain correlation act as the dominant production terms at the downstream edge of the curtain. Downstream of the curtain, there are no production mechanisms and viscous dissipation dominates.

The magnitude of the terms in the budget are observed to increase with increasing $\varPhi _p$ and decrease with increasing shock Mach number. This reduction at higher Mach number is not due to enhanced dilatational dissipation, but rather a reduction of all terms, similar to what has been observed in single-phase compressible shear layers (Sarkar 1995; Pantano & Sarkar 2002).

3.5. Energy spectra

Two-dimensional energy spectra of the phase-averaged streamwise velocity fluctuations are computed at different locations along the curtain. Special care is taken to account for the presence of particles. At each location along the $x$ -axis, the instantaneous energy spectrum is defined as

(3.14) \begin{equation} E_{uu}(x,t) = \widehat { \sqrt {\mathcal{I} \rho } u^{\prime\prime}} \widehat { \sqrt {\mathcal{I} \rho } u^{\prime\prime}}^*, \end{equation}

where the $\widehat {(\boldsymbol{\cdot })}$ notation denotes the two-dimensional Fourier transform and $^*$ indicates its complex conjugate. The integration of $E_{uu}$ at each streamwise location is taken over a circular shell in the $[\kappa _y \times \kappa _z]$ space, where $\kappa$ represents the wavenumber. This definition of the Fourier coefficient is consistent with the classic compressible turbulence literature (Kida & Orszag 1990; Lele 1992), extended to include the indicator function to account for particles.

Figure 9 shows the energy spectra for Case 5 ( $M_s=1.66$ , $\varPhi _p=0.2$ ) at various $x$ locations within the particle curtain at $t/\tau _L=2$ . The spectra for the initial 40 grid points ( $x-x_0\lt D$ ) are excluded because the turbulence is not fully developed in this region. It is evident that the spectra remain relatively consistent across the streamwise positions, exhibiting minimal variation from the ensemble average of all spectra. Thus, although the flow is inhomogeneous in $x$ , the turbulence is relatively homogeneous in the majority of the curtain. Consequently, subsequent figures will only display the ensemble average.

Figure 9. One-dimensional spectra of streamwise velocity fluctuations for $M_s=1.66$ and $\varPhi _p=0.2$ at $t/\tau _L=2$ . The colour bar corresponds to different locations in the particle curtain. Ensemble average of all the spectra within the curtain (thick cyan line). Vertical lines indicate relevant length scales in the flow. Solid and dashed lines correspond to slopes of $-5/3$ and $-3$ , respectively.

The inclusion of the discontinuous indicator function in (3.14) introduces oscillations throughout the spectrum, known as a ‘ringing’ artefact. While the ringing can be mitigated by applying a Butterworth filter or similar methods, such filtering was not employed to avoid the introduction of ad hoc user-defined parameters.

Most of the energy resides at length scales that coincide with the mean interparticle spacing, $\lambda$ . The interparticle spacing is found to differentiate the energy-containing range from the inertial subrange, indicating that wakes in the interstitial spaces between particles are responsible for the generation of PTKE. An inertial subrange is evident at scales smaller than $\lambda$ , characterized by an energy spectrum that follows a $-5/3$ power law before transitioning to a steeper $-3$ power law at higher wavenumbers. The energy diminishes rapidly at scales below $2\Delta x$ , which is attributed to numerical dissipation. Interestingly, part of the inertial subrange aligns with characteristics of homogeneous single-phase turbulence, displaying a $-5/3$ power law, while the smaller scales align with BIT, evidenced by a $-3$ power law. However, the presence of noise in the spectra makes it challenging to draw definitive conclusions.

In figure 10, the ensemble-averaged spectra are compared across different cases at $t/\tau _L = 2$ . A broadband reduction in $E_{uu}$ is observed with increasing $M_s$ , which is consistent with the observations made in the PTKE budget. As before, the turbulence levels are largely invariant with $\varPhi _p$ . For each case, the mean interparticle spacing is found to delineate the inertial subrange. Compensated spectra are also shown to better identify the power-law scaling, which appears consistent in each case. The spectrum decays with a $-5/3$ law for non-dimensional wavenumbers roughly between $5$ and $40$ , while at higher wavenumbers there is a steeper $-3$ decay. It remains unclear whether this steepening is due to gas-phase compressibility, interphase exchange with particles, or both. The following section decomposes the turbulent velocity field into solenoidal and dilatational components to gain further insight.

Figure 10. (a,c,e) Mean and (b,d,f) compensated energy spectra of streamwise velocity fluctuations within the particle curtain at $t/\tau _L=2$ for $(a,b)$ $\varPhi _p=0.1$ , $(c,d)$ $\varPhi _p=0.2$ and $(e,f)$ $\varPhi _p=0.3$ . Here $M_s=1.2$ (light blue, square), $M_s=1.66$ (lavender, circle) and $M_s=2.1$ (purple, triangle).

3.5.1. Helmholtz decomposition

A Helmholtz decomposition of the velocity field is performed to analyse the solenoidal (divergence-free) and dilatational (curl-free) components separately, according to (Kida & Orszag 1990; Yu, Xu & Pirozzoli 2019)

(3.15) \begin{equation} \boldsymbol{u} = \boldsymbol{u}_{sol} +\boldsymbol{u}_{dil}, \end{equation}

where $\boldsymbol{u}_{sol} = \boldsymbol{\nabla }\times \boldsymbol{A}$ and $\boldsymbol{u}_{dil} = \boldsymbol{\nabla }\varphi$ . Here $\boldsymbol{A}$ is the vector potential satisfying ${\nabla} ^2 \boldsymbol{A} = - \boldsymbol{\omega }$ , where $\boldsymbol{\omega }=\boldsymbol{\nabla }\times \boldsymbol{u}$ is the local vorticity. The velocity potential $\varphi$ satisfies ${\nabla} ^2 \varphi = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$ .

Figure 11 shows two-dimensional slices of the instantaneous streamwise velocity components. The solenoidal component exhibits significant fluctuations throughout the curtain, capturing particle wakes. In contrast, the dilatational velocity field remains relatively small within the curtain and increases sharply at the downstream edge, where the flow chokes. This indicates that the majority of PTKE is concentrated in the solenoidal portion, with compressibility playing a minor role except near large volume fraction gradients.

Figure 11. A two-dimensional slice of the $(a)$ solenoidal and $(b)$ dilatational streamwise velocity fields at $t/\tau _L=2$ for Case $5$ .

Figure 12 shows energy spectra of the streamwise solenoidal and dilatational velocity components at $t/\tau _L=2$ . The solenoidal energy spectrum is approximately two orders of magnitude larger than the dilatational component across all wavenumbers and tends to decrease with increasing $M_s$ , while the dilatational component increases with increasing Mach number. These findings align with observations from direct numerical simulations of compressible homogeneous isotropic turbulence (Donzis & Jagannathan 2013). Interestingly, only the solenoidal spectrum demonstrates a $-3$ power law decay, while the dilatational component maintains an approximate $-5/3$ scaling throughout the inertial subrange. Consequently, the $-3$ power law decay may be attributed to incompressible wakes rather than compressible effects.

Figure 12. (a,c,e) Mean and (b,d,f) compensated spectra of the streamwise velocity fluctuations computed using solenoidal ( ) and dilatational ( ) velocity fields at $t/\tau _L=2$ for $(a,b)$ $\varPhi _p=0.1$ , $(c,d)$ $\varPhi _p=0.2$ and $(e,f)$ $\varPhi _p=0.3$ . Colour scheme same as figure 10.

A similar decomposition of the dissipation rate was performed following Sarkar et al. (1991) and Donzis & John (2020), yielding solenoidal ( $\bar {\rho }\epsilon _s= \overline {\mu \omega _i^{\prime\prime} \omega _i^{\prime\prime}}$ ; $\boldsymbol{\omega }^{\prime\prime}=\boldsymbol{\nabla }\times \boldsymbol{u}^{\prime\prime}$ ) and dilatational ( $\bar {\rho } \epsilon _d =\overline {4/3\mu ( {\partial u_i^{\prime\prime}}/{\partial x_i} )^2}$ ) components. It was found that $\epsilon _s$ dominates the total dissipation of PTKE, exceeding $\epsilon _d$ by approximately two orders of magnitude.

4. Two-fluid turbulence model

In this section, we propose a two-equation model for PTKE and its dissipation. This turbulence model is integrated into a one-dimensional Eulerian-based two-fluid framework. The hyperbolic equations for particle-laden compressible flows include added mass and internal energy contributions, derived from kinetic theory based on the recent work of Fox (2019) and Fox, Laurent & Vié (2020). The section ends with an a posteriori analysis of the turbulence model and comparisons are made with the particle-resolved simulations.

4.1. A kinetic-based hyperbolic two-fluid model

Particle-resolved simulations require grid spacing significantly smaller than the particle diameter to adequately resolve boundary layers and capture relevant aerodynamic interactions. Eulerian-based two-fluid models are a widely used coarse-grained modelling approach that assume the properties of both solid and fluid phases can be expressed as interpenetrating continua interacting through interphase drag terms. Unlike particle-resolved simulations, the computational cost of modelling the particle phase scales with the number of grid cells rather than the number of particles, making it a more efficient option for simulating systems with a large number of particles.

The added mass is included in the mass, momentum and energy balances, augmented to account for particle wakes. These equations are fully hyperbolic and avoid the ill-posedness common in conventional compressible two-fluid models with two-way coupling (Fox et al. 2020). To match the conditions used in the particle-resolved simulations, stationary monodisperse particles are considered (i.e. the particle velocity $\boldsymbol{u}_p=0$ , granular temperature $\varTheta _p =0$ and $\alpha _p \rho _p$ is constant in the curtain, where $\alpha _p=1-\alpha _g$ is the particle volume fraction and $\rho _p$ is the particle density). Heat transfer between the phases is neglected. For brevity, brackets and tildes are omitted and it is implied that the equations are written in terms of Favre- and phase-averaged quantities.

The governing equation for mass balance (added mass, gas phase) in one spatial dimension are given by

(4.1) \begin{equation} \begin{aligned} \frac {\partial }{\partial t} (\alpha _a \rho _a ) &= S_a, \\ \frac {\partial }{\partial t} (\alpha _g^\star \rho ) + \frac {\partial }{\partial x} (\alpha _g^\star \rho u) & = -S_a. \end{aligned} \end{equation}

The gas-phase momentum balance is

(4.2) \begin{equation} \begin{aligned} \frac {\partial }{\partial t} (\alpha _g^\star \rho u) + \frac {\partial }{\partial x} (\alpha _g^\star \rho u^2 + \hat {p} + \alpha _p^\star \alpha _g^\star \rho u^2 ) &= -\frac {\alpha _p^\star \rho }{\tau _p} u + \alpha _p^\star \Big (\frac {\partial }{\partial x} \hat {p} + F_{pg} \Big ) - S_{gp}, \end{aligned} \end{equation}

and the gas-phase total energy balance is

(4.3) \begin{equation} \begin{gathered} \frac {\partial }{\partial t} (\alpha _g^\star \rho E ) + \frac {\partial }{\partial x} (\alpha _g^\star \rho u E + \alpha _g^\star u \hat {p} ) = - S_E. \end{gathered} \end{equation}

The added-mass internal energy balance is

(4.4) \begin{equation} \begin{gathered} \frac {\partial }{\partial t} (\alpha _a \rho _a e_a) = S_E. \end{gathered} \end{equation}

Here, $\alpha _a$ is the volume fraction of the added-mass phase and $\rho _a$ is its density. The gas-phase volume fraction is replaced by $\alpha _g^\star =\alpha _g-\alpha _a$ , $\alpha _p^\star = \alpha _p + \alpha _a$ , $\alpha _g^\star = 1 - \alpha _p^\star$ and $e_a$ is the specific internal energy of the added mass. The gas- and added-mass phases have the same pressure $p$ , but different temperatures $T$ and $T_a$ , found from $e$ and $e_a$ , respectively. Here $S_a$ represents mass exchange between the two phases through added mass, leading to momentum $S_{gp}$ and energy $S_E$ exchange. The particle response time scale $\tau _p=4\rho _p D^2/(3\mu C_D Re_p)$ depends on the drag coefficient $C_D$ modelled using the drag law from Osnes et al. (2023). This model takes into account the effects of local volume fraction, the particle Reynolds number $Re_p=\rho |u| D/\mu$ and particle Mach number $M_p=|u|/c$ based on slip velocity $|u|$ ( $|u| = |u - u_p|$ , $u_p=0$ ). The remaining parameters are provided in Appendix B. Note that PTKE contributes to the modified pressure $\hat {p}$ .

The equations are solved using a standard finite-volume method implemented in MATLAB. A HLLC (Harten–Lax–van Leer–Contact) scheme (Toro, Spruce & Speares 1994) is employed to solve the hyperbolic part of the system. Further details on the implementation and discretization of the one-dimensional two-fluid model can be found in Boniou & Fox (2023).

4.2. Two-equation model for PTKE

To capture PTKE in the Eulerian framework, a two-equation $k_g$ – $\epsilon$ model is proposed that retains only the significant source terms from the budget,

(4.5) \begin{equation} \frac {\partial }{\partial t}(\alpha _g^\star \rho k_g) + \frac {\partial }{\partial x}(\alpha _g^\star \rho k_g u \ ) = \mathcal{P}_s + \mathcal{P}_D - (1+M_t^2) {\alpha _g^\star \rho } \epsilon , \end{equation}

(4.6) \begin{equation} \frac {\partial }{\partial t}(\alpha _g^\star \rho \epsilon ) + \frac {\partial }{\partial x}(\alpha _g^\star \rho \epsilon u \ ) = C_{\epsilon ,1} \frac {\epsilon }{k_g} \mathcal{P}_S + \frac {C_{\epsilon ,D}}{\tau _D}\mathcal{P}_D - C_{\epsilon ,2} \ \alpha _g^\star \rho \frac {\epsilon ^2}{k_g}, \end{equation}

where $C_{\epsilon ,1}=1.44$ and $C_{\epsilon ,2}=1.92$ are constants from single-phase turbulence modelling (Pope 2000). The mean-shear production term is $\mathcal{P}_s = - \alpha _g^\star \rho \widetilde {u^{\prime\prime}u^{\prime\prime}} (\partial u / \partial x)$ . Drag production is $\mathcal{P}_D = \alpha _p^\star \rho _p u^2/\tau _p$ . Here $\tau _D= D /|u|$ is the rate of drag dissipation. Here $\epsilon$ represents the solenoidal component of dissipation, while the dilatational component is captured by the compressibility correction ( $1+M_t^2$ ) (Sarkar et al. 1991).

The mean-shear production term, $\mathcal{P}_s$ , includes the streamwise component of the Reynolds stress, $\widetilde {u^{\prime\prime}u^{\prime\prime}}$ . Based on the findings from § 3.3, the anisotropy was found to be relatively constant across the curtain and independent of volume fraction and shock Mach number (see figure 5). The streamwise component of Reynolds stress is therefore given by

(4.7) \begin{equation} \begin{aligned} \widetilde {u^{\prime\prime}u^{\prime\prime}} = 2 \left ( b_{11} + \frac {1}{3} \right ) k_g,\\ \widetilde {v^{\prime\prime}v^{\prime\prime}} = \widetilde {w^{\prime\prime}w^{\prime\prime}} = 2 \left ( b_{22} + \frac {1}{3} \right ) k_g, \end{aligned} \end{equation}

with $b_{11}=0.2$ and $b_{22}=-0.1$ .

In the two-equation model, the only remaining term requiring closure is $C_{\epsilon ,D}$ , a model coefficient that controls the portion of PTKE produced through drag that ultimately gets dissipated. In the limits of homogeneity and steady state with $M_t=0$ , the unsteady, convective and mean shear production terms go to zero, thus reducing (4.5) and (4.6) to

(4.8) \begin{equation} \mathcal{P}_D = \alpha _g^\star \rho \epsilon \implies \epsilon = \frac {\alpha _p^\star \rho _p}{\alpha _g^\star \rho } \frac {u^2}{\tau _p} \end{equation}

and

(4.9) \begin{equation} \frac {\mathcal{C}_{\epsilon ,D}}{\tau _D}\mathcal{P}_D = C_{\epsilon ,2} \alpha _g^\star \rho \frac {\epsilon ^2}{k_g} \implies C_{\epsilon ,D} = \frac {C_{\epsilon ,2} \tau _D \tau _p}{u^2} \frac {\alpha _g^\star \rho } {\alpha _p^\star \rho _p}\frac {\epsilon ^2}{k_g} , \end{equation}

respectively. Substituting (4.8) into (4.9) yields

(4.10) \begin{equation} \frac {k_g}{u^2} = \frac {C_{\epsilon ,2}}{C_{\epsilon ,D} }\frac {\alpha _p^\star \rho _p}{\alpha _g^\star \rho } \frac {\tau _D}{\tau _p} \implies C_{\epsilon ,D} = \frac {3}{4} \frac {\alpha _p^\star }{\alpha _g^\star }\frac {u^2}{k_g} C_{\epsilon ,2} C_D. \end{equation}

We now calibrate the algebraic model for PTKE proposed by Mehrabadi et al. (2015) for homogeneous particle suspensions valid for $ \alpha _gRe_p \lt 300$ , $M_p=0$ and $\alpha _p \geqslant 0.1$ as

(4.11) \begin{equation} \frac {k_g}{u^2} = 2\alpha _p \left (1+1.25\alpha _g^3\exp (-\alpha _p \sqrt {\alpha _gRe_p}) \right ). \end{equation}

Plugging (4.11) into (4.10) provides closure for $C_{\epsilon ,D}$ and ensures the model returns the correct level of PTKE in the limit of incompressible, homogeneous, steady flow. Because $\alpha _p^\star \to \alpha _p$ when $\alpha _p \to 0$ , $C_{\epsilon ,D}$ remains finite outside the particle curtain. With the expression for $C_{\epsilon ,D}$ , the two-equation transport model for PTKE and dissipation is now fully closed and is implemented in the two-fluid model from § 4.1.

4.3. A posteriori analysis

One-dimensional shock–particle interactions are simulated using the two-fluid model detailed above with the parameters used in the particle-resolved simulations. It should be noted that the results will depend significantly on the volume fraction profile. To ensure a fair comparison, one-dimensional volume fraction profiles are extracted from the particle-resolved simulations and used in the model (see figure 13).

Figure 13. One-dimensional particle volume fraction profiles obtained from the particle-resolved simulations for $\varPhi _p=0.1$ (light pink), $\varPhi _p=0.2$ (pink) and $\varPhi _p=0.3$ (purple).

Figure 14 shows comparisons of the mean streamwise velocity between the two-equation model and particle-resolved simulations. Overall excellent agreement is observed. The locations of the transmitted and reflected shocks are predicted correctly. The model can be seen to predict choked flow at the downstream edge resulting in supersonic expansion, closely matching the particle-resolved simulations.

Figure 14. Comparison of mean streamwise velocity from particle-resolved simulations ( ) with results from the two-equation model ( ). Here $(a,d)$ $M_s=1.2$ , $(b,e)$ $M_s=1.66$ , $(c,f)$ $M_s=2.1$ . Here (a–c) $t/\tau _L=1$ and (d–f) $t/\tau _L=2$ . The colour scheme for different volume fraction cases is the same as in figure 13.

Figure 15 shows comparisons of PTKE between the two-equation model and particle-resolved simulations at two time instances. Results show good agreement for all cases considered except for the cases with $M_s=1.2$ at higher volume fractions. The model predicts an increase in PTKE with $\varPhi _p$ , which is not observed in the particle-resolved simulations. Despite this, the model results show overall good agreement both within the curtain and downstream.

Figure 15. Comparison of PTKE between particle-resolved simulations ( ) with the two-equation model ( ). Here $(a,d)$ $M_s=1.2$ , $(b,e)$ $M_s=1.66$ , $(c,f)$ $M_s=2.1$ . Here (a–c) $t/\tau _L=1$ and (d–f) $t/\tau _L=2$ . The colour scheme for different volume fraction cases is the same as in figure 13.

The terms in the PTKE budget computed from the two-equation model are compared with particle-resolved simulation data to identify and explain the observed discrepancies in PTKE. Specifically, the dominant terms – drag production $\mathcal{P}_D$ , viscous dissipation $\alpha _g \rho \epsilon$ and mean-shear production $\mathcal{P}_s$ – are examined. Figure 16 presents the comparison for one case, with similar results observed across all cases. Overall, excellent agreement is found over the three time instances shown. The largest discrepancies occur at the upstream and downstream edges of the curtain, where the particle-resolved simulations predict a sharper increase in drag production at the upstream edge and greater dissipation at the downstream edge. These differences may be attributed to numerical diffusion in the coarse-grained model.

Figure 16. Comparison of terms in the PTKE budget between the two-equation model ( ) and particle-resolved simulations ( ) for $(a)$ $M_s=1.2$ and $\varPhi _p=0.2$ , $(b)$ $M_s=1.66$ and $\varPhi _p=0.2$ and $(c)$ $M_s=2.1$ and $\varPhi _p=0.3$ at $t/\tau _L=2$ . Here $\mathcal{P}_{D}^P$ (red), $\mathcal{P}_s$ (purple) and $\alpha _g \rho _g \epsilon _g$ (blue).

The streamwise and spanwise fluctuations are reconstructed using (4.7) and compared with particle-resolved simulations in figure 17. Overall, the results show good agreement, with Cases $2$ and $3$ exhibiting the most discrepancies. These discrepancies may arise from the drag model, the choice of $C_{\epsilon ,D}$ , or the omission of the viscous term in the two-fluid model. At the downstream edge, streamwise fluctuations are slightly underpredicted, likely due to an overprediction of viscous dissipation, as observed in the previous figure. This overprediction is ultimately linked to the choice of $C_{\epsilon ,D}$ or the drag model. Despite these issues, the two-equation model predicts the overall behaviour well, including the PTKE downstream of the curtain, in the pure gas.

Figure 17. Comparison of pseudo-turbulent Reynolds stresses between the particle-resolved simulations ( ) and the model ( ). Here $(a,d)$ $M_s=1.2$ , $(b,e)$ $M_s=1.66$ , $(c,f)$ $M_s=2.1$ . Colour scheme defined in figure 4.

The gas-phase turbulence downstream of the curtain lacks any production terms and, according to the budget, should only advect and diffuse. The cases considered so far extend only $6D$ from the downstream curtain edge to the right domain boundary. Here, we examine Case 10 from table 1, with $M_s=1.66$ , $\varPhi _p=0.3$ , and a domain extending $34D$ ( $2L$ ) downstream. Figure 18 shows PTKE comparisons after the flow reaches a statistically stationary state. The PTKE decay resembles grid-generated turbulence, and the model captures the turbulence transport and decay well.

Figure 18. Comparison of PTKE of the longer domain (Case 10) between particle-resolved simulations ( ) and the two-equation model ( ) at $(a)$ $t/\tau _L = 3$ , $(b)$ $t/\tau _L = 4$ and $(c)$ $t/\tau _L = 5$ .

5. Conclusions

When a shock wave interacts with a suspension of solid particles, momentum and energy exchanges between the phases give rise to complex flow. Particle wakes induced by the transmitted shock generate velocity fluctuations referred to as ‘pseudo-turbulence’. Phase-averaging the viscous compressible Navier–Stokes equations reveals a route for turbulence generation through drag production within the particle curtain and localized mean-shear production at the edge of the curtain. This turbulence generation is balanced by viscous and dilatational dissipation.

Three-dimensional particle-resolved simulations of planar shocks interacting with stationary spherical particles were used to analyse the characteristics of pseudo-turbulence for a range of shock Mach numbers and particle volume fractions. In each case, PTKE is generated through interphase drag coupling, contributing to $20\,\%{-}50\,\%$ of the post-shock kinetic energy. The abrupt change in volume fraction at the downstream edge of the curtain chokes the flow, resulting in supersonic expansion where PTKE is maximum. The pseudo-turbulent Reynolds stress is highly anisotropic but approximately constant throughout for the range of volume fractions and Mach numbers considered. The energy spectra of the streamwise gas-phase velocity fluctuations reveal an inertial subrange that begins at the mean interparticle spacing and decays with a $-5/3$ power law then steepens to $-3$ at smaller scales. This $-3$ scaling only exists in the solenoidal component of the velocity field and is attributed to particle wakes.

A one-dimensional two-equation turbulence model was formulated for PTKE and its dissipation and implemented within a hyperbolic two-fluid framework. Drag production is closed using a drag coefficient that takes into account local volume fraction, Reynolds number and Mach number. A new closure is proposed for drag dissipation that ensures the proper amount of PTKE is obtained in the limit of statistically stationary and homogeneous flow. An a posteriori analysis demonstrated the ability of the model to predict PTKE accurately during shock–particle interactions and capture flow-choking behaviour. Such a turbulence model can be adopted into Eulerian two-fluid models or Eulerian–Lagrangian frameworks.