Particle-resolved simulations of shock-induced flow through particle clouds at different Reynolds numbers

TL;DR

This paper uses particle-resolved large eddy simulations to analyze how shock-induced flows through particle clouds depend on Reynolds number, revealing flow behavior, force statistics, and implications for modeling dispersed flows.

Contribution

It introduces detailed LES-based analysis of shock-particle cloud interactions across Reynolds numbers, enhancing understanding and modeling of such flows.

Findings

Reflected shock strength increases as Reynolds number decreases.

Maximum particle drag aligns with inviscid predictions.

Flow velocity correlates with flow fluctuation magnitudes.

Abstract

This study investigates the Reynolds-number dependence of shock-induced flow through particle layers at 10\% volume fraction, using ensemble-averaged results from particle-resolved large eddy simulations. The advantage of using large eddy simulations to study this problem is that they capture the strong velocity shears and flow separation caused by the no-slip condition at the particle surfaces. The shock particle cloud interaction produces a reflected shock wave, whose strength increases with decreasing particle Reynolds number. This results in important changes to the flow field that enters the particle cloud. The results show an approximate proportionality between the mean flow velocity and the flow fluctuation magnitudes. Maximum particle drag forces are in excellent agreement with previous inviscid studies, and we complement these results with statistics of time-averaged particle forces as well as the variation of temporal oscillations. The results of this work provides a basis for development of improved simplified dispersed flow models.

Tables

Table 1. Table 1: Average ratio of viscous to pressure forces acting on the particles (force ratio) over the time frame where the shock wave interacts with each particle, the Mach numbers of the transmitted ( Ma TS subscript Ma TS \mathrm{Ma}_{\mathrm{TS}} ) and the reflected ( Ma RS subscript Ma RS \mathrm{Ma}_{\mathrm{RS}} ) shock waves, the particle Reynolds number based on the flow state behind the reflected shock wave ( Re p , RS subscript Re p RS \mathrm{Re}_{\mathrm{p,RS}} ), and the local Mach number of the flow behind the reflected shock wave ( Ma ∗ superscript Ma \mathrm{Ma}^{*} ).

${\mathrm{Re}}_p$	Force ratio	${\mathrm{Ma}}_{\mathrm{TS}}$	${\mathrm{Ma}}_{\mathrm{RS}}$	${\mathrm{Re}}_{p,\mathrm{RS}}$	${\mathrm{Ma}}^{\ast }$
$500$	$0.142$	$2.076$	$1.580$	$360$	$0.379$
$2000$	$0.059$	$2.139$	$1.558$	$1495$	$0.404$
$5000$	$0.028$	$2.174$	$1.530$	$3904$	$0.435$
$10000$	$0.015$	$2.192$	$1.516$	$7967$	$0.451$

Equations

$${\partial }_t\rho +{\partial }_k\left(\rho u_k\right)=0,$$

$${\partial }_t\left(\rho u_i\right)+{\partial }_k\left(\rho u_i{u}_k\right)=-{\partial }_{i}p+{\partial }_j{\sigma }_{ij},$$

$${\partial }_t\left(\rho E\right)+{\partial }_k\left(\rho E{u}_k+p{u}_k\right)={\partial }_j\left({\sigma }_{ij}{u}_i\right)-{\partial }_k\left(\lambda {\partial }_{k}T\right),$$

$$\begin{array}{rl}{\partial }_t\left(\alpha ⟨\rho ⟩{\tilde{u}}_1\right)+{\partial }_x\left(\alpha ⟨\rho ⟩{\tilde{u}}_1{\tilde{u}}_1+\alpha ⟨p⟩\right)={\partial }_x\left(\alpha ⟨\sigma {⟩}_{11}\right)& \\ -{\partial }_x\left(\alpha ⟨\rho ⟩{\tilde{R}}_{11}\right)+\frac{1}{V}{\int }_{S}p{n}_{1}dS-\frac{1}{V}{\int }_S{\sigma }_{1k}{n}_{k}dS& ,\end{array}$$

$${\tau }_L=L/u_{\mathrm{s}},$$

$$C_{\mathrm{D}}=\frac{{\int }_{S_{\mathrm{p}}}(p{n}_1-{\sigma }_{1k}{n}_k)dS}{0.5⟨\rho ⟩{\tilde{u}}_1^2{A}_{\mathrm{p}}},$$

$${\mathrm{Re}}_{\mathrm{p}}=\frac{\rho u{D}_{\mathrm{p}}}{\mu }.$$

$$⟨\rho ⟩(\mathbf{x},t)=\frac{1}{N_{\mathrm{sim}}}\sum _{i=0}^{N_{\mathrm{sim}}}\frac{1}{V^i(x)}\sum _{j=0}^{N_{\mathrm{CV}}^i(x)}{\rho }^j(t)V_{\mathrm{CV}}^j,$$

$$\widetilde{u_i^{\prime \prime }{u}_j^{\prime \prime }}=\widetilde{u_i{u}_j}-{\tilde{u}}_i{\tilde{u}}_j.$$

$$F_{\mathrm{p}}=\frac{1}{0.5N_{\mathrm{p}}(x)(\rho u^2{)}_{\mathrm{IS}}{A}_{\mathrm{p}}}\sum _{j=0}^{N_{\mathrm{p}}(x)}{\int }_{S_j}(p{n}_1-{\sigma }_{1k}{n}_k)dS,$$

$$⟨\rho ⟩{\tilde{R}}_{11}={\int }_{-\infty }^{\infty }f(u^{\prime \prime })d{u}^{\prime \prime }.$$

$$⟨u_1^{\prime \prime }{u}_1^{\prime \prime }⟩={\tilde{u}}_1^2\frac{{\alpha }_{\mathrm{sep}}}{\alpha -{\alpha }_{\mathrm{sep}}},$$

$$u_{\mathrm{free}}=\tilde{u}\frac{\alpha }{\alpha -{\alpha }_{\mathrm{sep}}},$$

Full text

I Introduction

Interactions between shock waves and dense suspensions of particles occur in many engineering applications and industrial processes. Some examples are liquid and solid fuel engines (Chang and Kailasanath,, 2003; Davis et al.,, 2013; Ren et al.,, 2018), heterogeneous explosives (Zhang et al.,, 2001), blast mitigation (Chaudhuri et al.,, 2013; Milne et al.,, 2014) and dust-explosion safety measures (Sichel et al.,, 1995). Natural occurrences of shock-wave particle cloud interactions include volcanic eruptions (Chojnicki et al.,, 2006) and supernovae (Inoue et al.,, 2009; Silvia et al.,, 2012).

The process of shock wave particle cloud interaction has recently received attention in both experimental (Wagner et al.,, 2012; Theofanous et al.,, 2016; DeMauro et al.,, 2017) and numerical studies (Ling et al.,, 2012; Houim and Oran,, 2016; Theofanous and Chang,, 2017; Shallcross and Capecelatro,, 2018; Sugiyama et al.,, 2019). Numerical studies have primarily utilized two different approaches. The first type employs Eulerian-Eulerian or Eulerian-Lagrangian methods to simulate flows with a large number of particles. The model equations are formulated to account for the volume and density of the dispersed phase. There are numerous issues with such simplified dispersed flow models, such as non-hyperbolic equation sets, as discussed in Lhuillier et al., (2013); Theofanous and Chang, (2017). The second approach circumvents these problems by utilizing particle-resolved simulations. This is the method utilized in this study. Particle-resolved simulations are valuable because they provide highly resolved data, in both space and time, which enables a comprehensive analysis. The primary limitation of such simulations is that only small-scale problems are computationally feasible.

There is a considerable body of work investigating shock wave particle cloud interaction using particle-resolved inviscid simulations (Regele et al.,, 2014; Theofanous et al.,, 2018; Mehta et al.,, 2018; Mehta et al., 2019a, ; Mehta et al., 2019b, ; Sen et al.,, 2018). The number of studies that have included viscous effects in such simulations is, however, limited. Particle-resolved studies that have included viscous effects are, e.g., the two-dimensional simulations of Naiman and Knight, (2007); Hosseinzadeh-Nik et al., (2018) and the three-dimensional simulations of Vartdal and Osnes, (2018); Osnes et al., 2019a . Viscous effects have not been analyzed directly in experimental studies either. There is, however, a moderate range of particle Reynolds numbers $\mathrm{Re}_{\mathrm{p}}$ where experimental results exist. Based on the shock wave strength, the particle Reynolds numbers in Wagner et al., (2012) are approximately $2650-4200$. The study of DeMauro et al., (2017) spans $\mathrm{Re}_{\mathrm{p}}\approx 900-5700$, whereas Theofanous et al., (2016) investigated configurations in the range $\mathrm{Re}_{\mathrm{p}}\approx 7000-70000$. Since these studies vary Reynolds number simultaneously with Mach number and ratio of layer thickness to particle diameter, isolating the Reynolds-number effects is difficult.

The Reynolds number is a measure of the relative importance of inertial to viscous forces within the flow. Viscosity has an effect on the flow field even at very high particle Reynolds numbers, as long as the flow remains in the continuum regime. In particular, viscous effects are responsible for flow separation behind particles. This is in fact an important phenomenon in shock wave particle cloud interaction, as will be shown in this study.

This study considers particle Reynolds numbers in the range $500-10000$. In incompressible, single-particle flows, this particle Reynolds number range spans three different flow regimes, by their definitions in Tiwari et al., ress . These flow regimes are commonly labeled the vortex shedding regime ($270\leq\mathrm{Re}_{\mathrm{p}}\leq 800$), the separating vortex regime ($800<\mathrm{Re}_{\mathrm{p}}\leq 3000$), and the sub-critical regime ($3000<\mathrm{Re}_{\mathrm{p}}<3.4\times 10^{5}$). Sakamoto and Haniu, (1990) discusses the flow properties over this range and draws the following conclusions. In the vortex shedding regime, laminar vortices are shed periodically. Above $\mathrm{Re}_{\mathrm{p}}=800$, vortex tubes and vortex loops are formed in the vortex sheet, and can interact with the large scale laminar vortices, introducing irregularity in the vortex shedding. Above $\mathrm{Re}_{\mathrm{p}}=3000$, the vortex sheet becomes turbulent, and the turbulence intensity increases with $\mathrm{Re}_{\mathrm{p}}$. As a result of this, the Strouhal number, i.e., the non-dimensional vortex shedding frequency, decreases with Reynolds number. At $\mathrm{Re}_{\mathrm{p}}=6000$ and above, the vortex sheet is fully turbulent and the Strouhal number is approximately constant.

Compressible flows over isolated particles in this Reynolds number range have not received nearly as much attention. Nagata et al., (2016, 2018) investigated $50\leq\mathrm{Re}_{\mathrm{p}}\leq 1000$. Notably, the drag coefficient increases by more than a factor of two when the flow goes from subsonic to supersonic in this $\mathrm{Re}_{\mathrm{p}}$ range. The Mach number, $\mathrm{Ma}$, also has a stabilizing effect on the flow, and Nagata et al., (2018) found that at $\mathrm{Ma}=2$ the flow is stable to at least $\mathrm{Re}_{\mathrm{p}}=1000$.

It is not straightforward to apply these findings to shock-wave particle cloud interaction. Importantly, the incident flow becomes subsonic due to the generation of a planar, upstream propagating, reflected shock (Boiko et al.,, 1997; Wagner et al.,, 2012; DeMauro et al.,, 2017). The strength of this shock determines which flow conditions the particles are exposed to, and these flow conditions define a particle Reynolds number and a local Mach number that is not known a priori. These parameters are likely more appropriate for classification of the flow through the particle cloud after the initial shock-induced transient than the ones based on flow properties behind the incident shock. Additionally, the flow Mach number increases with downstream distance within the particle cloud (Regele et al.,, 2014; Theofanous et al.,, 2018; Osnes et al., 2019a, ) in a manner similar to that observed in Fanno flows (Emanuel and Jones,, 1968). The presence of nearby particles also affects the flow around each particle by changing the direction of the incoming flow, exposing it to particle wakes, etc.

As will be shown in this study, the particle Reynolds number affects the level of particle-scale fluctuations within the particle cloud. These fluctuations are defined here as the deviations from the volume averaged velocity, and need not be turbulent fluctuations; a laminar particle wake generates significant fluctuations. When the particle concentration increases, the impact of the wake-induced fluctuations also increases. In particular, particle wakes start to comprise a significant part of the volume within the cloud. This causes a problem for models based on volume averaging. The volume averaged velocity is shifted towards the average particle velocity, and as a result, it represents neither the flow in immediate proximity of each particle, nor the freer flow between them. In the current work, we explore what this means for the fluctuation levels and particle drag coefficients.

The current study is concerned with how the particle Reynolds number affects shock-wave particle cloud interaction. To this end, we perform particle resolved, three-dimensional, viscous simulations of the passage of a shock wave through a particle layer with volume fraction 0.1. The particle positions are assumed to be fixed over the time frame considered. We examine shock wave attenuation and reflection, mean flow fields and fluctuation statistics. We also explore the Reynolds number dependence of forces acting on the fixed particles, including both maximal values, time-averages and distributions of the drag coefficients. The results from the current study are important as a baseline for development of closures for simplified dispersed flow models, as well as the formulation of particle force models for shock-wave particle cloud interaction.

This paper is structured as follows. In section II we describe the equations governing the flow through the particle cloud, the mathematical framework used in the flow analysis and the computational method. Section III specifies the problem under investigation. The simulation results are presented and discussed in section IV, and section V contains concluding remarks.

II Governing equations and computational approach

The motion of the fluid is governed by the compressible Navier-Stokes equations, which can be written

[TABLE]

[TABLE]

[TABLE]

where $\rho(\mathbf{x},t)$ is the mass density, in which $\mathbf{x}$ is the spatial coordinate vector and $t$ denotes time, $u(\mathbf{x},t)$ is the velocity, $p(\mathbf{x},t)$ is the pressure, $\sigma_{ij}(\mathbf{x},t)=\mu(\partial_{j}u_{i}+\partial_{i}u_{j}-2\partial_{k}u_{k}\delta_{ij}/3)$ is the viscous stress tensor, $\mu(\mathbf{x},t)$ is the dynamic viscosity, $E(\mathbf{x},t)=\rho e+0.5\rho u_{k}u_{k}$ is the total energy per unit volume, $e(\mathbf{x},t)$ is the internal energy per unit mass, $\lambda(\mathbf{x},t)$ is the thermal conductivity, and $T(\mathbf{x},t)$ is the temperature. Here, $\partial_{t}$ denotes partial differentiation with respect to time and $\partial_{i,j,k}$ denotes partial differentiation with respect to space. The subscripts denote tensor components, and repeated indices imply summation from 1-3. Internal energy, pressure and density are related by the ideal gas equation of state, with $\gamma=1.4$. Temperature and internal energy are related by a constant specific heat capacity. We assume that the viscosity varies with temperature as $\mu(\mathbf{x},t)=\mu^{0}(T/T^{0})^{0.76}$, where $\mu^{0}$ is the viscosity at $T=T^{0}$. We further relate the thermal diffusivity to the viscosity through a constant Prandtl number of $0.7$.

The analysis in this work is based on the volume averaged momentum conservation equations. A discussion of volume averaging for flows containing particles can be found in, e.g., Crowe et al., (2011); Schwarzkopf and Horwitz, (2015). Since the configuration is statistically one-dimensional in space, only the streamwise momentum equation is relevant. It can be expressed

[TABLE]

where $\alpha$ is the gas volume fraction, $\tilde{R}_{ij}(\mathbf{x},t)=\widetilde{u_{i}^{\prime\prime}u_{j}^{\prime\prime}}$ is the velocity fluctuation correlation tensor, $V$ denotes the averaging volume, $n$ is the normal vector pointing into the continuous phase, and $S$ denotes the continuous phase boundary. Here, $\langle\cdot\rangle$ denotes a phase average, $\widetilde{\cdot}$ denotes a Favre average, and fluctuations from the Favre average are denoted $\cdot^{\prime\prime}$. Furthermore, we denote the fluctuating kinetic energy per unit volume, $\langle\rho\rangle\tilde{R}_{ii}/2$, by $k$, and mean kinetic energy per unit volume by $K$. The particle volume fraction is denoted by $\alpha_{\mathrm{p}}$.

In the analysis of the computational results, the temporal evolution is discussed in terms of a time-scale related to the passage of the shock through the particle layer. This time-scale is

[TABLE]

where $u_{\mathrm{s}}$ is the speed of the incident shock wave and $L$ is the length of the particle cloud. Furthermore, we use $t_{0}$ to denote the time point when the shock wave is at $x=0$.

The particle drag coefficient will be used in the discussion of particle forces. It is

[TABLE]

where $S_{\mathrm{p}}$ is the particle surface and $A_{\mathrm{p}}$ is the projected area of the particle in the flow direction.

II.1 Computational method

The present simulations are conducted using the compressible flow solver ”CharLES” (Brès et al.,, 2017), developed by Cascade Technologies. The compressible, filtered, Navier-Stokes equations are solved on an unstructured Voronoi-based grid (Bres et al.,, 2018). CharLES is a finite volume large-eddy simulation code. It computes fluxes between control volumes based on an entropy-stable modified Lax-Friedrich flux. A more detailed discussion of the entropy stable scheme can be found in Masquelet et al., (2017). A third order explicit Runge-Kutta scheme is applied to advance the solution in time.

III Set-up and computational grid

We consider the interaction between a planar shock wave and a cloud of particles at different $\mathrm{Re}_{\mathrm{p}}$. Figure 1 shows a sketch of the computational set-up. The computational domain is $-0.9L\leq x\leq 2L$, $0\leq y\leq 8\sqrt[3]{4}D_{\mathrm{p}}$ and $0\leq z\leq 8\sqrt[3]{4}D_{\mathrm{p}}$, where $D_{\mathrm{p}}$ is the particle diameter. The particles are located within $0\leq x\leq L$, and the particle volume fraction is $0.1$. In all configurations, the length of the particle layer and the particle diameter are related by $L=12\sqrt[3]{16}D_{\mathrm{p}}$. The particles are drawn randomly within this range, and we require that no particle boundary extends outside $0\leq x\leq L$. The minimum distance between each particle is set to roughly $0.05D_{\mathrm{p}}$. The shock wave position is initially $-0.1L$, and is indicated by the thick line in the figure. Its direction of propagation is from left to right.

The computational grid is an unstructured Voronoi-based grid, which is body-fitted to the particles. We use three different refinement levels, as indicated by the roman numerals in fig. 1. Refinement level I has control volumes eight times larger (by volume), than level II, and similarly for levels II and III. The length scale of the smallest control volumes (level III) is roughly $0.036D_{\mathrm{p}}$. Due to the density of the particles, most of the region $0\leq x\leq L$ contains level III control volumes. The total number of control volumes is roughly $10^{8}$. The control volume sizes are similar to those used in Osnes et al., 2019a ; Osnes et al., 2019b . Those studies investigated grid sensitivity based on the grid-convergence of the particle forces, resolution of viscous shear length scales, and convergence of spanwise velocity fluctuation correlations.

This study investigates the effect of particle Reynolds number over the range $500-10000$ on the shock-induced flow through particle clouds. The particle Reynolds number is defined as

[TABLE]

Throughout this paper, the Reynolds number is based on flow properties behind the incident shock wave, unless otherwise specified. Subscript IS will be used to denote post-shock values for the other flow variables. The Reynolds number is systematically varied by altering the viscosity, rather than varying the particle diameter or the shock wave strength. This enables the use of the same computational grid for different Reynolds numbers.

The initial conditions consist of two homogeneous spatial regions, which are separated by the shock wave at $x=-0.1L$. The pre-shock conditions are chosen to correspond to air at atmospheric conditions, i.e., $\rho^{0}=1.2048$ $\mathrm{kg/m^{3}}$, $u^{0}=0$ m/s, $p^{0}=1.01325\times 10^{5}$ Pa. The viscosity $\mu^{0}$, is varied to obtain the desired particle Reynolds number. At $\mathrm{Re}_{\mathrm{p}}=5000$, this corresponds approximately to the viscosity of air at standard atmospheric conditions ($\mu^{0}\approx 1.8\times 10^{-5}\ \mathrm{kg/ms}$). The Mach number of the incident shock wave is $2.6$, and behind it is a homogeneous flow with the post-shock state, which is $\rho_{\mathrm{IS}}=4.16$ $\mathrm{kg/m^{3}}$, $u_{\mathrm{IS}}=633$ m/s, and $p_{\mathrm{IS}}=782$ kPa. At the $x=-0.9L$ boundary, we apply constant inflow conditions, corresponding to the post-shock state. At $x=2L$, we use a zero-gradient outlet. The spanwise boundaries are periodic. The spanwise domain lengths are both roughly 13 particle diameters, and this should be sufficient to hinder any periodic artifacts. The simulations are run for a total time of $3.75\tau_{L}$.

Flow field statistics are obtained by averaging flow variables over predetermined volumes. In the current work, we utilize volumes spanning the computational domain in the $y$ and $z$ directions, with a streamwise length of $L/60\approx 0.5D_{\mathrm{p}}$.

A large number of particles is required in the spanwise directions to achieve well resolved statistics of the time-dependent flow field. Since only a limited number of particles is feasible to include in a single simulation, we perform an ensemble of simulations at each $\mathrm{Re}_{\mathrm{p}}$. The gas-phase variables are presented in terms of volume averaged quantities. The volume averages from the simulations are averaged over the ensemble of simulations, so that for example

[TABLE]

where $N_{\mathrm{sim}}$ is the number of simulations (ten in this study), $N_{\mathrm{CV}}$ is the number of control volumes within the averaging volume $V(x)$ and superscripts denote quantities belonging to simulation $i$ or control volume $j$. Fluctuation correlations are computed implicitly as the difference between the products of the averages and the averages of the products. Thus, the velocity correlations are

[TABLE]

We consider only statistics up to second order moments, and we therefore use the convergence of the fluctuating kinetic energy, $k$, to estimate the required ensemble size. We find that the tenth simulation changes the statistics of $k$ by about $0.5\%$ on average. We therefore consider an ensemble of ten simulations at each $\mathrm{Re}_{\mathrm{p}}$ sufficient for the current analysis.

IV Results

This section presents the simulation results. Section IV.1 examines the attenuation of the main shock wave and the strength of the reflected shock wave. Section IV.2 discusses mean fields and fluctuations in the gas phase. Finally, Section IV.3 presents particle drag coefficient trends and distributions.

IV.1 Shock-wave attenuation and reflection

Figure 2 shows the main shock position as a function of time. The results are plotted as the deviation from the hypothetical trajectory of a shock wave with constant Mach number 2.6. This deviation is due to shock wave attenuation caused by the presence of the particle cloud. Both shock wave reflection and particle drag contributes to the attenuation. The arrival time at a given location is delayed for lower Reynolds numbers.

If there is to be a considerable effect of $\mathrm{Re}_{\mathrm{p}}$ on the shock wave attenuation, the impulse from viscous forces must be non-negligible on a time-scale comparable to the passage of the shock wave over a particle. Viscous forces are caused by the no-slip condition at the particle surface, which sets up a strong velocity gradient around the particle. An estimation of the importance of viscous forces on the shock wave attenuation can be obtained by comparing viscous and pressure forces acting on the particles during the time when the shock interacts with them. The average ratio of the viscous to pressure forces on the particles during the shock-particle interaction is given in table 1. It can be seen that for the lower Reynolds numbers, the viscous forces are appreciable even during the initial phase.

The Mach number of the transmitted shock wave, based on the average transmitted shock wave speed over $L\leq x\leq 2L$, is also given in table 1. The Mach number is reduced by $5.3\%$ as $\mathrm{Re}_{\mathrm{p}}$ decreases from 10000 to 500. For the particle layer lengths examined here, the additional shock wave attenuation by viscous effects is only moderate, but it might become significant if longer layers are considered.

The position of the reflected shock wave as a function of time is also shown in fig. 2. Since the incident flow is supersonic, the reflected shock wave moves slowly in this reference frame. Its strength increases with decreasing Reynolds number. This is expected, since the magnitude of the particle drag increases with reduced Reynolds number. There is a larger resistance for the flow to pass through the particle layer due to the viscous forces, which builds up pressure in the flow. Therefore the reflected shock strength must increase. Note that the strength of the reflected shock increases with time.

The average Mach number of the reflected shock is given in table 1. It is based on the average speed of the reflected shock wave over the time-frame where it is located upstream of $x=-0.1L$. The lower Reynolds numbers have higher reflected shock Mach numbers, and the difference between $\mathrm{Re}_{\mathrm{p}}=10000$ and $\mathrm{Re}_{\mathrm{p}}=500$ is $4\%$.

Table 1 also shows the particle Reynolds number and the local Mach number of the flow based on the properties behind the reflected shock, assuming a constant reflected shock Mach number. For classification of the late-time flow field in the particle cloud, these two parameters are more appropriate than the Mach number behind the incident shock wave and $\mathrm{Re}_{\mathrm{p}}$ based on $\rho_{\mathrm{IS}}$, $u_{\mathrm{IS}}$ and $\mu_{\mathrm{IS}}$. The change in Reynolds number is moderate, and does not imply any qualitative changes in the flow, based on the characterization in single-particle studies. The change in local Mach number is more important because the flow behind the incident shock is supersonic while the flow behind the reflected shock is subsonic with a relatively low Mach number.

IV.2 Flow fields and fluctuations

Figure 3 shows the normalized streamwise velocity, pressure and fluctuating kinetic energy at three different time points. The fluctuating kinetic energy is normalized by the inflow kinetic energy. Similarly, fig. 4 shows the Mach number, particle Reynolds number and normalized average streamwise particle forces, i.e.,

[TABLE]

where $N_{\mathrm{p}}(x)$ is the number of particles within the control volume centered at $x$, and $S_{j}$ is the surface of particle $j$.

The fluctuating kinetic energy is the only flow quantity in figs. 3 and 4 that has a non-monotone dependence on $\mathrm{Re}_{\mathrm{p}}$. The highest level of fluctuating kinetic energy is found for $\mathrm{Re}_{\mathrm{p}}=2000$. The three highest Reynolds numbers are very similar at early time points, but the difference increases with time.

The Mach number plots reveal that the flow becomes sonic just upstream of the downstream particle cloud edge. Similar observations have been made in other studies (Theofanous and Chang,, 2017; Theofanous et al.,, 2018; Osnes et al., 2019a, ), and this behaviour is very similar to that found in Fanno flows. As discussed in Emanuel and Jones, (1968), the Mach number increases with downstream distance in Fanno flows due to the decreasing pressure. If the duct is sufficiently long, the flow chokes at the downstream exit, just like we observe here.

In fig. 4 it can be seen that the particle forces increase drastically at the downstream layer edge for the two latest time points shown. The Reynolds number increases slightly as well, but is insufficient for explaining the increase in particle forces. Instead, the higher particle forces are due to the increased Mach number.

An important part of the high speed flow through particle clouds is the particle-scale fluctuations. These are commonly neglected in Eulerian-Lagrangian and Eulerian-Eulerian simulations of shock-wave particle cloud interactions, see , e.g., Ling et al., (2012); Houim and Oran, (2016); Theofanous and Chang, (2017); Sugiyama et al., (2019). The gradient of the streamwise velocity fluctuations is an integral part of the momentum balance around the particle cloud edges. Osnes et al., 2019a showed that at the upstream particle cloud edge, the streamwise Reynolds stress gradient is of equal importance as the particle forces after the initial shock-related transient. As can be seen in fig. 3, the fluctuating kinetic energy has a very sharp gradient, and its magnitude is considerable, especially at early times. Figure 5 shows the ratio of fluctuating kinetic energy to total kinetic energy, integrated over the interior of the particle cloud, as a function of time. It can be seen that the ratio exceeds $0.25$ for all Reynolds numbers around $(t-t_{0})/\tau_{L}=0.5$, i.e., when the shock wave is slightly less than halfway through the particle cloud. When the shock wave exits the cloud, the ratio decays slowly over time. The trend with $\mathrm{Re}_{\mathrm{p}}$ is not monotone, which is also the case for the fluctuating kinetic energy, c.f. fig. 3. Regele et al., (2014) showed that the pressure field computed in a one-dimensional Eulerian-Lagrangian simulation of shock-wave particle cloud interaction corresponded to the sum of the volume averaged pressure and the Reynolds stress in the their two-dimensional simulations. The overestimation of the pressure in their one-dimensional simulations is therefore the result of neglecting the energy deposited in velocity fluctuations. For the current simulations, the ratio $k/p$ takes values up to $0.25$ when the shock wave is inside the particle cloud. It remains non-negligible even at late times, where the value is between $0.02$ and $0.14$ within the particle cloud. The ratio is lower for $\mathrm{Re}_{\mathrm{p}}=500$ than the other Reynolds numbers, which are quite similar.

Currently, there are very few models for flow field fluctuations that have been developed specifically for shock-wave particle cloud interaction. It is therefore useful to accurately identify flow phenomena that contribute to the fluctuation correlations, so that model development can be based on a fundamental understanding of the flow. In the current work, we consider only the velocity fluctuations. Figure 6 shows the contribution to the Reynolds stress by velocity fluctuations of different magnitudes at $t=3.75\tau_{L}$. The figure shows the function $f$, which is defined by

[TABLE]

$f$ is a measure of the importance of velocity fluctuations of a given amplitude with regards to the velocity fluctuation correlations. For $\mathrm{Re}_{\mathrm{p}}=500$, the peak in $f$ is very sharp, and the contribution from negative fluctuations below the peak is negligible. As the Reynolds number increases, the peak decreases. It moves towards more negative fluctuations because the mean flow velocity increases. The peak also broadens, and the contribution from strong negative fluctuations increases. The total contribution by negative fluctuations increases with Reynolds number, comprising $62.9\%$, $67.1\%$, $69.0\%$ and $69.2\%$ for Reynolds numbers $500$, $2000$, $5000$ and $10000$, respectively. The contribution from positive fluctuations can be seen to have a non-monotone dependence on Reynolds number over the range studied here. The peak value is very similar for the two lowest Reynolds numbers. When the Reynolds number is increased further, the peak value decreases to about $80\%$ and $65\%$ of the $\mathrm{Re}_{\mathrm{p}}=2000$ peak. Similar to the negative fluctuations, the curve broadens, but the effect is small.

A similar plot of $f$ was used in Osnes et al., 2019a to argue that the bulk contribution to the Reynolds stress was due to the circulation region behind each particle. Those results were for a flow configuration very similar to the $\mathrm{Re}_{\mathrm{p}}=5000$ case considered here. Figure 7 shows contours of $f$ in space and confirms that the circulation regions are those with the highest contributions to the streamwise Reynolds stress. The flow in the circulation region has an average velocity of zero and therefore differs drastically from the mean flow speed. This implies that the volume of the circulation region is an important parameter for determining the streamwise Reynolds stress magnitude. It is also evident that most of the flow volume contains velocities that are higher than the volume averaged velocity. This effect arises because the separation regions shift the average velocity away from the velocity of the flow between the particles. At $\mathrm{Re}_{\mathrm{p}}\geq 2000$ there are very few areas within the particle cloud that have a velocity equal to the volume-averaged velocity.

Figure 6 also shows the contribution to the spanwise Reynolds stresses by spanwise velocity fluctuations of different magnitudes. The function $g(|w^{\prime\prime}|)$ is the analogue to $f(u^{\prime\prime})$ for spanwise fluctuations, and is averaged over the two spanwise velocity components. Additionally, since $g(|w^{\prime\prime}|)$ should be symmetric around $w^{\prime\prime}=0$, we plot the contribution as a function of the fluctuation magnitude. The shapes of the curves are quite similar for all Reynolds numbers, but the magnitude is significantly lower for $\mathrm{Re}_{\mathrm{p}}=500$. The location of the peak value is at $0.08$ for $\mathrm{Re}_{\mathrm{p}}=500$ and around $0.135$ for the other Reynolds numbers. The trend with $\mathrm{Re}_{\mathrm{p}}$ is non-monotone, and the highest magnitude is obtained at $\mathrm{Re}_{\mathrm{p}}=2000$. This is similar to the trend for positive streamwise fluctuations, but the reduction at higher $\mathrm{Re}_{\mathrm{p}}$ is less significant for the spanwise fluctuations.

The areas that contribute most to the spanwise Reynolds stress components are shown in fig. 8. It can be seen that at $\mathrm{Re}_{\mathrm{p}}=500$, the fluctuations are located in large coherent regions. These are significantly larger than a single particle, which indicates that it is the inter-particle spacing rather than the particle size that is important for the spanwise Reynolds stress components. The same phenomenon is seen for the higher Reynolds numbers, but to a lesser degree. At $\mathrm{Re}_{\mathrm{p}}=10000$, the scale of the coherent regions only appears slightly larger than a single particle. The regions contributing to the spanwise fluctuations are not the same as those that are important for the streamwise fluctuations. The streamwise fluctuations are primarily particle wakes, whereas the spanwise fluctuations are mostly due to flow deflection around particles and streaming through channels between them. The spanwise fluctuations are therefore located slightly in front and to the side of each particle, as well as in the open spaces within the cloud.

In Figure 9 and Figure 10, the streamwise and spanwise Reynolds stress contributions are plotted at different streamwise locations within the particle layer. $\mathrm{Re}_{\mathrm{p}}=2000$ and $10000$ have very similar trends as $\mathrm{Re}_{\mathrm{p}}=5000$ and are therefore not shown. It can be seen that negative streamwise fluctuations at $\mathrm{Re}_{\mathrm{p}}=500$ have very little variation with streamwise position. On the other hand, positive fluctuations shift towards higher magnitudes as distance increases, and the curve broadens. This is also the case for $\mathrm{Re}_{\mathrm{p}}=5000$, although the shifts are slightly lower. The shift in positive direction is related to the flow acceleration, which is stronger for $\mathrm{Re}_{\mathrm{p}}=500$ than for $\mathrm{Re}_{\mathrm{p}}=5000$. For $\mathrm{Re}_{\mathrm{p}}=5000$, the shape of the curve for negative fluctuations changes with downstream distance. Smaller magnitudes become more important, presumably due to increasingly turbulent flow.

For spanwise fluctuations, we find that at $\mathrm{Re}_{\mathrm{p}}=500$, the magnitude increases with downstream distance throughout the particle cloud. Additionally, the peak location shifts towards higher fluctuation magnitudes. At $\mathrm{Re}_{\mathrm{p}}=5000$, the trend is similar up to $x/L=0.5$, but further downstream distance has little effect on the fluctuations.

Osnes et al., 2019a found that the streamwise fluctuations in the interior of the particle cloud could be represented by

[TABLE]

where $\alpha_{\mathrm{sep}}$ is an estimate of the average volume of the circulation region behind each particle. This model approximates the flow as two constant states. One is the flow between the particles and the other is the separated flow behind a particle, where the velocity is assumed to be zero. We note that this model does not predict zero velocity fluctuations when $\alpha$ goes to one, nor is it Galilean invariant, and therefore needs a modification if it is to be used in simplified dispersed flow models. A way to correct this is to replace $\tilde{u}_{1}$ by the relative velocity between the particles and the gas, and replace $\alpha_{\mathrm{sep}}$ by $\alpha_{\mathrm{p}}C$ in eq. 12, where $C$ is a constant or a function of dimensionless flow variables. Figure 11 shows $\alpha_{\mathrm{sep}}$ at the end of each simulation, estimated from the simulation results using eq. 12. $\alpha_{\mathrm{sep}}$ changes rapidly close to the cloud edges, but over the interior of the cloud it has an approximately constant value. $\mathrm{Re}_{\mathrm{p}}=500$ deviates slightly more from a constant value than the higher Reynolds numbers. The trend with $\mathrm{Re}_{\mathrm{p}}$ is monotone over this range, and lower Reynolds numbers have larger separation volumes. The trends found here agree qualitatively with observations from single sphere studies in the incompressible regime, where the separation length, i.e., the average distance from the particle to the end of the separated flow region, was found to be roughly 30% shorter at $\mathrm{Re}_{\mathrm{p}}=10000$ than at $\mathrm{Re}_{\mathrm{p}}=3700$ (Rodríguez et al.,, 2013). Indications of a similar trend were also found for compressible flows, where $\mathrm{Re}_{\mathrm{p}}=1000$ at $\mathrm{Ma}=0.8$ and $1.2$ had a shorter separation length than $\mathrm{Re}_{\mathrm{p}}=750$ (Nagata et al.,, 2018).

Figure 11 also shows the analogous function to $\alpha_{\mathrm{sep}}$ for the spanwise fluctuations. It can be seen that this function has little variation over $0.2\leq x/l\leq 0.95$. It does not follow from the derivation of eq. 12 that the spanwise fluctuations should vary this way. The results should be interpreted as confirming that the velocity fluctuations vary proportionally to the mean streamwise flow velocity throughout most of the particle layer. We note that for the spanwise fluctuations, there is a longer distance at the start of the layer where the function increases with distance. We suspect that this is because spanwise fluctuations are more related to inter-particle distance than particle diameter. If this is the case, it is not surprising that the build-up distance is longer, since the inter-particle distance is larger than the particle diameter on average. In contrast to the monotonous behavior of $\alpha_{\mathrm{sep}}$ with $\mathrm{Re}_{\mathrm{p}}$, there is a non-monotone trend for the spanwise fluctuations. We also find that for $\mathrm{Re}_{\mathrm{p}}=10000$, the function begins to increase towards the end of the layer. We suspect that this is the result of increasingly turbulent flow.

The simple relationship between the mean flow velocity and the velocity fluctuations is surprisingly robust across different flow conditions. The flow field within the particle cloud features rapidly varying velocity-, pressure- and density-magnitudes as well as varying gradients. Additionally, the Mach numbers vary from 0.3 to 1. The insensitivity to these variations is encouraging for development of velocity fluctuation correlation models that can be used as closures for simplified dispersed flow models.

IV.3 Particle drag

Accurate description of particle forces is of key importance in dispersed flow simulations. In dense particle suspensions, the local particle configuration can lead to large variations in the directions and magnitudes of the particle forces. The maximal drag for each particle occurs during its interaction with the shock wave. The behaviour of maximal drag forces has been examined in inviscid simulations by Mehta et al., (2018); Mehta et al., 2019b . We compare the peak streamwise drag forces in the current simulations to their inviscid simulations in fig. 12. The figure shows the trend in average peak streamwise drag forces with downstream distance and the distribution of peak streamwise drag forces. Note that for the results shown in fig. 12, the drag forces are presented as a drag coefficient using $\rho_{\mathrm{IS}}$ and $u_{\mathrm{IS}}$. In figs. 12, 13, 14 and 15, subscripts $x$ and $y$ are used to indicate streamwise or spanwise drag coefficients, respectively. The shock wave Mach number is $2.6$ in the current simulations and $3$ in the inviscid simulations. We thus expect slightly higher drag coefficients in our simulations. The lower Reynolds number simulations show higher peak drag forces, and this is to be expected because the contribution from viscous forces increases. It can also be seen that the slope is steeper for the lower Reynolds numbers. This is due to the stronger shock wave attenuation occurring for lower Reynolds numbers. The distribution of the deviation of peak drag forces from the linear trend with $x$ is also shown in fig. 12. The results are in excellent agreement with the inviscid simulations. It can also be seen that the distribution only has an inconsequential dependence on Reynolds number. This means that the peak drag variation is controlled by processes that are almost unaffected by viscosity, such as shock wave focusing and diffraction.

After the initial shock-induced transient, the importance of viscosity increases. It is therefore of interest to investigate how the distributions of instantaneous and time-averaged drag coefficients depend on Reynolds number. Figure 13 shows the distribution of the drag coefficients averaged over $1.875\leq t/\tau_{L}\leq 3.75$. Starting at $t=1.875\tau_{L}$ ensures that any shock-related transients have decayed. The linear trend with $x$ has been subtracted from the drag coefficients and they are also normalized with that function. For the spanwise distributions, we have normalized the drag coefficients based on the linear trend of the streamwise drag coefficients. The highest Reynolds number has a wider distribution of streamwise drag coefficients than the lower Reynolds numbers. The other three are very similar, but it seems that the distribution for $\mathrm{Re}_{\mathrm{p}}=2000$ is slightly higher and narrower than those for $\mathrm{Re}_{\mathrm{p}}=500$ and 5000. The distributions are slightly skewed. This is most pronounced for $\mathrm{Re}_{\mathrm{p}}=10000$. The spanwise distributions are more similar for all $\mathrm{Re}_{\mathrm{p}}$. The highest Reynolds number has the widest distribution, and the distributions for the two middle $\mathrm{Re}_{\mathrm{p}}$ are almost equal.

Figure 14 shows the distribution of the deviation of the instantaneous drag coefficients from a linear trend in time. For each particle, the linear trend in time is also used for normalization. Since the trend in time has been subtracted, these distributions reflect the magnitude and frequency of temporal oscillations. The trend with $\mathrm{Re}_{\mathrm{p}}$ is similar for both the streamwise and spanwise drag-coefficients, where higher $\mathrm{Re}_{\mathrm{p}}$ result in wider distributions. Note that the spanwise drag coefficient distribution is about twice as wide as that of the streamwise drag coefficient. The trends are more pronounced for the distribution of instantaneous than mean drag coefficients.

Both mean flow and instantaneous drag coefficient distributions have a dependency on particle Reynolds number. There are separate implications for drag modeling from each of these distributions. The distributions of time-averaged drag coefficients imply that drag models should incorporate terms that depend on the local particle configuration. These models can be based on the relative position of nearby particles or on the particle density. The PIEP model developed in Akiki et al., (2017) is an example of the former model type. In the latter model type, the distributions shown in fig. 13 can be used directly by imposing random fluctuations drawn from these distributions. The results here also imply that the importance of including such distributions increases with Reynolds number since the distribution is widest at $\mathrm{Re}_{\mathrm{p}}=10000$.

The distribution of instantaneous drag coefficients is a result of temporal flow fluctuations, and these have strong Reynolds number dependencies. It is likely that a part of these distributions is related to the shedding of vortices from both the particle itself and the nearby particles. The variation imposed by vortex shedding from nearby particles could be modeled based on the local particle configuration, although this is likely quite challenging. In addition to the trends with Reynolds number observed here, the distributions are likely dependent on the particle volume fraction. Investigating this effect is a topic for future studies.

Figure 15 shows the streamwise drag coefficient at $t=3.75\tau_{L}$. The figure also shows the drag coefficients where a modified flow velocity, defined as

[TABLE]

is used to calculate the drag coefficient. The modified velocity was suggested in Osnes et al., 2019a and attempts to correct for the contributions of the circulation regions behind the particles, which shift the average velocity away from the ”free” velocity within the particle layer. Also shown in the figure are the single-particle drag coefficients obtained by Nagata et al., (2016, 2018) and the drag correlation by Clift and Gauvin, (1970). For reference, the average Mach numbers within the particle cloud are in the range $0.3-0.7$, c.f. fig. 4. Here, the lower Mach number is for $\mathrm{Re}_{\mathrm{p}}=500$ and the highest Mach number is for $\mathrm{Re}_{\mathrm{p}}=10000$, where the velocity correction has also been used for the Mach number. Even with the correction, which reduces the drag coefficients, the current simulations result in significantly higher drag coefficients than predicted by models based on single-particle results. The single-particle simulation results of Nagata et al., (2016, 2018) only cover the lower two Reynolds numbers considered here. However, if the trend from these single particle simulation results can be extrapolated, it seems that the deviation is larger for low particle Reynolds numbers than high.

Other studies have found that confinement has a strong effect on the drag on single particles (Achenbach,, 1974; Akutsu,, 1977; Yeung,, 2009; Krishnan and Kaman,, 2010). For $\mathrm{Re}_{\mathrm{p}}$ up to about 300, Akutsu, (1977) found that a $30\%$ area blockage increased the drag coefficient by more than a factor of two. At $\mathrm{Re}_{\mathrm{p}}$ above $3\times 10^{4}$, Achenbach, (1974) found a $10\%$ increase for $50\%$ blockage. We expect that a similar effect occurs due to the presence of nearby particles. If a correction for this is combined with the corrections due to Mach number and separation volume, it is likely that the average drag coefficients will approach the values seen in isolated particle studies.

V Concluding remarks

This study has investigated the propagation of a Mach 2.6 shock wave through particle clouds with particle volume fraction 0.1. Particle resolved large-eddy simulations were utilized, and the particle Reynolds number, based on the incident shock state, was varied between 500 and 10000.

The results show that the shock wave attenuation increases with decreasing particle Reynolds number, but the effect is small over this Reynolds number range. The strength of the reflected shock wave increases with decreasing Reynolds number. This results in considerable changes in the flow quantities upstream of and inside the particle layer, i.e., higher pressures, lower velocities and lower Mach numbers with reduced particle Reynolds numbers.

We examined the development of fluctuations from the volume averaged flow quantities within the particle layer. The velocity fluctuation intensity has a non-monotone dependency on particle Reynolds number. The most intense velocity fluctuations were found at $\mathrm{Re}_{\mathrm{p}}=2000$. Streamwise velocity fluctuations are located in the particle wakes, and most importantly in the circulation region behind each particle. Spanwise fluctuations are located in the inter-particle regions, and consist of coherent regions that extend further than the size of single particles. As the Reynolds number increases, the sizes of these regions decrease. We find an approximate proportionality between the mean flow velocity and the flow fluctuations at late times. This relationship holds for both the streamwise and the spanwise fluctuations. Interestingly, the proportionality factor varies monotonously with $\mathrm{Re}_{\mathrm{p}}$ for the streamwise fluctuations, but not for the spanwise fluctuations.

The distributions of particle forces were examined. Maximal drag forces and their spatial variation agree well with the inviscid simulations of Mehta et al., (2018). In the late-time flow field, we find that both the variation of time-averaged particle forces and the temporal variation increases with particle Reynolds number. The trend with particle Reynolds number is more pronounced for the temporal variation than the time-averaged distribution. These distributions provide a basis for development of particle force variation models, that can be used in simplified dispersed flow models.

The particle drag coefficients within the particle cloud were found to be higher than isolated particle drag coefficients on average, and the largest deviation occurred for the lowest Reynolds number. We demonstrated that utilizing a correction to the mean flow velocity yields drag coefficients closer to the single-particle results. A combination of this correction with drag correlations that incorporate Mach number and flow blockage effects, can possibly reproduce the observed average drag coefficients in shock-wave particle cloud interactions.