**Dynamics of a 2D Vibrated Model Granular Gas in Microgravity**

Yan Grasselli, Georges Bossis, Alain Meunier and Olga Volkova

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.68277

#### Abstract

We are reporting an experimental study performed on a granular gas enclosed into a 2D cell submitted to controlled external vibrations. Experiments are performed in microgravity during parabolic flights. High-speed optical tracking allows to obtain the kinematics of the particles and the determination of all inelastic parameters as well as the translational and rotational velocity distributions. The energy into the medium is injected by submitting the experimental cell to an external and controlled vibration. Two model gases are studied beads and disks; the latter being used to study the rotational part of the particle's dynamics. We report that the free cooling of a granular medium can be predicted if we consider the velocity dependence of the normal restitution coefficient and that the experimental ratio of translational versus rotational temperature decreases with the density of the medium but increases with the driving velocity of the cell. These experimental results are compared with existing theories. We also introduce a model that fairly predicts the equilibrium temperatures along the direction of vibration.

Keywords: granular, microgravity, translational temperature, rotational temperature

## 1. Introduction

Granular gases are suspensions in air of macroscopic particles whose dynamics is ruled by momentum transfer during collisions between the particles. Unlike molecular gases, these collisions are not elastic, and the dissipation resulting of each collision gives important qualitative differences, since without bringing energy to the system, the motion of the particles will quickly stop. The supply of energy can be natural, as it is the case with gravity forces during avalanches or due to a flow of fluid through a bed of particles, or artificial for instance by shaking a box containing the grains. Due to the importance of granular flows in many industries, they have been the subject of intensive research and numerous reviews [1–5].

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Gravity is, of course, a fundamental parameter, which governs the density distribution of particles with height in a sheared flow or in a vibrated container. The understanding of the complex physics of these granular systems is then complicated by the presence of the gravity. Besides, numerical simulations allow to obtain information on the dynamics of model systems of granular particles without the need to use sophisticated experiments in parabolic flights, drop tower or suborbital rocket flight but can't replace real experiments [6–9].

Two specific phenomena of the dynamics of a vibrated granular are the clustering instabilities, which occur due to the dissipation in multiple collisions between grains and the violation of the equipartition energy between each translational and rotational degree of freedom.

Recent results, obtained by molecular dynamics simulations of a box with two opposite vibrating walls and fixed side walls, have shown that in zero gravity the cluster of particles oscillate around its equilibrium position [10]. Experiments made in a parabolic flight in similar conditions with two opposite vibrating walls and with two different sizes of particles (diameter of 1 and 2 mm) were compared to simulation results. A phase diagram of clustering versus the volume fraction of each species obtained by numerical simulation was well agreeing with the experimental results and showed a segregation effect [11]. A different kind of cluster consisting of regular alignments of particles along the velocitylinesin a Couette flowwas also foundin parabolic flight experiments [12].

Concerning the temperature of a vibrated granular gas, it was shown to follow a power law: <sup>T</sup> <sup>¼</sup> CV<sup>α</sup> <sup>p</sup> , where Vp is the peak amplitude of the vibration velocity and 1 < α < 2. In a recent paper [13], it was demonstrated that the different values of the power α were related to the ratio, W, of the energy injected by the vibration mV<sup>2</sup> <sup>p</sup> to the gravitational potential energy. When W is large, that is, to say when the gravity is negligible, the value α ¼ 2 is recovered. A balance of the energy flux injected by the vibrated wall with the dissipation induced by particle-particle and particle-wall collisions allows to demonstrate that the temperature along the direction of vibration should be larger than the transverse temperature and that this ratio increases when the radial restitution coefficient decreases [43]. It was also found theoretically [38] that, even if the temperature of translational and rotational degrees of freedom were initially identical, the decrease of the translational temperature, after switching off the energy supply, was much faster than the rotational one but both of them are predicted to have a decrease in t �<sup>2</sup> as predicted by Haff's law, although at longer time, the system becomes inhomogeneous and follows a decrease with t �6=<sup>5</sup> [14]. The t �<sup>2</sup> decrease and the difference between translational and rotational energies were also obtained by numerical simulations for particles with a needle shape [15]. In simulations of the cooling of a gas of ellipsoidal grains with different aspect ratio [16], the authors have also found a t �<sup>2</sup> law in all cases, but when the aspect ratio (a/b) increases the difference between the translational and rotational temperature decreases and had totally disappeared for a/b ¼ 2, even becoming slightly larger than one for larger values of a/b. This is explained by the fact that the coupling between translational and rotational velocity is strongly increased by the shape's anisotropy. Also, they do not observe significant deviation from the Gaussian distribution for the velocity distribution. On the other hand, experiments in a rocket flight with a box having three moveable walls and containing needles of aspect ratio close to 10 [17] show a non-Gaussian distribution of the velocities in the excitation direction. In this experiment, the temperature perpendicular to the excitation direction as well as the rotational temperature was approximately two times less than the temperature in the excitation direction. In this last experiment, the positions of the particles were determined with the help of two video cameras at right angle, and the determination of the positions and orientations was done manually from each video frame. Although the physics of 2D and 3D system can differ in some aspects, the use of cells where particles are confined in 2D allows a much easier tracking of the particles and an automatic detection of their position. Such experiments made in low gravity with a determination of the trajectories of each particle are scarce but very useful to check the validity of models and of numerical simulations. This is the case for instance with the confirmation of the difference of the temperatures in the directions parallel and perpendicular to the external vibration [18] and for the comparison of the cooling time with the theoretical expression [19] where the experiment gave a time of the same order of magnitude as the theory: τexp ¼ 38 ms against τth ¼ 60 ms [20].

Two model systems made of inelastic hard spheres or disks will be used as reference models to study the dynamics. The general experimental situation is to have the particles enclosed in a vibrated box (for energy input) where the vibration parameters (amplitude and frequency) are monitored. Direct optical observations can lead to the dynamics of individual particle and to retrieve the physical data. The collisions between particles are leading the dynamics of the system through the inelastic interactions and momentum transfer. The normal and tangential restitution coefficients depend on the material of the particles but sometimes also on the impact velocity. Gravity being one of the main issues to overcome when studying a granular medium, the experimental results presented here have been performed in a low-gravity environment. Experiments were boarded in the Airbus Zero-G from Novespace (www.novespace.fr), and the results presented obtained during parabolic flights. 2D cells containing the granular particles were mounted on a vibrated device and high-speed video recording was used to register and track the motion of each individual particles.

We will first report an experimental study on the free cooling of a granular medium made of beads: focusing on the time relaxation of the energy of the medium, and then, we will present similar experiments realized with disks in order to get access to the granular temperatures for the translation and rotational part of the particle's energy.

## 2. Free cooling

Gravity is, of course, a fundamental parameter, which governs the density distribution of particles with height in a sheared flow or in a vibrated container. The understanding of the complex physics of these granular systems is then complicated by the presence of the gravity. Besides, numerical simulations allow to obtain information on the dynamics of model systems of granular particles without the need to use sophisticated experiments in parabolic flights,

Two specific phenomena of the dynamics of a vibrated granular are the clustering instabilities, which occur due to the dissipation in multiple collisions between grains and the violation of

Recent results, obtained by molecular dynamics simulations of a box with two opposite vibrating walls and fixed side walls, have shown that in zero gravity the cluster of particles oscillate around its equilibrium position [10]. Experiments made in a parabolic flight in similar conditions with two opposite vibrating walls and with two different sizes of particles (diameter of 1 and 2 mm) were compared to simulation results. A phase diagram of clustering versus the volume fraction of each species obtained by numerical simulation was well agreeing with the experimental results and showed a segregation effect [11]. A different kind of cluster consisting of regular alignments of particles along the velocitylinesin a Couette flowwas also foundin parabolic flight experiments [12]. Concerning the temperature of a vibrated granular gas, it was shown to follow a power law:

<sup>p</sup> , where Vp is the peak amplitude of the vibration velocity and 1 < α < 2. In a recent

�<sup>2</sup> as predicted by Haff's law, although at longer time, the system becomes

�6=<sup>5</sup> [14]. The t

<sup>p</sup> to the gravitational potential energy.

�<sup>2</sup> decrease and the difference

�<sup>2</sup> law in all cases, but when the

paper [13], it was demonstrated that the different values of the power α were related to the

When W is large, that is, to say when the gravity is negligible, the value α ¼ 2 is recovered. A balance of the energy flux injected by the vibrated wall with the dissipation induced by particle-particle and particle-wall collisions allows to demonstrate that the temperature along the direction of vibration should be larger than the transverse temperature and that this ratio increases when the radial restitution coefficient decreases [43]. It was also found theoretically [38] that, even if the temperature of translational and rotational degrees of freedom were initially identical, the decrease of the translational temperature, after switching off the energy supply, was much faster than the rotational one but both of them are predicted to have a

between translational and rotational energies were also obtained by numerical simulations for particles with a needle shape [15]. In simulations of the cooling of a gas of ellipsoidal grains

aspect ratio (a/b) increases the difference between the translational and rotational temperature decreases and had totally disappeared for a/b ¼ 2, even becoming slightly larger than one for larger values of a/b. This is explained by the fact that the coupling between translational and rotational velocity is strongly increased by the shape's anisotropy. Also, they do not observe significant deviation from the Gaussian distribution for the velocity distribution. On the other hand, experiments in a rocket flight with a box having three moveable walls and containing needles of aspect ratio close to 10 [17] show a non-Gaussian distribution of the velocities in the excitation direction. In this experiment, the temperature perpendicular to the excitation

ratio, W, of the energy injected by the vibration mV<sup>2</sup>

inhomogeneous and follows a decrease with t

with different aspect ratio [16], the authors have also found a t

the equipartition energy between each translational and rotational degree of freedom.

drop tower or suborbital rocket flight but can't replace real experiments [6–9].

<sup>T</sup> <sup>¼</sup> CV<sup>α</sup>

72 Granular Materials

decrease in t

The cooling of a granular gas can be experimentally investigated by considering a granular medium submitted to a continuous external energy input (generally done by submitting the medium to a controlled vibration), then removing it and observing how the medium goes back to rest. Experimental studies on granular have to deal with gravity effects and studying model particles (disks in general) on an air flow table can overcome gravitational effects. Our approach was to perform experiments in a low-gravity environment by boarding the experimental apparatus in the Airbus Zero-G from Novespace. The airplane undergoes successive parabolic flights allowing around 22 s of microgravity per parabola. The relative gravity is recorded during the flight (Figure 1) allowing to monitor the quality of microgravity environment.

Irons beads with radius a ¼ 1 mm enclosed in a 2D cell (Figure 2) have been used as model granular particles. The area fraction of the medium was φ<sup>i</sup> ¼ 19%. The cell was chosen with a circular geometry to ensure a homogeneous energy input into the system while submitted to the vibration. The walls of the cell are made of glass to cancel as much as possible electrostatic effects and to reduce undesired friction effects between the particles and the walls. The cell is mounted on a vibrating device ("Modal exciter, 100N, Bruel&Kjaer") to allow periodic (sine oscillations) external vibrations with different frequencies, ν, and amplitudes A. The maximum cell's velocities can be changed from 30 cm/s up to 250 cm/s. This experimental set up prevents us from density fluctuations found in fluidized beds or strong rolling contributions encountered when particles move over a horizontal vibrated plate. The motion of granular particles is recorded with a high-speed camera at 470 fps during 6 s in each experiment. About 3000 pictures (320 � 320 pixels) can be retrieved from each recording.

The high-speed video recording allows us to track each particle and then to get access to their positions inside the cell. From this knowledge, the dynamics of the medium can be retrieved through the velocities of each particle. The experimental processing is performed by image analysis [21]. Each particle p is tracked individually allowing to obtain the positions xp(t) and yp(t) as a function of time. It is interesting to note that from these sets of coordinates all experimental parameters required to describe the collective motion can be directly determined such as the velocities components, the normal restitution coefficient, e, but also the pair correlation function g(r) (Figure 3). The maximum of g(r) is found at the particle diameter

Figure 1. Typical behavior of the relative gravity during a parabola. These curves are used to check the quality of the microgravity environment.

Irons beads with radius a ¼ 1 mm enclosed in a 2D cell (Figure 2) have been used as model granular particles. The area fraction of the medium was φ<sup>i</sup> ¼ 19%. The cell was chosen with a circular geometry to ensure a homogeneous energy input into the system while submitted to the vibration. The walls of the cell are made of glass to cancel as much as possible electrostatic effects and to reduce undesired friction effects between the particles and the walls. The cell is mounted on a vibrating device ("Modal exciter, 100N, Bruel&Kjaer") to allow periodic (sine oscillations) external vibrations with different frequencies, ν, and amplitudes A. The maximum cell's velocities can be changed from 30 cm/s up to 250 cm/s. This experimental set up prevents us from density fluctuations found in fluidized beds or strong rolling contributions encountered when particles move over a horizontal vibrated plate. The motion of granular particles is recorded with a high-speed camera at 470 fps during 6 s in each experiment. About

The high-speed video recording allows us to track each particle and then to get access to their positions inside the cell. From this knowledge, the dynamics of the medium can be retrieved through the velocities of each particle. The experimental processing is performed by image analysis [21]. Each particle p is tracked individually allowing to obtain the positions xp(t) and yp(t) as a function of time. It is interesting to note that from these sets of coordinates all experimental parameters required to describe the collective motion can be directly determined such as the velocities components, the normal restitution coefficient, e, but also the pair correlation function g(r) (Figure 3). The maximum of g(r) is found at the particle diameter

Figure 1. Typical behavior of the relative gravity during a parabola. These curves are used to check the quality of the

microgravity environment.

74 Granular Materials

3000 pictures (320 � 320 pixels) can be retrieved from each recording.

Figure 2. Snapshots of the experimental cell. The external vibration is applied along the y-direction (direction of the normal gravity). High-speed video recording is used to track the motion of each individual particle. (a) Vibration is on: the central part of cell contains an almost constant density of particles (dashed region). (b) The external vibration has been stopped and the overall motion of particles stops due to the inelastic collisions.

Figure 3. Experimental pair correlation function g(r) retrieved from the positions of the particles. This curve is averaged over all pictures recorded and on the spatial configurations of particles in the central area of the cell.

which proves that electrostatic effects are negligible. The small non-null value of the pair correlation function observed "before" the particle diameter is due to the uncertainty in the particle's position by image treatment.

In order to study the free cooling, that is, to relate the loss of energy of the medium due to the inelastic collisions between particles, the external vibration is switched on prior the microgravity occurs. In zero-g, the particles will then fill the entire region of the cell, and the video recording is started. After few seconds, the external vibration is switched off, and we observe the return to equilibrium (particles at rest throughout the cell). It is worth noting that in the presence of the vibration, two types of different regions clearly appear in the cell: two hot (and dilute) regions at the top and bottom of the cell while a dense region exists in the center of the cell (dashed area shown in Figure 2a). This experimental configuration gives us the possibility to study a homogeneous bed of particles in contact with two hot regions responsible for the energy input. As the external vibration is cancelled, the particles continue to move freely throughout the cell and tend to come to rest rapidly because of inelastic collisions between particles inducing energy loss. For cooled granular media, the formation of dense clusters of particles is often reported in experiments but it is not clearly observed in our situation: we rather observe some alignments of particles along "wavy lines" but there is no evidence of high and low-density regions as the main part of the energy loss is supposed to occur along the normal direction between two particles. The relative low area fraction of particles is also a possible reason for this non-observation of this clustering effect. Moreover, g-jitter still exists could add a general motion of particles in a given direction. But a short

Figure 4. Volume fraction of particles in the central area of the cell as a function of the recording time (see Figure 2a). In this region, we will assume that the volume fraction of the granular medium remains constant.

time, after the vibration has been removed, we generally observe that the particles tend to stop in the center of the cell without evidence for clustering. During the recording, we have verified that the concentration in the central part of cell remains constant, and we have based all of other study on the dynamics of this area (Figure 4).

which proves that electrostatic effects are negligible. The small non-null value of the pair correlation function observed "before" the particle diameter is due to the uncertainty in the

In order to study the free cooling, that is, to relate the loss of energy of the medium due to the inelastic collisions between particles, the external vibration is switched on prior the microgravity occurs. In zero-g, the particles will then fill the entire region of the cell, and the video recording is started. After few seconds, the external vibration is switched off, and we observe the return to equilibrium (particles at rest throughout the cell). It is worth noting that in the presence of the vibration, two types of different regions clearly appear in the cell: two hot (and dilute) regions at the top and bottom of the cell while a dense region exists in the center of the cell (dashed area shown in Figure 2a). This experimental configuration gives us the possibility to study a homogeneous bed of particles in contact with two hot regions responsible for the energy input. As the external vibration is cancelled, the particles continue to move freely throughout the cell and tend to come to rest rapidly because of inelastic collisions between particles inducing energy loss. For cooled granular media, the formation of dense clusters of particles is often reported in experiments but it is not clearly observed in our situation: we rather observe some alignments of particles along "wavy lines" but there is no evidence of high and low-density regions as the main part of the energy loss is supposed to occur along the normal direction between two particles. The relative low area fraction of particles is also a possible reason for this non-observation of this clustering effect. Moreover, g-jitter still exists could add a general motion of particles in a given direction. But a short

Figure 4. Volume fraction of particles in the central area of the cell as a function of the recording time (see Figure 2a). In

this region, we will assume that the volume fraction of the granular medium remains constant.

particle's position by image treatment.

76 Granular Materials

As the behavior of the medium is governed by inelastic collisions, we have determined experimentally the normal restitution coefficient as a function of the relative normal velocities of two colliding particles. A systematic investigation of binary collisions of particles has been realized either in the presence of the external vibration or without it. From the optical tracking and the knowledge of the trajectories, we can compute the directions and the magnitudes of the velocities before, VR ! and after, V<sup>0</sup> R ! an impact between particles. By tracking these changes in the direction of motion of each particle when a nearest neighbor is present, we are able to precisely determine the binary collisions from the trajectories and so the exact position of the colliding particles. It is then possible to consider the positions of particles around the location of the collision (Figure 5) and insure that the trajectories before and after collision are linear to

Figure 5. Experimental trajectories recorded during a binary collision between particles. The circles represent the positions retrieved from optical tracking. For a better understanding, we have added on the experimental trajectories the direction of motion (arrows) before and after collision. We can precisely obtain the position of each particle at impact and the direction of the normal direction n.

qualify this collision for processing. The direction of the normal direction at contact is then possible, and the restitution coefficient is obtained from e ¼ �j n ! � <sup>V</sup><sup>0</sup> R �! j=j n ! � VR<sup>j</sup> �! .

The behavior of the restitution coefficient vs. the normal relative impact velocity is presented in Figure 6. For high relative velocities, we obtain a value of the restitution coefficient of 0.9 (typical value for steel beads). The most amazing observation is that the restitution coefficient shows a sharp decrease for "small" impact velocity. This is a situation encountered in the case of wet particles when e ¼ 0 for Stokes number St ¼ ðmVi=6πηa2Þ smaller than a critical value [22, 23]. This comes from the viscous dissipation but for dry particles, most experimental investigations report that the restitution coefficient increases for increasing impact velocities. Most of the experiments are made in labs (i.e., with gravity present) with impact velocities larger than 1 m/s (for a height h ¼ 5 cm the impact velocity of a bead on the plane: Vi <sup>¼</sup> ffiffiffiffiffiffiffi 2gh p is already 1 m/s). The restitution coefficient between two dry beads attached by strands in a pendulum device has also been studied [24], where it has also been found that a low value of the restitution coefficient was reported at low velocities (typically below 20 cm/s), and such behavior is well confirmed in our study without experimental drawbacks.

Figure 6. Experimental dependence of the normal restitution coefficient, e, as a function of the relative normal impact velocity, obtained from the experimental trajectories of the particles. A clear decrease at low impact velocities is observed.

To investigate the free cooling more precisely, a typical record on how the energy decreases once the external energy input has been cancelled is presented in Figure 7. This behavior is monitored through the average velocities of the particles in the central area of the cell. One can observe the rapid decay of the average velocity. The non-zero value measured for "long times" comes from the small gravity fluctuations occurring during the parabolic flight.

qualify this collision for processing. The direction of the normal direction at contact is then

The behavior of the restitution coefficient vs. the normal relative impact velocity is presented in Figure 6. For high relative velocities, we obtain a value of the restitution coefficient of 0.9 (typical value for steel beads). The most amazing observation is that the restitution coefficient shows a sharp decrease for "small" impact velocity. This is a situation encountered in the case of wet particles when e ¼ 0 for Stokes number St ¼ ðmVi=6πηa2Þ smaller than a critical value [22, 23]. This comes from the viscous dissipation but for dry particles, most experimental investigations report that the restitution coefficient increases for increasing impact velocities. Most of the experiments are made in labs (i.e., with gravity present) with impact velocities larger than 1 m/s (for a height h ¼ 5 cm the impact velocity of a bead on the

attached by strands in a pendulum device has also been studied [24], where it has also been found that a low value of the restitution coefficient was reported at low velocities (typically below 20 cm/s), and such behavior is well confirmed in our study without experimental

Figure 6. Experimental dependence of the normal restitution coefficient, e, as a function of the relative normal impact velocity, obtained from the experimental trajectories of the particles. A clear decrease at low impact velocities is observed.

2gh p is already 1 m/s). The restitution coefficient between two dry beads

! � <sup>V</sup><sup>0</sup> R �! j=j n ! � VR<sup>j</sup> �! .

possible, and the restitution coefficient is obtained from e ¼ �j n

plane: Vi <sup>¼</sup> ffiffiffiffiffiffiffi

drawbacks.

78 Granular Materials

We can first consider the energy decay assuming a constant restitution coefficient (typically e ¼ 0.9 for stainless steel beads). The time dependence of the energy is predicted to behave as EðτÞ ¼ 1=ð1 þ τÞ 2 , <sup>τ</sup> ¼ ð<sup>1</sup> � <sup>e</sup><sup>2</sup>Þt=tE, where tE is the Enskog time [25]: tE ¼ ð<sup>a</sup> ffiffiffi <sup>π</sup> <sup>p</sup> <sup>Þ</sup><sup>=</sup> � ffiffiffi 2 <sup>p</sup> <sup>φ</sup>V0gðr<sup>Þ</sup> � . V<sup>0</sup> is the initial average velocity in the medium. With our experimental set up, we can determine experimentally all the parameters involved. A quantitative comparison with experiments is presented in Figure 8 (squared symbols) for a cell velocity of 75 cm/s and with the following experimental values: φ ¼ 0:297 � 0:027, V<sup>0</sup> ¼ ð0:11 � 0:01Þ m=s and gðr ¼ 2aÞ ¼ 2:23 � 0:02. The drop of energy found from experiments is much faster than the theoretical one while considering a constant restitution coefficient. A possible explanation may result from the friction of the particles on walls of the cell, introducing an additional loss of energy. Nevertheless, a precise and systematic analysis of the trajectory of a single particle after the vibration has been cut off shows a linear motion at constant speed between two collisions of particles; we may then reject this possibility. The Enskog collision time is tE ¼ 0:0172 s. In order to check this value, we

Figure 7. Average translational velocity as a function of time obtained in the central region of the cell. The vibration is removed during the microgravity period. A clear decrease of the energy can be observed (max cell velocity of 74.6cm/s).

Figure 8. Experiments (plain curve) and theory of the energy decrease. Squares: theory including a constant restitution coefficient. Circles: theory considering the rotational energy. Dashed line: theory focusing only the translational energy but including the velocity dependence of the restitution coefficient (see Figure 6).

have performed a large statistics on our experiments to obtain the average time interval separating two consecutive collisions in the central part of the cell. We found an average time interval of (0.0127 � 0.0021) s by direct measurements in rather good agreement with the theoretical value.

So, the possible discrepancy between energy decays observed experimentally and the predicted one by theory may come from the rotational kinetic energy which also dissipates a part of the energy through the surface roughness of the particles [26]. To get a complete description of the binary collision, we have introduced a tangential restitution coefficient, β, in order to characterize the contribution of the rotation of the particles. The time dependence of the translational and rotational energy is obtained from coupled differential equations (Eq. (15) in Ref. [26]): note that the parameters considered in this description are all retrieved from experiments, except β. This system of equations was solved numerically. We have introduced our experimental results for the inelastic parameters of particles and setting β ¼ 0.1 (Figure 8, black circle—if we cancel the rotation, that is, β ¼ �1, we recover the situation of a constant normal restitution coefficient). We see that the energy decreases more rapidly but it seems that the rotational kinetic energy has limited impact whatever the value of the tangential restitution coefficient and it is still not enough to represent the experimental behavior.

To improve the agreement between theory and experiments, we may consider the velocity dependence of the restitution coefficient. We can express the rate of decrease of the translational kinetic energy T like dT=dt ¼ �ncð1 � e2ÞT with nc the rate of binary collisions. In 2D, nc ¼ � 2VφgðrÞ � =ðπaÞ with V, the average velocity. Moreover, introducing the normalized energy E ¼ T=T<sup>0</sup> and the velocity ratio V=V<sup>0</sup> ¼ T=T0, the rate of decrease of energy can be rewritten in the form:

$$\frac{dE}{dt} = -\frac{g(r)\phi V\_0}{a} \sqrt{\frac{2}{\pi}} \left(1 - e(E)^2\right) E^{3/2} \tag{1}$$

But now e is assumed to depend on the normal relative velocity. We may assume that the average relative velocity is of the same order as the average velocity; then from Figure 6, we can obtain the following trend <sup>e</sup>ðEÞ ¼ <sup>0</sup>:<sup>82</sup> � <sup>0</sup>:5e�2:5<sup>E</sup>. Substituting this last relation in Eq. (1) and solving it numerically gives the behavior presented in Figure 8 (dashed line). Compared to the case including the rotation, we observe a more pronounced decrease of the energy with time. This is understandable since we observed that during cooling the restitution coefficient decreases, increasing the loss of energy. Equation (1) is obtained from a rough approximation based on the average velocity, while the probability distribution may be considered. Moreover, the tangential coefficient may probably also be dependent on the relative angular velocities of colliding particles.

This first approach on the dynamics of a granular medium shows interesting results but as stated before, the analysis is not complete due to the lack of consideration on rotational effects. With beads and our experimental setup, accessing these data is not possible. We will then introduce in the next part recent experimental investigations based on the same principle but replacing beads by disks in order to obtain a complete description of the dynamical behavior.

## 3. Translational and rotational temperatures

have performed a large statistics on our experiments to obtain the average time interval separating two consecutive collisions in the central part of the cell. We found an average time interval of (0.0127 � 0.0021) s by direct measurements in rather good agreement with the theoretical value. So, the possible discrepancy between energy decays observed experimentally and the predicted one by theory may come from the rotational kinetic energy which also dissipates a part of the energy through the surface roughness of the particles [26]. To get a complete description of the binary collision, we have introduced a tangential restitution coefficient, β, in order to characterize the contribution of the rotation of the particles. The time dependence of the translational and rotational energy is obtained from coupled differential equations (Eq. (15) in Ref. [26]): note that the parameters considered in this description are all retrieved from experiments, except β. This system of equations was solved numerically. We have introduced our experimental results for the inelastic parameters of particles and setting β ¼ 0.1 (Figure 8, black circle—if we cancel the rotation, that is, β ¼ �1, we recover the situation of a constant normal restitution coefficient). We see that the energy decreases more rapidly but it seems that the rotational kinetic energy has limited impact whatever the value of the tangential restitution coefficient and it is still not enough to represent the experimental

Figure 8. Experiments (plain curve) and theory of the energy decrease. Squares: theory including a constant restitution coefficient. Circles: theory considering the rotational energy. Dashed line: theory focusing only the translational energy

but including the velocity dependence of the restitution coefficient (see Figure 6).

To improve the agreement between theory and experiments, we may consider the velocity dependence of the restitution coefficient. We can express the rate of decrease of the

behavior.

80 Granular Materials

In this part, our aim is to provide experimental data both for the normal and tangential restitution coefficients and for the different quantities related to the rotational and translational degrees of freedom such as the distribution functions and the rotational and translational temperatures. As introduced previously, the kinematics of granular particles submitted to a vertical vibration will still be used in a low-gravity environment. We shall particularly focus on the ratio between rotational and translational temperatures. Several other groups have already presented experimental results on granular flow under such conditions [27–30], but to our knowledge, this is the first experimental study giving access to rotational and translational velocities and so the corresponding temperatures.

We have used the same 2D cell from the previous part by now with a rectangular shape made in Duralumin. The cell has a height Ly ¼ 6.8 cm and a width Lx ¼ 6 cm. The particles studied were brass disks having a diameter <sup>σ</sup> <sup>¼</sup> 6 mm (radius a <sup>¼</sup> 3 mm) and mass m <sup>¼</sup> 4.6�10�<sup>4</sup> kg. The initial area fraction φ of the medium is obtained from the number of disks N (12 or 24, i.e., area fractions of 8.3 or 16.6%). For this experimental study, we chose this rectangular shape in order to easily monitor the energy input into the medium. The external vibration is still periodic (sine oscillations) with different frequencies, ν, and amplitudes A, and it is still applied along the y-direction (which is the direction of normal gravity). In order to increase the precisions of experimental data, the video recording is performed during the whole parabola with a higher frame rate (i.e., 900 fps) and higher image resolution (720 � 720 Pixels); each record gives access to around 22,000 pictures per parabola. To reduce the effect of friction between the disks and the glass plates of the cell, we have added three small steel beads on each side of a disk. In addition, it also reduces the tilting of the disks when the external vibration is on. The key question being now the rotational aspects, each disk is pierced with two small holes, symmetric about the center of the disk and video observations are realized by light transmission (Figure 9).

This set up grants us with images having a high contrast and quality. The position of the disk is retrieved from the tracking of the two holes for each disk as a function of time. The barycenter of the holes gives access to the x- and y-position of the disk and by following the variations of these positions as a function of time to the components of the velocity vx(t) and vy(t). Nevertheless, by computing the time dependence of the angle θðtÞ (Figure 10) (obtained from the angular position of the holes about the horizontal), allows to retrieved the angular velocity <sup>θ</sup>\_ðtÞ.

Figure 9. Picture of the medium recorded in microgravity when being submitted to the external vibration (along the ydirection). Optical observations are performed from light transmission. The two holes can be clearly identified. A side and top sketch of one disk is also shown.

Dynamics of a 2D Vibrated Model Granular Gas in Microgravity http://dx.doi.org/10.5772/intechopen.68277 83

applied along the y-direction (which is the direction of normal gravity). In order to increase the precisions of experimental data, the video recording is performed during the whole parabola with a higher frame rate (i.e., 900 fps) and higher image resolution (720 � 720 Pixels); each record gives access to around 22,000 pictures per parabola. To reduce the effect of friction between the disks and the glass plates of the cell, we have added three small steel beads on each side of a disk. In addition, it also reduces the tilting of the disks when the external vibration is on. The key question being now the rotational aspects, each disk is pierced with two small holes, symmetric about the center of the disk and video observations are realized by

This set up grants us with images having a high contrast and quality. The position of the disk is retrieved from the tracking of the two holes for each disk as a function of time. The barycenter of the holes gives access to the x- and y-position of the disk and by following the variations of these positions as a function of time to the components of the velocity vx(t) and vy(t). Nevertheless, by computing the time dependence of the angle θðtÞ (Figure 10) (obtained from the angular position of the holes about the horizontal), allows to retrieved the angular velocity

Figure 9. Picture of the medium recorded in microgravity when being submitted to the external vibration (along the ydirection). Optical observations are performed from light transmission. The two holes can be clearly identified. A side and

light transmission (Figure 9).

top sketch of one disk is also shown.

<sup>θ</sup>\_ðtÞ.

82 Granular Materials

Figure 10. Disks used as the granular particles. Light transmission allows very high contrast pictures. The optical tracking of the two holes of a single disk permits to compute the orientation angle of the disk as a function of time.

The orientation angles of the disks can be fully determined from 0 to 360 degrees. A typical experimental record of θ(t) is presented in Figure 11. On such record, a sharp change in the direction of rotation (positive or negative slope) or a significant change in the slope is the proof that a collision occurs with another particle. On the contrary, when the particle experiences no collision (e.g., time larger than 5 s in Figure 11), the angular velocity remains quite constant, indicating the absence of friction with the lateral walls. During a parabolic flight, the aircraft is subjected to g-jitter along the three directions (Figure 1). Experiments were submitted to g-jitter with period of fluctuations of about 1 s and amplitudes of about 0.01 g (Figure 1). Although these fluctuations may play a role during the collision of the particles with the moving walls of the cell, they have limited impact on the motion of particles located in the central region of the cell where experimental data were retrieved. Moreover, a systematic

Figure 11. Experimental record of the angle of orientation, θ(t) of a disk when both microgravity and external vibration are present. Collisions can be clearly identified from a change of rotation or value of the angular velocity (i.e., The slope).

analysis of inelastic parameters (normal e, and tangential β, restitution coefficients) was achieved by analyzing the trajectory of each disk.

The collision between disks is processed as we did for the beads in the previous part except that now, the relative velocity includes the rotational part VR ! ¼ v<sup>1</sup> ! � v<sup>2</sup> ! �aðθ\_ 1 ! <sup>þ</sup> <sup>θ</sup>\_ 2 ! Þ� n ! where the subscripts 1 and 2 stand for the two colliding particles at a given time. Again, the normal restitution coefficient is obtained through e ¼ �j n ! � <sup>V</sup><sup>0</sup> R ! j=j n ! � VR<sup>j</sup> ! and the tangential restitution coefficient by β ¼ �j n ! � <sup>V</sup><sup>0</sup> R ! j=j n ! � VR<sup>j</sup> ! . Another way to express the tangential restitution coefficient is to introduce the angle γ between n ! and VR ! (Figure 12), we have for disks the relation [31]: 1 þ β ¼ �3ð1 þ eÞμ cotðγÞ. The initial slope of β versus cot (γ) allows the computation of μ, the friction coefficient.

We have obtained experimentally an average value of e ¼ 0:64 � 0:03. Although it is sometimes noticed in such situation [24, 32], we did not observe for the disks we used any clear dependence of e on the relative impact velocity. The experimental determination of the restitution coefficient α, between a particle and the walls of the cell reports a value α ¼ 0:71 � 0:04. We were also able to determine the behavior of the experimental tangential restitution coefficient as a function of cotðγÞ.The results are presented in Figure 13. From the initial slope, one can compute an average value for the friction coefficient when particles are at contact: μ ¼ 0:14 � 0:01. As most of the binary collisions are head-on collisions (due to shape of the experimental cell), we have decided to take an average value of the tangential restitution β ¼ 0:7 � 0:05.

The density and local velocity profiles of particles within the cell can be determined again from the positions of particles. In Figure 14, we have plotted the profiles of the x- and y- components of the disks' velocity. The area fraction of particles in the center of the cell is almost twice the initial one while close to the top and bottom walls, the value found is smaller. This is directly related to the inelastic nature of collisions which form clusters of particles [33, 34]. For this

Figure 12. Sketch of a binary collision. v<sup>1</sup> ! and v<sup>2</sup> ! , and θ\_ <sup>1</sup> ! and θ\_ <sup>2</sup> ! represent, respectively, the linear and rotational velocities of the particles before and after impact. VR ! is the relative velocity and n the normal direction at collision. The impact angle γ is defined from n to VR ! .

analysis of inelastic parameters (normal e, and tangential β, restitution coefficients) was

The collision between disks is processed as we did for the beads in the previous part except

where the subscripts 1 and 2 stand for the two colliding particles at a given time. Again, the

disks the relation [31]: 1 þ β ¼ �3ð1 þ eÞμ cotðγÞ. The initial slope of β versus cot (γ) allows

We have obtained experimentally an average value of e ¼ 0:64 � 0:03. Although it is sometimes noticed in such situation [24, 32], we did not observe for the disks we used any clear dependence of e on the relative impact velocity. The experimental determination of the restitution coefficient α, between a particle and the walls of the cell reports a value α ¼ 0:71 � 0:04. We were also able to determine the behavior of the experimental tangential restitution coefficient as a function of cotðγÞ.The results are presented in Figure 13. From the initial slope, one can compute an average value for the friction coefficient when particles are at contact: μ ¼ 0:14 � 0:01. As most of the binary collisions are head-on collisions (due to shape of the experimental cell), we have decided to take an average value of the tangential restitution

The density and local velocity profiles of particles within the cell can be determined again from the positions of particles. In Figure 14, we have plotted the profiles of the x- and y- components of the disks' velocity. The area fraction of particles in the center of the cell is almost twice the initial one while close to the top and bottom walls, the value found is smaller. This is directly related to the inelastic nature of collisions which form clusters of particles [33, 34]. For this

! � VR<sup>j</sup> !

! ¼ v<sup>1</sup> !

! and VR !

! � <sup>V</sup><sup>0</sup> R ! j=j n ! � VR<sup>j</sup> !

� v<sup>2</sup> !

. Another way to express the tangential

represent, respectively, the linear and rotational

is the relative velocity and n the normal direction at collision. The

�aðθ\_ 1 ! <sup>þ</sup> <sup>θ</sup>\_ 2 ! Þ� n !

and the tangential

(Figure 12), we have for

achieved by analyzing the trajectory of each disk.

restitution coefficient by β ¼ �j n

Figure 12. Sketch of a binary collision. v<sup>1</sup>

impact angle γ is defined from n to VR

velocities of the particles before and after impact. VR

! and v<sup>2</sup> ! , and θ\_ <sup>1</sup> ! and θ\_ <sup>2</sup> !

! . !

β ¼ 0:7 � 0:05.

84 Granular Materials

the computation of μ, the friction coefficient.

that now, the relative velocity includes the rotational part VR

! � <sup>V</sup><sup>0</sup> R ! j=j n

normal restitution coefficient is obtained through e ¼ �j n

restitution coefficient is to introduce the angle γ between n

Figure 13. Experimental behavior of the tangential restitution β as a function of cotðγÞ, where γ is the angle between the normal direction at contact and the direction of the relative velocity. The plain line is a linear regression used to compute the value of the friction coefficient when particles are at contact.

experimental study, the cell may again be divided into two different and well-identified regions [35]: a central one that we will refer as the "cold" area and the ones close to the top and bottom walls (where the energy is injected into the medium), referred to "hot" areas. In the following, subscripts H and C will be used, respectively, to identify the "hot" and "cold" regions of the cell. Considering all experiments, we have noticed first that the values found to characterize the "hot" zone by the height hH were not related to the amplitude of vibration as one could expect, and second that an average value hH ¼ 9 mm ≈ 3.5a was acceptable in our experimental situations. Moreover, computation of the mean free path of particles in the "cold" zone gives a distance of about 12a (when using 12 disks) and 6a (for situations with 24 disks). These statements are found by considering the density of the "cold" region where typical values are found to be 13% (12 disks) and 30% (24 disks). We may conclude that the behavior of particles in the "cold" area is mainly governed by particle-particle collisions. We do not meet situations in which particles are moving through the bulk without being struck by another particle. Last, checking experimentally how the density profiles evolve with time does not present low-frequency oscillations like the situation reported in Ref. [36].

The temperatures of the granular medium are retrieved from the velocities of particles and from both contributions: the translation Ttr <sup>¼</sup> mV<sup>2</sup> <sup>=</sup>2, and the rotation Trot <sup>¼</sup> <sup>I</sup>ω<sup>2</sup>=2 (I the moment of inertia). In steady state, the temperature is calculated from a balance between two opposite fluxes: the energy brought to the medium by the particles in the "hot" areas of the cell and the dissipation in the bulk (i.e., the "cold" region). All experimental data presented below were retrieved from the "cold" region. To obtain reliable values, the whole set of the 20,000 images recorded are treated for each experimental value reported in this paper. We may also

Figure 14. Velocity distributions of the component along the direction of vibration (y-direction) and transverse to it (xdirection). The experimental curves are drawn with plain lines. The dashed lines correspond to a Gaussian plot with the average velocity determined experimentally.

consider an undesired effect of g-jitter arising in parabolic flights but we will take it into account when comparing experimental results with theories. An average number NH of particles present in the "hot" regions can be obtained directly from the density profiles at any time.

Figure 15. Typical angular velocity distribution of the particles (experiment: plain curve). The dashed line corresponds to the mathematical plotting of a Maxwell distribution, which includes the average angular velocity determined experimentally.

Typical experimental distributions for translation and rotational velocities are shown, respectively, in Figures 14 and 15: A Gaussian behavior can be observed. The dashed line on the figures represents the plots of the Gaussian distribution in which the experimental values of the squared velocities have been introduced, and one can notice the good agreement.

Due to the rectangular shape of the experimental cell used and to the relatively low area fraction, the main contribution to the temperature was expected to be found along the direction of the external vibration (the y� direction). The temperature ratios Ty=Tx and Ttr=Trot with Ttr ¼ ðTx þ TyÞ=2, in terms of the maximum cell's velocity Aω (ω ¼ 2πn) for the two area fractions used, have been calculated. One may note that the ratio Ttr=Trot is not drastically affected if one considers only Ty as the only contribution to the energy. For the fraction area of 16.6%, the ratio Ttr=Trot is ranging from about 4 to 10, respectively, for maximum cell velocity from 20 to 40 cm/s, while Ty=Tx goes from 2 to 4 at maximum, while for the lowest fraction area (8.3%), Ttr=Trot is ranging from 11 up to 24 and Ty=Tx from 5 to 7 only under the same conditions of maximum velocities. We shall analyze these experimental results by focusing first on the ratio Ty=Tx which is clearly dependent on the area fraction of the medium (the larger the ratio, the smaller the area fraction). Without any surprise, the temperatures obtained

consider an undesired effect of g-jitter arising in parabolic flights but we will take it into account when comparing experimental results with theories. An average number NH of particles present in the "hot" regions can be obtained directly from the density profiles at any time.

Figure 14. Velocity distributions of the component along the direction of vibration (y-direction) and transverse to it (xdirection). The experimental curves are drawn with plain lines. The dashed lines correspond to a Gaussian plot with the

average velocity determined experimentally.

86 Granular Materials

along the direction of the vibration are always larger than the ones in the transverse direction. This comes from the fact that the main part of energy injected is along the y-direction, and the relatively low area fraction does not allow the redistribution of this energy toward the perpendicular direction. At low area fraction, the particles can move easily and the y-direction drives the general motion. On the other hand, we also observe a net increase of the ratio Ty=Tx with the driving velocity of the cell, but less pronounced for the lower area fraction. However, the driving velocity is not the only parameter of the problem and the amplitude can also play a role. For example, the ratio Ty=Tx ¼ 1:98 found was obtained for the smallest amplitude (A ¼ 0:556 mm) and the largest frequency (60 Hz). Under the conditions of high frequency but small amplitude, we clearly see a concentration of particles in the central area of the cell, and consequently, the corresponding energy input is small. It then explains the low ratio found in such experiment, to be compared to a similar value of Aω ¼ 0:22 m/s but with much larger amplitude (A ¼ 2:3 mm). We have set the frequency scale between 10 Hz and 30 Hz, and the corresponding amplitudes of vibration used are large enough to avoid the aggregation of particles in the center of the cell. Last, we can clearly identify the "cold" area straight from the density profiles. The second result is related to the ratio Ttr=Trot which obviously increases with Aω and which is also dependent on the area fraction. The translational temperature is quite one order of magnitude larger than the rotational temperature. Because almost all collisions between particles are quite, head-on as reflected by the high value of Ty=Tx, can explain why the transfer from translational to rotational energy is rather weak, mainly at the lowest area fraction.

As a first step to describe the experimental behavior on granular temperatures, we can use existing theoretical models using a mean-field theory [37]. In this description, the rate of change of the temperature of a granular medium is determined through two coupled equations

$$\begin{cases} \frac{dT\_{tr}}{dt} = J\_{dr} + G[-AT\_{tr}^{3/2} + BT\_{tr}^{1/2}T\_{rot}]\\ \qquad \frac{dT\_{rot}}{dt} = 2G[B'T\_{tr}^{3/2} - CT\_{tr}^{1/2}T\_{rot}] \end{cases} \tag{2}$$

where Ttr and Trot represent, respectively, the translational and rotational temperatures, and <sup>G</sup> <sup>¼</sup> <sup>16</sup> <sup>σ</sup> ffiffiffiffiffi <sup>π</sup><sup>m</sup> <sup>p</sup> ϕg2ðϕÞ is related to the collision rate between particles; g2ðϕÞ being the pair correlation function at contact. In two dimensions, g2ðϕÞ¼ð1 � 7ϕ=16Þ=ð1 � ϕÞ2. A, B, B<sup>0</sup> , and C are constants, which depend only on the inelastic properties of the particles [38]. Jdr is the energy flux input into the medium and a homogeneous input of energy into the medium is assumed. These four constants are positive so that the minus signs are related to the loss of energy during the collision of particles. We state that the driving energy is acting on the translational temperature because of the preferred collisions with normal incidence. We may note that when the rotation mainly governs the behavior of the granular, Jdr is included in the second equation of (2) [39, 40].

Several inelastic modelizations were proposed by Herbst et al. [38] going from a simple consideration of a constant tangential restitution coefficient up to more complex ones where the tangential restitution depending on γ<sup>12</sup> (the contact angle obtained neglecting the rotational velocities) or on the real contact angle γ. From the second equation of Eq. (2), the energy ratio Ttr=Trot can be obtained considering the medium in steady state dTrot=dt ¼ 0, allowing to get the relation Ttr=Trot ¼ C=B<sup>0</sup> . In this equilibrium regime, dTtr=dt ¼ 0 and replacing Trot in the first equation of Eq. (1) also give Ttr ¼ ½CJdr=GðAC � BB<sup>0</sup> Þ�<sup>2</sup>=<sup>3</sup> . Depending on the model used, the expressions of the constant C and B<sup>0</sup> can be evaluated and only depend on the inelastic properties of particles and to their inertia; however, the area fraction or the driving energy flux Jdr is never considered. Using our experimental values for the normal and tangential restitution and friction coefficient, we can numerically solve the models. Results show values of the ratio of Ttr=Trot at maximum of 5, and more importantly, the results are not dependent on the maximum velocity of the cell. These models do not deal with an anisotropic temperature because their predictions are usually compared to numerical simulations where the energy is supposed to be added into the medium isotropically. This is why these models cannot represent our results.

along the direction of the vibration are always larger than the ones in the transverse direction. This comes from the fact that the main part of energy injected is along the y-direction, and the relatively low area fraction does not allow the redistribution of this energy toward the perpendicular direction. At low area fraction, the particles can move easily and the y-direction drives the general motion. On the other hand, we also observe a net increase of the ratio Ty=Tx with the driving velocity of the cell, but less pronounced for the lower area fraction. However, the driving velocity is not the only parameter of the problem and the amplitude can also play a role. For example, the ratio Ty=Tx ¼ 1:98 found was obtained for the smallest amplitude (A ¼ 0:556 mm) and the largest frequency (60 Hz). Under the conditions of high frequency but small amplitude, we clearly see a concentration of particles in the central area of the cell, and consequently, the corresponding energy input is small. It then explains the low ratio found in such experiment, to be compared to a similar value of Aω ¼ 0:22 m/s but with much larger amplitude (A ¼ 2:3 mm). We have set the frequency scale between 10 Hz and 30 Hz, and the corresponding amplitudes of vibration used are large enough to avoid the aggregation of particles in the center of the cell. Last, we can clearly identify the "cold" area straight from the density profiles. The second result is related to the ratio Ttr=Trot which obviously increases with Aω and which is also dependent on the area fraction. The translational temperature is quite one order of magnitude larger than the rotational temperature. Because almost all collisions between particles are quite, head-on as reflected by the high value of Ty=Tx, can explain why the transfer from translational to rotational

As a first step to describe the experimental behavior on granular temperatures, we can use existing theoretical models using a mean-field theory [37]. In this description, the rate of change

tr <sup>þ</sup> BT<sup>1</sup>=<sup>2</sup>

tr � CT<sup>1</sup>=<sup>2</sup>

tr Trot�

ð2Þ

, and C are con-

tr Trot�

of the temperature of a granular medium is determined through two coupled equations

dt <sup>¼</sup> Jdr <sup>þ</sup> <sup>G</sup>½�AT<sup>3</sup>=<sup>2</sup>

0 T<sup>3</sup>=<sup>2</sup>

where Ttr and Trot represent, respectively, the translational and rotational temperatures, and

stants, which depend only on the inelastic properties of the particles [38]. Jdr is the energy flux input into the medium and a homogeneous input of energy into the medium is assumed. These four constants are positive so that the minus signs are related to the loss of energy during the collision of particles. We state that the driving energy is acting on the translational temperature because of the preferred collisions with normal incidence. We may note that when the rotation mainly governs the

Several inelastic modelizations were proposed by Herbst et al. [38] going from a simple consideration of a constant tangential restitution coefficient up to more complex ones where the tangential restitution depending on γ<sup>12</sup> (the contact angle obtained neglecting the rotational velocities) or on the real contact angle γ. From the second equation of Eq. (2), the energy ratio Ttr=Trot can be obtained considering the medium in steady state dTrot=dt ¼ 0, allowing to get the

<sup>π</sup><sup>m</sup> <sup>p</sup> ϕg2ðϕÞ is related to the collision rate between particles; g2ðϕÞ being the pair correlation

dt <sup>¼</sup> <sup>2</sup>G½<sup>B</sup>

function at contact. In two dimensions, g2ðϕÞ¼ð1 � 7ϕ=16Þ=ð1 � ϕÞ2. A, B, B<sup>0</sup>

behavior of the granular, Jdr is included in the second equation of (2) [39, 40].

energy is rather weak, mainly at the lowest area fraction.

dTtr

8 >><

>>:

<sup>G</sup> <sup>¼</sup> <sup>16</sup> <sup>σ</sup> ffiffiffiffiffi

88 Granular Materials

dTrot

When an external vibration is acting on the granular, it can be viewed as a medium which dissipates energy while energy is added into it through the vibration per unit time. The equilibrium temperature (and state) can be found from the equilibrium equation Jdr þ Qd ¼ 0, where Jdr is the energy flux injected in the medium by the collisions of particles with the walls of the cell and Qd, the energy flux dissipated during the binary collisions between particles. Jdr is found to act in the "hot" areas of the cell, while Qd is computed in the bulk. The results obtained from experiment and geometry clearly show that the main energy input on particles occurs along the direction of the vibration. From our observations, we have considered the areas of energy input by defining two layers of thickness hH and having the same width Lx. Then, the particles' density is much smaller than the one of the medium, and we have introduced NH as the average number of particles present at any time. Thus, the bulk of the medium (i.e., the "cold" zone) reduces to dimensions hC ¼ Ly � 2hH where only NC ¼ N � 2NH particles are present at any time; the surface of this zone is then SC ¼ hCLx. In the "cold" zone, the dissipated energy depends on the collision frequency f <sup>E</sup>ðTÞ which in turns depends on the temperature <sup>T</sup> of the medium <sup>T</sup> <sup>¼</sup> <sup>m</sup>〈v<sup>2</sup> <sup>x</sup> <sup>þ</sup> <sup>v</sup><sup>2</sup> <sup>y</sup>〉=2. If we neglect the loss of energy coming from tangential restitution coefficient, the energy dissipated per collision is given by:

$$
\Delta E\_{pp} = m \frac{(e2 - 1)}{4} \langle [(\overrightarrow{v\_1} - \overrightarrow{v\_2}) \cdot \overrightarrow{n}]^2 \rangle = \frac{(e^2 - 1)}{2} T \tag{3}
$$

The frequency collision which is the inverse of the Enskog time is given in 2D by <sup>f</sup> <sup>E</sup> <sup>¼</sup> ffiffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> NC SC <sup>σ</sup>g2ðϕÞ〈v〉 <sup>¼</sup> <sup>2</sup> NC f N <sup>E</sup> where NC=SC represents the number density in the "cold" region and f N <sup>E</sup> is the number of collisions between N particles per unit time. Finally, the dissipated energy flux will be expressed as follows:

$$\begin{split} Q\_d &= f\_E^N \Delta E\_{pp} = \frac{N\_\odot^2}{h\_\odot L\_\chi} \frac{1 - e^2}{2} \sigma g\_2(\rho) \sqrt{\frac{\pi}{m}} T^{3/2} \\ &= \frac{N\_\odot^2}{h\_\odot L\_\chi} \frac{1 - e^2}{4} \sigma g\_2(\rho) \sqrt{\frac{\pi}{m}} T\_y^{3/2} \left( 1 + \frac{T\_x}{T\_y} \right)^{3/2} \end{split} \tag{4}$$

To express the energy input into the medium, we must now take into account the flux coming from the collisions between the particles and the walls of the cell. When a collision occurs, the kinetic energy change for one particle is: ΔEpw ¼ mðv0 2 <sup>y</sup> � <sup>v</sup><sup>2</sup> <sup>y</sup>Þ=2, where v<sup>0</sup> <sup>y</sup>2 and v<sup>2</sup> <sup>y</sup>, respectively, are the velocities of the particle after and before collision with the cell's wall. The cell is assumed to move with a velocity Vdr. The relative velocity equation gives v<sup>0</sup> � Vdr ¼ αðVdr � vÞ, where α is the normal restitution coefficient between the particle and the wall. The change in kinetic energy of one particle may then be rewritten as <sup>Δ</sup>Eðvy, VdrÞ ¼ <sup>m</sup> <sup>2</sup> ½ð1 þ αÞ 2 V2 dr � 2ð1 þ αÞVdrvy� v2 <sup>y</sup>ð1 � α2Þ� and the energy flux j dr, associated with particles going toward the wall, can be expressed as j dr <sup>¼</sup> NH <sup>2</sup>hH vyΔEðvy, VdrÞ, where we have assumed that NH=2 particles are going toward the wall. The net energy flux for a given wall velocity is then obtained by averaging the flux of the incoming particles with the velocity distribution function, fðvyÞ associated with the "cold" region and integrating on the velocities directed towards the wall:

$$J\_{dr}(V\_{dr}) = \bigcap\_{0}^{\prime\prime} j\_{dr} f(v\_y) dv\_y \tag{5}$$

where the distribution function of the velocity is the Gaussian one retrieved from experiments (Figure 14).

The integral (5) over the velocities gives the following result

$$J\_{dr}(V\_{dr}) = \frac{m}{4} \frac{N\_H}{h\_H} \left[ (1+a)^2 V\_{dr}^2 I\_1 - 2(1+a)V\_{dr}I\_2 - (1-a2)I\_3 \right] \tag{6}$$

<sup>I</sup>1, <sup>I</sup>2, and <sup>I</sup><sup>3</sup> are the integrals <sup>ð</sup><sup>∞</sup> 0 vi <sup>y</sup>fðvyÞ dvy (i ¼ 1::3) which are, respectively, given by:

$$I\_1 = \sqrt{\frac{T\_y}{2\pi m}} I\_2 = \frac{T\_y}{2m} I\_3 = \left(\frac{T\_y}{m}\right)^{\frac{3}{2}} \sqrt{\frac{2}{\pi}}\tag{7}$$

We may assume that the particles go from the bulk towards the "hot" areas (double collisions are neglected) so that with the average velocity found experimentally and that we are using in the comparisons. The last thing to do is to compute the average on the wall velocity: the linear term in Vdr cancels while the term related to V<sup>2</sup> dr averages to ðAωÞ2=2. Multiplying by 2 (because of 2 moving walls) gives the following expression for the injected flux of energy:

$$J\_{dr} = m \frac{N\_H}{2\hbar\_H} \left[ (1+\alpha)2(V\_{dr})2\sqrt{\frac{T\_y}{2\pi m}} - (1-\alpha2)\left(\frac{T\_y}{m}\right)^{\frac{3}{2}}\sqrt{\frac{2}{\pi}} \right] \tag{8}$$

For perfectly elastic walls (α ¼ 1), the expression proposed by Soto [35] for a sinusoidal vibration taking for their function q � T=mðAωÞ2 � is recovered; the constant <sup>q</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffi <sup>2</sup>=<sup>π</sup> <sup>p</sup> <sup>¼</sup> <sup>0</sup>:8 is a very good approximation for our experimental conditions and our values of T=mðAωÞ2. The additional contribution to the energy input from g-jitter (even if it has quite no impact on the velocities of the free floating particles in microgravity) can create an additional contribution to the wall velocity and so to the energy injected into the medium. This can be estimated as 〈δV<sup>2</sup> <sup>¼</sup> <sup>0</sup>:<sup>005</sup> <sup>m</sup><sup>2</sup>=s<sup>2</sup>〉 [41]. Since it does not exceed 15% of <sup>V</sup><sup>2</sup> dr in the worst case, it was neglected. The equilibrium between injection (Eq. (8)) and dissipation (Eq. (4)) gives:

kinetic energy change for one particle is: ΔEpw ¼ mðv0

<sup>y</sup>ð1 � α2Þ� and the energy flux j

dr <sup>¼</sup> NH

v2

expressed as j

90 Granular Materials

(Figure 14).

energy of one particle may then be rewritten as <sup>Δ</sup>Eðvy, VdrÞ ¼ <sup>m</sup>

The integral (5) over the velocities gives the following result

4 NH hH

> ð∞ 0 vi

> > I<sup>1</sup> ¼

JdrðVdrÞ ¼ <sup>m</sup>

term in Vdr cancels while the term related to V<sup>2</sup>

NH 2hH

Jdr ¼ m

vibration taking for their function q

I1, I2, and I<sup>3</sup> are the integrals

"cold" region and integrating on the velocities directed towards the wall:

JdrðVdrÞ ¼

½ð1 þ αÞ 2 V2

> ffiffiffiffiffiffiffiffiffi Ty 2πm

ð1 þ αÞ2ðVdrÞ2

�

<sup>I</sup><sup>2</sup> <sup>¼</sup> Ty 2m

(because of 2 moving walls) gives the following expression for the injected flux of energy:

ffiffiffiffiffiffiffiffiffi Ty 2πm

" # r

r

For perfectly elastic walls (α ¼ 1), the expression proposed by Soto [35] for a sinusoidal

a very good approximation for our experimental conditions and our values of T=mðAωÞ2. The additional contribution to the energy input from g-jitter (even if it has quite no impact on the velocities of the free floating particles in microgravity) can create an additional contribution to

�

T=mðAωÞ2

We may assume that the particles go from the bulk towards the "hot" areas (double collisions are neglected) so that with the average velocity found experimentally and that we are using in the comparisons. The last thing to do is to compute the average on the wall velocity: the linear

r

2 <sup>y</sup> � <sup>v</sup><sup>2</sup>

<sup>2</sup>hH vyΔEðvy, VdrÞ, where we have assumed that NH=2 particles are going

are the velocities of the particle after and before collision with the cell's wall. The cell is assumed to move with a velocity Vdr. The relative velocity equation gives v<sup>0</sup> � Vdr ¼ αðVdr � vÞ, where α is the normal restitution coefficient between the particle and the wall. The change in kinetic

toward the wall. The net energy flux for a given wall velocity is then obtained by averaging the flux of the incoming particles with the velocity distribution function, fðvyÞ associated with the

> ð ∞

0 j

where the distribution function of the velocity is the Gaussian one retrieved from experiments

<sup>y</sup>Þ=2, where v<sup>0</sup>

<sup>2</sup> ½ð1 þ αÞ

dr, associated with particles going toward the wall, can be

2 V2

drfðvyÞdvy ð5Þ

drI<sup>1</sup> � 2ð1 þ αÞVdrI<sup>2</sup> � ð1 � α2ÞI3� ð6Þ

dr averages to ðAωÞ2=2. Multiplying by 2

<sup>y</sup>fðvyÞ dvy (i ¼ 1::3) which are, respectively, given by:

� ð<sup>1</sup> � <sup>α</sup>2<sup>Þ</sup> Ty

m � �<sup>3</sup> 2 ffiffiffi 2 π

is recovered; the constant <sup>q</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffi

<sup>I</sup><sup>3</sup> <sup>¼</sup> Ty m � �<sup>3</sup> 2 ffiffiffi 2 π r

<sup>y</sup>2 and v<sup>2</sup>

<sup>y</sup>, respectively,

dr � 2ð1 þ αÞVdrvy�

ð7Þ

ð8Þ

<sup>2</sup>=<sup>π</sup> <sup>p</sup> <sup>¼</sup> <sup>0</sup>:8 is

$$T\_y = \frac{\frac{N\_H}{2h\_H}(1+a)2}{\frac{N\_\subset^2}{h\_\subset l\_x} \pi \sigma g\_2(\rho) \frac{1-\alpha}{2} \left(1 + \frac{1}{R\_\Gamma}\right)^{3/2} + 2\frac{N\_H}{h\_H}(1-a2)} m(A\omega)2} m(A\omega)2\tag{9}$$

The temperature is proportional to the square of the amplitude of the driving velocity as it should have been shown in Ref. [13]. The densities in the "cold" and "hot" areas are known so that we can compare the predictions of Eq. (9) with our experimental values of Ty calculated. To consider the dissipation due to the tangential restitution coefficient β, we introduce an effective restitution coefficient re proposed by McNamara and Luding [41] instead of e in Eq. (9): re ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e<sup>2</sup> � � <sup>q</sup>ð<sup>1</sup> � <sup>β</sup><sup>2</sup><sup>Þ</sup> � =ð1 þ 2q � βÞ r . Using q ¼ 0:5 for a disk and β ¼ 0:7;, we obtain re ¼ 0:462 instead of e ¼ 0:64. The comparison between the theoretical temperatures Ty obtained from Eq. (9) with the experimental ones calculated in the "cold" region is presented in Figure 16.

Figure 16. Comparison of the equilibrium temperature computed from Eq. (9) as a function of the driving velocity of the cell (Aω).

With the two area fractions, we have used in this study, one can see that the agreement is good. Of course, to improve the theoretical prediction of the temperature with the driving velocity, being able to predict the density NH=hH close to the top and bottom walls instead of considering the experimental value obtained from the profiles.

The anisotropy found for the temperatures created by a vibrating wall is scarcely reviewed in the literature. A recent experiment with a setup including a 3D-cylindrical [42] where the anisotropic behavior of the ratio RT ¼ Ty=Tx is reported as a function of the volume fraction of particles shows a strong increase for values below 10%, but it still remains smaller than our results. Moreover, a theoretical study including two different Maxwellian distributions for parallel and perpendicular directions about the vibration and a density along the vibration axis is presented in Ref. [43]. A balance between energy fluxes along the direction of vibration and perpendicular to it gives the ratio RT and predicts that, for perfectly elastic walls, this ratio would only depend on the restitution coefficient. This is not our experimental situation where the ratio RT is much larger for the lowest density. Thus, it impossible to relate this theoretical approach with our study since our density profile is much different for a one driven by gravity. Nevertheless, we may predict values of RT if we assume a constant density in the "cold" area.

## 4. Conclusion

We have reported experimental investigations on the dynamics of a model granular medium. Experiments have been performed in a low-gravity environment. The cell containing the medium is subjected to external vibration which drives the collective motion of the particles. As the dynamical behavior of the medium is driven by the kinematics of the particles, highspeed video recording coupled to an individual particle tracking technique allows to obtain the trajectory of each particle. From these raw data, the inelastic parameters of the particles which are at the origin of the dynamics of the whole medium can be retrieved as well as a direct measurement of the energy (or temperatures). We have found that depending on the type of particles used, the normal restitution coefficient can be dependent on the relative impact velocity between two particles but not always. One way to characterize the inelastic nature of the collisions is to look to the energy decay once the medium is freely evolving. We have obtained smaller experimental relaxation times of this energy than the ones predicted by theories at least if we do not take into account the velocity dependent of this restitution coefficient. It is also interesting to note that the effect played by the rotation of the particles can significantly affect the whole behavior of the medium. In particular, we have reported the translational temperatures along and perpendicular to the direction of vibration as well as the rotational temperatures. When compared to existing theories, it appears that there are significant differences which also depend on the driving velocity and on the concentration of the medium. Two major points on the comparison can be raised: First the density is not homogeneous in the cell and second the translational velocities are much higher in the direction of vibration than perpendicular to it (versus a homogeneous input of energy as considered in theories). We report that the balance of the energy fluxes along the vibration can correctly represent the behavior of the granular temperature with the driving velocity of the cell and with the area fraction. In this balance, the contribution of the tangential velocities to the dissipation must be considered. At least the distinction between the dissipation due to the collisions between the particles which is proportional to the average temperature T ¼ ðTx þ TyÞ=2 and the driving flux, which depends only on Ty, was introduced, but on the basis of the experimental ratio Ty=Tx. This ratio increases when the volume fraction decreases and it also depends on the driving velocity. A theoretical determination of Ty=Ty which could reproduce these behaviors should involve the non-elastic collisions with the lateral walls, but is let for a future developments.

## Acknowledgements

With the two area fractions, we have used in this study, one can see that the agreement is good. Of course, to improve the theoretical prediction of the temperature with the driving velocity, being able to predict the density NH=hH close to the top and bottom walls instead of consider-

The anisotropy found for the temperatures created by a vibrating wall is scarcely reviewed in the literature. A recent experiment with a setup including a 3D-cylindrical [42] where the anisotropic behavior of the ratio RT ¼ Ty=Tx is reported as a function of the volume fraction of particles shows a strong increase for values below 10%, but it still remains smaller than our results. Moreover, a theoretical study including two different Maxwellian distributions for parallel and perpendicular directions about the vibration and a density along the vibration axis is presented in Ref. [43]. A balance between energy fluxes along the direction of vibration and perpendicular to it gives the ratio RT and predicts that, for perfectly elastic walls, this ratio would only depend on the restitution coefficient. This is not our experimental situation where the ratio RT is much larger for the lowest density. Thus, it impossible to relate this theoretical approach with our study since our density profile is much different for a one driven by gravity. Nevertheless, we may predict values of RT if we assume a constant density in the "cold" area.

We have reported experimental investigations on the dynamics of a model granular medium. Experiments have been performed in a low-gravity environment. The cell containing the medium is subjected to external vibration which drives the collective motion of the particles. As the dynamical behavior of the medium is driven by the kinematics of the particles, highspeed video recording coupled to an individual particle tracking technique allows to obtain the trajectory of each particle. From these raw data, the inelastic parameters of the particles which are at the origin of the dynamics of the whole medium can be retrieved as well as a direct measurement of the energy (or temperatures). We have found that depending on the type of particles used, the normal restitution coefficient can be dependent on the relative impact velocity between two particles but not always. One way to characterize the inelastic nature of the collisions is to look to the energy decay once the medium is freely evolving. We have obtained smaller experimental relaxation times of this energy than the ones predicted by theories at least if we do not take into account the velocity dependent of this restitution coefficient. It is also interesting to note that the effect played by the rotation of the particles can significantly affect the whole behavior of the medium. In particular, we have reported the translational temperatures along and perpendicular to the direction of vibration as well as the rotational temperatures. When compared to existing theories, it appears that there are significant differences which also depend on the driving velocity and on the concentration of the medium. Two major points on the comparison can be raised: First the density is not homogeneous in the cell and second the translational velocities are much higher in the direction of vibration than perpendicular to it (versus a homogeneous input of energy as considered in theories). We report that the balance of the energy fluxes along the vibration can correctly represent the behavior of the granular temperature with the driving velocity of the cell and

ing the experimental value obtained from the profiles.

4. Conclusion

92 Granular Materials

The authors like to thank NOVESPACE and the CNES for giving them the possibility to board the A300-zero G in order to perform the experimental study.

## Author details

Yan Grasselli<sup>2</sup> \*, Georges Bossis<sup>1</sup> , Alain Meunier<sup>1</sup> and Olga Volkova1

\*Address all correspondence to: yan.grasselli@skema.edu

1 Laboratoire de la Matière Condensée, University of Nice Sophia Antipolis, Parc Valrose, France

2 SKEMA Business School, University of Côte d'Azur, Sophia-Antipolis, France

## References


[7] Opsomer E, Ludewig F, Vandewalle N. Phase transitions in vibrated granular systems in

[8] Falcon E, Wunenburger R, Evesque P, Fauve S, Chabot C, Garrabos Y, Beysens D. Cluster formation in a granular medium fluidized by vibrations in low gravity. Physical Review

[9] Falcon E, Aumaıtre S, Evesque P, Palencia F, Lecoutre-Chabot C, Fauve S, Beysens D, Garrabos Y. Collision statistics in a dilute granular gas fluidized by vibrations in low

[10] Noirhomme M, Opsomer E, Vandewalle N, Ludewig F. Granular transport in driven

[11] Opsomer E, Noirhomme M, Vandewalle N, Falcon E, Merminod S. Segregation and pattern formation in dilute granular media under microgravity conditions. Nature. Article number:

[12] Bossis G, Grasselli Y, Volkova O. Granular rheology in zero gravity. Journal of Physics:

[13] Bhateja A, Sharma I, Singh JK. Scaling of granular temperature in vibro-fluidized grains.

[14] Pathak SN, Jabeen Z, Das D, Rajesh R. Energy decay in three-dimensional freely cooling

[15] Huthmann M, Aspelmeier T, Zippelius A. Granular cooling of hard needles. Physical

[16] Villemot F, Talbot J. Homogeneous cooling of hard ellipsoids. Granular Matter. 2012;14:91 [17] Harth K, Kornek U, Trittel T, Strachauer U, Home S, Will K, Stannarius R. Granular gases of rod-shaped grains in microgravity. Physical Review Letters. 2013;110:144102

[18] Yan-Pei C, Evesque P, Mei-Ying Chin H. Breakdown of energy equipartition in vibro-

[19] Brilliantov NV, Pöschel. Kinetic Theory of Granular Gases. Oxford University Press; 2004 [20] Tatsumi S, Murayma Y, Hayakawa H, Sano M. Experimental study on the kinetics of granular gases under microgravity. Journal of Fluid Mechanics. 2009;641:521

[21] Grasselli Y, Bossis G. Three-dimensional particle tracking for the characterization of micrometer-size colloidal particles. Journal of Colloid and Interface Science. 1995;1:269

[22] Gondret P, Lance M, Petit L. Bouncing motion of spherical particles in fluids. Physics of

[23] KantakAdvait Ashok. Wet particles collisions [thesis]. University of Colorado at Boulder:

fluidized granular media in micro-gravity. Physics Letters. 2012;29:074501

microgravity. Physical Review E. 2011;84:051306

gravity. Europhysics Letters. 2006;74:830

Condensed Matter. 2004;16:3279

Physics of Fluids. 2016;28:043301

Review E. 1999;60:654

Fluids. 2002;14:268

2005. DOI: AAT 3190381

granular gas. The European Physical Journal E. 2015;38:9

granular gas. Physical Review Letters. 2014;112:038001

Letters. 1999;83:440

94 Granular Materials

1, 2017;3


## **Particle Jetting Induced by the Impulsive Loadings**

Kun Xue, Xiaoliang Shi, Kaiyuan Du and Haoran Cui

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.68795

### **Abstract**

[41] McNamara S, Luding S. Energy flows in vibrated granular media. Physical Review E.

[42] Windows-Yule CRK, Parker DJ. Boltzmann statistics in a three-dimensional vibrofluidized granular bed: Idealizing the experimental system. Physical Review E. 2013;87:022211 [43] van der Meer D, Reimann P. Temperature anisotropy in a driven granular gas. Europhysics

1998;58:813

96 Granular Materials

Letters. 2006;74:384

Particle rings/shells/cylinders dispersed by the radial impulsive loadings ranging from strong blast waves to moderate shock waves form a dual coherent jetting structure consisting of particle jets which have different dimensions. In both circumstances, the primary jets are found to initiate from the inner surface of particle layers and propagate through the thickness of particle layers, which are superimposed by a large number of much smaller secondary jets initiating from the outer surface of particle layers upon the reflection of the shock wave. This chapter first presents a summary of the experimental observations of the hierarchical particle jetting mainly via the cinematographic techniques, focusing on the characteristics of the primary particle jet structure. Due to the distinct behaviors of particles subjected to the strong blast and moderate shock waves, specifically solid-like and fluid-like responses, respectively, the explosive and shock-induced particle jetting should be attributed to distinct mechanisms. A dual particle jetting model from the perspective of continuum is proposed to account for the explosive-induced particle jetting. By contrast the shock-induced particle jetting arises from the localized particle shear flows around the inner surface of particle layers which result from the heterogeneous network of force chains.

**Keywords:** particle jetting, blast wave, shock wave, force chains, discrete element method, multiphase flows

## **1. Introduction**

When particles are dispersed by an impulsive pressure loading, the expanding particle cloud typically forms a nonuniform structure that takes the form of particle jets whose leading edges are agglomerates of constituent grains [1–12]. A host of experimental evidence from a wide range of sources shows that the expanding cloud of explosively disseminated material comprises of

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

"particles" or fragments that have significantly different dimensions from those associated with the original material as shown in **Figure 1(a)** and **(b)** [1–7, 10–13]. Photographic evidence shows characteristic jets or fingers behind these expanding fragments. These coherent conical particle jets travel ballistically as shedding mass along the trajectories with increasingly diffuse edges.

Particle jetting has been widely observed in volcanic eruptions, supernovae, explosion of landmines, thermobaric explosion (TBX), fuel-are explosion (FAE), and dense inert metal explosive (DIME) [14–17]. The formation of particle jetting has also been observed during the impact of solid projectile on granular media [18]. The structure of particle jets in terms of the jet number of size is important to the viability of many applications. For instance, the strength of TBX and FAE needs to be enhanced by the after-burning of the reactive particles dispersed in the payload cloud. The detonation performance of the particle laden cloud depends on both the shape and concentration of the cloud which in turn is a result of the jet mixing [19]. In contrast with the large-scale injury radius of TBX and FAE, DIME utilizes the high-speed dense metal–particle jets to hit the targets in close range. Thus the momentum of particle jets determines the lethal radius. Another opposite application is mitigation of the blast pressure

**Figure 1.** Explosive dispersal of dry (a) and wetted glass beads (b) using cylindrically stratified configurations [12]. Shock dispersal of flour (c) and quartz sand (d) particles using semi-two-dimensional configurations.

(both prompt and quasi-static) associated with the detonation, since a commonly used technique to reduce effects of blast from explosives is to surround the explosive with a layer of liquid, powder, or a slurry mixture of the two. Drag is seen as a potential mechanism to transfer energy from the blast wave to the disseminated particles or droplets so the size of particles or formation of jets is important in determining the efficiency of this mechanism. Frost and Zhang have reviewed many of the processes occurring in heterogeneous blast including jet formation [15, 17, 20, 21].

"particles" or fragments that have significantly different dimensions from those associated with the original material as shown in **Figure 1(a)** and **(b)** [1–7, 10–13]. Photographic evidence shows characteristic jets or fingers behind these expanding fragments. These coherent conical particle jets travel ballistically as shedding mass along the trajectories with increasingly diffuse edges.

98 Granular Materials

Particle jetting has been widely observed in volcanic eruptions, supernovae, explosion of landmines, thermobaric explosion (TBX), fuel-are explosion (FAE), and dense inert metal explosive (DIME) [14–17]. The formation of particle jetting has also been observed during the impact of solid projectile on granular media [18]. The structure of particle jets in terms of the jet number of size is important to the viability of many applications. For instance, the strength of TBX and FAE needs to be enhanced by the after-burning of the reactive particles dispersed in the payload cloud. The detonation performance of the particle laden cloud depends on both the shape and concentration of the cloud which in turn is a result of the jet mixing [19]. In contrast with the large-scale injury radius of TBX and FAE, DIME utilizes the high-speed dense metal–particle jets to hit the targets in close range. Thus the momentum of particle jets determines the lethal radius. Another opposite application is mitigation of the blast pressure

**Figure 1.** Explosive dispersal of dry (a) and wetted glass beads (b) using cylindrically stratified configurations [12].

Shock dispersal of flour (c) and quartz sand (d) particles using semi-two-dimensional configurations.

Extensive experimental investigations of the explosive- or shock-induced particle jetting mainly using cinematographic techniques provide fundamental data regarding both structure and evolution of particle jets. Closer look into the high-speed photos of either explosive- or shock-induced particle jetting reveals a dual structure (see **Figure 1(c)**) [1, 2, 12]. Primary jets initiated on the inner surface of the particle layers take shape during the first dozens of microseconds after the detonation of the central explosive evidenced by the light stripes detected from the radiographs of the explosive dispersal of particle shells [5, 6]. Upon the reflection of the shock wave on the outer surface of particle layers, a large number of smaller jets begin to burgeon from the outer surface and quickly develop into a full bloom [3]. The dominant primary jets are expelled from the outer surface and overtake the smaller secondary jets, merging of secondary jets occurring through the aerodynamic interaction. The respective evolutions of the primary and secondary jets are not so distinguishable from the radiographs and high-speed photos of the explosive dispersal of particles (see **Figure 2**). But the statistic distribution of jet size unravels two distinctive peaks representing the primary and secondary jets, respectively [13]. In order to overcome the difficulties in distinguishing the primary and secondary jets, a semi-two-dimensional configuration based on the Hele-Shaw cell that will be discussed in Section 3.1 was employed to access the evolution of both sets of jets subjected to the radial shock loading. Although the overpressure of weak shock waves is several orders of magnitude lower than that of blast waves, the formation and evolution of the primary and secondary jets as shown in **Figure 3** have astonishingly similar characteristics in terms of the initiation sequence and the signature structure [2]. Whereas whether or not the jetting process in these two extreme conditions follow the same path is still debatable.

Great efforts have been devoted to investigate the dependence of the jet number on a variety of parameters, including the configuration of charge, the mass ratio of the payload and the explosive (M/C), the inner and outer radius of particle layers, the particle material and size, and the moisture content, etc., mainly in the case of explosive dispersal of particles [2, 4, 6, 8, 13]. Specifically, Zhang et al. found that the numbers of primary and secondary jets dispersed by the 44 mm diameter of central explosive cylinder are 1.8 and 1.5 times those with the 10 mm diameter of explosive [14]. Frost et al. found that the jetting phenomenon is much more visible in cases of explosive dispersal of brittle or ductile powders, such as quartz sand, glass beads, SiC powders, aluminum powders, copper powders, compared with rigid and hard powders, like stain steel particles that are dispersed into the particle cloud rather than particle jets [10]. Frost and Xue both found that the addition of the interstitial water/oil significantly increases the jet number [12, 13].

**Figure 2.** Radiographs and high-speed photos of explosive dispersal of glass beads (a) [5], dry (above panel of (b)) and wet (bottom panel of (b)) quartz sand grains [13]. (c): the statistic distribution of size of the explosive induced dry and wet sand jets at t = 2ms [13].

Some fundamental problems need to be addressed in this regard. First, several variables are correlated rather than independent so that it is impossible to single out the effect of the individual variable. For instance, changing the inner or outer radius of the particle layers or particle materials would inevitably alter the M/C that proves to be key factor determining the jet number. Rodriguez et al. proposed an alternative way to measure the effects of pertinent factors [2]. Acceleration of the outer surface of particle ring, which is a function of a variety of parameters, is found to be key determining factor. Therefore, choosing some proper dynamic

**Figure 3.** High-speed photos of semi-two-dimensional shock-induced particle jetting [2].

variables instead of structural parameters may well provide a new perspective in this regard, but entailing a thorough understanding of the physics underlying the particle jetting.

Second, distinguishing the primary and secondary jets from the radiographs or high-speed photos of the explosive dispersal of particles using either spherical or cylindrical stratified configurations is so difficult if not impossible that the validity of the experimental results is questionable thanks to the superimposition of two sets of jets on the timescale of microseconds. Adopting a semi-two-dimensional configuration in which particle rings are dispersed by the radial propagating shock waves seems to be promising approach to this problem. Besides, no detonation product gases obscuring the particle jets and substantially prolonged duration of jets facilitate the observation of particle jetting. But to what extent the shock-induced particle jetting can mimic that driven by the central explosion is quite questionable taking into account that the overpressure of shock waves is several orders lower than that of blast waves.

Predicting the jet number entails the knowledge of the mechanisms governing the primary and secondary jets, respectively. Several theories have been put forward, but understanding the origin of particle jetting still remains a significant challenge [3, 8, 9, 13, 22, 23]. The timescale for the formation of primary jets predicted by the Rayleigh-Taylor instability is much slower than the experimental observation [6]. Another interface instability theory involves the perturbation on the inner and outer surfaces of the particle layers that act as the microjets precipitating the macrojets propagating into the bulk. Riple et al. demonstrated the evolution of the initial perturbation (see **Figure 4**) into well-developed jets and argued that the casing fragments and other imperfections may provide the initial perturbation [3]. However, particle jetting occurs regardless of the presence of the inner and outer casings and shows similar structure. Certain intrinsic imperfections with the length scale similar to the jets should exist if this theory holds. An increasing number of investigators have focused their attention on the bulk fracture of powder bed. Frost et al. postulated that the breakup of a layer of particles at high strain rates was governed by a balance of expansion inertia effects tending to fracture the layer versus viscous dissipation that tends to maintain the stability of the layer [24]. Along this line, Xue et al. developed a theoretical model account for the instability onset of the expanding powder shell [13]. Milne et al. conjectured that the powder is explosively compacted into

Some fundamental problems need to be addressed in this regard. First, several variables are correlated rather than independent so that it is impossible to single out the effect of the individual variable. For instance, changing the inner or outer radius of the particle layers or particle materials would inevitably alter the M/C that proves to be key factor determining the jet number. Rodriguez et al. proposed an alternative way to measure the effects of pertinent factors [2]. Acceleration of the outer surface of particle ring, which is a function of a variety of parameters, is found to be key determining factor. Therefore, choosing some proper dynamic

**Figure 2.** Radiographs and high-speed photos of explosive dispersal of glass beads (a) [5], dry (above panel of (b)) and wet (bottom panel of (b)) quartz sand grains [13]. (c): the statistic distribution of size of the explosive induced dry and

wet sand jets at t = 2ms [13].

100 Granular Materials

**Figure 4.** Schematics of the formation of primary and secondary jets caused by the inner and outer cases, respectively. (a): fragmentation of the inner and outer cases; (b): formation of primary and secondary jets arising from the gaps between fragments of inner and outer cases, respectively; (c): primary jets overtake the secondary jets. Inset: the corresponding snapshots from the hydrodynamic simulations of Ripley et al. [3,19].

a brittle solid which then forms cracks as the shell expands [5]. This conjecture is consistent with the observations that the primary jetting occurs during the first wave transit times. The major obstacle of this argument is that the compacted powder cannot sustain the tension or the surface energy, both among the essential components comprising the brittle fragmentation of solids. A few attempts try to understand the secondary jets, and the earlier works of Ripley et al. focused on the Richtmyer-Meshkov instability (RMI), which showed well-defined persistent jetting structures matching the number of prescribed outer surface perturbations [3]. However, the timescale for formation was slow and the surface instability did not propagate into the bulk [3]. Xue et al. modified the hollow sphere expansion model that originally accounts for the spallation of shocked solids so that the external particle jetting can be seen as parallel to the solid spallation [22].

Despite the resembling phenomenal features sheared by the explosive- and shock-induced particle jetting, the shock interaction with particles in the explosive dispersal is substantially stronger than that in the weak shock dispersal. In the former case, particles are compressed into solids with the density almost same as that of the constituent materials when the particle jetting commences. It suggests that a continuum approach is appropriate to model the explosion-driven particle jetting. By contrast, the weak shock wave only initiates the homogeneous or localized unsteady flows on the particle scale. The shocked particles behave more like fluids rather than solids. Unsteady and heterogeneous particle flows occurring during the weak shock interaction with particles entail a particle scale approach. Xue et al. described the particle scale formation and evolution of particle jets via the discrete element method (DEM), shedding some lights on the distinctive origins of the shock-induced particle jetting [25].

This chapter first reviews the up-to-date understanding of the phenomenology and physics of the particle jetting in both explosion-driven and shock-induced cases. Special attention is focused on theoretical progresses in unraveling the mechanism behind the respective particle jetting and establishing models account for the onset of jetting, which is elaborated in Sections 2 and 3. Further work and possible breakthrough in this regard would be discussed in Section 4. The conclusion is presented in Section 5.

## **2. Explosion-driven particle jetting**

## **2.1. Strong shock interaction with particles**

a brittle solid which then forms cracks as the shell expands [5]. This conjecture is consistent with the observations that the primary jetting occurs during the first wave transit times. The major obstacle of this argument is that the compacted powder cannot sustain the tension or the surface energy, both among the essential components comprising the brittle fragmentation of solids. A few attempts try to understand the secondary jets, and the earlier works of Ripley et al. focused on the Richtmyer-Meshkov instability (RMI), which showed well-defined persistent jetting structures matching the number of prescribed outer surface perturbations [3]. However, the timescale for formation was slow and the surface instability did not propagate into the bulk [3]. Xue et al. modified the hollow sphere expansion model that originally accounts for the spallation of shocked solids so that the external particle jetting can be seen as

**Figure 4.** Schematics of the formation of primary and secondary jets caused by the inner and outer cases, respectively. (a): fragmentation of the inner and outer cases; (b): formation of primary and secondary jets arising from the gaps between fragments of inner and outer cases, respectively; (c): primary jets overtake the secondary jets. Inset: the corresponding

Despite the resembling phenomenal features sheared by the explosive- and shock-induced particle jetting, the shock interaction with particles in the explosive dispersal is substantially stronger than that in the weak shock dispersal. In the former case, particles are compressed into solids with the density almost same as that of the constituent materials when the particle jetting commences. It suggests that a continuum approach is appropriate to model the explosion-driven particle jetting. By contrast, the weak shock wave only initiates the homogeneous or localized unsteady flows on the particle scale. The shocked particles behave more like fluids rather than solids. Unsteady and heterogeneous particle flows occurring during the weak shock interaction with particles entail a particle scale approach. Xue et al. described the particle scale formation and evolution of particle jets via the discrete element method (DEM), shedding some lights on the distinctive origins of the shock-induced particle jetting [25].

This chapter first reviews the up-to-date understanding of the phenomenology and physics of the particle jetting in both explosion-driven and shock-induced cases. Special attention is focused on theoretical progresses in unraveling the mechanism behind the respective

parallel to the solid spallation [22].

102 Granular Materials

snapshots from the hydrodynamic simulations of Ripley et al. [3,19].

One generally accepted fact of the explosive-driven particle jetting is that particle instabilities occur during the first dozens of the microseconds after the detonation of the central explosive. It is thus necessary to elucidate the interactions between particles, shock waves, and detonation product gases. Hydrodynamic simulations [22] have been performed to reveal the evolution of dry and saturated sand layers surrounding the spherical central explosive (TNT or HXM), the configuration illustrated in **Figure 5(a)**. In order to accurately describe the dynamic responses of wet sand with different degrees of saturation β, we adopted a modified version of Laine and Sandvik model developed by Grujicic et al. [26] to account for the effect of moisture content via explicitly incorporating the degree of saturation in the equation of state (EOS) and the strength model. Given the relative incompressibility of the water phase, the compressibility of the wet sand is increasingly reduced with the degree of saturation as illustrated by the EOS of the wet particles with varying saturation (see **Figure 5(b)**). Besides, the wet sand's yield stress is reduced due to the moisture-induced interparticle lubrication effects leading to a reduced effective friction coefficient (see **Figure 5(c)**. For details of the modified compaction model, readers can be referred to Refs. [26, 27] (see **Figure 6**).

The evolvement of the sand shell upon the blast wave can be well embodied by the variations in its radial density profile as shown in **Figure 3**. The sequence of events basically resembles those occurring in the shock-loaded water shell described by Milne et al. [5, 6]. When the shock front reflects upon the outer surface of the particle shell, the rarefaction wave travels back into particles and pulls away a thin spall layer moving forward into air. The compressive

**Figure 5.** (a) Schematic of the spherical stratified configuration used in the hydrodynamic simulations. (b) EOS curves of the sand with varying degree of saturation. (c) Variations in dependence of the sand's yield stress on the pressure with increasing moisture contents [22].

**Figure 6.** Evolutions of density profiles in dry sand (a) and saturated sand (b) after the detonation of the central explosive (TNT) [22].

stresses in the compacted particles are relaxed in the wake of the rarefaction wave accompanied by the rapid decrease of the packing density. The expansion of detonation product gases sends a shock wave into the particles, which arrests the rarefaction wave in its path in the case of dry sand or recompact particles diluted by rarefaction wave in the case of saturated sand. As a result, besides the outmost thin spall layer, the particle shell evolves into two distinct layers, namely the inner compact layer and outer dilute layers. The inner compact layer retains the maximum density almost as that of pure quartz and expands as an incompressible shell during a relatively long time, at least during the first hundred of microseconds after the detonation of central explosive. The hypothesis is supported by the consistent velocity across the thickness of the inner compact layer (see **Figure 7(a)** and **(b)**). Opposedly, particles inside the outer dilute layer lose the persistent contacts in the wake of the rarefaction wave. The

**Figure 7.** Evolutions of velocities of the inner and outer surfaces of the compact dry/saturated sand layer driven by the detonation of central TNT (a) and HXM (b).

mass ratio of the compact and dilute layers depends on the geometry and the composition of granular shell, as well as the strength of central explosive.

Due to the trivial compressibility of saturated sand, the acceleration of the compact saturated sand layer is much stronger than that of dry sand since less shock energy is dissipated among the compaction. The expanding velocity of the compact saturated sand layer is much larger than that of the dry sand (see **Figure 7(a)**).

## **2.2. Dual particle jetting model**

stresses in the compacted particles are relaxed in the wake of the rarefaction wave accompanied by the rapid decrease of the packing density. The expansion of detonation product gases sends a shock wave into the particles, which arrests the rarefaction wave in its path in the case of dry sand or recompact particles diluted by rarefaction wave in the case of saturated sand. As a result, besides the outmost thin spall layer, the particle shell evolves into two distinct layers, namely the inner compact layer and outer dilute layers. The inner compact layer retains the maximum density almost as that of pure quartz and expands as an incompressible shell during a relatively long time, at least during the first hundred of microseconds after the detonation of central explosive. The hypothesis is supported by the consistent velocity across the thickness of the inner compact layer (see **Figure 7(a)** and **(b)**). Opposedly, particles inside the outer dilute layer lose the persistent contacts in the wake of the rarefaction wave. The

**Figure 6.** Evolutions of density profiles in dry sand (a) and saturated sand (b) after the detonation of the central explosive

**Figure 7.** Evolutions of velocities of the inner and outer surfaces of the compact dry/saturated sand layer driven by the

detonation of central TNT (a) and HXM (b).

(TNT) [22].

104 Granular Materials

The decomposition of the particle shell into the inner compact and outer dilute layers as a result of shock interaction prompts us to speculate that the fragmentation of the inner and outer layers correspond to the primary and secondary particle jetting, respectively. This speculation satisfies some fundamental facts that (1) the primary and secondary particle jets initiate from the inner and outer surface of particle shells, respectively; (2) the secondary particle jetting occurs upon the reflection of the shock wave on the outer surface; (3) the primary jets overtake the primary jets in later times. Therefore, a dual particle jetting model illustrated in **Figure 8** has been put forward to account for the formation of the primary and secondary jets. The following task is to elaborate the proper models describing the respective fragmentation of the inner and outer particle layers. These models should be based on the underlying mechanisms and validated against the experimental results, the onsets of primary/secondary

**Figure 8.** Illustration of the dual particle jetting model, which consists of the formation of the inner compact and outer dilute layers, and the breakup of these two distinct layers [22].(a):the initial annular configuration; (b): expansion of detonation gases issues the compression wave; (c): the reflected rarefaction wave causes the spallation of outermost layer; (d):fragmentation of the inner compact layer; (e): protrusion of secondary jets; (f): overtake of secondary jets by primary jets.

jetting, the size of primary/secondary jets, and the dependence of the jet number on a variety of factors as well.

## **2.3. Primary particle jetting model: the destabilization of the expanding shell**

The consistent density and velocity across the thickness of the inner compact layer indicate that the compacted layer expands as the incompressible shell. Under this premise, we consider a sphere shell characterized by an inner radius *R*<sup>1</sup> and outer radius *R*<sup>2</sup> as shown in **Figure 9**, which can be determined by the hydrodynamic simulations (see **Figure 6**). The thickness of the shell is *R*<sup>2</sup> –*R*<sup>1</sup> . Adopting the spherical coordinate system associated with the frame (*e***<sup>r</sup>** , *e***θ**, *e***φ**), the outward divergent motion of the continuous sand shell demonstrated in experiments is modeled by applying a uniform velocity *Vr e***<sup>r</sup>** at the inner surface (*R* = *R*<sup>1</sup> ), which can also be derived from the hydrodynamic simulation (see **Figure 7**).

Applying the continuity and momentum equations to the incompressible granular shells that can be described as viscoplastic materials, the analytical circumferential stress can be derived as follows (details of formulation can be referred to Ref. [13]).

$$\sigma\_{\theta\theta} = \frac{\rho}{2} \frac{V\_1^2 R\_1^4}{2} \left(\frac{1}{R^4} - \frac{1}{R\_1^4}\right) + \frac{\eta}{\sqrt{2}} \frac{V\_1 R\_1^2}{\left(R^3 + \frac{2}{R\_2^3}\right)}$$

$$-2\rho \{V\_1^2 R\_1 + A \ R\_1^2\} \left(\frac{1}{R} - \frac{1}{R\_2}\right) + \sqrt{3} \ \tau\_c \left(2\ln\frac{R}{R\_2} + 1\right) \tag{1}$$

where *ρ* is the mass density of the sand shell, *τ<sup>c</sup>* is the yield stress, and *η* is dynamic viscosity. Bearing in mind that the yield stress, *τ<sup>c</sup>* , is a function of both saturation degree and the pressure applied on the inner surface which is in the order of O(10<sup>0</sup> –101 ) Mpa (see **Figure 5(c)**), the yield stress of saturated sand (~1 MPa) is much lower than that that of dry sand (~13.7 MPa) due to the lubrication effect assuming average pressure *P*‾‾~<sup>10</sup> <sup>M</sup>Pa. The dynamic viscosity, *η*, is in the order of O(10−1).

To predict the instability onset of the expanding sand shell, we will invoke a criterion for instability that has been shown to reasonably emulate more rigorous stability analysis [28]. This method can be viewed as an application of Le Chatelier's principle that states that for a system to be stable any deviation from equilibrium must bring about forces that tend to restore equilibrium. In general, the loss of stability is assumed to take place when an increment in strain occurs with no simultaneous increase in pressure or in load.

To obtain the circumferential pressure in the expanding shell, the circumferential stress from Eq. (1) is integrated through the thickness *h* of the shell,

$$\begin{split} T &= \int\_{\tilde{\mathcal{R}}\_{i}}^{\tilde{\mathcal{R}}\_{i}} \sigma\_{\theta\theta} \, d\mathcal{R} \\ &= \frac{\rho}{2} \frac{V\_{r}^{2} R\_{i}}{2} \Big( \frac{1}{3} - \frac{4}{3} \frac{R\_{i}^{3}}{R\_{\mathcal{Z}}^{3}} + \frac{R\_{i}^{4}}{R\_{\mathcal{Z}}^{4}} \Big) + \eta \, V\_{r} \sqrt{2} \Big[ \frac{1}{4} + \frac{3}{4} \frac{R\_{i}^{2}}{R\_{\mathcal{Z}}^{2}} - \frac{R\_{i}^{3}}{R\_{\mathcal{Z}}^{3}} \Big] \\ &- 2\rho \Big( V\_{r}^{2} R\_{1} \Big) \Big( \ln \frac{R\_{\sharp}}{R\_{\mathcal{I}}} - 1 + \frac{R\_{i}}{R\_{\mathcal{Z}}} \Big) + \sqrt{3} \, \tau\_{c} \, R\_{1} \Big( 1 - \frac{R\_{\sharp}}{R\_{\mathcal{I}}} - 2 \ln \frac{R\_{i}}{R\_{\mathcal{Z}}} \Big). \end{split} \tag{2}$$

**Figure 9.** Configuration of an expanding spherical shell with inner radius *R*<sup>1</sup> and outer radius *R*<sup>2</sup> [13].

jetting, the size of primary/secondary jets, and the dependence of the jet number on a variety

The consistent density and velocity across the thickness of the inner compact layer indicate that the compacted layer expands as the incompressible shell. Under this premise, we con-

**Figure 9**, which can be determined by the hydrodynamic simulations (see **Figure 6**). The

Applying the continuity and momentum equations to the incompressible granular shells that can be described as viscoplastic materials, the analytical circumferential stress can be derived

yield stress of saturated sand (~1 MPa) is much lower than that that of dry sand (~13.7 MPa) due to the lubrication effect assuming average pressure *P*‾‾~<sup>10</sup> <sup>M</sup>Pa. The dynamic viscosity, *η*, is

To predict the instability onset of the expanding sand shell, we will invoke a criterion for instability that has been shown to reasonably emulate more rigorous stability analysis [28]. This method can be viewed as an application of Le Chatelier's principle that states that for a system to be stable any deviation from equilibrium must bring about forces that tend to restore equilibrium. In general, the loss of stability is assumed to take place when an incre-

To obtain the circumferential pressure in the expanding shell, the circumferential stress from

<sup>4</sup>) + *η Vr* √

\_\_ 2[ \_1 4 + 3 *R*<sup>1</sup> 2 \_ 4 *R*<sup>2</sup> <sup>2</sup> <sup>−</sup> *<sup>R</sup>*<sup>1</sup> 3 \_ *R*2

<sup>3</sup> *<sup>τ</sup><sup>c</sup> <sup>R</sup>*1(<sup>1</sup> <sup>−</sup> *<sup>R</sup>*\_2

*R*1 − 2 ln *R*\_1 *<sup>R</sup>*2).

, *e***θ**, *e***φ**), the outward divergent motion of the continuous sand shell demonstrated

) <sup>+</sup> *<sup>η</sup> <sup>V</sup>*<sup>1</sup> *<sup>R</sup>*<sup>1</sup>

2 \_\_\_\_\_ √ \_\_ 2 ( \_\_1 *<sup>R</sup>*<sup>3</sup> <sup>+</sup> \_\_2 *R*2 3 ) 

and outer radius *R*<sup>2</sup>

. Adopting the spherical coordinate system associated with the

<sup>3</sup> *<sup>τ</sup>c*(<sup>2</sup> ln\_\_*<sup>R</sup> R*2 <sup>+</sup> <sup>1</sup>),

, is a function of both saturation degree and the pres-

3]

–101

is the yield stress, and *η* is dynamic viscosity.

as shown in

),

(1)

(2)

at the inner surface (*R* = *R*<sup>1</sup>

) Mpa (see **Figure 5(c)**), the

**2.3. Primary particle jetting model: the destabilization of the expanding shell**

which can also be derived from the hydrodynamic simulation (see **Figure 7**).

<sup>2</sup> *R*<sup>1</sup> 4 \_\_\_\_\_ <sup>2</sup> (

ment in strain occurs with no simultaneous increase in pressure or in load.

*R*\_2 *R*1 − 1 + *R*\_1 *<sup>R</sup>*2) <sup>+</sup> <sup>√</sup> \_\_

\_\_1 *<sup>R</sup>*<sup>4</sup> <sup>−</sup> \_\_1 *R*2 4

<sup>2</sup> *R*<sup>1</sup> + *A R*<sup>1</sup> 2 )( \_\_1 *<sup>R</sup>* <sup>−</sup> \_\_1 *<sup>R</sup>*2) <sup>+</sup> <sup>√</sup> \_\_

sider a sphere shell characterized by an inner radius *R*<sup>1</sup>

–*R*<sup>1</sup>

in experiments is modeled by applying a uniform velocity *Vr e***<sup>r</sup>**

as follows (details of formulation can be referred to Ref. [13]).

*σθθ* <sup>=</sup> *<sup>ρ</sup> <sup>V</sup>*<sup>1</sup>

sure applied on the inner surface which is in the order of O(10<sup>0</sup>

−2*ρ*(*V*<sup>1</sup>

Eq. (1) is integrated through the thickness *h* of the shell,

*T* = ∫ *R*1 *R*2 *σθθ dR* = *<sup>ρ</sup> Vr* <sup>2</sup> *<sup>R</sup>* \_\_\_\_\_1 2 ( \_1 <sup>3</sup> <sup>−</sup> <sup>4</sup> *<sup>R</sup>*<sup>1</sup> 3 \_ 3 *R*<sup>2</sup> 3 + *R*1 4 \_ *R*2

− 2*ρ*(*Vr*

<sup>2</sup> *R*1)(ln

where *ρ* is the mass density of the sand shell, *τ<sup>c</sup>*

Bearing in mind that the yield stress, *τ<sup>c</sup>*

in the order of O(10−1).

of factors as well.

106 Granular Materials

thickness of the shell is *R*<sup>2</sup>

frame (*e***<sup>r</sup>**

**Figure 10(a)** and **(b)** plot the variations of circumferential tension in dry and saturated sand shells with the expansion of the shell driven by the detonation of central TNT or HMX, respectively. The parameters are chosen as follows: *ρ* = 2.1 × 10<sup>3</sup> kg/m3 , *Vr*,dry,tnt = 220 m/s, *Vr*,dry,hmx = 330 m/s, *Vr*,saturated,tnt = 380 m/s, *Vr*,saturated,hmx = 420 m/s, *τc*,dry = 15 MPa, *τc*,saturated = 0.5 MPa. The terms with the coefficient involving *η* can reasonably be ignored as a result of the dimensional analysis. The instability onset is identified as the point at which *dT*/*dR*<sup>2</sup> = 0, beyond which the increase of strain does not render the corresponding increase of the pressure or loads. Specifically, the critical radius of dry and saturated sand shells corresponding to the destabilization onset driven by detonation of TNT or HMX are *Rc*,dry,tnt = 75 mm, *Rc*,dry,hmx = 80 mm, *Rc*,saturted,tnt = 98 mm, *Rc*,saturted,hmx = 105 mm, respectively. Clearly, faster detonation velocity of explosive and addition of interstitial fluids can effectively delay the destabilization onset of the inner compact layer, equivalently the initiation of the primary jetting, consistent with the experimental observations. Likewise, we can predict the destabilization onsets of expanding sand shells with varying moisture contents as plotted in **Figure 11**, which agree well with those derived from the experimental observations. Note that the observed destabilization onsets of particle shells were determined from the high-speed photos that show the visible patterns in the surface of charge, which actually occurs after the destabilization onset.

The fragment size following breakup is substantially determined by the wavelength of the most unstable disturbance that has the greatest growth rate. Determination of a dominant unstable wave length is difficult due to the time-varying nature of the mean flow. Louis suggested that for a small value of *Γ*, the most disturbances are in a range of wavelengths between *O*(1) and *O*(1/*Γ*) times the instant thickness of the shell, where *Γ* is the dimensionless number as follows

$$
\Gamma^2(t) = \frac{6\rho \, d\_{e0}^2(t) \, h^2(t)}{\tau\_c}.\tag{3}
$$

**Figure 10.** The variations of circumferential tension in dry and saturated sand shells with the expansion of the shell driven by the detonation of central TNT (a) or HMX (b).

**Figure 11.** Theoretically predicted (curve line) and experimentally observed (red circles) critical radii of expanding sand shells with varying saturation degree driven by the detonation of central TNT.

In Eq. (3), *d*‾‾*θθ* and *h* are average circumferential strain and the instant shell thickness, respectively. Along this line, we give the first order of the estimation of the range of fragment mass of dry sand shells, or equally jet mass *mjet*, as a function of the yield stress as shown in **Figure 12**. The experimental determined jet mass for the dry and moderately wetted sand falls well into the range predicted by the aforementioned model.

**Figure 12.** Primary jet mass vs. yield stress of sand shells [13].

In Eq. (3), *d*‾‾*θθ*

108 Granular Materials

the range predicted by the aforementioned model.

shells with varying saturation degree driven by the detonation of central TNT.

driven by the detonation of central TNT (a) or HMX (b).

and *h* are average circumferential strain and the instant shell thickness, respec-

tively. Along this line, we give the first order of the estimation of the range of fragment mass of dry sand shells, or equally jet mass *mjet*, as a function of the yield stress as shown in **Figure 12**. The experimental determined jet mass for the dry and moderately wetted sand falls well into

**Figure 11.** Theoretically predicted (curve line) and experimentally observed (red circles) critical radii of expanding sand

**Figure 10.** The variations of circumferential tension in dry and saturated sand shells with the expansion of the shell

## **2.4. Secondary particle jetting model: cavitation model based on the expansion of hollow spheres**

A micromechanical approach describing the cavitation process originally applied to ductile damage in solids has been proposed to account for the spallation in a liquid (or melt metal) subjected to a pulsed tensile load [29]. Xue et al. [22] adapted this cavitation-based spallation model to account for the disintegration of the outer particle layer, or equivalently, the formation of the secondary particle jetting, which is initiated by the unloading wave opposed to the tensile loading.

The incipient spallation of the outer particle layer takes the form of the macroscopic dilation in the wake of rarefaction waves. The dependence of the volumetric variation on the pressure is schematically plotted in **Figure 13(a)**. Within the frame of cavitation model, the bulk of the sample is seen as a collection of adjacent hollow spheres of internal and external radii *a*(*t*) and *b*(*t*) (see **Figure 13(b)**), respectively. The initial outer radius of the sphere *b*<sup>0</sup> can be interpreted as the mean half-length between two neighboring nucleation sites as depicted in **Figure 13(b)**. As *b*<sup>0</sup> defines the mass volume involved in the cavitation pattern, the "microscopic" pressure invoked by the cavitation varies with *b*<sup>0</sup> . Since the "microscopic" pressure should agree with the "macroscopic" pressure dictated by the volumetric variation, this compatibility provides a criterion for the determination of *b*<sup>0</sup> . To ensure the expansion of the microscopic hollow sphere is compatible with the dilation of the macroscopic outer particle layer, the microscopic expansion rate of the sphere, *3b/b*, should remain consistent with the macroscopic dilation rate of the particle layer, *V/V*, where *V* is the volume of the outer particle layer. The dilation rates at these two length scales are thereafter detonated by a single parameter *D*.

The spallation or, equally, the dilatation process of the outer layer consists of three stages. The first so-called hollow sphere expansion stage is prescribed by the relaxation of the accumulated pressure when the volumetric increase is dictated by the dilatation rate *D.* During the phase I, voids hardly begin to grow due to the inertial resistance. The end of the phase I of cavitation coincides with the full relaxation of the pressure marked by the restoration of the initial packing density. Afterward, the rapidly expanding matrix progressively becomes gaseous so that the particles interact by collision and the continuous displacement/stress field does not exist. Thus, the matrix and the void of the hollow sphere undergo the independent inertial expansion. The gaseous regime of the matrix is hereinafter detonated as the phase II of cavitation, which sustains as long as the matrix remains diluter than the initial packing state. Examining the packing density of the matrix in the dry sand suggests that the gaseous state of matrix maintains even when the fragmentation starts. By contrast, the gaseous saturated sand is soon transformed to the dense granular flows when the loose particles get recompressed

**Figure 13.** (a) Pressure relaxation experienced by the dry (dashed line) and saturated (solid line) particle layer accompanied by the dilation; (b) schematic of the hollow sphere pattern; (c) the expansion of the individual hollow sphere [22].

by the unconstrained outward expansion of the void. The subsequent expansion of the void, detonated as the phase III, is conditioned by the dense granular flow in the incompressible matrix.

with the macroscopic dilation rate of the particle layer, *V/V*, where *V* is the volume of the outer particle layer. The dilation rates at these two length scales are thereafter detonated by

The spallation or, equally, the dilatation process of the outer layer consists of three stages. The first so-called hollow sphere expansion stage is prescribed by the relaxation of the accumulated pressure when the volumetric increase is dictated by the dilatation rate *D.* During the phase I, voids hardly begin to grow due to the inertial resistance. The end of the phase I of cavitation coincides with the full relaxation of the pressure marked by the restoration of the initial packing density. Afterward, the rapidly expanding matrix progressively becomes gaseous so that the particles interact by collision and the continuous displacement/stress field does not exist. Thus, the matrix and the void of the hollow sphere undergo the independent inertial expansion. The gaseous regime of the matrix is hereinafter detonated as the phase II of cavitation, which sustains as long as the matrix remains diluter than the initial packing state. Examining the packing density of the matrix in the dry sand suggests that the gaseous state of matrix maintains even when the fragmentation starts. By contrast, the gaseous saturated sand is soon transformed to the dense granular flows when the loose particles get recompressed

**Figure 13.** (a) Pressure relaxation experienced by the dry (dashed line) and saturated (solid line) particle layer accompanied by the dilation; (b) schematic of the hollow sphere pattern; (c) the expansion of the individual hollow

a single parameter *D*.

110 Granular Materials

sphere [22].

Analytical modeling of these three sequent phases can be referred to Ref. [22]. This cavitation model estimates that the fragment size or equally the secondary jet size for dry and saturated sand ranges from 4 to 6 mm and 1.6 to 3.3 mm, respectively. Applying the proper fragmentation criterion, the predicted onset of secondary particle jetting occurs at 200–300 μs for the dry sand and 50–100 μs for the saturated sand after the detonation, respectively. The cavitation model is capable of predicting the fragmentation onset and the fragment size consistent with the experimental results. Therefore, cavitation is inferred here to be the most probable spallation mechanism of the outer particle layer.

The size of the secondary jets represented by twice the length between two activated nucleation sites, 2*b*<sup>0</sup> , is dictated by the compatibility of the "microscopic" and macroscopic pressures during the unloading of the compacted particles. Mathematically, smaller *b*<sup>0</sup> in saturated sand is rendered by the significantly elevated dilation rate due to the larger elastic energy and faster moving release waves in the saturated sand. Micromechanically, it is the results of the competition between two neighboring cavities. Analogous to the scenario involving the Mott waves traveling between fractures (see **Figure 14(a)**), the expansion of cavity emanates the compressive waves into the neighborhood so as to suppress the potential cavitation nucleation in the encompassed area. The combined travel length of the compression waves emanating from the neighboring nucleation sites can be taken as the upper limit of the spacing between nucleation sites, namely 2*b*<sup>0</sup> (see **Figure 14(b)**). The unloading duration in saturated sand is almost one order shorter than that in dry sand, leading to the significantly shortened distance between two neighboring cavities.

**Figure 14.** (a) A schematic of the Mott cylinder model with regard to the dynamic fragmentation of the solid cylinder (ring); (b) in particles, compression waves propagate away from an activated nucleation sites (above) retarding any activation of the nucleation sites within the travel radius and collide with those emanating from the adjacent nucleation sites [22].

## **3. Shock-induced particle jetting**

## **3.1. Quasi-two-dimensional particle jetting under moderate impulsive loads: phenomenal description**

It is difficult to visualize the particle jet spread in the spherical or cylindrical experiments of the explosive-driven particle jetting due to the superimposition of jets, obscured by detonation gases, and the very short timescale as well. To overcome these disadvantages, Rodriguez et al. [1, 2] studied the particle jetting in quasi-two-dimensional configurations using moderate pressure loads induced by shock-tube-type facilities connected to a Hele-Shaw cell. With this convenient experimental setup, it is possible to conduct repetitive reliable experiments using a ring of particles in radial expansion trapped in a Hele-Shaw cell as shown in **Figure 15(a)**. More importantly, it is much easier to visualize and distinguish the primary and secondary jets. Xue et al. carried out the experiments of the shock-induced particle jetting using the apparatus similar to that devised by Rodriguez and reported similar observations of the particle jetting process.

**Figure 15(b)** shows the evolutions of dual particle jets of flour ring dispersed by the shock wave with the overpressure of 3.33 bar. The perturbation of the inner surface of ring can be detected at *t* = 1 ms. The primary jets cutting through the inner surface are well defined in

**Figure 15.** (a) Schematic of the quasi-two-dimensional experimental setup for shock induced particle jetting. Insets: photo of the four ring sample (left) and the overpressure histories at the exit of shock tube (right). (b) High speed photos of particle induced particle jetting.

the first several milliseconds. A large number of secondary jets burst out of the outer surface of ring 1.5 ms after the shock front reaches the outer surface. Afterward, the needle-like secondary jets undergo dramatic growth during the following one millisecond, while the tips of primary jets seem to be arrested at the bottom of secondary jets. It takes another several milliseconds that the primary jets overtake the secondary jets.

## **3.2. Particle-scale evolution of shock-induced particle jetting: DEM investigation**

**3. Shock-induced particle jetting**

**description**

112 Granular Materials

of particle induced particle jetting.

**3.1. Quasi-two-dimensional particle jetting under moderate impulsive loads: phenomenal** 

It is difficult to visualize the particle jet spread in the spherical or cylindrical experiments of the explosive-driven particle jetting due to the superimposition of jets, obscured by detonation gases, and the very short timescale as well. To overcome these disadvantages, Rodriguez et al. [1, 2] studied the particle jetting in quasi-two-dimensional configurations using moderate pressure loads induced by shock-tube-type facilities connected to a Hele-Shaw cell. With this convenient experimental setup, it is possible to conduct repetitive reliable experiments using a ring of particles in radial expansion trapped in a Hele-Shaw cell as shown in **Figure 15(a)**. More importantly, it is much easier to visualize and distinguish the primary and secondary jets. Xue et al. carried out the experiments of the shock-induced particle jetting using the apparatus similar to that devised by Rodriguez and reported similar observations of the particle jetting process. **Figure 15(b)** shows the evolutions of dual particle jets of flour ring dispersed by the shock wave with the overpressure of 3.33 bar. The perturbation of the inner surface of ring can be detected at *t* = 1 ms. The primary jets cutting through the inner surface are well defined in

**Figure 15.** (a) Schematic of the quasi-two-dimensional experimental setup for shock induced particle jetting. Insets: photo of the four ring sample (left) and the overpressure histories at the exit of shock tube (right). (b) High speed photos Experimental observations can only provide the configurational evolution of particle ring having no access to the particle-scale information, such as the particle velocities and forces. DEM has proven to be an effective tool to investigate the particle-scale velocity and stress fields in particles subjected to the static or dynamic loadings. Xue et al. performed the DEM simulations of the shock-induced particle jetting using the same geometrical configuration as in the experimental. Parametric studies were carried out to quantify the effect of a variety of variables, including the overpressure of shock loading (*p*<sup>0</sup> ), the inner and outer radii of ring (*Rin* and *Rout*), the packing density (*χ*), and particle size (*dp* ). Details of the simulation can be found in Ref. [25].

**Figure 16** shows the shock dispersal of particle rings in terms of variations in velocity profiles. The shock-loaded particle rings with different initial parameters develop into the resembling jet structures with distinctive features as demonstrated in **Figure 16**. The formation and evolution of the primary jets in all cases, which are barely accurately described using experimental techniques, undergo two distinctive phases, namely the nucleation of the incipient jets and the competitive growth of the incipient jets. Here, the incipient jets are referred to as the localized shear flows or, equivalently, the fast moving particle clusters as shown in the innermost frame in each subfigure of **Figure 16**. The inner surface of ring remains smooth without visible dents or ripples so that the first phase is almost impossible to identify from the experimental observations.

The azimuthal velocity profiles of particle ring in early times shown in **Figure 17** demonstrate the nucleation of the incipient jets. No consistent pattern persists during the first millisecond, the spikes in the azimuthal velocity profile being transient and irregular. The flows behind the shock front are largely homogeneous around the perimeter. The following several milliseconds saw the dramatic transition of azimuthal velocity profile from irregular oscillations to regular fluctuations that are consistent throughout. This transition is clearly manifested by the substantial jump around *t* = 0.5–1 ms in the variations of correlation coefficient of the two sequential azimuthal velocity profile (see **Figure 17(b)**). The peaks indicated in **Figure 17(a)** correspond to the localized shear flows, or equivalently the incipient jets identified in **Figure 17(d)**.

The radial growth of incipient jets in terms of the penetration depth into the bulk and the cross-sectional width is strongly uneven, the strong jets mushrooming outwards opposed to the retarded weak jets. As a result, the substantial elimination and the coalescence of weak jets prevail throughout the second phase. By contrast, the mushroom-like strong jet occasionally would split into multiple subjets, which is more likely to occur in rings with low packing density (see **Figure 16(c)** and **(d)**). Interestingly, the multiplication of strong jets can take place multiple times. The evolutional characteristics of incipient jets revealed by the DEM simulations are substantiated by the experimental observations (see **Figure 18**).

**Figure 16.** Evolutions of the velocity profiles in particle rings with different parameters. Particles are shaded according to the magnitude of velocities. (a) *dp* = 2 mm, *p*<sup>0</sup> = 5 bar, *Rin* = 20 cm, *χ* = 0.55; (b) *dp* = 2 mm, *p*<sup>0</sup> = 5 bar, *Rin* = 35 cm, *χ* = 0.55; (c) *dp* = 2 mm, *p*<sup>0</sup> = 5 bar, *Rin* = 20 cm, *χ* = 0.42; (d) *dp* = 1 mm, *p*<sup>0</sup> = 5 bar, *Rin* = 20 cm, *χ* = 0.45.

The elimination of the weak jets significantly influences the temporal variations of the jet number as shown in **Figure 19**. After the chaotic initiation of incipient jets during the first several milliseconds evidenced by the strong oscillation of jet number, the jet number plummets dramatically in the following 5–10 ms. Afterward, the jet number undergoes much more gradual decrease until the jets are expelled from the outer surface of ring. Taking into account these fundamentals demonstrated in **Figure 19**, a physics-based equation as follows can be derived to describe the temporal variation of jet number, *Njet*.

$$N\_{\rm jet} = N\_{\rm jet}(R\_{\rm nr}, d\_{p'} \chi) - V\_{\rm jet}(p\_{o'} \chi) \Delta t \{h, p\_{o'} \chi\}.\tag{4}$$

In Eq. (4), *Njet,i* represents the number of initial activated incipient jets; *Vjet* represents the decline rate of jet number during phase II (number per unit time); Δ*t* is the duration of phase II. Surprisingly, the overpressure of shock waves does not have the discernible effect on the number of initial jet, *Njet,i*, which instead is a function of the inner radius of ring, *Rin*, the particle diameter, *dp* , and the packing density, *χ*. It suggests that *Njet,i* is indicative of some intrinsic characteristics of particles, analogous to the intrinsic flaws of solids. The decline rate of jet number, *Vjet*, is clearly elevated by stronger shock loadings. Besides, lower packing density seems to hinder the elimination of jets. The duration of phase II, Δ*t*, is among the most important factors governing the jet

**Figure 17.** (a) Azimuthal velocity profiles of particle ring in early times and (b) variations of correlation coefficient of the two sequential azimuthal velocity profile of particle ring. Snapshots of particle ring at *t* = 0.2 ms (c) and *t* = 3 ms (d).

**Figure 18.** High speed photos of shock dispersal of corn quartz sand ring.

The elimination of the weak jets significantly influences the temporal variations of the jet number as shown in **Figure 19**. After the chaotic initiation of incipient jets during the first several milliseconds evidenced by the strong oscillation of jet number, the jet number plummets dramatically in the following 5–10 ms. Afterward, the jet number undergoes much more gradual decrease until the jets are expelled from the outer surface of ring. Taking into account these fundamentals demonstrated in **Figure 19**, a physics-based equation as follows can be

**Figure 16.** Evolutions of the velocity profiles in particle rings with different parameters. Particles are shaded according

= 1 mm, *p*<sup>0</sup>

= 5 bar, *Rin* = 20 cm, *χ* = 0.55; (b) *dp*

, *χ*) − *V*jet(*p*<sup>0</sup>

In Eq. (4), *Njet,i* represents the number of initial activated incipient jets; *Vjet* represents the decline rate of jet number during phase II (number per unit time); Δ*t* is the duration of phase II. Surprisingly, the overpressure of shock waves does not have the discernible effect on the number of initial jet, *Njet,i*, which instead is a function of the inner radius of ring, *Rin*, the particle diameter,

, and the packing density, *χ*. It suggests that *Njet,i* is indicative of some intrinsic characteristics of particles, analogous to the intrinsic flaws of solids. The decline rate of jet number, *Vjet*, is clearly elevated by stronger shock loadings. Besides, lower packing density seems to hinder the elimination of jets. The duration of phase II, Δ*t*, is among the most important factors governing the jet

, *χ*)*Δt*(*h*, *p*<sup>0</sup>

= 5 bar, *Rin* = 20 cm, *χ* = 0.45.

= 2 mm, *p*<sup>0</sup>

, *χ*). (4)

= 5 bar, *Rin* = 35 cm, *χ* = 0.55;

derived to describe the temporal variation of jet number, *Njet*.

= 2 mm, *p*<sup>0</sup>

= 5 bar, *Rin* = 20 cm, *χ* = 0.42; (d) *dp*

*N*jet = *N*jet,*<sup>i</sup>*(*R*in, *dp*

to the magnitude of velocities. (a) *dp*

*dp*

(c) *dp*

114 Granular Materials

= 2 mm, *p*<sup>0</sup>

**Figure 19.** Temporal variations of the jet number with varying overpressure of shock loadings (a) and inner radius of ring (b). The inner radius of ring remains constant, Rin = 20 cm in (a). The overpressure peak in (b) is 5 bar.

number, since the significant increase of jet number either due to the stronger shock loadings or larger inner radius of ring is dominantly caused by the truncated phase II. In another way, there is not enough time for the elimination of jets to fully unfold. The thickness of ring, *h*, the overpressure of shock loadings, *p*<sup>0</sup> , and the packing density, *χ*, are among the parameters influencing Δ*t*.

#### **3.3. Mechanisms governing the shock-induced particle jetting**

The analytical formulation of Eq. (4) entails a thorough understanding of the underlying mechanisms, specifically the formation and elimination mechanisms of incipient jets. With regard to the formation of incipient jets, it is necessary to unlock the transition of the homogeneous flows to the localized shear flows. Unlike solids or liquids, the stress waves in particles travel through particle contact points and are primarily transmitted by the "force chains" that carry most of load in the granular materials [18, 30]. Meanwhile, the shock energy is dissipated by the random particle collisions. Because of the strong energy dissipation and nonlinear characteristics of granular systems, the inter-particle forces are transmitted through heterogeneous architecture of force chains such as shown in **Figure 20**, where the inter-particle contact forces are represented by inter-particle lines scaled with the magnitude of the contact forces. The initial contact network of particles (see the top panel in **Figure 20**) appears to be homogeneous in general with particle-scale heterogeneities. The cylindrical shock loading activates the contacts aligning with the local radial directions. Besides the intricate contact network in the innermost particle layers, a handful of long linear force chains extend radially from the inner surface toward the outer surface (shaded red in the second panel in **Figure 20**). These long linear force chains act as the arteries from which a growing number of short force chains are initiated, forming distinguishable clusters of force chains at *t* = 1 ms with the dimensions much larger than that of constituent particles.

**Figure 20.** Snapshots of the network of force chains in the bottom section of the particle ring subjected to the shock loading of P0 = 20 bar in early times. Force chains (denoted by the thick dashed red lines )at t = 0. 6 ms indicate the long linear force chains acting as the nuclei of the force chain clustering [25].

number, since the significant increase of jet number either due to the stronger shock loadings or larger inner radius of ring is dominantly caused by the truncated phase II. In another way, there is not enough time for the elimination of jets to fully unfold. The thickness of ring, *h*, the overpres-

**Figure 19.** Temporal variations of the jet number with varying overpressure of shock loadings (a) and inner radius of ring

(b). The inner radius of ring remains constant, Rin = 20 cm in (a). The overpressure peak in (b) is 5 bar.

The analytical formulation of Eq. (4) entails a thorough understanding of the underlying mechanisms, specifically the formation and elimination mechanisms of incipient jets. With regard to the formation of incipient jets, it is necessary to unlock the transition of the homogeneous flows to the localized shear flows. Unlike solids or liquids, the stress waves in particles travel through particle contact points and are primarily transmitted by the "force chains" that carry most of load in the granular materials [18, 30]. Meanwhile, the shock energy is dissipated by the random particle collisions. Because of the strong energy dissipation and nonlinear characteristics of granular systems, the inter-particle forces are transmitted through heterogeneous architecture of force chains such as shown in **Figure 20**, where the inter-particle contact forces are represented by inter-particle lines scaled with the magnitude of the contact forces. The initial contact network of particles (see the top panel in **Figure 20**) appears to be homogeneous in general with particle-scale heterogeneities. The cylindrical shock loading activates the contacts aligning with the local radial directions. Besides the intricate contact network in the innermost particle layers, a handful of long linear force chains extend radially from the inner surface toward the outer surface (shaded red in the second panel in **Figure 20**). These long linear force chains act as the arteries from which a growing number of short force chains are initiated, forming distinguishable clusters of force chains at *t* = 1 ms with the dimensions

**3.3. Mechanisms governing the shock-induced particle jetting**

much larger than that of constituent particles.

, and the packing density, *χ*, are among the parameters influencing Δ*t*.

sure of shock loadings, *p*<sup>0</sup>

116 Granular Materials

The variations in the circumferential distributions of strong contact density, *ρ*contact, in early times (see **Figure 21**), demonstrate how the particle-scaled heterogeneities evolve into the macroscale clusters of strong contacts indicated by the contact force peaks with width much larger than the particle size. Note that the agglomeration of force chains is well ahead of the formation of the nonuniform velocity profile that signifies the beginning of the particle clustering. Since the momentum alongside the stresses is being transmitted along the force chains, leaving the particles disconnected from the force chains, there are few chances to obtain the momentum. Particles connected by the strong force chains are supposed to move faster than those cut off from the contact network. Force chains thus act as the main channels of momentum at least in early times as suggested by the strong correlation between the Azimuthal distribution of contact density *ρ*contact and radial velocity *Vr* in the first millisecond as shown in **Figure 21(a)**.

Force chains also play a major role in the elimination of weak jets caused by the dilating strong jets as demonstrated in **Figure 22**. With the incipient jets (composed of the red circles in **Figure 22**) moving ahead of the slow-moving particles (denoted by the greendashed circles in **Figure 22**), velocity differences across the edges of the incipient jets retard any sustained contacts, leading to the weakened lateral confinement imposed on the jets. Therefore, nontrivial transverse flows occur along the edges of jet, the jet front flaring out significantly (see the middle panel in **Figure 22**). The lateral expansion of adjacent jets, especially the jet heads, squeezes the slow-moving particles in between (denoted by blue-dotted circles in **Figure 22**) establishing an intricate network of force chains therein

**Figure 21.** (a) Azimuthal distribution of contact density *ρ*contact and radial velocity *Vr* of particle ring at different times and (b) temporal evolution of correlation coefficient between *ρ*contact and *Vr* .

(see the middle and bottom panels in **Figure 22**). The newly constructed force chains with the dominant transverse orientation hinder the radial transport of the momentum that is instead channeled along the transversely aligned force chains (see the middle panel in **Figure 22**). The growth of the burgeoning minor jets between two major jets is thus likely to be suppressed or even retarded. The minor jets composed of particles indicated by the dotted circle in the middle panel of **Figure 22** are degraded to the slow-moving cluster. With the slow-moving particles increasingly lagging behind, more spaces are left outside the edges of jets, resulting in the intensified transverse flows along the edges. By contrast, the radial compaction leads to the enhanced radial resistance restraining the radial advance of the jet front such as illustrated in the bottom panel of **Figure 22**. At some point, the transverse flows along the edges of jets are expected to overwhelm the radial propagation. The edges of major jets curl outward toward the opposite directions so that the major jet splits into several subjets (indicated by the circles in the bottom panel in **Figure 22**). The subjets with the propagation direction deviating from that of the parental jet would undergo the same development described above until they are expelled from the outer surface.

Given that the suppression of weak jets by the strong jets is mainly responsible for the decrease of jet number, the decline rate of jet number, *Vjet*, decidedly depends on the spatial density of

**Figure 22.** Illustrations of the evolution of the jetting pattern as well as the contact network. The red circles, dashed-line filled circles and blue filled circles represent the fast-moving particles connected by force chains, slow-moving particles without effective contacts among them, and slow-moving particles connected by transversely oriented force chains, respectively [25].

(see the middle and bottom panels in **Figure 22**). The newly constructed force chains with the dominant transverse orientation hinder the radial transport of the momentum that is instead channeled along the transversely aligned force chains (see the middle panel in **Figure 22**). The growth of the burgeoning minor jets between two major jets is thus likely to be suppressed or even retarded. The minor jets composed of particles indicated by the dotted circle in the middle panel of **Figure 22** are degraded to the slow-moving cluster. With the slow-moving particles increasingly lagging behind, more spaces are left outside

.

of particle ring at different times and

**Figure 21.** (a) Azimuthal distribution of contact density *ρ*contact and radial velocity *Vr*

(b) temporal evolution of correlation coefficient between *ρ*contact and *Vr*

118 Granular Materials

incipient jets, the perimeter of ring, and the transverse expansion of strong jets. The average spacing of initial jets varies little with *Rin* and *p*<sup>0</sup> , whereas decreases with decreasing packing density and particle size. The transverse expansion of strong jets strongly correlates with the radial propagation of jets that are driven by the impulsive loadings. Accordingly, stronger shock loading intensifies both the radial and transverse expansion of jets, hastening the suppression of the adjacent weak jets.

**Figure 23** highlights the key events characterizing the formation and competitive growth of incipient jets. An excessive large number of strong force chains extruding into the bulk serves as the nuclei of incipient jets. The jets born earlier or showing stronger shear flows undergo considerable transverse flare up, annihilating the burgeoning weak jets. A substantial portion of initial incipient jets cannot survive the first instants of the phase II.

**Figure 23.** Illustration of key events dominating the formation (left) and elimination (right) of incipient jets.

## **4. Discussion**

Despite the resembling jetting pattern driven by the central explosion and radial shock loadings, the underlying mechanisms are fundamentally different as required by the distinct behaviors of particles subjected to strong blast waves and modest shock waves. In the former case, particle layers are compacted so tightly that they expand like the solids of the constituent materials. Thus, the (primary) particle jetting may well be understood from the continuum perspective since the hydrodynamic instability of interface, such as RT instability, fails to predict the jetting timescale comparable with the experimental data. Bulk fracture of compacted expanding particle layers becomes the promising candidate. The dynamic fragmentation theories of solids may well be applicable to the theoretical model of the explosion-driven particle jetting. But some major alterations need to be made to adapt these theories to the fragmentation of particle assemble. Since particles cannot sustain the tensile stresses nor have the surface energy, the fracture criterion of solids involving these two pivotal variables does not hold in particles. Experimental results suggest that the inception of particle jetting initiates shortly after the propagation of the rarefaction wave. This observation implies that the unloading of particles triggers the particle jetting. The fractures of solids mainly nucleate at the intrinsic flaws that determine the statistics of fragment size. By contrast, the dimension of flaws in particle system, namely the inter-grain pores, contradicts with that of particle jets. It is plausible to assume that the nuclei of the particle jets may well be brought in by a strong shock interaction. The implosion of particles causes dozens of shear bands across the thickness of particle ring with attrited grains [31]. Recent experiments also collected the sintered clumps of aluminum powders after the explosive dispersal of powders [32]. The heterogeneous thermodynamic activities occurring in the blast loaded particles, such as shear banding, should be the focus of the future study. A thorough understanding in this regard needs the rigorous examination of previous experimental data and development of adequate experimental and numerical techniques providing more direct evidences.

Shock-induced particle jetting opens a fundamentally different domain but attracts relatively less attention compared with the explosion-driven particle jetting. This scenario offers an ideal opportunity to look into the transient particle flows. This chapter presents some preliminary investigations into this problem via both experimental and numerical methods. Much more work is needed to clarify the origin of the nuclei of incipient jets and the interplays between jets with varying strengths. In this regard, force chains play an essential role via introducing the inhomogeneity and modulating the jetting pattern.

## **5. Conclusion**

incipient jets, the perimeter of ring, and the transverse expansion of strong jets. The average

density and particle size. The transverse expansion of strong jets strongly correlates with the radial propagation of jets that are driven by the impulsive loadings. Accordingly, stronger shock loading intensifies both the radial and transverse expansion of jets, hastening the sup-

**Figure 23** highlights the key events characterizing the formation and competitive growth of incipient jets. An excessive large number of strong force chains extruding into the bulk serves as the nuclei of incipient jets. The jets born earlier or showing stronger shear flows undergo considerable transverse flare up, annihilating the burgeoning weak jets. A substantial portion

Despite the resembling jetting pattern driven by the central explosion and radial shock loadings, the underlying mechanisms are fundamentally different as required by the distinct behaviors of particles subjected to strong blast waves and modest shock waves. In the former case, particle layers are compacted so tightly that they expand like the solids of the constituent materials. Thus, the (primary) particle jetting may well be understood from the continuum perspective since the hydrodynamic instability of interface, such as RT instability, fails to predict the jetting timescale comparable with the experimental data. Bulk fracture of compacted expanding particle layers becomes the promising candidate. The dynamic fragmentation theories of solids may well be applicable to the theoretical model of the explosion-driven particle jetting. But some major alterations need to be made to adapt these theories to the fragmentation of particle assemble. Since particles cannot sustain the tensile stresses nor have the surface energy, the fracture criterion of solids involving these two pivotal variables does

**Figure 23.** Illustration of key events dominating the formation (left) and elimination (right) of incipient jets.

of initial incipient jets cannot survive the first instants of the phase II.

, whereas decreases with decreasing packing

spacing of initial jets varies little with *Rin* and *p*<sup>0</sup>

pression of the adjacent weak jets.

120 Granular Materials

**4. Discussion**

Both explosion-driven and shock-induced particle jetting exhibit the dual jetting structure, namely, the primary jets initiating from the inner surface and the secondary jets initiating from the outer surface of particle rings/cylinders/shells. The primary and secondary jets have fundamentally different size and occur in different times so that respective mechanisms are required. More importantly, distinct behaviors of particles subjected to strong blast waves and weak shock waves dictate different mechanisms underpinning the particle jetting in both cases. Accordingly, we adopt a continuum approach to model the explosion-driven particle jetting. Specifically, a destabilization model of expanding shell is proposed to account for the onset of the primary jetting. The secondary jetting can be described by a cavitation spallation model based on the expansion of hollow spheres. The timescale and characteristic size of primary/secondary jets predicted by theoretical models agree well with the experimental data. By contrast, the shock-induced particle jetting is studied via the DEM method, which can access the particle-scale information, such as particle velocities and contact forces. The investigation reveals a two-staged evolution of particle (primary) jets, the formation and competitive growth of incipient jets. The formation of incipient jets is characterized by the transition from the homogeneous flows to the localized shear flows. The ensuing evolution of incipient jets is accompanied by the substantial annihilation of weak jets and the multiplication of strong jets. The mechanisms underlying these two phases are found to be closed related with the network of force chains.

## **Author details**

Kun Xue\*, Xiaoliang Shi, Kaiyuan Du and Haoran Cui

\*Address all correspondence to: xuekun@bit.edu.cn

Key State Laboratory of Explosive Science and Technology, Beijing Institute of Technology, Beijing, China

## **References**


[13] Xue K, Li F, Bai C. Explosively driven fragmentation of granular materials. The European Physical Journal E. 2013;**36**(8):1-16

**Author details**

122 Granular Materials

Beijing, China

**References**

2010;**20**(1):41-51

2014;**24**(4):393-402

Kun Xue\*, Xiaoliang Shi, Kaiyuan Du and Haoran Cui

shocked particle ring. Physical Review E. 2014;**90**:043013

wave acceleration. Physical Review E. 2013;**88**(6):063011

Journal of Physics: Conference Series. 2014;**500**(15):152012

AIP Conference Proceedings. 2011;**1426**(1):1615-1618

spherical geometry. Shock Waves. 2014;**24**(5):501-513

ticles. AIP Conference Proceedings. 2011;**1426**(1):1623-1626

persal of solid particles. Physics of Fluids. 2012;**24**(9):091109

Key State Laboratory of Explosive Science and Technology, Beijing Institute of Technology,

[1] Rodriguez V, Saurel R, Jourdan G, Houas L. External front instabilities induced by a

[2] Rodriguez V, Saurel R, Jourdan G, Houas L. Solid-particle jet formation under shock-

[3] Ripley RC, Zhang F. Jetting instability mechanisms of particles from explosive dispersal.

[4] Ripley R, Donahue L, Zhang F. Jetting instabilities of particles from explosive dispersal.

[5] Milne AM, Floyd E, Longbottom AW, Taylor P. Dynamic fragmentation of powders in

[6] Milne A, Parrish C, Worland I. Dynamic fragmentation of blast mitigants. Shock Waves.

[7] Grégoire Y, Sturtzer M-O, Khasainov BA, Veyssière B. Cinematographic investigations of the explosively driven dispersion and ignition of solid particles. Shock Waves.

[8] Gregoire Y, David F, Oren P. Development of instabilities in explosively dispersed par-

[9] Fue-Sang L, Tao X, Fan Z. The role of vorticity and turbulence on the instability of a dense solid particle flow. AIP Conference Proceedings. 2011;**1426**(1):1619-1622

[10] Frost DL, Loiseau J, Marr BJ, Goroshin S. Particle segregation during explosive dispersal of binary particle mixtures. AIP Conference Proceedings. 2017;**1793**(1):120020

[11] Frost, D.L., Gregoire, Y., Goroshin, S., and Zhang, F., "Interfacial instabilities in explosive gas-particle flows," Proceedings of 23rd International Colloquium on the Dynamics

[12] David LF, Yann G, Oren P, Samuel G, Fan Z. Particle jet formation during explosive dis-

of Explosions and Reactive Systems, Univ. of California, Irvine, July 24-29, 2011

\*Address all correspondence to: xuekun@bit.edu.cn


**Solid State of Granular Materials**

[28] Borg JP, Grady D, Cogar JR. Instability and fragmentation of expanding liquid systems.

[29] Huang X, Ling Z, Dai LH. Cavitation instabilities in bulk metallic glasses. International

[30] Tai Q, Sadd MH. A discrete element study of the relationship of fabric to wave propagational behaviours in granular materials. International Journal for Numerical and

[31] Nesterenko VF, Meyers MA, Chen HC. Shear localization in high-strain-rate deforma-

[32] Frost DL, Loiseau J, Goroshin S, Zhang F, Milne A, Longbottom A. Fracture of explosively compacted aluminum particles in a cylinder. AIP Conference Proceedings. 2017;**1793**(1):

International Journal of Impact Engineering. 2001;**26**(1-10):65-76

Journal of Solids and Structures. 2013;**50**(9):1364-1372

Analytical Methods in Geomechanics. 1997;**21**:295-311

120019

124 Granular Materials

tion of granular alumina. Acta Materialia. 1996;**44**(5):2017-2026
