Instabilities in a Shock Interaction with a Perturbed Curtain of Particles

Hydrodynamic instabilities at the interface of two materials of different densities are a critical issue in high energy density physics (HEDP). The Rayleigh–Taylor instability (RTI) occurs when a fluid accelerates another fluid of high density [1,2]. The RTI is ubiquitous in HEDP phenomena, such as high Mach number shocks and jets, radiative blast waves and radioactively driven molecular clouds, gamma-ray bursts and accreting black holes. Due to the importance in physics mentioned above, there have been many studies related to RTI from both physical and numerical simulation points of view in the literature.

It was shown [3] that when two immiscible fluids of different densities are accelerated in a direction perpendicular to their interface, this interface is stable or unstable according to whether the acceleration is directed from the heavier to the lighter fluid or vice versa. The case with gravity (with acceleration g pointing vertically downwards) is equivalent to the two fluids being accelerated vertically upwards with the same acceleration g.

A wide variety of fluid motions can be generated following the interaction of a shock wave with an interface separating two fluids of different properties. Any perturbation initially present on the interface will, in most cases, be amplified following the refraction of the shock. This class of problems is generally referred to as the Richtmyer–Meshkov instability (RMI). The basic mechanism for the amplification of perturbations at the interface is baroclinic vorticity generation resulting from misalignment of the pressure gradient of the shock and the local density gradient across the interface. The growth of perturbations soon enters into a nonlinear regime with the appearance of bubbles of light fluid rising into heavy fluid and spikes of heavy fluid falling into the light fluid. As the interface between the two fluids becomes more distorted, secondary instabilities, such as the Kelvin–Helmholtz shearing instability, develop, and a region of turbulence and mixing (usually referred to as mixing zone) ultimately forms.

The first analytical study on impulsive acceleration of an interface by a shock wave was performed by Markstein [4], who investigated the interaction of a shock wave with a density interface (flame front). His analytical result was the same as Taylor’s result for constant acceleration [2]. But the first real treatment of the impulsive acceleration of the interface by a shock wave was given by Richtmyer [3]. He studied the growth of a sinusoidally perturbed interface in the linear regime following the shock wave impingement. In developing his impulsive model, he assumed that the shock wave is not strong enough to cause perturbation velocities comparable to the speed of sound. Ten years later, Meshkov [5] confirmed Richtmyer’s prediction experimentally.

In studies of RMI an interface between two species with different density is usually considered. In detonations and explosion, a discrete phase of particles (also called dispersed phase) is often a part of the physical setting. In dispersed phase flows, particles that are materially not connected to each other form a separate phase. These include gas-particle and liquid-particle flows in which the particles constitute the dispersed phase. Few studies report on the instabilities that occur on the interface of a curtain of particles interacting with a shock. Balakrishnan et al. [6] investigate the flow field subsequent to the detonation of a spherical charge of TNT with an ambient distribution of dilute aluminum particles. Jacobs et al. [7] have studied the interaction between blast waves and clouds of particles.

In this paper, we study the interaction between a moving, normal shock and a perturbed curtain of particles. We show that the curtain exhibits instabilities that are very similar to RMI and RTI, however, as earlier research [8, 9] indicates, these instabilities are non-baroclinic, but rather driven by interphase velocities and changes in number density in the particle phase. Here and further on, ‘interphase velocity’ refers to a velocity difference between the gaseous and particle phases.

In the next section, we present the physical and numerical model. We then discuss the qualitative features of a shock interaction with a sinusoidal perturbed curtain of particles. The growth rates are compared with RMI theory in the next section. Conclusions are reserved for the last section.

We model the particle-laden flow with the particle-source-in-cell (PSIC) method [10]. The conservation equations are solved for the carrier flow in the Eulerian frame, while particles are traced in the Lagrangian frame. In the following, we shall denote the subscript p for the particle variables and f for the gas variables at the particle position.

The governing equations for the carrier flow are the inviscid two-dimensional Euler equations in Cartesian coordinates given by:

where

are the advective fluxes, and the pressure is determined by

The equation of state closes the system

Here $M=U_{r e f} / \sqrt{\gamma R T_{r e f}}$ is a reference Mach number determined with the reference velocity Uref and reference temperature $T_{\text {ref }}$. The source term $\vec{S}$ accounts for the effect of the particles on the carrier gas and will be discussed in more detail below.

The particle phase is solved by tracking the so-called point particles individually using equations for an empirically corrected Stokes flow around a spherical particle. The kinematic equation describing the particle position x_p is given by

The particle acceleration is governed by Newton’s second law forced by the drag on the particle. With particles assumed spherical, we take the drag as a combination of the Stokes drag with corrections for high Reynolds and Mach number and the pressure drag leading to the following equations governing the particle velocity:

Here, $\vec{v}_f$ is the velocity of the gas at the particle position and $\rho_{\mathrm{p}}$ is the particle density. The first term on the right hand side describes the particle acceleration resulting from the velocity difference between the particle and the carrier gas flow. $f_1$ is an empirical correction factor [11] that yields an accurate determination within 10% of measured particle acceleration for higher relative particle Reynolds numbers up to $R_f=10,000$ and relative particle Mach numbers up to $M f=\left|\vec{v}_f\right| / \sqrt{T_f}=1.2$:

The second term is the particle acceleration induced by the pressure gradient in the carrier flow at the particle position. The particle time constant is $\tau_p=R d_p^2 \rho_p / 18$, where $d_p$ is the particle diameter. This time constant is a measure for the reaction time of the particle to the changes in the carrier gas. $\mathrm{R}=\mathrm{UL} / v$ is the Reynolds number of the carrier gas flow, where L is the reference length and is the dynamic viscosity. Because R is large, we do not model viscous effects in the governing Euler equations for the gas flow. The particle temperature is mostly affected by convection. From the first law of thermodynamics and Fourier law for heat transfer, the equation for temperature is derived as

where $\operatorname{Pr}=1.4$ is the Prandtl number, taken at its typical value for air in this paper. $N u=2+$$\sqrt{R}_f P_r^{0.33}$ is the Nusselt number corrected for high Reynolds number.

Each particle generates momentum and energy that affect the carrier flow. The volume averaged summation of all these contributions gives a continuum source contribution on the momentum and energy equation in eqn 1 as

where $\vec{K}(x, y)=\vec{K}(|x-y|) / V$ is a normalized distribution function that distributes the influence of each particle onto the carrier flow. $W_m$ and $W_e$ are weight functions describing the momentum and energy contribution of one particle and are

$m_p$ is the mass of one spherical particle, which can be derived from $\mathrm{T}_p$ and $\rho_p$. $N_p$ is the total number of particles in a finite volume V.

The Langrangian algorithm traces each particle’s position, velocity and temperature over time. The carrier phase properties such as velocities and temperature required at each particle’s position are obtained via interpolation from the grid points that surround the particle to the position of the particle. To account for the back-way coupling effect of particles on the carrier phase, the momentum and energy of the particle phase forcing is averaged and coupled to the carrier phase equations through a source term. Crowe et al. [10] first proposed this Particle-Source-in-Cell (PSIC) model for gas-droplet flows. The PSIC method is a particle-mesh type algorithm where the carrier phase is resolved in a static mesh, while the particle dynamics are traced along their path in the Lagrangian frame.

Here we use the high-resolution PSIC algorithm developed by Jacobs and Don and described in [7]. The Jacobs algorithm is a high-order, Eulerian–Lagrangian algorithm based on the weighted essentially non-oscillatory (WENO) finite difference method on a structured grid for the computation of the gas phase. Higher-order ENO interpolation between the Eulerian gas PDEs and the Lagrangian particle ODES ensure that the high-order nature of the WENO scheme is preserved in the EL-WENO method. The high-order resolution scheme captures shocks sharply without diffusing essential interactions between shocks, turbulence and particles.

We simulate the interaction of a perturbed curtain with a 4% volume fraction loading of 135,000 particles with a moving planar shock as schematically summarized in Fig. 1. The particle response time and density are $\tau_p=1.7845 \times 10^3$ and $\rho_p=1200$, respectively.

Particles are initially evenly distributed in a perturbed curtain stretching a width of 0.2941 (in x-direction) and a height of 0.2 (in y-direction) with 150 x 900 particles. The function of the initial perturbation follows a cosine profile allocating the particles uniformly. The x and y locations of each particle are assigned as follows:

Here $x_{\mathrm{p}}$ and $y_{\mathrm{p}}$ are the position of the particle in x and y directions, respectively, $x_{\mathrm{o}}$ is the initial location of the curtain, $\delta x_p=\frac{a}{n_{p a r t x}}$ is the distance between particles in x-direction, where a is the width of the curtain and $n_{\text {part-x }}$ is the number of particles in x-direction, $\delta y_p=\frac{b}{n_{\text {party }}}$ the distance between particles in y-direction where b is the height of the curtain and $\eta_{\text {part-y }}$ is the number of particles in y-direction. The variable k is a counter for each row of particles ranging from 1 to $n_{\text {part-x }}$ while j is a counter for each column of particles ranging from 1 to nparty. wo is the initial amplitude of the curtain. Tcloud is the wavelength of the perturbation. $\delta y_p=\frac{b}{n_{\text {party }}}$ also identifies the coordinate of a particle row in y direction.

For simulation of the gas phase we use a fifth-order WENO scheme with a mesh of 800 by 200 uniformly spaced grid points in the x and y directions, respectively.

The computational domain spans an area of –0.5 to 6.5 in the x-direction, which enables the tracing of relevant instabilities in the carrier and dispersed phase. The computational domain’s height ranges from –0.1 to 0.1 in the y-direction. A uniform inflow boundary condition is specified at x = –0.15 and a periodic boundary condition is specified on the spanwise boundaries of the domain. The right boundary of the domain specifies outflow boundary conditions, as shown in Fig. 2. For calculation of the viscous forcing on the particles the Reynolds number is set to $\mathrm{R}_f=9.6503 \cdot 10^5$. The shock Mach number is M = 2.8.

Tables 1 and 2 summarize the parameters and initial conditions of the particle and carrierphase respectively. The initial configuration of the simulation is depicted in Fig. 3 that shows the initial density of the carrier and the initial position of the particles in the cloud. Throughout this paper the color coding of the particles (Fig. 3) is used to identify the initial location of columns of particles in the curtain.

While there are no baroclinic effects in the dispersed phase, the growth of the initial perturbation in the curtain’s geometry is very similar to that of a baroclinically driven instability. The mechanism is of course very different and we will show is closely connected to changes in the local number density of the dispersed phase and related changes in the interphase velocity between the particles and the fluid. If the number density is higher, the carrier gas flow is decelerated by the denser particle phase. A lowered carrier phase velocity reduces the forcing of the particle motion by the gas phase and hence leads to a reduction in the particle phase velocity.

When the shock interacts with the (left) farthest upstream particles at an early time of t = 1.1, the middle part of the curtain is compressed, as observed in the number density $\Phi=\sum_i^{N_p}\left(K\left(x-x_{p, t}\right)\right)$ contours in Fig. 4. The middle part of the curtain remains more compressed as compared to the side parts throughout the flow development (e.g. at t = 2.2 and t = 3.3 in Figs. 5 and 6, respectively) because of this earlier interaction with the accelerated flow behind the shock. The resulting lower particle phase velocity in streamwise direction, $u=\frac{1}{\Phi} \sum_i^{N_p}\left(K\left(x-x_{p, i}\right)\right) u_{p, i}$ in the center as compared to the sides (Fig. 4) leads to growth of the initial perturbation.

The difference in velocity magnitude in the middle and the side creates a vortex pair in the gas phase as observed by the vorticity, $\Omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ contours in Fig. 7. The vortex pair transports particles from the sides to the middle at the downstream location of the curtain. The positive particle velocity in y-direction, $\overrightarrow{\tilde{v}}$and corresponding symmetric negative velocity on the bottom (Fig. 8), leads to an increase in the number density along the centerline (Fig. 6).

In the upstream part of the curtain the $\overrightarrow{\tilde{v}}$ velocity is positive and negative on the top and bottom half of the domain, respectively. Particles are hence transported towards the outer boundaries of the domain. Because we have specified a periodic boundary condition in spanwise direction, the particles that leave the domain through the bottom boundary enter through the top boundary and vice versa at later times (t > 3), as highlighted in Fig. 6. The number density trends are directly connected to the structures in the velocity contours in Fig. 9. We observe that the centerline velocity is lower in areas where the particle number density is high. The periodic particle streaks at the sides of the domain redirect the flow and create a visible change in the velocity at those locations.

By assuming a potential flow, one can perform analysis of the acceleration driven normal instability between two fluids of different densities. Richtmyer [12] found that the growth of the perturbed interface is linear according to

where $\omega_o$ is the initial amplitude, $[u]=\dot{\mathrm{u}}_p / d t$ is the initial change in the velocity of the interface imparted by the shock wave, $A=\frac{\rho_2-\rho_1}{\rho_1+\rho_2}$ is the Atwood number and $k=\frac{2 \pi}{\lambda}$ is the wave number of the initial perturbation.

The amplitude $\omega$ grows with the time. Both light-heavy (A > 0) and heavy-light (A < 0) cases are unstable. For A = 0 and A >0 the initial amplitude is monotonically increasing. For A <0, the amplitude is initially zero, then reverses sign and grows according to equation 13.

Richtmyer linearized the compressible computations using finite difference techniques and accounted for an initial compression of the interface using post-shock values A and $u$ proposed by Richtmyer was

To apply Richtmyer’s formulas to a particle-laden flow, we determine the post-shock Atwood number based on density of the mixture, $p_{\text {mix }}$, as

The mixture density is determined as

where the volume concentration of the particles is V=0.04 and p₁ is the density of the gas before the shock wave.

Using the computations of the shock-curtain interaction, we compare the computational growth rate against this theoretical growth rate for three different values of particle density, $\rho_{\text {part }}=1200$, 2400 and 8800. Note that if the density of the particle is changed, then the particle’s diameter changes also according to

where $d_p$ is the diameter of the particle. $\tau_{\mathrm{p}}$ is the particle response time and $R_f$ is the Reynolds Number around the particle. For a fixed number of particles the volume concentration reduces with an increased particle density and a reduced particle diameter. In Table 3, we summarize the parameters for the three cases considered.

In Fig. 10, growth of the amplitude of the curtain’s perturbation is plotted in time. All cases show an initial compression of the curtain followed by a linear growth. The compression of is smallest for p_p = 1200 and increases with increasing particle density as predicted by theory.

The growth rates for the three cases compare remarkably well with theory, considering that the instability in the particle phase is non-baroclinic. The growth rates are m=0.3828, 0.3357 and 0.2991 for p_p=1200, 2400 and 8800, respectively. Comparison of these growth rates (Table 3) to the theoretical growth rates of $\frac{d \omega}{d t}=0.4527$, 0.3766 and 0.2381, respectively, shows that the theory and computation are within at most 15% The deviation is a largest for the lowest particle density that was considered here.

A study is conducted on the growth of a perturbation in a curtain of particles when a planar shock moves through it. Two-dimensional computations show that the perturbation of the particle curtain grows much like a baroclinically driven instability in a gas phase. However, the growth mechanism is different: because of an increased particle number density at the farthest upstream part of the perturbation, an interphase velocity difference develops in spanwise direction. This in turn yields a spanwise velocity difference in the particle phase that drives the growth of the perturbation.

An important finding is that the growth rate of this instability initially follows the theoretically linear baroclinic growth. In a parametric study on the effect of the particle density, we found an agreement between the theoretical equation and the simulation within 15%. This observation is also consistent with recent experimental results [9].

Our current efforts focus on a parametric study with different dependencies to further explore instability growth in the particle phase. We are also working on an experiment based on this simulation. We expect to find similarities, not only in the behavior of the instability, but also in the growth rate.