1. Introduction
High-speed flows through particulate media occur in diverse applications, such as detonation blasts (Zhang et al. Reference Zhang, Frost, Thibault and Murray2001), volcanic eruptions (Chojnicki, Clarke & Phillips Reference Chojnicki, Clarke and Phillips2006; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020), coal-dust explosions (Sapko et al. Reference Sapko, Weiss, Cashdollar and Zlochower2000; Zheng et al. Reference Zheng, Feng, Jing, Qian, Li, Liu and Huang2009), pulsed-detonation engines (Chang & Kailasanath Reference Chang and Kailasanath2003; Roy et al. Reference Roy, Frolov, Borisov and Netzer2004) and plume–surface interactions during interplanetary landings (Plemmons et al. Reference Plemmons, Mehta, Clark, Kounaves, Peach, Renno, Tamppari and Young2009; Morris et al. Reference Morris, Goldstein, Varghese and Trafton2011; Capecelatro Reference Capecelatro2022). 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 Reference Sagaut and Cambon2008; Jagannathan & Donzis Reference Jagannathan and Donzis2016). 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 Reference Brown and Roshko1974; Bradshaw Reference Bradshaw1977; Sarkar et al. Reference Sarkar, Erlebacher, Hussaini and Kreiss1991). Early work by Zeman (Reference Zeman1990) and Sarkar et al. (Reference Sarkar, Erlebacher, Hussaini and Kreiss1991) 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 (Reference Sarkar1995) 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 (Reference Vreman, Sandham and Luo1996) and Pantano & Sarkar (Reference Pantano and Sarkar2002) 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 (Reference Kida and Orszag1990) 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 (Reference Donzis and Jagannathan2013) 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, Reference Kolmogorov1941b
). 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 Reference Kida and Orszag1990; Donzis & Jagannathan Reference Donzis and Jagannathan2013; Wang et al. Reference Wang, Shi, Wang, Xiao, He and Chen2011, Reference Wang, Shi, Wang, Xiao, He and Chen2012; San & Maulik Reference San and Maulik2018). Compressibility effects are typically attributed to
$\boldsymbol{u}_d$
, and both Wang et al. (Reference Wang, Shi, Wang, Xiao, He and Chen2011) and Donzis & Jagannathan (Reference Donzis and Jagannathan2013) 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 Reference Lance and Bataille1991; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015), a term also applied to bubble-induced turbulence (BIT) in liquid flows (Risso Reference Risso2018). Lance & Bataille (Reference Lance and Bataille1991) 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. Reference Mercado, Gomez, Van Gils, Sun and Lohse2010; Risso Reference Risso2018; Mezui et al. Reference Mezui, Obligado and Cartellier2022, Reference Mezui, Cartellier and Obligado2023). 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 Reference Risso2018). 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 Reference Capecelatro, Desjardins and Fox2018).
Over the past few decades, turbulence models have evolved to incorporate the effects of particles (Troshko & Hassan Reference Troshko and Hassan2001; Fox Reference Fox2014; Ma et al. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017). 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 (Reference Elghobashi and Abou-Arab1983), 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 Reference Ma and Ahmadi1990). Crowe, Troutt & Chung (Reference Crowe, Troutt and Chung1996) 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. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015) 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 Reference Shallcross, Fox and Capecelatro2020).
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 (Reference Khalloufi and Capecelatro2023) 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 Reference Sommerfeld1994; Patel et al. Reference Patel, Rubio, Shekhtman, Parziale, Rabinovitch, Ni and Capecelatro2024). Two-dimensional particle-resolved simulations of shock–particle curtain interactions revealed PTKE magnitudes comparable to the resolved kinetic energy (Regele et al. Reference Regele, Rabinovitch, Colonius and Blanquart2014; Hosseinzadeh-Nik, Subramaniam & Regele Reference Hosseinzadeh-Nik, Subramaniam and Regele2018). In three-dimensional inviscid simulations, Mehta, Jackson & Balachandar (Reference Mehta, Jackson and Balachandar2020) 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. Reference McFarland, Black, Dahal and Morgan2016; Middlebrooks et al. Reference Middlebrooks, Avgoustopoulos, Black, Allen and McFarland2018). 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. (Reference Osnes, Vartdal, Omang and Reif2019) 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. (Reference Shallcross, Fox and Capecelatro2020) 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. Reference Ma, Santarelli, Ziegenhein, Lucas and Fröhlich2017) – 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 Reference Ling, Haselbacher and Balachandar2011). The simulations are designed to emulate the multiphase shock-tube experiments of Wagner et al. (Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012). 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.
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.
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
\begin{equation} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}) = 0, \end{equation}
\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
\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
\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
\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
\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 Reference Capecelatro2023). Spatial derivatives are approximated using narrow-stencil finite-difference operators that satisfy the summation by parts (SBP) property (Strand Reference Strand1994; Svärd et al. Reference Svärd, Carpenter and Nordström2007). 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 Reference Pirozzoli2011), 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. Reference Svärd, Carpenter and Nordström2007). High-order SBP dissipation operators (Mattsson, Svärd & Nordström Reference Mattsson, Svärd and Nordström2004) 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 (Reference Kawai, Shankar and Lele2010). To limit the artificial diffusivity to regions of high compression (shocks), we employ the sensor originally proposed by Ducros et al. (Reference Ducros, Ferrand, Nicoud, Weber, Darracq, Gacherieu and Poinsot1999) and later improved by Hendrickson, Kartha & Candler (Reference Hendrickson, Kartha and Candler2018). More details can be found in Khalloufi & Capecelatro (Reference Khalloufi and Capecelatro2023). 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 (Reference Mohd-Yusof1997) and later extended to compressible flows by Chaudhuri, Hadjadj & Chinnayya (Reference Chaudhuri, Hadjadj and Chinnayya2011) 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 Reference Khalloufi and Capecelatro2023), 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.
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
\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
\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
\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. Reference Mehta, Neal, Salari, Jackson, Balachandar and Thakur2018; Theofanous, Mitkin & Chang Reference Theofanous, Mitkin and Chang2018; Osnes et al. Reference Osnes, Vartdal, Omang and Reif2019). Mehta et al. (Reference Mehta, Neal, Salari, Jackson, Balachandar and Thakur2018) 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.
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.
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. (Reference Mehta, Jackson and Balachandar2020) 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
\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. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015). 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. (Reference Osnes, Vartdal, Omang and Reif2019). 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}}$
.
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 Reference Batchelor and Townsend1948; Mohamed & Larue Reference Mohamed and Larue1990; Kurian & Fransson Reference Kurian and Fransson2009). According to Batchelor & Townsend (Reference Batchelor and Townsend1948), the decay of turbulence intensity downstream of a grid (or screen) with mesh width
$\Delta M$
follows a power law, given by
\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
\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.
$(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 Reference Mohamed and Larue1990; Kurian & Fransson Reference Kurian and Fransson2009). 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.
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 (Reference Sagaut and Cambon2008), 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. (Reference Sarkar, Erlebacher, Hussaini and Kreiss1991), 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 (Reference Vartdal and Osnes2018)), given by
\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
\begin{align} \mathcal{P}_S &= - \alpha _g \overline {\rho } \widetilde {u^{\prime\prime} u^{\prime\prime}} \frac {\partial \widetilde {u}}{\partial x}, \end{align}
\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}
\begin{align} {\varPi } &= \alpha _g \ \overline { p^{\prime} \frac {\partial u_i^{\prime\prime}}{\partial x_i} }, \end{align}
\begin{align} \alpha _g \overline {\rho } \epsilon & = - \alpha _g \ \overline { \sigma _{ik}^{\prime} \frac {\partial u_i^{\prime\prime}}{\partial x_k} }, \end{align}
\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}
\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.
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 Reference Sarkar1995; Pantano & Sarkar Reference Pantano and Sarkar2002).
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
\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 Reference Kida and Orszag1990; Lele Reference Lele1992), 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.
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.
(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 Reference Kida and Orszag1990; Yu, Xu & Pirozzoli Reference Yu, Xu and Pirozzoli2019)
\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.
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 Reference Donzis and Jagannathan2013). 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.
(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. (Reference Sarkar, Erlebacher, Hussaini and Kreiss1991) and Donzis & John (Reference Donzis and John2020), 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 (Reference Fox2019) and Fox, Laurent & Vié (Reference Fox, Laurent and 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. Reference Fox, Laurent and Vié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
\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
\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
\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
\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. (Reference Osnes, Vartdal, Khalloufi, Capecelatro and Balachandar2023). 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 Reference Toro, Spruce and Speares1994) 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 (Reference Boniou and Fox2023).
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,
\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}
\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 Reference Pope2000). 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. Reference Sarkar, Erlebacher, Hussaini and Kreiss1991).
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
\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
\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
\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
\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. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015) for homogeneous particle suspensions valid for
$ \alpha _gRe_p \lt 300$
,
$M_p=0$
and
$\alpha _p \geqslant 0.1$
as
\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).
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.
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.
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.
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.
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.
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.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10831.
Acknowledgements
The authors would like to acknowledge the computing resources and assistance provided by Advanced Research Computing at the University of Michigan, Ann Arbor. The authors would also like to thank the resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.
Funding
A portion of this work was supported by the National Aeronautics and Space Administration (NASA) SBIR contract no. 80NSSC20C0243 and SBIR contract no. HDTRA125P0006.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available in GitHub at https://github.com/jessecaps/jCODE (Capecelatro Reference Capecelatro2023).
Appendix A. Convergence studies
This section quantifies the effects of varying domain size and particle configurations within the curtain in particle-resolved simulations. A grid refinement study of the numerical solver for periodic compressible flow over a homogeneous suspension is detailed in our previous work (Khalloufi & Capecelatro Reference Khalloufi and Capecelatro2023).
A.1. Effect of domain size
In this section, we examine the effects of varying the domain size in the periodic spanwise (
$y$
and
$z$
) directions. A series of three-dimensional simulations were performed with
$M_s=1.66$
and
$\varPhi _p=0.2$
to evaluate the impact of domain size on the individual terms in the PTKE budget. Table 2 summarizes the cases considered. The streamwise domain length
$L_x$
is kept constant, while the spanwise dimensions
$L_y$
and
$L_z$
are varied. Uniform grid spacing is maintained at
$\Delta =D/40$
.
Parameters used for the domain size study. For each case,
${M}_s=1.66$
and
$\varPhi _p=0.2$
.
Effect of domain size on the PTKE budget for
$M_s=1.66$
and
$\varPhi _p=0.2$
at
$t/\tau _L=1.5$
. Case A (
), Case B (
), Case C (
). Colours correspond to figure 8.
Figure 19 presents comparisons of the individual PTKE budget terms. The results indicate that variations in the periodic domain lengths have minimal influence on the budget terms, suggesting that volume-averaging over two-dimensional
$y{-}z$
slices can be performed without significantly affecting the one-dimensional statistics. Consequently, for the case studies presented in the main paper, we adopt
$L_y=L_z=12D$
.
A.2. Effect of varying random particle distributions
Particles are randomly distributed within the curtain while avoiding overlap. The drag force on individual particles is known to depend on the arrangement of their neighbours (Akiki, Moore & Balachandar Reference Akiki, Moore and Balachandar2017; Lattanzi et al. Reference Lattanzi, Tavanashad, Subramaniam and Capecelatro2022; Osnes et al. Reference Osnes, Vartdal, Khalloufi, Capecelatro and Balachandar2023). This section investigates the impact of different random particle configurations within the curtain on PTKE for
$M_s=1.2$
and
$\varPhi _p=0.3$
. Three distinct realizations are considered, keeping all parameters constant except for the random arrangement of particles.
Figure 20 shows the PTKE budget terms for each realization at
$t/\tau _L=1$
, when the shock has just passed the downstream edge of the curtain. All realizations exhibit similar trends with negligible discrepancies. Therefore, we conclude that the random distribution of particles does not significantly affect the statistics.
Effect of random particle placement on the PTKE budget for
$M_s=1.2$
and
$\varPhi _p=0.3$
at
$t/\tau _L=1$
. Realization 1 (
), realization 2 (
), realization 3 (
). Colours correspond to figure 8.
Appendix B. One-dimensional two-fluid model parameters
Starting from the conserved variables
$X_1 = \alpha _a \rho _a$
and
$X_2 = \alpha _g^\star \rho$
with known
$\alpha _p$
, the primitive variables are found using the following formulae:
\begin{equation} \hat {\kappa } = \frac {X_1}{X_2}; \quad \kappa = \frac {T}{T_a}; \quad T = \frac {\gamma e}{C_p}; \quad T_a = \frac {\gamma e_a}{C_p}; \quad e = E - \frac {1}{2}u^2 - k_g ; \end{equation}
\begin{equation} \alpha _g = 1- \alpha _p; \quad \alpha _a = \frac {\hat {\kappa }}{\hat {\kappa }+\kappa } \alpha _g; \quad \alpha _p^\star = \alpha _p+ \alpha _a; \quad \alpha _g^\star = \alpha _g - \alpha _a; \quad \rho = \frac {X_2}{\alpha _g^\star }; \end{equation}
\begin{equation} p = (\gamma -1)\rho e ; \quad \hat {p} = p + 2/3 \rho k_g. \end{equation}
The remaining model parameters are defined as follows:
\begin{equation} P_{pfp} = \rho ( \alpha _p^\star u )^2 ; \quad F_{pg} = u^2 \partial _x \rho + 2/3 \rho (\partial \alpha _g^\star u / \partial x)u; \end{equation}
\begin{equation} S_a = \frac {\rho }{\tau _a}(c_m^\star \alpha _p \alpha _g - \alpha _a); \quad S_{gp}=\max (S_a,0)u ; \quad S_E = \max (S_a,0)E + \min (S_a,0)e_a; \end{equation}
\begin{equation} Re_p = \frac {\rho D u}{\mu }; \quad Pr = \frac {C_p \mu }{k}; \end{equation}
\begin{equation} c_m^\star = \frac {1}{2}\min (1+2\alpha _p,2); \quad \tau _a = 0.001\tau _p ; \quad \tau _p = \frac {4\rho _p D^2}{3 \mu C_DRe_p}. \end{equation}
The drag coefficient
$C_D$
is given by Osnes et al. (Reference Osnes, Vartdal, Khalloufi, Capecelatro and Balachandar2023).
Open access
