Drag force scaling for penetration into granular media
DDrag force scaling for penetration into granular media
Hiroaki Katsuragi , and Douglas J. Durian Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA Department of Earth and Environmental Sciences, Nagoya University, Nagoya 464-8601, Japan (Dated: September 20, 2018)Impact dynamics is measured for spherical and cylindrical projectiles of many different densitiesdropped onto a variety non-cohesive granular media. The results are analyzed in terms of thematerial-dependent scaling of the inertial and frictional drag contributions to the total stoppingforce. The inertial drag force scales similar to that in fluids, except that it depends on the internalfriction coefficient. The frictional drag force scales as the square-root of the density of granularmedium and projectile, and hence cannot be explained by the combination of granular hydrostaticpressure and Coulomb friction law. The combined results provide an explanation for the previously-observed penetration depth scaling.
PACS numbers: 45.70.Ht, 45.70.Cc, 83.80.Fg, 89.75.Da
Collision is one of the most fundamental processes innature, and can be exploited to uncover the basic physicsof systems ranging from planets to elementary particles.Impact of projectiles into a pack of grains has been of in-creasing interest for highlighting and elucidating the un-usual mechanical properties of granular materials. Thisincludes studies of crater morphology [1–3], penetrationdepth [4–6], dynamics [7–18], boundary effects [16, 19–23] and packing-fraction effects [15, 23, 24]. One findingis that the penetration depth scales as d ∼ D p / H / where D p is projectile diameter and H is the total dropdistance [4, 5]. Granular impact is also important formilitary and industrial applications [25, 26], and conepenetration tests are used for in-situ soil characteriza-tion [27, 28].In a previous study [13] we measured the impact dy-namics of a 2.54 cm steel ball onto a packing of glassbeads. The stopping force was found to be the sumof an inertial drag, proportional to the square of thespeed v , and a frictional drag, proportional to depth z .This has since been supported by several other experi-ments [14, 16–18, 21–23, 29]. The equation of motionduring impact is thus ma = − mg + mv /d + k | z | , (1)where m is projectile mass, g = 980 cm/s , and d and k are materials parameters expressed with units ofa length and a spring constant, respectively. By using ma = d K/ d z and an integrating factor, as in Ref. [10],this can be solved for speed versus depth: v v = e − zd − kd zmv + (cid:18) gd v + kd mv (cid:19) (cid:16) − e − zd (cid:17) , (2)where v is the impact speed at z = 0. The final pene-tration depth d is then given by the limit of v → dd = 1 + 2 mgkd + W (cid:18) mv − mgd − kd kd e mg/kd (cid:19) , (3) where W ( x ) is the Lambert W -function. An additionalconstant stopping force f [6, 18] can be included in theseexpressions by replacing g with g − f /m .Some open questions are how d and k scale with thematerials properties of the projectile and the granularpacking, and how this conspires to give d ∼ D / p H / .For inertial drag, exactly as for hydrodynamic drag athigh Reynolds number Re, momentum transfer givesthe expectation mv /d ∼ Aρ g v where A is the pro-jected projectile area and ρ g is the mass density of thegranular medium. For frictional drag, the combina-tion of hydrostatic pressure and Coulomb friction gives k | z | ∼ µgρ g A | z | where the internal frictional coefficientis µ = tan θ r and θ r is the angle of repose. The scalingof d and k would thus be d /D p ∼ ρ p /ρ g , (4) kD p /mg ∼ µρ g /ρ p , (5)where D p and ρ p correspond to diameter and densityof projectile, respectively. Here, we test the speed versusdepth prediction of Eq. (2), and we compare fitted valuesof d and k with the above expectations. While neitherturns out to be quite correct, we connect the results tothe observed penetration depth scaling.Our basic experimental setup is identical with previouswork [13]. The velocity of the projectile v ( t ) at time t iscomputed by particle-image velocimetry applied to finestripes on a vertical rod glued to the top of the projectile.The system has 20 µ s temporal resolution and 100 nmspatial resolution, which is fine enough to compute posi-tion and acceleration from v ( t ). The primary differencefrom Ref. [13] is that now we vary both the projectiles andthe granular medium. We begin with 0.35 mm diameterglass beads, prepared to a reproducible random packingstate by slowly turning off a fluidizing up-flow of air. Intothis we drop a wide variety of projectiles, listed in Ta-ble I, from free-fall heights between 0 to 85 cm. The firsttype is spheres. The second type is aluminum cylinders,which are dropped lengthwise with the axis horizontal a r X i v : . [ c ond - m a t . s o f t ] M a y Projectile ρ p (g/cm ) D p (cm) L p (cm)Hollow PP ball 0.51 2.54 -Wood ball 0.95 2.54 -Delrin ball 1.65 2.54 -PTFE ball 2.46 2.54 -Steel ball 7.96 - 159 0.3175 - 5.08 -Tungsten carbide ball 15.3 2.54 -Aluminum cylinder 2.89 - 4.26 0.635 -1.27 5.08 - 15.24TABLE I: Projectile properties. The steel sphere diametersare D p = 5 .
08, 2 .
54, 1 .
27, 0 . . D p = 1 .
27 and 0 .
635 cm; eachhave lengths L p = 5 .
08, 10 .
16, and 15 .
24 cm. The density ρ p is projectile plus rod mass divided by projectile volume. and parallel to the surface of the granular medium. Inboth cases the effective density is given by projectile plusrod mass divided by projectile volume. The length of thecylinders is varied, but we find that it does not affect thedynamics or final penetration depths.Example speed versus position data are plotted inFig. 1 for D p = 2 .
54 cm diameter projectiles of four differ-ent densities, dropped onto the glass bead packing, withinitial impact speeds ranging from zero to 400 cm/s. Forslow impacts, the speed first increases and then decreaseswith depth. For faster impacts, the speed versus depthcurves gradually change from concave down to concaveup. Generally, there is a rapid decrease of speed to zero atthe final penetration depth. We obtain good fits to thesecomplex trajectories by adjusting k and d , as shown inFig. 1. The displayed level of agreement is typical for allprojectiles, including the cylinders.Unfortunately a good simultaneous fit for a given pro-jectile to a single pair of k and d values at all impactspeeds can be obtained only for denser projectiles. Forthe less dense projectiles, the fitting parameters becomeconstant only at high impact speeds. Then the samevalues for k and d are obtained as from the analysismethod of Ref. [13]. This holds roughly for d > D p / ρ p > ρ g as shown in Fig. 2. We speculate that forlow density projectiles and small impact speeds, the pen-etration can be shallow enough that the detailed shapeof projectile must be taken into account [30]. In thisregime surface flows and surface roughness could playsa role, though identical penetration behavior was foundearlier for slick and tacky projectiles of the same size anddensity [4]. We note that including a constant force f as a third free parameter does not noticeably change thefits; the largest fit value is f /m = 2 . , which issmall compared to g . We note, too, that f cannot bechosen such that k and d become constant. For the restof the paper, we restrict attention to conditions where thededuced k and d values do not depend on impact speedand hence can be considered as materials parameters.We now compare the fitting parameters with Eqs (4-5).The expectation for d is based on momentum transfer, -400-2000 v [ c m / s ] -8-6-4-20 z [cm] -400-2000 v [ c m / s ] -400-2000 v [ c m / s ] -400-2000 v [ c m / s ] (a)(b)(c)(d)Wood DelrinPTFE Steel FIG. 1: (color online) Example velocity versus depth datafor D p = 2 .
54 cm diameter spheres. The impact begins at z = 0, and proceeds downward in the − z direction. Theblack dashed curves are a simultaneous fit to Eq. (2), wherea single pair of k and d values is found for different initialimpact speeds. The grey dotted curves are also fits to Eq. (2),but where k and d are adjusted for each impact speed andhence do not represent well-defined materials parameters. so that the inertial drag force is mv /d ∼ Aρ g v just likean object moving in a fluid at high Reynolds number. Forspheres and cylinders, the characteristic length is thus d ∼ m/Aρ g and can be written as d ∼ D (cid:48) p ρ p /ρ g if wetake D (cid:48) p to be 1 times diameter for spheres and 3 π/ m/A ). For a unified analysis, we thereforedivide the fitted value of d by D (cid:48) p and plot versus ρ p inFig. 3a. We find that all data, including the cylinderdata, collapse nicely onto a power-law with the expected d / D p (cid:108) p / (cid:108) g FIG. 2: (color online) Scatter plot of penetration depth versusprojectile density, scaled respectively by projectile diameterand grain density. An open circle is plotted for conditionssuch that the k and d fitting parameters are constant; across is plotted otherwise. Blue is for glass beads; green is forrice; pink is for beach sand; yellow is for sugar. slope of one; i.e. d/D (cid:48) p ∼ ρ p holds as per expectation.For the fitting parameter k that sets the quasi-staticfriction force, k | z | , we now make a similar compari-son with expectation by plotting kD (cid:48) p /mg versus ρ p inFig. 3b. All the data, including the cylinder, again col-lapse nicely to a power law in projectile density. However,the expected power-law kD p /mg ∼ /ρ p is clearly wrong.The data are instead consistent with kD p /mg ∼ / √ ρ p .Therefore, the nature of the quasistatic frictional drag isdifferent from the simple combination of Coulomb fric-tion and hydrostatic pressure.To investigate this further, we perform a second seriesof experiments where D p =2.54 cm diameter steel spheresare dropped into rice, beach sand, and sugar (with ma-terials properties listed in Table II). As in Fig. 1, speedvs position data for a range of impact speeds are fit toEq. (2) to obtain values of d and k . Based on Fig. 3,and assuming that the x-axis of this figure is correctlymade dimensionless by dividing out the bulk density ρ g of the granular medium, the observed scalings so far are d /D (cid:48) p ∼ ( ρ p /ρ g ) and kD (cid:48) p /mg ∼ ( ρ p /ρ g ) − / . There-fore, we divide out this density dependence and plot it inFig. 4 versus the only remaining material property of thegrains: the internal friction coefficient µ = tan θ r where θ r is the draining angle of repose. In this figure the datafor glass beads, from Fig. 3, all lie at µ = 0 .
45. The rangeof µ values is less than a factor of two, but to withinuncertainty the scaled d and k parameters collapse topower-laws in µ . For the quasistatic frictional drag co-efficient we find k ∼ µ , which is expected for Coulombfriction. For the speed-squared inertial drag coefficient,we find d ∼ /µ . Therefore the inertial drag force isproportional to µ , which could correspond to an added-mass effect whereby the volume of grains boosted to theprojectile speed grows in proportion to µ . This is unlikethe case of simple fluids, where the fluid flow and the k D p ' / m g (cid:108) p [g/cm ] -1/2-1 d / D p ' (a)(b) FIG. 3: (color online) (a) Dimensionless inertial drag lengthscale d /D (cid:48) p , and (b) dimensionless quastistatic drag coeffi-cient kD (cid:48) p /mg , versus projectile density, for impact into glassbeads. Here D (cid:48) p is the effective projectile diameter; m is itstotal mass, including rod; ρ p is m divided by projectile vol-ume; g = 980 cm/s . Open symbols are for spheres made ofhollow PP ( (cid:66) ), wood ( (cid:67) ), PTFE ( (cid:53) ), delrin ( (cid:52) ), steel ( (cid:13) ),and tungsten carbide ( (cid:3) ); closed diamonds are for aluminumcylinders. Symbol sizes increase with D (cid:48) p , with values givenin Table I. The solid gray lines denote power-laws as labeled.The dashed line in (b) is the expected scaling, kD (cid:48) p /mg ∼ ρ − p .Granular Material ρ g (g/cm ) θ r Grain size (mm)Glass beads 1.45 24 ◦ ◦ × ◦ ◦ ρ g is bulk density; θ r is drainage angle of repose; size is the range of grain diameters,except for rice where it is the length of short and long axes. inertial drag force at high Re depend only on the densityof the fluid.As an alternative analysis for the µ dependence of d ,one could imagine that an inertial drag force ∼ ρ g v A loads the contacts and gives an additional friction forceof µ times this loading. Then the total velocity-squaredforce is mv /d ∼ (1+ αµ ) ρ g v A , which gives d ∝ / (1+ αµ ). This reasonably fits the data, as shown in Fig. 4awith α = 2 . ± .
6. The residuals are smaller for thepower-law form.Combining the power-laws in Figs. 3-4, and the ac-tual numerical prefactors, the inertial and frictional drag ( k D p ' / m g )( (cid:108) p / (cid:108) g ) / µ ( d / D p ' )( (cid:108) g / (cid:108) p ) -1 (a)(b) FIG. 4: (color online). (a) ( d /D (cid:48) p )( ρ g /ρ p ) and (b)( kD (cid:48) p /mg )( ρ p /ρ g ) / , versus µ = tan θ r , where θ r is the re-pose angle of the granular medium. The symbols at µ = 0 . µ are for rice (green), beachsand (pink), and sugar (yellow). The gray lines indicatepower-law behavior. The dashed curve in (a) is a fit to d ∝ / (1 + αµ ), giving α = 2 . ± . coefficients are altogether found by measurements for arange of projectiles and grains to be consistent with d /D (cid:48) p = (0 . /µ )( ρ p /ρ g ) , (6) kD (cid:48) p /mg = 12 µ ( ρ g /ρ p ) / . (7)The two main differences from the simple expectation ofEqs. (4-5), are the factor of 1 /µ in d and the densityratio exponent of 1/2 rather than 1 in k . These resultsmay be inspected differently by re-writing the equationof motion as ma = − mg + 2 . µρ g v A + 8 . µ ( ρ p ρ g ) / g | z | A. (8)Note that the numerical prefactor for the depth-dependent frictional drag is significantly larger thanunity, and that √ ρ p ρ g is larger than ρ g for dense pro-jectiles. For both reasons, frictional drag exceeds thevalue expected from hydrostatic pressure and Coulombfriction. One might have expected the opposite effect,either by a decrease in contact area between projectileand grains due to ejection of grains or by fluidization ofthe grains from the motion of the projectile. Our re-sults instead appear to indicate that frictional contactsare loaded by the motion of the projectile, so that themedium is stronger than set by gravity alone. Such be-havior is not seen for the fast horizontal rotation of bars d [ c m ] µ -1 ( (cid:108) p / (cid:108) g ) D p ' H [cm] H D p (cid:108) p (cid:108) g µ Geometric mean 1/3 scaling
FIG. 5: (color online). Final penetration depth scaling plot.All trials data are shown and they reasobaly agree with theempirical scaling d = 0 . µ − ( ρ p /ρ g ) / D (cid:48) / p H / (dottedline) [4, 5]. Colored curves show the force law predictions.They are also close to data and the empirical scaling. [29], where the depth- and speed-dependent forces areeasily disentangled, and warrants further attention.As a final test we now compare data for the final pen-etration depth d with Eq. (8). Prior observations [4, 5]are consistent with the empirical form d = 0 . µ − ( ρ p /ρ g ) / D (cid:48) / p H / . (9)Thus in Fig. 5 we plot penetration depth data for alltrials versus the quantity ( ρ p /ρ g ) / D (cid:48) / p H / from theright-hand side of this expression. This collapses our newdata to within the experimental uncertainty, includingthat for the cylinder, to the line y = 0 . x and henceshows agreement with prior work. However Eqs. (3,8)do not predict perfect power-law behavior. To comparewith data, we first calculate the geometric mean of eachof the five variables H , D p , ρ p , ρ g , and µ over the rangeof experimental conditions employed here. The predictedpenetration depth for these mean values is shown by asingle red open circle in Fig. 5. Then we vary each of vari-ables, one at a time, keeping all others fixed at their meanvalue. The resulting five predicted penetration depthcurves are included in Fig. 5. They are not identical,or even perfect power laws, but are all close together andin fair agreement with the data. Thus the empirical pen-etration depth data scaling is satisfactorily understoodin terms of the nature of the stopping force exerted bythe medium onto the projectile.In conclusion, we developed an exact solution of Eq. (1)for the dynamics of penetration, and we conducted awide range of granular impact experiments to elucidatethe materials dependence of the inertial and frictionalcontributions to the stopping force. The final equationof motion, Eq. (8), is roughly consistent with the em-pirical scaling of penetration depth versus drop heightand materials parameters. However it is not consistentwith the apparently simplistic expectation of Eqs. (4-5).So there is unanticipated physics, yet to be understood,which could possibly arise from motion-loading of fric-tional contacts or from granular flow fields that dependon the internal friction coefficient.This work was supported by the JSPS PostdoctoralFellowships for Research Abroad (HK) and the NationalScience Foundation through Grant Nos. DMR-0704147and DMR-1305199. [1] J. C. Amato and R. E. Williams, Am. J. Phys. , 141(1998).[2] A. M. Walsh, K. E. Holloway, P. Habdas, and J. R.de Bruyn, Phys. Rev. Lett. , 104301 (2003).[3] X. J. Zheng, Z. T. Wang, and Z. G. Qiu, Eur. Phys. J. E , 321 (2004).[4] J. S. Uehara, M. A. Ambroso, R. P. Ojha, and D. J.Durian, Phys. Rev. Lett. , 194301 (2003).[5] M. A. Ambroso, C. R. Santore, A. R. Abate, and D. J.Durian, Phys. Rev. E , 051305 (2005).[6] J. R. de Rruyn and A. M. Walsh, Can. J. Phys. , 439(2004).[7] K. E. Daniels, J. E. Coppock, and R. P. Behringer, Chaos , S4 (2004).[8] M. P. Ciamarra, A. H. Lara, A. T. Lee, D. I. Goldman,I. Vishik, and H. L. Swinney, Phys. Rev. Lett. , 194301(2004).[9] D. Lohse, R. Rauh´e, R. Bergmann, and D. van der Meer,Nature , 689 (2004).[10] M. A. Ambroso, R. D. Kamien, and D. J. Durian, Phys.Rev. E , 041305 (2005).[11] M. Hou, Z. Peng, R. Liu, K. Lu, and C. K. Chan, Phys.Rev. E , 062301 (2005).[12] L. S. Tsimring and D. Volfson, in Powders and Grains2005 , edited by R. Garcia-Rojo, H. Herrmann, and S. Mc-Namara (Balkema, Rotterdam, 2005), p. 1215.[13] H. Katsuragi and D. J. Durian, Nature Phys. , 420 (2007).[14] D. I. Goldman and P. B. Umbanhowar, Phys. Rev. E ,021308 (2008).[15] P. B. Umbanhowar and D. I. Goldman, Phys. Rev. E ,010301 (2010).[16] F. Pacheco-V´azquez, G. A. Caballero-Robledo, J. M.Solano-Altamirano, E. Altshuler, A. J. Batista-Leyva,and J. C. Ruiz-Su´arez, Phys. Rev. Lett. , 218001(2011).[17] A. H. Clark, L. Kondic, and R. P. Behringer, Phys. Rev.Lett. , 238302 (2012).[18] A. H. Clark and R. P. Behringer, Europhys. Lett. ,64001 (2013).[19] J. F. Boudet, Y. Amarouchene, and H. Kellay, Phys. Rev.Lett. , 158001 (2006).[20] E. L. Nelson, H. Katsuragi, P. Mayor, and D. J. Durian,Phys. Rev. Lett. , 068001 (2008).[21] A. Seguin, Y. Bertho, and P. Gondret, Phys. Rev. E ,010301 (2008).[22] S. von Kann, S. Joubaud, G. A. Caballero-Robledo,D. Lohse, and D. van der Meer, Phys. Rev. E , 041306(2010).[23] H. Torres, A. Gonzalez, G. Sanchez-Colina, J. C. Drake,and E. Altshuler, Rev. Cub. Fis. , 1E45 (2012).[24] J. R. Royer, B. Conyers, E. I. Corwin, P. J. Eng, andH. M. Jaeger, Europhys. Lett. , 28008 (2011).[25] D. J. Roddy, R. O. Pepin, and R. B. Merril, eds., Impactand explosion cratering (Pergamon Press, NY, 1978).[26] J. A. Zukas, ed.,
High velocity impact dynamics (Wiley,NY, 1990).[27] H. S. Yu and J. K. Mitchell, J. Geotech. Eng. , 140(1998).[28] M. Toiya, J. Hettinga, and W. Losert, Granular Matter , 323 (2007).[29] T. A. Brzinski and D. J. Durian, Soft Matter , 3038(2010).[30] K. A. Newhall and D. J. Durian, Phys. Rev. E68