Abstract
Understanding the response of granular matter to intrusion of solid objects is key to modelling many aspects of behaviour of granular matter, including plastic flow. Here we report a general model for such a quasistatic process. Using a range of experiments, we first show that the relation between the penetration depth and the force resisting it, transiently nonlinear and then linear, is scalable to a universal form. We show that the gradient of the steadystate part, K_{ ϕ }, depends only on the medium’s internal friction angle, ϕ, and that it is nonlinear in μ = tan ϕ, in contrast to an existing conjecture. We further show that the intrusion of any convex solid shape satisfies a modified Archimedes’ law and use this to: relate the zerodepth intercept of the linear part to K_{ ϕ } and the intruder’s crosssection; explain the curve’s nonlinear part in terms of the stagnant zone’s development.
Introduction
The ubiquity of dense granular matter in nature and the important role it plays in human society cannot be overestimated. Its significance focused much attention for centuries, but research into the fundamental science of granular matter underlying its rich behaviour exploded in the past couple of decades. In spite of the intensive activities, both the statics and dynamics of granular matter are not well understood with many remaining open problems^{1,2}.
A canonical problem in the field is the modelling of the penetration dynamics of a large object within a granular material made of much smaller, but macroscopic particles. This problem has been receiving growing attention in recent years for its multiple applications, e.g., locomotion of terrestrial animals^{3,4}, robots working on granular substrates^{5,6,7}, and crater formation in geological and astrophysical fields^{8,9,10}. Studies on this problem also play significant roles in understanding static properties of granular materials^{11,12}, characterising the dynamics of granular flows^{13,14}, and understanding the response of perturbed granular materials, such as shear banding^{15,16}, sinking effect^{17}, and the jamming transition^{18,19,20,21}.
The penetration dynamics is modelled usually using macroscopic laws of drag force exerted on the intruder^{22,23,24}, which describe the drag force as a combination of a hydrostaticlike force F_{ z } and a viscous force F_{v}. The hydrostaticlike force F_{ z } results from the frictional plasticity of the granular matter^{22,25}, and the viscous force F_{v} results from momentum transfer between the grains^{26,27}. To bridge between the developing knowledge of static granular matter and its dynamic flow, we must understand and construct a predictive model for the penetration process under quasistatic conditions. Advancing such an understanding is the main aim of this paper.
A number of experiments demonstrated that the force, F_{ z }, is hydrostaticlike, varying linearly with penetration depth h for submerged or flatbottom intruders^{5,25}. Based on this observation, an empirical resistive force theory (RFT)^{3} has been proposed for the resistance force on objects intruding granular media quasistatically at different orientations. This approach proved adequate to model the kinematics of slowmoving locomotors^{3,6,28}. It also agrees with a local friction force model (LFFM), recently proposed in ref. ^{23}
In this relation, dS is an infinitesimal area element pointing normal to the intruder surface, P = ρ_{s}gh is the hydrostaticlike pressure proportional to the packing density ρ_{s} of the granular matter, K_{ ϕ } is a coefficient (>1) proportional to the internal friction coefficient μ via a fitting parameter k.
Both RFT and LFFM are based on the assumption that the friction of the granular tangential flow against the intruder’s surface contribute negligibly to the resistance force. Askari and Kamrin^{29} then showed by FEM simulations, in which the granular medium is modelled as a continuum of bulk density ρ_{s} and internal friction coefficient μ, that these features are a consequence of two properties: cohesionless and a friction yield criterion.
Yet, existing models for the quasistatic resistance force have two main deficiencies: they are mainly phenomenological, thus providing little understanding of the physics of the penetration process, and the parameters governing the force magnitude need to be determined experimentally, undermining the model’s predictive power. It is essential to develop a physicsbased quantitative model for the quasistatic resistance force in cohesionless dry granular matter.
Here, we address this issue both theoretically and experimentally. First, we carry out experiments to measure the resistance force on cylinders penetrating vertically into five different granular materials. By varying the bottom area and shape of the intruders, we find a unified dimensionless pressure–depth curve. This curve consists of an initial nonlinear segment, extending to a similar depth (when properly scaled) for all experiments, corresponding to a compression of the granular material ahead of the advancing intruder, which gives rise to formation of a stagnant zone (SZ). This transient is followed by an extended linear region, whose gradient depends only on the medium’s internal friction angle ϕ. This region suggests that the medium yields in a fluidlike manner. We then analyse the process, treating the granular matter as a continuum at a critical state, characterised by parameters ρ_{s} and ϕ. We conjecture that the initial penetration process gives rise to a conical SZ below the cylinder, caused by shear jamming^{5,29}. Using the Mohr–Coulomb (MC) yield criterion and the method of characteristics, we model the steady state of the conical SZ and obtain theoretically an explicit expression for the curve’s gradient. This gradient is found to be independent of the intruder’s geometry and it provides the coefficient K_{ ϕ } in Eq. (1). However, we find that this coefficient is not linear in μ = tan ϕ, casting doubt on Eq. (1). This leads us to conclude that the quasistatic force on any cylinder, whether its bottom is flat, conical or half spherical, can be described universally by Archimedes’ law, albeit scaled by the factor K_{ ϕ } > 1, which we calculate. We then show that the Archimedes’ law can explain the force–depth relation, including the dependence of the constant term on medium nature and intruder crosssection area, as well as the effective steadystate size of the SZ.
Results
Experiments
The experimental apparatus is sketched in Fig. 1a. A cylindrical intruder was connected to a servocontrolled beam through a force sensor that records the total axial force opposing the downward motion. The intruder is pushed into a barrel of diameter 45 cm, filled with granular material up to a height of 21 cm. Throughout the process, the cylinder was kept well away from the barrel’s boundaries, eliminating boundary effects^{30,31}. For generality, several granular materials were used: three kinds of glass beads, made of different particle size distributions, dry quartz sand, and millet. Unlike the glass beads, the quartz sand grains are much more angular and irregular, while the millet grains are softer and less spherical. The physical parameters of these materials are shown in Fig. 2f, where the internal friction angle ϕ refers to the angle of repose, accurate to ±2°.
In the first set of experiments, we used four aluminium cylinders of different radii, R = 15, 20, 25, and 35 mm, but the same length L = 50 mm. The penetration depth, h, was limited to be below L to avoid both end effects at the top of the cylinder and jamming effects close to the bottom of the barrel. To ensure velocity independent behaviour, the penetration velocity was restricted to below a critical value^{32}, \(v_{\rm c} = \sqrt {2gd_{\rm g}} {{\rm /}}10\), where g is gravity acceleration and d_{g} is the mean diameter of granules. This was tested by measuring the force on a R = 20 mm cylinder penetrating into sand at velocities ranging from 10 mm min^{−1} to 300 mm min^{−1}. Indeed, Fig. 1b verifies that the resultant force–depth curves are hardly different. The velocities were then limited to below 30 mm min^{−1} in all experiments.
The experiments reported in refs. ^{5,23} demonstrated that the friction between the intruder surface and grains is negligible relative to the force required to push material below the cylinder out of the way. This allowed us to define a mean pressure, p_{ u } = F/S, where F is the measured vertical resistance force and S is the cylinder’s crosssection area. To nondimensionalise our results, as well as for reasons to be understood below, we scaled the mean pressure, \(\tilde p_u = p_u{{\rm /}}\left( {\rho _{\rm s}gR} \right)\), and the penetration depth, \(\tilde h = h{{\rm /}}R\).
Figure 2a–e shows the dependence of \(\tilde p_u\) on \(\tilde h\) for the five granular materials, when penetrated by the four different cylinders. For each medium, the curves collapse onto a master curve for all the different intruders. This result agrees with the surfacelevel superposition rules postulated by RFT^{3} and LFFM^{23}.
All the \(\tilde p_u  \tilde h\) curves consist of an initial nonlinear regime, \(\tilde p_u = f_1\left( {\tilde h \le \tilde h_0} \right)\), corresponding to material compression ahead of the intruder, which culminates in a fully formed rigid SZ (see below), followed by a much longer linear regime, in which \(\tilde p_u\sim \tilde h\). The nonlineartolinear crossover was determined, by inspection, to be at \(\tilde h_0 = 0.15 \pm 0.06\). Identifying \(\tilde p_u\left( {\tilde h_0} \right)\) as the crossover pressure to the steady state, we define the gradient of the linear region beyond this point as K_{ ϕ }. The dimensionless pressure–depth relation is then
Similar behaviours were also observed in previous resistance force measurements^{5,24,27,33,34,35}. Particle image velocimetry measurements^{5,36,37} and discrete element method simulations^{38} showed that, during the transient nonlinear regime, a conical SZ forms under the flat bottom of intruder segment, driven by a local shear jamming process^{19}. This region then advances as a rigid cone ahead of the intruder^{39,40}. While the cone formation process is complex and not fully understood^{5,40}, this analysis provides a way to model it in terms of the intruder’s crosssection area and the location of the crossover depth to steady state (see below).
The conical SZ advances ahead of the intruder at the same downward speed, parting the medium and wedging matter sideways^{5,8,34}. This means that the quasistatic resistance force on the intruder is governed by the intergranular contact friction at the cone advancing surface^{23}. It follows that the slope of the linear part, K_{ ϕ }, should depend on ϕ. In our experiments, the value of K_{ ϕ } ranges very widely for different granular materials: from about 20 to more than 100.
To test whether K_{ ϕ } depends on intruder properties or is a pure material constitutive parameter, we carried out a suite of penetration experiments with flatbottom prisms of four crosssection shapes: square, equilateral triangle, rectangle, and rightangled triangle, shown in Fig. 3a. The prisms were made of aluminium and L = 120 mm long. The granular medium was the type 1 glass beads. The depth and pressure scaling depend on the intruder’s crosssection area, S and we define an equivalent radius \(R_{\rm e} = \sqrt {S{{\rm /}}\pi }\). This scaling indeed leads to a collapse of the steady state of all the pressure–depth curves onto a unique master curve, regardless of crosssection shape, as shown in Fig. 3b. Moreover, this curve is identical to the one for the cylindrical intruder. This establishes that K_{ ϕ } is indeed independent of the shape and size of the intruder and, therefore, that it is a constitutive characteristic of the granular material^{23,34}. In the next section, we present a theoretical model to predict this parameter. Moreover, there is also a good collapse in the nonlinear regime when the aspect ratio of the crosssections is not too different than 1 (see Discussion below).
Theoretical analysis
Quasistatically yielding cohesionless dry granular media can be modelled as a continuum of bulk density ρ_{s} and internal friction angle ϕ^{29}. To model the stress field, the mechanical equilibrium conditions are closed by a yield criterion, for which there are several models. We choose here to use the MC criterion^{41}, \(\left {\tau {{\rm /}}\sigma _{\rm n}} \right = \mu \equiv {{\rm tan}}{\kern 1pt} \phi\), where τ is the local shear stress and σ_{n} is the corresponding normal stress. Several experiments established that the penetration process gives rise to a conical SZ ahead of the intruder, advancing as a rigid body^{36,37,38,42}. Taking this into consideration, the boundary conditions are sketched in Fig. 4a and include the downward pressure at depth h and the resistive force, F_{ z }, on the intruder. We assume that the free surface, left behind, is flat, which is a good assumption as long as any deviation from flatness is appreciably smaller than the size of the SZ. Supplementary Fig. 1 supports this assumption.
Expressing the stress components (Supplementary Fig. 2) in terms of the mean of the major and minor principal stresses, σ_{0}, and the angle between the major principal stress and the radial direction, α, gives a hyperbolic set of equations in these two variables. The equations can be solved by characteristics (Supplementary Note 1 and Fig. 3), which are paths in the plane along which they reduce to two ordinary differential equations^{43}. The solution yields two families of characteristics, whose local gradients are
Conveniently, the angle between two + and − characteristic curves, passing at any one point in the rz plane, is always 2β. Focusing on the steady state, when the SZ is fully formed, and noting that α = −π/2 everywhere on the cone surface, we can solve for σ_{0} along a characteristic path:
where 0 ≤ η ≤ 1 parameterises the cone surface and runs from its apex to its base (see Fig. 4b) and A(η, ϕ) is derived in the Methods section. This analysis is in good agreement with the observations in refs. ^{36,37,38,42}.
Obtaining σ_{ z } from σ_{0} and integrating it along the cone surface to obtain F_{ z } (Supplementary Note 2), we finally identify the coefficient of h:
This result is significant. It shows that: (i) K_{ ϕ } depends only on ϕ and on no property of the intruder; (ii) K_{ ϕ } is not linear in tan ϕ, in contrast to the proposed relation (1)^{23}. In Fig. 5a, we plot the values of K_{ ϕ }, computed numerically from Eq. (5), as a function of ϕ and compare them with the experimental measurements for the five granular media. The agreement between the model prediction and the experimental values is excellent. From the data reported in ref. ^{23} for glass beads with ϕ = 22°, we calculate K_{ ϕ } = 13.36. This also agrees very well with their measured result, which is equivalent to K_{ ϕ } = 14 ± 2. There is also a good agreement with ref. ^{5}, in which force–depth data are reported for intrusion of ‘foot’ objects into assemblies of poppy seeds at various volume fractions. We can compare their results with our data for cylindrical intruders into millet, the internal angle of which should be close to that of poppy seeds. From their data, we can estimate our \(\tilde h\), \(\tilde h_0\), and \(\tilde p_u\). Translating our \(\tilde h_0\) to their notations, using the different normalisation methods, we find that it corresponds to their δ = 3.8 mm, which is well within the range shown in the inset in their Fig. 3b. Specifically, using their fit in that inset, we get that this value corresponds to a packing fraction 0.600 ± 0.002. Comparing then the slope of their linear steadystate region in that packing fraction with our millet data, we find that it corresponds to our K_{ ϕ } = 60 ± 2. Using now our relations (5) and (2), we obtain that internal friction angle of their medium should be 32.3° ± 0.3°, in excellent agreement both with the literature value and with the value we measured for our millet, 32° ± 2°. The wide range of values of ϕ in our tests covers most common granular materials, implying that our theoretical model has a broad application.
Equation (5) makes it also possible to derive the relation between K_{ ϕ } and μ, which is shown in the inset of Fig. 5a. This relation is clearly superlinear, raising questions about the assumption of its linearity in the LFFM. The quality of the agreement between model and experiments also supports the prediction of Eq. (5) that K_{ ϕ } is a constitutive property of the granular medium alone, an issue explored extensively in penetration experiments^{5,23,34,37}.
Modified Archimedes’ law in granular matter
Within the LFFM^{23}, the resistance force arises from local forces acting normal to the intruder surface. Defining \( \widehat {\bf{z}}\) as the direction of gravity, with \(\widehat {\bf{z}}\) a unit vector in the z direction, we determine the total resistance force on a convex intruder, advancing slowly in this direction, by using Eqs. (1) and (5), as well as Gauss law for inversion from a surface to a volume integral:
with V the volume of the displaced granular materials. This expression is equivalent to Archimedes’ law^{44} in fluid mechanics, in which K_{ ϕ } = 1. The orderofmagnitude larger value of K_{ ϕ } in granular media is due to an effective interaction of the intruder with the entire cyan region shown in Fig. 4b.
For flatbottom intruders, Eq. (2) includes \(\tilde p_0\), a constant term also observed by others^{23,27}. As we show below, this is because the volume V in relation (6) should include the developing volume of the SZ, which gives rise to the initial nonlinear regime. We calculate this volume explicitly below.
For the nonflat intrusions (see Fig. 6), the quasistatic resistance force may be estimated directly from Eq. (6)^{23}. We test four aluminium cones with head angles 2θ of 52°, 60°, 62.3°, and 67.8°, all to ±0.1°, and four aluminium spheres with radii of 15, 20, 25, and 50 mm, all to 0.01 mm, in two granular materials (glass beads 2–3). The experimental measurements are plotted in Fig. 7 and compared with the corresponding theoretical calculation of the resistance force from Eq. (6). The excellent agreement between all the theoretical and experimental results has two significant implications. One is that the quasistatic resistance force on a moving intruder satisfies Eq. (6). The other is that the agreement for all intruder shapes within a given medium with the same value of K_{ ϕ } establishes that this parameter is a pure constitutive property of the granular medium. These conclusions can be used to derive theoretically the constant term in Eq. (2), which we proceed to do next.
Prediction of \(\tilde p_0\) and the SZ growth
Establishing the validity of the modified Archimedes’ law (6) for any intruder convex shape, allows us to derive the value of \(\tilde p_0\), as well as the dynamic development of the SZ ahead of flatbottom intruders. Consider an intruding prism of crosssection S, with a convexshaped bottom of volume V_{0} and extension h_{0} in the penetration direction. For example, for a cylinder of radius \(R{{\rm = }}\sqrt {S{{\rm /}}\pi }\) with either a halfspherical bottom or a cone of head angle 2θ, (h_{0}, V_{0}) = (R, 2πR^{3}/3) or (R/tanθ, πR^{3}/(3 tanθ)), respectively.
The intrusion volume is
where h is measured from the medium’s surface to the bottom of the linear part of the prism or cylinder. Substituting into the modified Archimedes’ law and rewriting it in terms of the normalised depth, \(\tilde h = \left( {\sqrt {\pi {{\rm /}}S} } \right)h\), and pressure, \(\tilde p_u = \sqrt \pi F_z{{\rm /}}\left( {\rho _{\rm s}gS^{3/2}} \right)\), we obtain
Comparing to Eq. (2), the first term on the righthand side is exactly \(\tilde p_0\), which appears to depend on both the medium, through K_{ ϕ }, and the intruder shape parameters. Testing this result against measurements of conical and spherical intruders, we find an excellent agreement with the theory, as shown in Fig. 7.
Next, we apply this result to model the dynamic development of the SZ ahead of flatbottom intruders. As discussed above, such intruders develop SZs that propagate ahead of them, parting the medium. The SZ buildup starts when the intruder enters the medium and ends once the SZ has reached an established steadystate shape. This dynamic process, illustrated in Fig. 6c, is the cause for the initial nonlinear dependence of \(\tilde p_u\) on \(\tilde h\).
To understand better this dynamic process, we model the developing SZ as an effective cone or polyhedron, whose apex is at \(\tilde H = \left( {\sqrt {\pi {{\rm /}}S} } \right)H\) ahead of the intruder. Initially, H = 0 and it reaches a steadystate value, H = H_{ss} at the end of the nonlinear regime. The volume of this zone is V_{0} = SH/3 and, using Eq. (8), we have
Fitting now a functional form to f_{1}(h) in Eq. (2), we can express \(\tilde p_u\) in terms of \(\tilde h\) in this regime, which provides the dependence of H on depth. Since the insertion of the intruder is at constant speed, this also provides the time evolution of the SZ.
In particular, we can find its steadystate size, \(\tilde H_{{\rm ss}}\), for each granular material by substituting the crossover point values: \(\tilde h = 0.15\) and \(\tilde p_u\left( {\tilde h = 0.15} \right)\) (see Fig. 2). We find \(\tilde H_{{\rm ss}}\) = 1.0 ± 0.1, 1.6 ± 0.3, 1.5 ± 0.2, 0.7 ± 0.3, and 3.1 ± 0.8 for glass beads types 1–3, sand, and millet, respectively. The error bars are the standard deviation of the measurements of \(\tilde p_0\) in the different experiments.
Discussion
In summary, we have proposed and tested a general model for the quasistatic penetration dynamics of convex solid objects into granular media. This is a key problem in granular science, relevant to diverse phenomena, such as animal and robotic locomotion^{3,4,5,6,7}, crater formation^{8,9,10}, granular flow dynamics^{13,14}, and granular response to external loading^{15,16,17,18,19,20,21}. A cone penetration test is also the standard for characterising soils.
First, by measuring experimentally the resistance forces on flatbottom cylindrical and prismatic intruders, we find that all the force–depth curves can be generally collapsed onto a dimensionless pressure–depth master curve. This curve crosses over from a short initial nonlinear regime to a long steadystate linear one. Next, assuming a solid granular SZ, propagating ahead of the intruder, we construct a parameterfree model for the steadystate rate of the linear increase and show that this rate, K_{ ϕ }, depends only on the internal friction angle, ϕ, and is therefore a constitutive property of the granular medium. Our calculation fits all our observations, as well as all the reports we could find in the literature, which provided data we could compare with. Our results also call into question a conjecture in the literature that the gradient of the linear part is proportional to the internal friction coefficient^{23}.
Then, using noncylindrical intruders, cones and spheres, we show that the force–depth curves of all convex intruders satisfy a modified Archimedes’ law in that the resistance force is proportional to the hydrostaticlike pressure, the volume of the intruding object, and a proportionality coefficient K_{ ϕ }. The value of the latter is the only real difference between Archimedes’ law for ordinary fluids and dense granular matter. Significantly, the value of K_{ ϕ } can be calculated theoretically, which renders our derived relation parameterfree! The validation of the model for the noncylindrical intruders allowed us to derive the force–depth relation explicitly in terms of the medium’s constitutive property and the intruder’s geometry.
Specifically, in our model, the initial nonlinearity is a result of the growth process of the SZ, while the linear part corresponds to the SZ having reached a steadystate size, after which the effective volume of the intruder and the SZ does not change. This also explains the dependence of the gradient of the linear part only on the medium, since the resistance force is dominated by the normal force between the plastically flowing medium and the jammed SZ, both of which are made from the same material. This we have validated successfully by measurements on materials with internal frictional angles ranging from 22° to 36.5°.
Using this picture, we solve for the stress in the plastic zone around the SZ, using, for simplicity, the MC yield criterion model. This solution made it possible to calculate the total force on the intruder and the gradient K_{ ϕ } explicitly as a function of ϕ. The solution agrees excellently with all the available experimental observations, ours and others’. Indeed, resistance force measurements on prisms of different crosssection shapes, confirm that K_{ ϕ } is a constitutive property of the granular material alone.
The penetration may be affected strongly by the presence of container boundaries^{30,31,33,34}, in particular by the side walls parallel to the intruder’s motion. Our theoretical model indicates that these effects arise from the interference between the stress field in the plastic zone and the walls. Theoretically, for avoiding such effects, the minimum container radius (MCR) is the horizontal extent of the longest characteristic, which is the green line in Fig. 4b, the length of which depends on R and ϕ. In Supplementary Table 1, we list the calculated MCR values for the five granular materials and four cylinder radii in our experiment. For some of the thick cylinders, it is larger than the container size we used. Nevertheless, the excellent collapse for all cases suggest no boundary effects. We conclude that the theoretical calculation, based on the MC criterion, provides an upper bound for the MCR. This is discussed in detail in Supplementary Note 3.
The modified Archimedes’ law also allows us to model quantitatively the SZ growth and predict its effective steadystate size. As a case study, we considered a flatbottom cylinder, intruding at a constant speed below the critical value. On intruding, the volumes of both the immersed cylinder and a growing SZ, V_{sz}, increase with depth. This gives rise to the initial nonlinear regime, in which the total volume increase is V(h) = Sh + V_{sz}(h). V_{sz} increases slower than linearly until it reaches a steadystate value. Generally, assuming that the SZ is a cone throughout its growth (or pyramidlike for noncylinders) of height H(h), we have V_{sz}(h) = SH(h)/3, independent of the intruder’s crosssection shape. The inset in Fig. 3b displays a magnification of the nonlinear regimes for the different crosssections. These all appear to collapse to a master curve, except when the crosssection aspect ratio is considerably different than one (see below).
Based on this insight, we used the nonlineartolinear crossover to estimate the steadystate height of the SZ. The errors of the estimates, for the different media, ranged from 10% to 40% and could stem from: (i) measurement errors of ϕ, which gives rise to errors in K_{ ϕ } and \(\tilde p_0\); (ii) error in estimating the nonlineartolinear transition point \(\tilde h_0\); (iii) an inaccurate yield criterion in the initial nonlinear regime. The latter may also be the reason for an overestimate in the horizontal extent of the response region in the granular medium (cyan in Fig. 4). We also note that the SZ may be blunted at its edges and apex, making our estimates of its steadystate volume an upper bound.
To conclude, we have derived a parameterfree model for the penetration of solid objects into granular media, showing that the behaviour can be fully understood in terms of a modified Archimedes’ law. This understanding should shed more light on a wide range of practical issues, including the standard cone penetration test in soil mechanics in civil engineering, locomotion of animals or robots in sand, digging on Earth and extraterrestrially, and crater formation, to name a few. It is also useful for furthering fundamental understanding of plastic flow of granular matter, e.g., in shear banding, the jamming transition, as well as bridging between the current different models for low and highvelocity penetration regimes.
Methods
Experimental methods
Samples preparation: Before each penetration experiment, the granular medium was stirred by hand vigorously for 60 s. Then the free surface was flattened carefully by passing a trowel across it. Finally, two bubble levels were placed perpendicularly to one another to determine flatness of the surface. The same person prepared all samples to minimised personal variability. The penetration speed was controlled by a servo system to an accuracy of ±0.5% and displacement precision 0.001 mm. The resistance force was measured by a sensor to accuracy of ±0.5% and recorded by a software at sampling frequency of 20 Hz. The intruder was positioned about 3 mm above the flat surface and accelerated to reach a constant speed before it touched the granular surface. For flatbottomed intruder, recording of the vertical resistance force and displacement started the instant the axial force exceeded 0.03 N, which was defined as the beginning of penetration. After checking that the data was independent of velocity (see Fig. 1), all the data was collected at an intrusion velocity of 30 mm min^{−1}. Each experiment was carried out on a number of freshly prepared samples. The sample sizes for the experiments shown in Fig. 2, were: 13, 26, and 27, for the glass beads types 1–3, respectively, 26 for the sand, and 21 for the millet. For Fig. 3, we measured the intrusion force on three samples for each crosssection. The error bars represent the standard deviation over the samples and the fits qualities were: R^{2} >0.997, 0.992, 0.998, and 0.993, respectively, for the square, rectangle, equilateral triangle, and rightangle triangle. Similarly, the measurements on the noncylindrical intruders, shown in Fig. 7, were taken on three samples for each of the four experiments, with the error bars in the figure also representing the standard deviation over the samples. The fits qualities in these cases were: R^{2} >0.991, 0.993, 0.993, and 0.995, respectively, for: (i) cones with four different head angles penetrating into glass beads type 2; (ii) same as (i) but into glass beads type 3; (iii) spheres of four different radii penetrating into glass beads type 2; (iv) same as (iii) but into glass beads type 3.
Data processing
The experimental value of K_{ ϕ } was determined by scaling the force–depth curves into a dimensionless form first and then using the firstorder polyfit command in MATLAB to leastsquares fit the linear slope. Figure 5 shows the results, with the points and error bars representing the mean and spread values of the slopes for the specific granular materials.
Computation of A(η, ϕ) of Eq. (5)
The characteristic paths, shown in Fig. 4b, consist of a curved, DE, and linear, CD, parts, with the latter’s slope being tanβ. Defining the distance from the cone apex along the path OB as η ∈ [0, 1], we have \(l_0(\eta ) \equiv \overline {BE} = \frac{{(1  \eta )R}}{{{{\rm sin}}{\kern 1pt} \beta }}\) and \(l_\psi (\eta ) \equiv \overline {BJ} = l_0(\eta ){{\rm e}}^{\psi {\kern 1pt} {{\rm tan}}{\kern 1pt} \phi }\) (see Fig. 4). The positions of points C, D, E are:
Focusing, for example, on the “+” family of Eq. (3), we have
and integrating Eq. (11) along the characteristic curve CDE from C to E gives
The rightmost term in this relation can be calculated noting that: (i) the position of any point along the curve DE can be expressed in terms of variable ψ:
and (ii) ψ ∈ [0, π/2] between points D and E:
Substituting this into Eq. (12), we obtain finally for the mean stress at point E
where \(A(\eta ,\phi ) \equiv \left( {\frac{{r_C^{1 + {{\rm tan}}^{{\rm 2}}\beta }}}{{r_D^{{{\rm tan}}^{{\rm 2}}\beta }r_E}}} \right)^{{{\rm sin}}{\kern 1pt} \phi }{{\rm e}}^{{{\rm sin}}{\kern 1pt} \phi {\kern 1pt} {{\rm tan}}{\kern 1pt} \beta Z}\) and \(\sigma _{0,C} = \frac{{\rho _{\rm s}gh}}{{1  {{\rm sin}}{\kern 1pt} \phi }}\), corresponding to the stress boundary condition at point C.
Data availability
The data are available from the corresponding author on request.
References
 1.
Nagel, S. R. Experimental softmatter science. Rev. Mod. Phys. 89, 025002 (2017).
 2.
Collins, A. L. et al. The effect of rod nose shape on the internal flow fields during the ballistic penetration of sand. Int. J. Impact Eng. 38, 951–963 (2011).
 3.
Li, C., Zhang, T. & Goldman, D. I. A terradynamics of legged locomotion on granular media. Science 339, 1408–1412 (2013).
 4.
Dickinson, M. H. et al. How animals move: an integrative view. Science 288, 100–106 (2000).
 5.
Aguilar, J. & Goldman, D. I. Robophysical study of jumping dynamics on granular media. Nat. Phys. 12, 278–283 (2016).
 6.
Maladen, R. D., Ding, Y., Li, C. & Goldman, D. I. Undulatory swimming in sand: subsurface locomotion of the sandfish lizard. Science 325, 314–318 (2009).
 7.
Hunt, M. L. Robotic walking in the real world. Science 339, 1389–1390 (2013).
 8.
Goldman, D. I. & Umbanhowar, P. Scaling and dynamics of sphere and disk impact into granular media. Phys. Rev. E 77, 021308 (2008).
 9.
Housen, K. R. & Holsapple, K. A. Impact cratering on porous asteroids. Icarus 163, 102–119 (2003).
 10.
Zhao, R., Zhang, Q., Tjugito, H. & Cheng, X. Granular impact cratering by liquid drops: understanding raindrop imprints through an analogy to asteroid strikes. Proc. Natl Acad. Sci. USA 112, 342–347 (2015).
 11.
Kamrin, K. & Koval, G. Nonlocal constitutive relation for steady granular flow. Phys. Rev. Lett. 108, 178301 (2012).
 12.
Henann, D. L. & Kamrin, K. Continuum modeling of secondary rheology in dense granular materials. Phys. Rev. Lett. 113, 178001 (2014).
 13.
Jop, P., Forterre, Y. & Pouliquen, O. A constitutive law for dense granular flows. Nature 441, 727–730 (2006).
 14.
Gnoli, A., Lasanta, A., Sarracino, A. & Puglisi, A. Unified rheology of vibrofluidized dry granular media: from slow dense flows to fast gaslike regimes. Sci. Rep. 6, 38604 (2016).
 15.
Singh, A., Magnanimo, V., Saitoh, K. & Luding, S. Effect of cohesion on shear banding in quasistatic granular materials. Phys. Rev. E 90, 022202 (2014).
 16.
Luding, S. Constitutive relations for the shear band evolution in granular matter under large strain. Particuology 6, 501–505 (2008).
 17.
Huang, L., Ran, X. & Blumenfeld, R. Vertical dynamics of a horizontally oscillating active object in a twodimensional granular medium. Phys. Rev. E 94, 062906 (2016).
 18.
Liu, A. J. & Nagel, S. R. Nonlinear dynamics: Jamming is not just cool any more. Nature 396, 21–22 (1998).
 19.
Bi, D., Zhang, J., Chakraborty, B. & Behringer, R. Jamming by shear. Nature 480, 355–358 (2011).
 20.
Peters, I. R., Majumdar, S. & Jaeger, H. M. Direct observation of dynamic shear jamming in dense suspensions. Nature 532, 214–217 (2016).
 21.
Vinutha, H. & Sastry, S. Disentangling the role of structure and friction in shear jamming. Nat. Phys. 12, 578–583 (2016).
 22.
Katsuragi, H. & Durian, D. J. Unified force law for granular impact cratering. Nat. Phys. 3, 420–423 (2007).
 23.
Brzinski, T. III, Mayor, P. & Durian, D. Depthdependent resistance of granular media to vertical penetration. Phys. Rev. Lett. 111, 168002 (2013).
 24.
Le Bouil, A., Amon, A., McNamara, S. & Crassous, J. Emergence of cooperativity in plasticity of soft glassy materials. Phys. Rev. Lett. 112, 246001 (2014).
 25.
Tsimring, L. & Volfson, D. Modeling of impact cratering in granular media. Powders Grains 2, 1215–1223 (2005).
 26.
Bester, C. S. & Behringer, R. P. Collisional model of energy dissipation in threedimensional granular impact. Phys. Rev. E 95, 032906 (2017).
 27.
Clark, A. H., Kondic, L. & Behringer, R. P. Particle scale dynamics in granular impact. Phys. Rev. Lett. 109, 238302 (2012).
 28.
Ding, Y., Gravish, N. & Goldman, D. I. Drag induced lift in granular media. Phys. Rev. Lett. 106, 028001 (2011).
 29.
Askari, H. & Kamrin, K. Intrusion rheology in grains and other flowable materials. Nat. Mater. 15, 1274–1279 (2016).
 30.
Nelson, E., Katsuragi, H., Mayor, P. & Durian, D. J. Projectile interactions in granular impact cratering. Phys. Rev. Lett. 101, 068001 (2008).
 31.
Seguin, A., Bertho, Y. & Gondret, P. Influence of confinement on granular penetration by impact. Phys. Rev. E 78, 010301 (2008).
 32.
Albert, R., Pfeifer, M., Barabási, A.L. & Schiffer, P. Slow drag in a granular medium. Phys. Rev. Lett. 82, 205 (1999).
 33.
Stone, M. et al. Local jamming via penetration of a granular medium. Phys. Rev. E 70, 041301 (2004).
 34.
Stone, M. B. et al. Stress propagation: getting to the bottom of a granular medium. Nature 427, 503–504 (2004).
 35.
Pang, Y. & Liu, C. Continuum description for the characteristic resistance sensed by a cylinder colliding against granular medium. Sci. China Phys. Mech. Astron. 56, 1428–1436 (2013).
 36.
Hamm, E., Tapia, F. & Melo, F. Dynamics of shear bands in a dense granular material forced by a slowly moving rigid body. Phys. Rev. E 84, 041304 (2011).
 37.
Murthy, T. G., Gnanamanickam, E. & Chandrasekar, S. Deformation field in indentation of a granular ensemble. Phys. Rev. E 85, 061306 (2012).
 38.
Viswanathan, K., Mahato, A., Murthy, T. G., Koziara, T. & Chandrasekar, S. Kinematic flow patterns in slow deformation of a dense granular material. Granul. Matter 17, 553–565 (2015).
 39.
Waitukaitis, S. R. & Jaeger, H. M. Impactactivated solidification of dense suspensions via dynamic jamming fronts. Nature 487, 205–209 (2012).
 40.
Han, E., Peters, I. R. & Jaeger, H. M. Highspeed ultrasound imaging in dense suspensions reveals impactactivated solidification due to dynamic shear jamming. Nat. Commun. 7, 12243 (2016).
 41.
Michalowski, R. An estimate of the influence of soil weight on bearing capacity using limit analysis. Soils Found. 37, 57–64 (1997).
 42.
Tapia, F., Espndola, D., Hamm, E. & Melo, F. Effect of packing fraction on shear band formation in a granular material forced by a penetrometer. Phys. Rev. E 87, 014201 (2013).
 43.
Sokolovski, V. Statics of Soil Media (Butterworths Scientific Publications, London, 1960).
 44.
Munson, B. R., Young, D. F. & Okiishi, T. H. Fundamentals of Fluid Mechanics, Vol. 3 (John Wiley & Sons, New York, 1990)
Acknowledgements
C.L., W.K. and Y.F. acknowledge gratefully the support from National Science Foundation of China under project 11472011. R.B. is grateful for the hospitality of the State Key Laboratory of Turbulence and Complex System, College of Engineering, Peking University, during the work on this project.
Author information
Affiliations
Contributions
C.L. conceived the study and supervised the project. W.K. and Y.F. conceived the experiments and the numerical modelling. C.L., R.B. and W.K. carried out the analysis and wrote the paper jointly.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Kang, W., Feng, Y., Liu, C. et al. Archimedes’ law explains penetration of solids into granular media. Nat Commun 9, 1101 (2018). https://doi.org/10.1038/s41467018033443
Received:
Accepted:
Published:
Further reading

Constant speed penetration into granular materials: drag forces from the quasistatic to inertial regime
Granular Matter (2021)

Projectile oblique impact on granular media: penetration depth and dynamic process
Granular Matter (2021)

The role of initial speed in projectile impacts into light granular media
Scientific Reports (2020)

Stagnant zone formation in a 2D bed of circular and elongated grains under penetration
Granular Matter (2020)

Discrete Element Analyses of a Realisticshaped Rock Block Impacting Against a Soil Buffering Layer
Rock Mechanics and Rock Engineering (2020)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.