Packing-dependent granular friction exerted on a rod withdrawn from a granular layer: the role of shear jamming

Resistance force acting on a rod vertically withdrawn from a granular layer is studied by experiments and numerical simulations. The initial packing fraction of the granular layer is controlled to evaluate the packing dependence of the resistance force. In both experiments and numerical simulations, the frictional resistance force in the steady slipping regime increases as the initial packing fraction is increased. According to the numerical results, the principal reason for the increase of frictional resistance is the increase of normal force acting on the rod. This large normal force stems from the effective solidification of the granular layer which is caused by the shearing induced by the withdrawn rod. Based on this shear-induced solidification, the experimental data are analyzed using a simple ‘cylindrical shear-jammed zone’ model. Then, the critical divergence of the size of solidified (shear-jammed) zone is observed by increasing the initial packing fraction towards the jamming packing fraction.


Introduction
Granular matter is usually defined by a collection of dissipative solid particles [1,2]. The principal sources of dissipation are inelastic collisions and friction. The former becomes dominant in dilute regime while the latter governs the dynamics in dense flow regime. The constituent particles of dense granular layer are basically solid but the bulk of dense granular matter sometimes behaves like dissipative fluid. Frictional property of bulk granular matter determines its flowability. Therefore, granular friction relates to a lot of natural phenomena. Examples include snow avalanches, seismic slip of a fault gouge layer, landsliding, and so on. Since the fault gouge and snow particles can be regarded as granular matter, granular friction has to be properly understood to model their slipping behaviors. Namely, the reasonable understanding of granular friction could provide useful information for understanding physical mechanisims governing these disasters that are often very fatal for us. In order to mitigate the disasters by these natural hazards, granular friction is one of the most crucial key factors to be understood. However, a lot of fundamental features on granular friction have been left unsolved.
The friction is usually characterized by two kinds of frictional coefficients: μ s and μ. Let us consider the frictional surface at which the tangential force F and normal force N are applied. When F/N exceeds the static frictional coefficient μ s , the slipping is triggered. On the other hand, the force balance F=μN is satisfied in the steady slipping state. The values of μ s and μ are basically assumed as constant and independent of the area of frictional surface and sliding speed. Although this Coulomb-Amonton frictional constitutive law is simple and useful for a rough estimate of the frictional effect, its applicable range is actually limited. When considering the friction in bulk granular layer, its complex frictional nature cannot be simply explained by Coulomb-Amonton law. To characterize granular friction, useful control parameters for quantifying the state of sheared granular matter have to be introduced.
In the previous studies, the state of sheared granular layer is characterized using a dimensionless number called inertial number(see e.g. [3]). By investigating the relation between the inertial number and frictional coefficient, a constitutive law for granular friction has been obtained. The inertial number I is defined as, , where, D P , , g ġ , and ρ g are the shear strain rate, particle diameter, confining pressure, and density of granular matter, respectively. Using I, some types of frictional constitutive laws for dense granular flow have been proposed [3][4][5][6].
The inertial number I is undoubtfully the most important parameter to characterize granular friction. However, only I is not sufficient to understand the complex behaviors of granular friction. For example, geometry dependence of granular friction has been reported in [7]. Packing fraction f of granular matter also plays an important role in granular friction [8][9][10]. In general, mechanical properties of dense granular matter strongly depends on f [11][12][13][14][15][16]. However, the detail relation between granular friction and packing fraction has not yet been solved, especially in the vicinity of the jamming point. The concept of jamming [17] could be a crucial factor to characterize the dense granular friction. Jamming is a very useful idea to unify the physics of random solidification of not only granular matter but also various soft materials such as gels, glasses, and foams, etc. Particularly, the shear-induced jamming [18,19] must be a key concept to discuss the granular friction because the frictional resistance occurs under the shearing. Nevertheless, the relation between granular friction and shear-induced jamming has not yet been studied so far. In addition, mechanical perturbations such as vibration and air fluidization could affect the granular mechanical behaviors [20][21][22][23][24]. The effects of these factors must be investigated to reveal the general features of granular friction.
In this study, f-dependent granular friction near the jamming point is mainly examined by experiments and numerical simulations. Based on the obtained data, we will discuss the role of shear jamming in the behavior of granular friction. Finally, an empirical constitutive law for granular friction exerting on a withdrawn rod is derived.

Experiments
In the experiment, a stainless-steel rod buried in a granular layer is slowly withdrawn. The packing fraction of the granular layer is carefully controlled and the pulling force is precisely measured. The experimental setup we used in this study is similar to that used in [7]. A schematic drawing of the experimental setup is shown in figure 1. The rod-withdrawing system is built on the arm and base unit of a microscope system (OLYMPUS BX41). A cylindrical vessel of inner diameter 51mm and height 130mm is mounted on the system. Diameter of the withdrawn rod D is varied from D=0.3 to 3.2mm. As for granular matter, we use roughly monodisperse glass beads of diameter of D g =0.4, 0.6, 0.8, or 2.0mm (AS-ONE, BZ04, BZ06, BZ08, BZ20). The size ratio D D g is a crucial parameter to discuss the granular friction [7]. In the experiment, glass beads are poured into the cylindrical vessel, in which the rod is placed in advance. Weight of the granular bed is always set M=0.32 kg, and the initial penetration (buried) depth of the rod is set h 0 ; 45 mm. The force acting on a rod vertically withdrawn from a granular layer is measured by a load cell (KYOWA ELECTRONIC INSTRUMENT LVS-500GA). During the withdrawing, vertical displacement of the rod δz is controlled by an electric linear ball screw actuator (NTN DMX3104SA-02+D 1-KY) and measured by a transducer (SEIYU ELECTRONIC DTD-3). The vertical direction corresponds to z axis and z=0 is the initial surface level of the granular layer. The actuator is controlled by an in-house developed controlling system, and the displacement δz and pulling force acting on the rod F are simultaneously acquired and recorded on a PC. Using this system called PCM-nano (nano-version Programmable Creep Meter), we can control the maximum pulling velocity v max in the range of v max =0.4-16 μm s −1 . In such slow-slip regime, the effect of inertia can be neglected [7]. Therefore, the inertial number I is irrelevant in this study.
In this experiment, the initial packing fraction f of the granular layer has to be precisely controlled before the withdrawing. To control f of the granular layer, we use both air fluidization and vibration motors. An air pump (FUJI MEDICAL INSTRUMENTS OSP-90W) is connected to an air reservoir to produce a uniform air flow. Twenty vibration motors (SHICOH C1034) are put on the sidewall of the vessel. A camera is placed to take magnified photos of the top surface of the granular layer to estimate its volume. After the initial preparation, the air flow and vibration are applied to control the packing fraction of the layer before withdrawing the rod. Note that air flow and vibration are used only for the sample preparation. They are turned off before the measurement. Namely, the rod is withdrawn without any perturbations by air flow or vibration.
Before withdrawing the rod, we take three steps to make a reproducible initial granular state. First, the granular layer is fluidized by using the air pump so that the granular layer is sufficiently boiled to erase the memory of prior experimental runs and to make a uniform layer. Next, the mechanical vibration is added by the motors to increase f. The value of f can be controlled by vibration strength and duration. Finally, after the above-mentioned operations, we take a photo of the granular surface to measure the volume of granular layer V gra . From the value of true density of glass beads ρ=2500 kg m −3 , f is computed as f=M/(ρ V gra ).
Right after the initial sample preparation, the rod is withdrawn. In this study, the force-control condition is employed. Specifically, the increasing rate of pulling force is controlled as R=1 mN s −1 . From the initial state (t=0) defined by the beginning of pulling, the pulling force F is monitored with a rate of 2000 Hz. If the measured force F is smaller than the control-signal value; F<Rt, the rod is withdrawn by the liner actuator. The minimum withdrawing length set by the actuator is 4 nm. On the other hand, when the measured F is equal to or larger than Rt, the rod is not withdrawn. During the experiment, force F and vertical displacement of the rod δz are recorded on a PC through a DAQ system (NI 9215) with a sampling rate of 10Samples s −1 .

Numerical simulations
We employ discrete element method in our numerical simulations. In the simulation, the particles are considered to be perfect spheres with each of these particles having a center at a point in space. The forces and torque on each of these particles are computed at each time t and their locations and velocities (both linear and angular) are updated at every t + dt based on the Newtonian equations of motion where dt is a small increment in time, i.e. the timestep used in the simulation. Supposing a particle located at position vector r i having a mass m i , diameter D i , and velocity v i . Then we calculate the normal force F ij n and tangential force F ij t of collision between two particles i and j based on the Hertzian contact model proposed by Brilliantov et al [25] as follows, In these equations, δ ij is the overlap between the two particles and is computed as  and s ij D is the tangential displacement vector between the two particles. Here K n and K t represent the normal and tangential nonlinear spring constants, respectively, and γ n and γ t are the normal and tangential dissipation coefficients. The magnitude of F ij t is limited to the maximum of F ij n c m | | to satisfy the slipping criteria where μ c is the frictional coefficient among particles. Table 1 describes the values of these simulation parameters.
In the simulations, we use a similar system to the one used in the experiments. We consider a cylindrical container with a rod of diameter D buried down to a certain depth. The diameter of cylindrical granular layer is ) . In the simulations, we pull the rod at a constant velocity of gd 10 3 in the positive z direction which is 5-10 times higher than the velocities in the experiments. We use this velocity because the computational time for such low velocity becomes extremely long. Decreasing the pulling velocity even further is infeasible. We believe that the simulations are still in the same regime as in experiments because we observe the sliding behaviors that are similar to the experiments. The total number of timesteps is 50 million which may take about few days of computational time for larger D g while several weeks for smaller D g on 42CPU. Additionally, the range of D D g in the simulation study is limited in D g /D1.4.

Experimental results
Typical examples of the experimentally acquired raw data of δz(t) and F(t) are shown in figure 3 (D g /D=1.45). In the early stage of withdrawing, static granular friction sufficiently resists the pulling force. In this regime, F increases roughly by R=1 mN s −1 and the displacement δz is almost constant (δz;0). When the pulling force exceeds a certain limit (at t;73 s; shown by vertical dotted line in figure 3), yielding (sliding) at the frictional interface among the rod and particles occurs. In the late stage, the steady sliding with a constant speed v max can be observed. In this regime, the granular dynamic friction balances with the pulling force. However, immediately after t=73 s, pulling force is not steady and it gradually approaches to the steady state. While the fluctuation of F in this transient state could be interesting, here we focus only on the steady state as a first-step approach.
To characterize the steady granular friction force F á ñ, we identify the steady state of F as follows. The force at the final stage, which is restricted by the limit of displacement of PCM-nano(;4 mm), is defined as F end . Then, the standard deviation in the last 5 s is measured and quantified as σ end . The average value of resistance force F á ñ are computed by the average in the range satisfying the condition F F 3 end end The relation between steady resistance force F á ñ and initial packing fraction f is shown in figure 4(a). For all D D g cases, one can confirm that F á ñ is an increasing function of f. This positive correlation between F á ñ and f is natural and has been observed in the similar previous experiments [8,9]. The measured F á ñ can be transformed into the frictional coefficient by the relation F N m = á ñ . The simplest assumption for the normal , where ρ g = ρf. Using these quantities, F N hydro m = á ñ can easily be computed as shown in figure 4(b). One peculiar feature observed in figure 4(b) is the extremely large μ value in the large f regime. The value of μ exceeds unity and reaches even around 6. Such a large μ value is unexpected because we use cohesionless particles (macroscopic glass beads). For cohesionless particle system, μ should be less than unity.
One might imagine the complex particle motions during the withdrawing for the origin of this peculiar behavior. However, it is difficult to detect the clear particle motions and/or variation of f before and after the withdrawing. This implies that any significant macroscopic deformation of the granular layer (rearrangements of particles) is not induced by the withdrawing. Therefore, the resistance force must originate from the frictional slipping at the interface among particles and rod. However, the estimated value of frictional coefficient μ is much greater than the expected value.

Jammed cylinder
From the above observations, we consider that the simple hydrostatic-pressure approximation does not work well in this experimental setup. Instead, here we consider that the force chain structures developed by the shearing [30] could form the solidified zone with radius nD g ( figure 5). Here, n indicates the normalized radius of solidified cylindrical zone. Then, this solidified region could cause the extremely large normal force on the surface of the rod because the shear-induced force chains are not parallel to the gravitational direction as schematically shown in figure 5. This oblique force-chain structures result in the increase of normal force acting on the rod. However, it is difficult to provide a direct evidence of the solidification based only on the experimental data. Thus, we use numerical simulations to access the internal force structures within the sheared granular layer.

Numerical results
In this subsection, we present the results obtained by the numerical simulations. We report the average withdrawal force F á ñ as the average of the instantaneous withdrawal force once the steady state is achieved just like experiments. The average withdrawal force F á ñ is shown in figure 6(a) in which the data trend is qualitatively similar to that observed in the experiments. The force exerting on the rod can be computed by using shear stress S and pressure P h ,rod ¢ where h¢ is the depth of the rod: the withdrawal force F S A d ò á ñ = and the normal force where dA is a small area on the surface of the rod and h d ¢ is a small length of the rod along z direction. Although we calculate F á ñ by averaging the instantaneous force acting on the rod in the steady state regime, N jam should be measured using this scheme instead of simple hydrostatic form. In addition, we compute the stress tensor for each particle including the ones conforming rod as follows (by assuming ith particle has a diameter D i and is in contact with jth particle) [31], Here, f ij is the contact force and r ij is the position vector from the center of ith particle to the point of contact.
We calculate the pressure Pi on each particle as P Tr Thus, the f-dependent increase of N jam is a plausible reason for the pseudo large friction found in the experiments. Note that any macroscopic compaction of the granular layer cannot be observed in the numerical simulations as well. Namely, tiny rearrangement of the packing structures induced by the rod pulling must develop the effective solidification. Furthermore, as mentioned in section 3.1, hydrostatic pressure cannot   explain the physical reason for large F á ñ value. Therefore, we will analyze the details of the difference between N hydro and N jam and the geometry of shear-jammed region, in the next subsection.

Evidence for solidification by shear-jamming
In order to understand the physical mechanism governing the growth of N jam , we firstly compute, the average pressure acting on particles as a function of depth from the surface z when the rod is not being pulled, i.e. in an initial static state. The static pressure P z hydro á ñ ( ) is calculated by the average of pressure on particles that are located at a depth of z d 2  . Measured P z gd hydro r á ñ ( ) with D g /D=1.422 and f=0.621 is plotted in figure 7(a). The pressure seems to saturate after a certain depth due to Janssen effect(see e.g. [1]). However, the rod is inserted down to the depth 35d where the pressure linearly increases as a function of depth (vertical line in figure 7(a)); hydrostatic regime in a static state. Except the region close to the wall, P(z) is constant and independent of radial distance from the rod when the rod is not being pulled out. P z hydro á ñ ( ) denotes the average at the fixed depth z. This result insists that the hydrostatic pressure is indeed the most crucial factor determining the pressure in the static (not sheared) granular layer. However, when the rod is being pulled out, this is not the case. Pressure on a particle becomes a function of the radial distance from the rod r¢ due to the shearing induced by the withdrawn rod. Here, we measure the average pressure on a particle in the range of a radial distance r d 2.55 ¢  and a depth z d 2 ;  P r z , á ¢ ñ ( ) . We observe that once the steady state is achieved, the particles near the rod experience relatively higher pressure which also increases with depth. To characterize the degree of pressure increase, we compare P r z , á ¢ ñ ( ) with a hydrostatic pressure P z hydro á ñ ( ) by using a form, which becomes zero when the depth-dependent pressure is simply identical to the hydrostatic one. In figure 7(b), we plot R z jam ( ) with D g /D=1.422 and f=0.621. It can be seen that a high pressure region is formed near the rod. This is a clear evidence of shear-induced pressure increase. We consider this pressure increase causes the effective solidification (shear jamming). In the large r¢ regime, P r z , á ¢ ñ ( ) approaches to P z hydro á ñ ( ) . This implies that the region of shear jamming is limited in the vicinity of the rod, as schematically illustrated in figure 5. One can observe that there exists a unique distance from the rod where P r z P z , hydro á ¢ ñ á ñ  ( ) ( ) is fulfilled. This proves that the shape of the shear zone is cylindrical. Through the simulations, it is difficult to verify the details of D g dependence of this crossover ( P r z P z , hydro á ¢ ñ á ñ  ( ) ( ) ) point, due to the limitation of computational time for small D D g cases or the limited quality of statistics for large D D g case. In figure 7(c), the variation of the radius of shear cylinder nD g against f with D g /D=1.422 is displayed. It can be clearly seen that n increases with f. Now, we are sure about the physical reason for the large F á ñ in the large f regime. The effective solidification due to the shear jamming is caused in the vicinity of the withdrawn rod. Then, the next problem to be studied is D D g and f dependence of the shear-jammed region. To this end, the systematic analysis of the experimental data is much more feasible than performing the complete set of numerical simulations. Therefore, we will go back to the analysis of the experimental data in the next section.

Size of shear-jammed (solidified) cylinder
As mentioned so far, the origin of large F á ñ is the growth of shear-jammed (solidified) zone. Thus, we discuss the relation between shear jamming by the rod withdrawing and initial packing fraction. When the rod is withdrawn, particles around the rod are sheared. As a consequence, if the initial packing fraction is large enough, the shear-induced jamming could be developed in the granular layer. By this shear-jam effect, a certain amount of particles close to the rod would be effectively solidified.
To understand fand D D g -dependent behaviors of F á ñ, we revisit the solidified cylinder model introduced in [7]. In the model, we assume that the frictional slipping occurs at the interface between the rod and particles. We also assume that the frictional coefficient of the slip can be regarded as a constant. We additionally assume that the normal force acting on the interface among the rod and particles is proportional to the mass of the solidified cylindrical region since a part of gravitational loading is redirected to the surface of the rod due to the inclined force chain structures. To specifically estimate the normal force acting on the rod surface N jam , the redirection ratio of the vertical force to the horizontal normal force, κ (Janssen parameter [1]), should be introduced. However, the precise determination of κ value is not very easy in general(see e.g. section 3.7 in [32]). Actually, the specific value of κ is not very important because we later discuss the critical divergence of the characteristic length scale (solidified cylinder radius) which is independent of the value of κ. Thus, in this study, the effective frictional coefficient jam * m m k = is used in the following calculation. Then, similarly to Coulomb-Amonton's frictional constitutive law, the force is expressed as F M g jam * m á ñ = , where M jam is the mass of solidified region and * m can be regarded as a constant. Using this solidified cylinder model, the increase of F á ñ in large f regime can be transformed into the expansion of the effective radius of solidified (jammed) region, n, in units of particle diameter, i.e. the actual radius of the solidified zone is nD g . By considering the cylindrical geometry, * m is written as [7], is employed. This value is reasonable by considering the typical frictional coefficient between glass and steel (μ=0.18 [33]) and the practical Janssen parameter value (κ;0.5). As can be seen in figure 8(a), n has a divergent tendency against f. However, the data show systematic variation depending on the size ratio D D g .

Divergence of the size of shear-jammed cylinder
To understand the underlying physics governing the behavior of n shown in figure 8(a), we analyze the D D g dependence of n. As a result, we find that the relation between n and D D g for various f is scaled by the powerlaw relation n D D g a -( ) with a characteristic exponent α = 0.56. We have confirmed that this relation is actually consistent with the previous study [7], in which the steady resistance force is a decreasing function of D D g . Using the relation, the scaled size of shear-jammed cylinder, n(D g /D) α , is plotted as a function of f in figure 8(b).
From the data trend observed in figure 8(b), we consider that n shows a critical divergence as f approaches the critical (jamming) point f c . To characterize the critical divergence of the length scale n, the data are fitted to the form, c 0 where β is a critical exponent characterizing the divergence and n 0 is a fitting parameter. The solid gray curve shown in figure 8(b) shows the least square fitting by equation (7). The already estimated value α=0.56 is used as a fixed parameter. Then, the fitting-parameter values are obtained as n 0 =0.35, f c =0.62, and β=1.3. In the inset of figure 8(b), the relation between n D D g a ( ) and f c -f is displayed in the log-log style. The data behavior is consistent with the critical divergence modeling.

Constitutive law for withdrawing granular friction
Altogether, now we can conclude the effective constitutive law for the granular withdrawal frictional force as follows,  * m pr á ñ = + - where n obeys the critical divergence form (equation (7)). This form is equivalent to equation ( ) , and β=1.3. This constitutive law includes the effective frictional resistance force governed by the shear-induced jamming and geometrical conditions (size ratio and cylindrical shear-jammed zone). Although this form is quite different from the conventional granular frictional constitutive law using the inertial number(e.g. [3][4][5][6]), the current experiment actually focuses only on the very small inertial number regime to reveal the effects of D D g and f. As discussed in [7], the rough estimate of I can be I v gh max = . By substituting characteristic values (v max ;10 −5 m s −1 and h;5 mm), we can estimate I is less than 10 −4 . In addition, we assume that the interfacial friction between the rod and particles is almost constant. While the current analysis relies on these assumptions, we find that the idea of critical divergence of the characteristic length scale by the jamming can be applied also to the vertical withdrawing system, just like other granular systems(e.g. [15,34,35]). Moreover, the form of critical divergence of the characteristic length scale written in equation (7) might be universal and independent of the details of experimental setup. However, the physical origin of the scaling exponents (α=0.56 and β=1.3) has not been revealed in this study. This is still an open problem. The computed F á ñ values for D g /D=0.13, 0.33, and 2.66 are shown as dotted curves in figure 4(a) by each corresponding color. Although these curves basically reproduce the experimental data, the quality of agreement is slightly limited. For instance, the constitutive law underestimates F á ñ in large D g /D(=2.66). Since the constitutive law is derived on the basis of simple scaling analysis, the law can be used to roughly estimate the value of F á ñ. In this study, the static granular friction without any perturbation has been investigated. As a next step, the effect of perturbations (air flow and/or vibration) applied during the withdrawing is also an interesting topic for future study to elucidate the physical mechanism of frictional constitutive law for granular matter.

Conclusion
In this study, granular frictional constitutive law was experimentally obtained from the relation among the steady resistance force, packing fraction, and size ratio between particles and withdrawn object. The increase of steady frictional resistance F á ñ by increasing f can be understood by the shear-induced solidification of the granular layer. The evidence of this shear-jammed (solidified) zone was obtained by numerical simulation. According to the numerical results, the cylindrical region around the withdrawn rod is effectively solidified by strong contact forces. The size of solidified region has a positive correlation to the initial packing fraction. From these numerical evidences, we conclude that the simple model of a cylindrical solidified zone can explain the pseudo increase of granular friction. Then, the size of cylindrical solidified zone was systematically estimated using the experimental data. Using these data, we built an empirical granular frictional constitutive law. In the modeling, we considered the elongation of the force chains developed by shearing. Due to this force chain developments, the cylindrical unit around the rod seems to be solidified by shear jamming. The normalized radius of this solidified cylindrical region n is estimated by the experimental result. The relation among n, D g /D, and f was analyzed and the empirical form for estimating n was obtained as equation (7). Although the obtained empirical form is consistent with previous studies, the physical meaning of the estimated critical exponents α and β in equation (7) has not yet been understood well. To establish a more comprehensive and universal granular frictional constitutive law, we have to pay attention also to other parameter dependencies as well (e.g. sliding velocity, particle shape etc).