Photon Doppler Velocimetry and Simulation of Ejection of Particles from the Surface of Shock-Loaded Samples

Abstract

The results of photon Doppler velocimetry of ejecta from shock-loaded metal samples are reported. The experiments have been performed with tin and lead samples of a given thickness and a given surface roughness. The direct numerical simulation of the process of mass ejection from the surface of shock-loaded samples is performed for conditions close to experimental by the smoothed particle hydrodynamics method. The areal density and initial velocity distribution of the volume density of ejecta are determined. Using these results, we calculate the time dependence of the profile of the volume density at the expansion of the formed dust cloud to air. Applying an approach based on the transport equation for the correlation function of the scattered field, the main parameters of the velocity distribution of ejecta, areal density of ejecta, etc. are reconstructed from spectral photon Doppler velocimetry data. The experimentally observed temporal dynamics of spectra, which is caused by the drag of dust in air, is described at an appropriately chosen size dispersion of dust particles. The masses of ejecta reconstructed from experimental data are in agreement with the smoothed particle hydrodynamics results.

Appendices

Simulation of the Initial Stage of Ejection by the Contact Smoothed Particle Hydrodynamics Method

The smoothed particle hydrodynamics (SPH) method is a method for the interpolation of a function f(r) whose values are specified on a discrete set of N points of space r₁, …, r_N. Using these known values and the smoothing function (kernel) W(r, h), one can approximate the following function at a certain point r of space:

$$f({\mathbf{r}}) = \sum\limits_{j = 1}^N {{{f}_{j}}\Delta {{V}_{j}}W({\text{|}}{\mathbf{r}} - {{{\mathbf{r}}}_{j}}{\text{|}},h),} $$

(A.1)

where ΔV_j is the corresponding volume element and h is the smoothing distance. The smoothing kernel satisfies the conditions

$$\int\limits_{|{\mathbf{r}} - {\mathbf{R}}| < h}^{} {W({\text{|}}{\mathbf{r}} - {\mathbf{R}}{\text{|}},h){{d}^{3}}r} = 1,$$

(A.2)

$$\mathop {\lim }\limits_{h \to 0} W({\text{|}}{\mathbf{r}} - {\mathbf{R}}{\text{|}},h) = \delta ({\mathbf{R}}),$$

(A.3)

where δ(R) is the Dirac delta function. This representation is convenient because the gradient of the function f(r) is expressed in terms of the smoothing kernel W(r, h) calculated for different points of space:

$$\nabla f({\mathbf{r}}) = \sum\limits_{j = 1}^N {{{f}_{j}}\Delta {{V}_{j}}\nabla W({\text{|}}{\mathbf{r}} - {{{\mathbf{r}}}_{j}}{\text{|}},h)} ,$$

(A.4)

$$\nabla W({\text{|}}{\mathbf{r}} - {{{\mathbf{r}}}_{j}}{\text{|}},h) = W{\kern 1pt} '({\text{|}}{\mathbf{r}} - {{{\mathbf{r}}}_{j}}{\text{|}},h)\frac{{{\mathbf{r}} - {{{\mathbf{r}}}_{j}}}}{{{\text{|}}{\mathbf{r}} - {{{\mathbf{r}}}_{j}}{\text{|}}}}.$$

(A.5)

Thus, the derivative of the function at a given point of space is approximated by simple summation using its nearest environment in the smoothing distance.

Using this representation of the spatial dependence of the function, one can construct a set of finite-difference SPH approximations of Eqs. (5)–(7). In this work, we use a scheme with the solution of the Riemann problem at the contact boundary between a pair of particles [59], which ensures the monotonicity of functions near the discontinuity point. In this case, the system of Eqs. (5)–(7) has the form

$$\frac{{d{{\rho }_{i}}}}{{dr}} = 2{{\rho }_{i}}\sum\limits_{j = 1}^N {\frac{{{{m}_{j}}}}{{{{\rho }_{j}}}}({{{\mathbf{U}}}_{i}} - {\mathbf{U}}_{{ij}}^{*})} \cdot \nabla {{W}_{{ij}}},$$

(A.6)

$$\frac{{d{{{\mathbf{U}}}_{i}}}}{{dt}} = - \frac{2}{{{{\rho }_{i}}}}\sum\limits_{j = 1}^N {\frac{{{{m}_{j}}}}{{{{\rho }_{j}}}}P_{{ij}}^{*}\nabla {{W}_{{ij}}},} $$

(A.7)

$$\frac{{d{{e}_{i}}}}{{dr}} = \frac{2}{{{{\rho }_{i}}}}\sum\limits_{j = 1}^N {\frac{{{{m}_{j}}}}{{{{\rho }_{j}}}}P_{{ij}}^{*}} ({{{\mathbf{U}}}_{i}} - {\mathbf{U}}_{{ij}}^{*}) \cdot \nabla {{W}_{{ij}}},$$

(A.8)

where ${\mathbf{U}}_{{ij}}^{*}$ and ${\mathbf{P}}_{{ij}}^{*}$ are the contact velocity and pressure for the ith and jth particles, respectively, which are obtained by solving the one-dimensional Riemann problem along the lines connecting their centers. The solution presented in [60] is used to calculate the contact values. The Wendland-C2 polynomial [61] is used as the smoothing function.

A large number of particles (high spatial resolution) is required for SPH calculations with a necessary accuracy. The parallel computation algorithm used in this case should be optimized in the computational time. In this work, we use the dynamic division of the medium into Voronoi cells [62]; the automated redistribution of particles between these cells ensures the balance of computational load. The essence of the approach is illustrated in Fig. 17. The simulated sample consists of a large number of particles (several millions), which are processed by parallel-operating processors. Each processor is applied only to a local set of data on particles (usually, 50 000–100 000 particles per processor core), which correspond to one Voronoi cell in the simulated space. The Voronoi diagram is unambiguously determined by the positions of the centers of cells. Consequently, the variation of the positions of these centers in order to improve the balance of the computational load leads to a change in the diagram and to the corresponding redistribution of particles between operating processors. Thus, parallel decomposition can be rapidly adapted to varying conditions. The algorithm was described in detail in [62].

The representation of the medium in the form of a set of SPH particles has a number of advantages over mesh methods. Eulerian mesh schemes require control of contact surfaces and free boundaries, which is the main source of errors at simulation by these methods. Special procedures of the rearrangement of the mesh in the region of formation of a jet at ejection are needed for Lagrangian mesh schemes; consequently, a significant numerical diffusion of the momentum and energy occurs in the rearrangement region and, as a result, small-scale processes accompanying jet formation are lost. The simulation by the meshless SPH method does not require overcoming the indicated technical difficulties. For this reason, it is optimal for the simulation of various physical phenomena such as jet formation, rotational and shear flows of compressible materials, and fragmentation of the computation region into parts because of effects of destruction.

Distribution of Dust Particles Decelerating in a Gas

The initial velocity and coordinate distribution of the ejected mass varies with time because of the deceleration of dust particles in air. The distribution function of particles with a certain size d at a given time t_m can be written in the form

$$\begin{gathered} N({v},z,{{t}_{{\text{m}}}}) = \frac{{{{N}_{{\text{a}}}}}}{{{{{v}}_{{{\text{fs}}}}}}}\int {d{{{v}}_{0}}F\left( {\frac{{{{{v}}_{0}}}}{{{{{v}}_{{{\text{fs}}}}}}}} \right)} \\ \, \times \delta (z - z({{{v}}_{0}},{{t}_{{\text{m}}}}))\delta ({v} - {v}({{{v}}_{0}},{{t}_{{\text{m}}}})), \\ \end{gathered} $$

(B.1)

where N_a is the number of particles ejected from unit area, (1/${{{v}}_{{{\text{fs}}}}}$)F(${{{v}}_{0}}$/${{{v}}_{{{\text{fs}}}}}$) is the initial normalized velocity distribution of dust particles (see Eq. (2)), ${v}$(${{{v}}_{0}}$, t_m) and z(${{{v}}_{0}}$, t_m) = $\int_0^{{{t}_{{\text{m}}}}} {dt{v}}$(${{{v}}_{0}}$, t) are the velocity and coordinate of a particle with the initial velocity ${{{v}}_{0}}$ at the time t_m. The dependence ${v}$(${{{v}}_{0}}$, t_m) is determined by means of the numerical integration of Eq. (1) taking into account a change in the drag law behind the shockwave front in air [36, 37]. The distribution given by (B.1) is obtained from the solution of the one-dimensional kinetic equation, where the drag force is taken into account and collisions between particles are neglected.

After integration over ${{{v}}_{0}}$, Eq. (B.1) can be represented in the form

$$\begin{gathered} N({v},z,{{t}_{{\text{m}}}}) = \frac{{{{N}_{{\text{a}}}}}}{{{{{v}}_{{{\text{fs}}}}}}}\sum\limits_i^{} {F\left( {\frac{{{v}_{0}^{{(i)}}}}{{{{{v}}_{{{\text{fs}}}}}}}} \right)} \\ \, \times {{\left| {\frac{{dz}}{{d{{{v}}_{0}}}}({v}_{0}^{{(i)}},{{t}_{{\text{m}}}})} \right|}^{{ - 1}}}\delta ({v} - {v}({v}_{0}^{{(i)}},{{t}_{{\text{m}}}})), \\ \end{gathered} $$

(B.2)

where summation is performed over the roots ${v}_{0}^{{(i)}}$ = ${v}_{0}^{{(i)}}$(z, t_m) of the equation z = z(${{{v}}_{0}}$, t_m). The distribution N(${v}$, z, t_m) on the (${v}$, z) plane is nonzero on the ${v}$ = ${v}$(z) curve (see Fig. 11), which is obtained by excluding ${{{v}}_{0}}$ from the system of equations ${v}$ = ${v}$(${{{v}}_{0}}$, t_m) and z = z(${{{v}}_{0}}$, t_m). Each point on this curve corresponds to a certain initial velocity ${{{v}}_{0}}$. Summation in Eq. (B.2) is equivalent to summation over groups of particles corresponding to monotonic segments in the dependence ${v}$(z).

Formulas (B.1) and (B.2) describe the evolution of the distribution of particles with a given diameter. To calculate the volume density distribution of the ejected mass ρ(${v}$, z, t_m) taking into account the dispersion of sizes d of dust particles, it is necessary to multiply Eq. (B.1) or (B.2) by the mass of a single particle m_d = ρ₀πd³/6 and to integrate over d with the distribution n(d). The result can be expressed in terms of the areal density of ejecta, using the relation N_a = ρ_a/$\overline {{{m}_{{\text{d}}}}}$, where $\overline {{{m}_{{\text{d}}}}}$ is the average mass of the particle. In this work, it is assumed that the particle size distribution n(d) is independent of their initial velocity (according to [47, 49], this is valid at ${{{v}}_{0}}$/${{{v}}_{{{\text{fs}}}}}$ > 1.1, i.e., for a larger fraction of particles involved in the multiple scattering of the probe signal).

If integration is replaced by discrete summation over sizes, the expression for the volume density ρ(${v}$, z, t_m) is reduced to summation over groups of particles each corresponding to a monotonic segment of the dependence ${v}$(z) for particles with a given size. The calculations in this work were performed in the 16-group approximation. The distribution n(d) was modeled by an ensemble of particles with four different sizes; four monotonic segments can be separated on the ${v}$(z) curve for each of four groups of particles (see Fig. 11). We approximated the log-normal distribution given by Eq. (14) with the width σ = 0.5 by a set of four discrete sizes d/d_m = 0.8, 1.3, 1.8, 2.3 with respective weight factors of 0.50, 0.32, 0.13, 0.05. An additional increase in the number of discrete values d/d_m hardly affects the results of calculation of the volume density and transport scattering coefficient.

Solution of the Transport Equation for the Multigroup Description of the Ejecta Cloud

To solve the transport equation (19) in the case of scattering of waves from the ensemble of particles, which is characterized an ambiguous dependence of the velocity of particles on the z coordinate, we use the multigroup representation of the scattering coefficient

$$\langle {{\sigma }_{{{\text{tr}}}}}\exp [ - i{{k}_{0}}(\mu - \mu {\kern 1pt} '{\kern 1pt} ){v}t]\rangle $$

in the form of expansion (21) into partial contributions. The substitution of Eq. (21) transforms the transport equation (19) to the form

$$\begin{gathered} \left( {\mu \frac{\partial }{{\partial z}} + {{\sigma }_{{{\text{tr}}}}}(z) + \kappa (z)} \right)I(z,\mu ,t) \\ = \frac{1}{{4\pi }}\sum\limits_k^{} {\sigma _{{{\text{tr}}}}^{{(k)}}(z)\exp [ - i{{k}_{0}}\mu {{{v}}_{k}}(z)t]{{\Phi }_{k}}(z,t),} \\ \end{gathered} $$

(C.1)

Here, summation is performed over groups each corresponding to a monotonic segment of the z dependence of the velocity ${{{v}}_{k}}$(z) for particles with a given diameter,

$${{\Phi }_{k}}(z,t) = 2\pi \int\limits_{ - 1}^1 {d\mu \exp [i{{k}_{0}}\mu {{{v}}_{k}}(z)t]I(z,\mu ,t).} $$

(C.2)

is the partial density, and σ_tr(z), κ(z), $\sigma _{{{\text{tr}}}}^{{(k)}}$(z) and ${{{v}}_{k}}$(z) parametrically depend on the motion time t_m.

Transforming Eq. (C.1) to the integral form and integrating it over the angular variable μ, we obtain the following equation for the partial density Φ_i:

$${{\Phi }_{i}}(\zeta ,t)\, = \,\sum\limits_k^{} {\int\limits_0^{{{\tau }_{{{\text{tot}}}}}} {d\zeta {\kern 1pt} '{\kern 1pt} {{h}_{{ik}}}(\zeta ,\zeta {\kern 1pt} '{\kern 1pt} ,t){{\Phi }_{k}}(\zeta {\kern 1pt} '{\kern 1pt} ,\zeta )\, + \,\Phi _{i}^{{(0)}}(\zeta {\kern 1pt} ,t)} ,} $$

(C.3)

where the optical depth

$$\zeta = \int\limits_z^\infty {dz{\kern 1pt} '[{{\sigma }_{{{\text{tr}}}}}(z{\kern 1pt} ') + \kappa (z{\kern 1pt} ')].} $$

(C.4)

is introduced instead of the variable z. The transport scattering coefficient σ_tr(z) in Eq. (C.4) is determined by the sum of partial contributions:

$${{\sigma }_{{{\text{tr}}}}}(z) = \sum\limits_k^{} {\sigma _{{{\text{tr}}}}^{{(k)}}(z).} $$

(C.5)

The total optical thickness of the ejecta cloud is τ_tot = ζ(z = ${{{v}}_{{{\text{fs}}}}}$t_m). In Eq. (C.3), $\Phi _{i}^{{(0)}}$(ζ, t) corresponds to the partial density from a source in the form of a monochromatic plane wave incident on the medium,

$$\Phi _{i}^{{(0)}}(\zeta ,t) = \exp [ - \zeta - i{{k}_{0}}{{{v}}_{i}}(\zeta )\,t].$$

(C.6)

The kernels h_ik(ζ, ζ′, t) in the system of integral equations (C.3) are given by the expressions

$$\begin{gathered} {{h}_{{ik}}}(\zeta ,\zeta {\kern 1pt} ',t) = \frac{\Lambda }{2}\frac{{\sigma _{{{\text{tr}}}}^{{(k)}}(\zeta {\kern 1pt} ')}}{{\sigma _{{{\text{tr}}}}^{{}}(\zeta {\kern 1pt} ')}}\int\limits_{ - 1}^1 {\frac{{d\mu }}{{{\text{|}}\mu {\text{|}}}}\eta \left( {\frac{{\zeta {\kern 1pt} '\, - \zeta }}{\mu }} \right)} \\ \times \exp \left( { - \frac{{\zeta {\kern 1pt} '\, - \zeta }}{\mu } - i{{k}_{0}}\mu [{{{v}}_{k}}(\zeta {\kern 1pt} ') - {{{v}}_{i}}(\zeta )]\,t} \right), \\ \end{gathered} $$

(C.7)

where it is assumed that the single scattering albedo Λ = σ_tr/(σ_tr + κ) is the same for all groups of particles and is independent of the coordinate. The Λ values is determined only by the material of the sample and the wavelength of radiation. Calculations within the Mie theory [54] for the wavelength λ = 1.55 μm, optical constants from [55], and particle sizes 1–5 μm give Λ = 0.89 and 0.83 for Pb and Sn, respectively.

The system of integral equations (C.3) is a multigroup generalization of the well-known Milne equation [53] to the case of moving particles. Equation (C.1) was previously considered in [30] in the single-group approximation (i.e., for the case of a single-valued function ${v}$(z)). In this work, we use the 16-group approximation (the ensemble of particles with four different sizes for each of which four monotonic segments can be distinguished on the ${v}$(z) curve (see Appendix B)).

The system of equations (C.3) is numerically integrated on the grid of ζ_n and μ_m values, which are roots of Legendre polynomials. For each time t entering into the Doppler phase factors in Eqs. (C.6) and (C.7), we solve the system of corresponding algebraic equations and, thus, form a set of density values Φ_k(ζ_n, t).

The backscattered signal intensity I(μ = 1, ζ = 0, t) is calculated by the formula

$$I(t)\, = \,\frac{\Lambda }{{4\pi }}\sum\limits_k^{} {\int\limits_0^{{{\tau }_{{{\text{tot}}}}}} {d\zeta \frac{{\sigma _{{{\text{tr}}}}^{{(k)}}(\zeta )}}{{\sigma _{{{\text{tr}}}}^{{}}(\zeta )}}} } \exp [ - \zeta {\kern 1pt} - \,i{{k}_{0}}{{{v}}_{k}}(\zeta )t]{{\Phi }_{k}}(\zeta ,t).$$

(C.8)

The Doppler spectrum of the backscattered signal I(ω) is then calculated by means of the discrete Fourier transform of Eq. (C.8) in the variable t.

Result (C.8) concerns the case of the backscattering of waves from the free layer (i.e., it is assumed that the reflection of waves from the free surface is absent).

Formula (C.8) can be generalized to the case of a reflecting free surface by introducing the corresponding boundary condition (see [30]) and can be represented in the form of an expansion of the solution in successive reflections from the boundary. This expansion can be summed for the case of a diffusely reflecting surface [30]. The simplest result is obtained under the assumption that the reflection coefficient of the surface is unity and reflection is specular [28, 63]. In this case, the backscattering intensity can be expressed in terms of the solution for a free layer with a doubled optical thickness. The intensity is given by Eq. (C.8) after the substitution of Φ_k(ζ, t) calculated for the layer with the thickness 2τ_tot, the substitution of 2τ_tot for the upper limit of the integral, and the replacement of one exponential factor by the sum

$$\begin{gathered} \exp [ - \zeta - i{{k}_{0}}{{{v}}_{k}}(\zeta )t] \\ + \exp [ - 2{{\tau }_{{{\text{tot}}}}} + \zeta - i{{k}_{0}}(2{{{v}}_{{{\text{fs}}}}} - {{{v}}_{k}}(\zeta ))t]. \\ \end{gathered} $$

This result makes it possible to estimate the maximum effect of reflection of waves from the surface on the backscattered signal.

The assumption of specular reflection from the surface [28, 63] is valid for two reasons. First, the character of reflection (specular or diffuse) is insignificant for multiply scattered radiation [52, 53]. Second, the reflection coefficient of radiation from metallic samples is close to unity (in the case of normal incidence, r = 0.92 and 0.87 for Pb and Sn, respectively [55], and the angular-averaged coefficients are 0.9 and 0.85, respectively).

For experimentally implemented optical thicknesses of the ejecta cloud τ_tot ≫ 1 and Λ values corresponding to Pb and Sn particles, the effect of reflection of waves from the free surface on the backscattered signal spectrum is insignificant (even for channel 2 in experiment 1).