The role of initial speed in projectile impacts into light granular media

Projectile impact into a light granular material composed of expanded polypropylene (EPP) particles is investigated systematically with various impact velocities. Experimentally, the trajectory of an intruder moving inside the granular material is monitored with a recently developed non-invasive microwave radar system. Numerically, discrete element simulations together with coarse-graining techniques are employed to address both dynamics of the intruder and response of the granular bed. Our experimental and numerical results of the intruder dynamics agree with each other quantitatively and are in congruent with existing phenomenological model on granular drag. Stepping further, we explore the ‘microscopic’ origin of granular drag through characterizing the response of granular bed, including density, velocity and kinetic stress fields at the mean-field level. In addition, we find that the dynamics of cavity collapse behind the intruder changes significantly when increasing the initial speed . Moreover, the kinetic pressure ahead of the intruder decays exponentially in the co-moving system of the intruder. Its scaling gives rise to a characteristic length scale, which is in the order of intruder size. This finding is in perfect agreement with the long-scale inertial dissipation type that we find in all cases.

We feel the drag force from air while walking, jogging or riding a vehicle. We also feel the drag force from sand while walking on the beach. The question is to which extent the drag force laws known for Newtonian fluid can be applied to granular materials, such as sand, powders and grains. Such an extension can not only shed light on a fundamental understanding of drag force in non-Newtonian fluids [1][2][3] , but also on widespread applications of granular materials including pile driving in soil mechanics, impact cratering in geo-and astrophysics, and last but not least, robotic locomotion 4,5 .
Drag induced by projectile penetration is an interesting topic with a long history, particularly in association with ballistic impact 6,7 . In the past decades, the interest has been broadened in connection with crater formation [8][9][10][11] , partly motivated by space exploration 12,13 as well as deciphering the history of planet formation 4 . Moreover, it has been extended to the influence of air cavity collapsing as well as granular jet formation 14 , the projectile shape 11,15,16 and the role of gravity on the final penetration depth and time 17 . As an extreme case, infinite penetration of a projectile into light-weighted EPS (expanded polystyrene) particles was observed, when the weight of the projectile exceeds a critical value 18 .
Based on experimental and numerical investigations, an empirical law based on Poncelet model 2, 19 is typically used to describe the penetration curve of a projectile. It includes a 'hydrostatic' term dependent on the penetration depth z 10 and a drag term ∝v 2 with v the instantaneous speed of the projectile with respect to the granular medium. The inertial nature of this drag force has been proposed in the past 10,20,21 . Nevertheless, that is not the complete story, as recent investigations also reveal a variation of the scaling with v for vibrofluidized granular materials 22 . Moreover, collisional damping ∝v 2 can also be understood, as a 'viscous' scale damping process, governed by a speed dependent collisional rate 23,24 and the fluctuations of force chains with a stochastic nature 25,26 . In short, a first-principle theory for describing granular drag is still lacking, partly owing to the challenges in predicting the rheological response of granular media, particularly on the brink of solid-liquid-like transition 27,28 .
From the perspective of applications, particles of expanded polypropylene (EPP) or expanded polystyrene (EPS) are widely used for construction, packaging, damping, as well as raw materials for molded forms, due to its light yet rigid, and good thermal isolation properties. For instance, the drop tower of the European Space Agency in Bremen Germany, uses EPS particles regularly as damping materials to capture drop capsules that weight typically over 400 kg. As such, a better understanding of granular drag, particularly the physical mechanisms behind the empirical laws, is desirable.
Here, we employ a recently developed microwave radar system that 'sees' through granular particles and DEM (discrete element method) simulations 17 to investigate the influence of the impact velocity on the drag force induced by light weighted EPP particles (Neopolen, P9255) with a mean diameter of 3 mm on a spherical projectile with a diameter = D 1 cm. Experimentally, the position of the projectile is determined from the phase shift of 10 GHz electromagnetic waves traveling from the transmission to the receiving antennae with a triangulation process 29 , as sketched in Fig. 1. The experimental conditions, including the dimensions of the intruder, granular particles, the container, as well as the material properties of the particles, are followed in the numerical simulations, which allows us to explore the origin of granular drag through characterizing the response of granular particles at the mean-field level. The results suggest that, despite the order-of-magnitude difference in the kinetic stress and granular energy levels, the profiles of scaled kinetic stress viewed from the perspective of the projectile is surprisingly independent on the initial speed of the intruder v (0) i . In another word, the kinetics of the projectile inside of a granular medium is pre-determined by the initial speed of the intruder v (0) i .

Intruder Dynamics
We perform a systematic study of the intruder penetration dynamics at various initial impact speed (v i = [0-2.3] m/s). Figures 2 and 3 show the vertical position z t ( ) and speed v t ( ) i of the intruder, obtained experimentally and numerically, respectively. For the experimental cases, the initial speed is estimated from a parabolic fit of the intruder trajectory before impact together with the initial falling height H. As noted, in all cases the intruder Figure 1. Schematic of the radar tracking set-up viewing from the side. The radar system compares the relative distance between the transmission and receiving antennae (only one sketched here) to obtain the projectile trajectory inside the EPP particles, which are transparent to EM waves.

Figure 2.
Sample penetration curve z(t); experiment (a) and simulation (b). In both cases, the data is compared with the empirical model Eq. (1). dynamics predicted numerically is in very good quantitative agreement with the experimental results and in line with previous findings 20, 30-32 . In general, the intruder dynamics is predominately conditioned by the initial conditions. For all cases, the intruder penetrates into the medium, slowly decelerates until stops at a certain depth. For low impact speeds, there exists an initial accelerating phase as gravity exceeds granular drag at impact, instead of continuous deceleration as in the cases with high impact velocities. Quantitatively, the intruder movement along the gravity direction can be described using its translational equation of motion.
which includes gravity mg, drag force and depth-dependent frictional terms with γ, κ and λ phenomenological parameters. Following previous findings 20 , the characteristic length scale λ is chosen to be the radius of the container cm 15 , which is sufficiently large to avoid boundary effects. The frictional term is linear and scale-free for small penetration depths and it saturates to κz for large scales 20,31 . In the past, several experimental works showed that the quadratic dependence in speed is sufficient to describe the mechanical energy dissipation of the intruder 20,31 . The energy dissipation from granular drag resembles an object moving in a fluid at high Reynolds number. The numerical solutions of Eq. (1) is also included in Figs. 2 and 3, for comparison in each case.
Nevertheless, we find the coefficient γ (see Fig. 4) is not a material constant, but decays with increasing v (0) i until it saturates at high impact speed. It suggests different mechanical response at low v i , which arises presumably from different perturbation distance in comparison to impact at high v i . Assuming that the coefficient γ can be estimated from the momentum transfer, γ ρ ϕ ∝ v A v i 2 p i 2 with A projected projectile area 31 , there outcomes would suggest that the initial impact reduce notably the packing fraction ϕ of the granular bed. Indeed, this trend was found earlier, examining intruders in fluidized beds, using an up-flow air supply 11 . In our experiment, however, we do not detect significant changes in volume fraction and, accordingly, our outcomes can not be explained following the same arguments. Moreover, we find the dependency of κ on the initial velocity is notably weaker 11  www.nature.com/scientificreports www.nature.com/scientificreports/ (see inset in Fig. 4). It is worth mention, that the intruder's trajectory Eq. (1) has also been described, in terms of its initial kinetic energy and introducing a depth-dependent γ 33,34 .
In Figure 5, we focus on the relation between the final penetration depth z max and v (0) i . As expected, the final depth increases with increasing initial speed (see Fig. 5a). In both experiment and simulations, we find a slightly non-linear response, which is more noticeable for low initial speed. It is important to mention that a linear dependency would be congruent with ∝ z h max 1/2 , which was earlier reported and justified as a liquid-like viscous damping 9 . On the other hand, the scaling ∝ z H max 1/3 was also reported in the past 10,31 and explained as a complex inertial response of granular beds to impact. For a better comparison, in Fig. 5b,c, we examine in log-log scale z max vs h and z max vs = + H z h max , respectively. As noted, our outcome does not match accurately on any of the two options, but both experiments and simulations are closer to the ∝ z H max 1/3 dependency, at least within the uncertainties.

Response of Granular Bed
Coarse-grained fields. Continuum approaches are also useful, addressing the intruder dynamics and the macroscopic response of the granular bed 35 . Here, we use a coarse-graining method [36][37][38][39] , which is a micro-macro mapping technique, to examine the macroscopic responses of the granular bed, as the intruder moves through it. Thus, from the DEM data of each individual particle, we compute the macroscopic fields: volume fraction , ) and kinetic stress σ r z t ( , , ) k . We took advantage of the cylindrical symmetry of the system because the intruder gets into the media by the axis of the cylindrical container. Moreover, no significant deviations from this trajectory were detected. Thus, we average the vertical and radial components for the studied magnitudes within an azimuthal representative volume element with the same size, obtaining the macroscopic fields in terms of the cylindrical coordinates ρ and z.
In Figs. 6 and 7, outcomes are illustrated for two initial speeds of the intruder, , respectively. First, it is noticeable that the penetration process does not perturb significantly the macroscopic density, as the impact is more noticeable in the neighborhood of the intruder. Moreover, the fields allow visualizing the time evolution of the void cavity at different depths, as the intruder moves through the system. As noted, the impacting sphere creates a well-defined cavity, which first increases its size and suddenly collapses at a given time t c . Remarkably, the cavity shapes obtained in both cases are comparable with previous experiments using X-ray radio-graph 14,40,41 .
The second row in Figs. 6 and 7 illustrates the vertical velocity field ρ v z ( , ) z during the penetration process, while the third row in both figures represents the radial velocity field ρ ρ v z ( , ). In general, the intruder notably perturbs the velocity fields of the granular bed, transmitting its mechanical energy through contacting particles over a certain distance. Interestingly, the field ρ v z ( , ) z indicates that the existence of a shear band, which correlates with the development of the void cavity. While the intruder penetrates, it pushes away the particles below it towards the bottom and the radial direction of the container. When the intruder travels deep enough, the particles move back to fill the cavity. Interesting, the color-maps of ρ v z ( , ) z ( Fig. 6 (row II) and Fig. 7 (row II)) also show that, even from the beginning of the process, the intruder notably perturbs the region that is later accessed.  . These findings indicate that the initial impact perturbs the system at large distances, enhancing its local plasticity and favoring further reordering events. As a result, the granular bed changes its macroscopic response, when varying the initial speed. The latter might explain that the coefficient γ (see Eq. (1) intruder's dynamics) is not a material constant, despite the weak changes in macroscopic volume fraction.
Additionally, we also connect the particle local activity to the intruder deceleration, exploring the behavior of the kinetic stress σ ρ z ( , ) k , which is the stress associated with velocity fluctuations (see the fourth row in Figs. 6 and 7). Note that momentum transfer as well as energy loss from inelastic collisions is closely associated with the deceleration trajectory of the intruder. This process is very heterogeneous and happens at a length scale comparable to the size of the intruder. The energy consumption rate due particle-particle collisions correlates with the strength of kinetic pressure Cavity collapse. The computed volume fraction fields allow us to explore nature of void collapse through a quantitative analysis on its dynamics. With that aim, we compute a sequence of φ ρ z ( , ) fields every millisecond. Subsequently, void profile R z t ( , ) is estimated using the definition φ ρ φ > z ( , ) / 2 RCP for granular bulk. This definition is plausible given the value of coarse-graining scale = w r p . Following this procedure, the void profile In the past, the dynamics of the cavity collapse was experimentally accessed using X-Ray Tomography 14,41 . Thus, it was found that the evolution of the cavity radius h t ( ) was consistent with a power law where t c is the collapsing time, and an exponent α = .
0 66 was reported 41 . Moreover, several groups have studied theoretically and experimentally the collapsing of air bubbles in liquids [42][43][44] . Furthermore, in void collapse due hydrostatic pressure, there are theoretical evidences supporting that the cavity radius should diminish as = − h t t t ( ) ( ) c 1/2 , resembling a Rayleigh-type singularity 44,45 . In Fig. 8, the evolution of cavity radius h t ( ) is illustrated for two different impactor speeds. In general, we find that h t ( ) smoothly diminishes in the neighborhood of t c , exhibiting a power law behavior = − α h t t t ( ) ( ) c (see the insets in Fig. 8a,b). This behavior is very similar to the cavity collapse in the wake of an impacting disk on a fluid 45 . Moreover, this power law controlled pitch-off was also encountered when a ball is dropped in sand 41 . We find that the exponent characterizing the collapsing process differs notably depending on the initial speed impact. For the case, 0 42 is estimated, while a noticeable lower value α = . 0 32 is found for = . v (0) 2 3 i m/s. Thus, the collapsing process gets slower as the initial speed gets larger, which might correlates with the amount of momentum transferred, the size of the perturbed zone and its corresponding dissipation time. It is important to remark, that in the case = . v (0) 2 3 i m/s, the last accessible data point always deviates from this trend because it is the result of the movement of individual particles. A more systematic and quantitative analysis of the collapsing dynamics will be a focus of future investigations.

Characteristic length and time scales.
In order to quantify the characteristic length and time scales of collisional energy dissipation, we examine the kinetic pressure profile in the vertical direction. The averaged values p z ( ) k correspond to an average on a radial sector of ∆ = r D, centered at the intruder location. Figure 9 illustrates several results p z ( ) k corresponding to two initial intruder speeds and different times; Fig. 9a,b, = v (0) 0 i and = . v (0) 2 3 i m/s, respectively. For a better comparison, the spatial coordinates correspond to the mobile reference system at the intruder ′ = − z z vt ( ) i and the values of kinetic pressure are re-scaled with their maximum values p k i (see insets of Fig. 9a,b), which is located in the neighborhood of the intruder. In general, the kinetic stress decays exponentially with the distance respect to the intruder location, regardless of the initial intruder velocity and time. Astonishingly, though the values of p z ( ) k differ in one order of magnitude, the characteristic length scale is practically the same in all cases and it is of order of D. In the past, a similar analysis was www.nature.com/scientificreports www.nature.com/scientificreports/ done for a 2D case 46 , examining a scaled steady-state velocity field, which was found to decay exponentially, with a correlation length in the order of the particle size.
During the penetration process of an intruder in a granular bed with macroscopic mass density ρ, dissipation occurs due to inelastic collisions between grains. In order to build a continuous theoretical approach of this process 47 , the starting point is the energy balance equation in terms of the granular temperature T, which written in the moving reference system of the intruder (2) m m T i In Eq. 2, one can use the Fourier's law λ = − ∇ q T T for the heat flux with diffusivity λ T ; while α T T represents the collisional granular dissipation and ε  is the macroscopic shear rate. Moreover, ρ ∇ v T m i accounts for the convective heat transmission.
Interestingly, the data shown in Fig. 9a,b suggest the presence of a very slow diffusive-convective heat transmission scenario. In that conditions ε ≈  0 and Eq. 2 in D 1 reads as, is a solution of Eq. (3), where T o is the granular temperature at the surface of the intruder. Remark that τ and ξ are the time and space characteristic values of the penetration process. Naively, τ is in the order of the penetration time, and ξ , is the positive root of the characteristic equation . From our data and assuming ′ ≅ ′ T z p z ( ) ( ) k , we find that the length scale of the energy dissipation ξ is practically constant ξ ≈ D 2 5 . Interesting, this finding is in perfect agreement with the long-scale inertial dissipation type γ = − F v d i 2 that we find in all cases. It is important to remark, however, that we find that the dissipation parameter γ decreases with increasing initial speed of the intruder. This fact correlates with the initial kinetic pressure that the intruder induces at impact, when mechanical energy is quickly transmitted through the granular bed by means of, e.g., shock waves 48,49 . Thus, right from the beginning the intruder perturbs the strong particle contacts on the force network, and enhances the plasticity of the system. Consequently, the system responses distinctively depending on the initial speed of the intruder. Finally, we note that the existence of a critical length scale ahead of a projectile was also observed in previous investigations on projectile impact on dust agglomerates for understanding the evolution of protoplanetesimal 50,51 . The authors also www.nature.com/scientificreports www.nature.com/scientificreports/ showed, based on dimensional analysis, that the characteristic length scale in the system should be directly related to the projectile size. Stepping further, our results show the length scale remains the same, when changing the initial speed and; i.e., the amount of energy transmitted by projectile to the surrounding particles.

Conclusions and Outlook
To summarize, we study experimentally and numerically the penetration of a spherical intruder in a granular bed at various impact velocities using radar particle tracking in experiments and coarse-graining techniques in simulations. The dynamics of the intruder is numerically reproduced and the obtained trajectories are in very good agreement with a well-accepted phenomenological model. However, the obtained coefficient of the inertial drag term, which is ∝v i 2 , is not a material property, but deviates from a constant value as the initial speed decreases to zero. The latter raises the question of what determines the nature of granular drag. To explain this, we compute the relevant macroscopic fields of the granular bed (volume fraction, macroscopic velocity and the kinetic stress tensor) by means of coarse-graining techniques. Our results show that the intruder perturbs the granular bed even at very long distances, although notable changes in volume fraction are not detected. Further analysis of the cavity collapse behind the intruder reveals a power law diminishes of the neck radius h t ( ), with an exponent depending on the initial impact speed. Moreover, the kinetic pressure profiles ′ p z ( ) k in the co-moving frame of the intruder decays always exponentially with the distance to the intruder ′ z , and the characteristic length is practically independent of the initial speed of the intruder. Analytic arguments show that this result indicates a very slow diffusive-convective collisional-energy transmission within the granular bed. In the future, further investigations on how gravity influences granular drag shall be carried out by means of oblique impacts as well as microgravity experiments.

Method
Experimental set-up. As illustrated in Fig. 1, the transmission horn antenna emits a 10 GHz electromagnetic (EM) wave into the free space, while a spherical projectile moves following an unknown trajectory. As the dielectric constant of styrofoam particles is very close to air, the granular bed is practically transparent to EM waves. Thus, the mobility of the target with respect to each receiving (Rx) antenna is obtained from the phase shift between the received signal and the emitted one. Finally, a coordinate transformation allows reconstructing the trajectory of the projectile accurately. Moreover, given the analog nature of the radar system output, an extremely high temporal resolution is also achieved. The limitation of the temporal signal only arises from the analog-digital (AD) converter. For the normal penetration processes studied in the current investigation, we use an AD converter (NI-DAQ6015) up to 200 kHz. Furthermore, the non-invasive nature and the high temporal resolution, make this technique suitable for monitoring the trajectory of spherical and non-spherical projectiles in three dimensions (3D). More details on the radar system as well as its advantages and disadvantages in comparison to other particle imaging techniques can be found in 29,52 and 53 , respectively.
In the experimental set-up, a steel ball with a radius of = .
R 0 5 cm is held by a styrofoam holder at various initial falling heights H. By gently pulling a thin thread that wrapped around the sphere, we allow the sphere to fall freely with minimized horizontal and angular velocities. The granular sample is composed of EPP particles (Neopolen, P9255) with a slightly ellipsoidal shape and a mean radius of = .
r 1 5 p mm, filled into a cylindrical container made of Styrofoam, which is also transparent to EM waves. The radius of the container is = R 15 c cm and the filling height is 23 cm, sufficiently deep to avoid the influence from the side wall. The radar system starts shortly before the projectile release until its final resting state. It operates with a power of 1 W with the capability of tracking the projectile at a distance of .1 5 m with a field of view of ~30 cm. DEM and coarse-graining methodology. We use discrete element modeling (DEM) 54 , describing a spherical intruder sinking into a granular bed composed of spherical particles. The implementation is a hybrid CPU GPU / algorithm that allows us to efficiently evaluate the dynamics of several hundred thousand particles 55,56 . We initiate each simulation by generating a random packing of polydisperse spheres with radius = . r 1 5 p mm, confined in a cylindrical container of radius R c with packing fraction φ = . ± .
0 63 0 02. The spherical intruder is released from the free granular surface with a given initial velocity, the range of which matches that in the experiments.
The DEM algorithm computes the movement of each particle = … i N 1 for the three translational degrees of freedom, and the rotational movement is described by a quaternion formalism. The interaction force between particle i and particle j, is composed of normal and tangential components Given the particles normal overlap δ n and the tangential relative displacement ξ → , it reads as, ij n n n r n t tr t Here we used a Hertz-Mindlin model 54  www.nature.com/scientificreports www.nature.com/scientificreports/ In all the simulations presented here, the system is composed of = N 180864 particles, the used contact parameters correspond approximately to particles with a Young's modulus = Y 800 kPa ( = G Y/10), density ρ = .
92 0 kg/m p 3 and normal restitution coefficient = . e 0 4 n . The intruder radius is = . × − R m 5 0 10 3 and density ρ = 7800 kg/m p 3 . The particle-wall interaction is modeled using the the same collision parameters used for particle-particle interaction. The parameters are chosen to match the experimental conditions.
Besides the time evolution of the intruder, the simulation data also provides the position and velocities of all the particles of the bed, during the entire impact process. Thus, in order to explore the macroscopic response of the system a coarse-graining methodology is also used [36][37][38] . As a result, the local response of the granular bed, as the intruder moves through it is addressed in detail, through the macroscopic fields of volume fraction φ ρ z ( , ), particle velocity ρ → v z ( , ) and kinetic stress σ ρ z ( , ) k . According to [36][37][38] , the macroscopic mass density of a granular flow, ρ → r ( ), at time t is defined by i N i i i 1 where the → v i represent the velocity of particle i. The macroscopic velocity field → → V r t ( , ) is then defined as the ratio of momentum and density fields, m Additionally, the propagation of the kinetic activity in the particle bed is quantified by values of the kinetic stress tensor, which reads, where → ′ v i is fluctuation of the velocity of particle i, respect to the mean field → → = → − → → ′ v t r v t V r t ( , ) ( ) ( , ) i i . It is generally accepted that the trace of kinetic stress tensor (kinetic pressure) is proportional to the granular temperature. Thus, it represents a measurement of the grain fluctuations respect to the mean velocity field. That's why it can be used to examine the propagation of the kinetic activity as the intruder moves through it.