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

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}_{{\rm{g}}}/D$ 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₀ ≃ 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⁻¹. In such slow-slip regime, the effect of inertia can be neglected [7]. Therefore, the inertial number I is irrelevant in this study.

Zoom In Zoom Out Reset image size

In this experiment, the initial packing fraction ϕ of the granular layer has to be precisely controlled before the withdrawing. To control ϕ 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 ϕ. The value of ϕ 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⁻³, ϕ is computed as ϕ = 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⁻¹. 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⁻¹.

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 ${{\boldsymbol{r}}}_{{\boldsymbol{i}}}$ having a mass m_i, diameter D_i, and velocity ${{\boldsymbol{v}}}_{{\boldsymbol{i}}}$. Then we calculate the normal force ${{\boldsymbol{F}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{n}}}$ and tangential force ${{\boldsymbol{F}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{t}}}$ of collision between two particles i and j based on the Hertzian contact model proposed by Brilliantov et al [25] as follows,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{n}}}=\sqrt{\displaystyle \frac{{D}_{i}{D}_{j}}{2({D}_{i}+{D}_{j})}}\sqrt{{\delta }_{{ij}}}({K}_{n}{\delta }_{{ij}}{\hat{{\boldsymbol{n}}}}_{{\boldsymbol{ij}}}-{m}_{{\rm{eff}}}{\gamma }_{n}{{\boldsymbol{v}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{n}}}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{t}}}=-\sqrt{\displaystyle \frac{{D}_{i}{D}_{j}}{2({D}_{i}+{D}_{j})}}\sqrt{{\delta }_{{ij}}}({K}_{t}{\boldsymbol{\Delta }}{{\boldsymbol{s}}}_{{\boldsymbol{ij}}}+{m}_{{\rm{eff}}}{\gamma }_{t}{{\boldsymbol{v}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{t}}}).\end{eqnarray} \tag{ 2 }$$

In these equations, δ_ij is the overlap between the two particles and is computed as ${\delta }_{{ij}}=({D}_{i}+{D}_{j})/2-| {{\boldsymbol{r}}}_{{\boldsymbol{i}}}-{{\boldsymbol{r}}}_{{\boldsymbol{j}}}|$, while the unit vector joining the centers of two particles ${{\boldsymbol{n}}}_{{\boldsymbol{ij}}}=\tfrac{{{\boldsymbol{r}}}_{{\boldsymbol{i}}}-{{\boldsymbol{r}}}_{{\boldsymbol{j}}}}{| {{\boldsymbol{r}}}_{{\boldsymbol{i}}}-{{\boldsymbol{r}}}_{{\boldsymbol{j}}}| }$. The terms ${{\boldsymbol{v}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{n}}}$ and ${{\boldsymbol{v}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{t}}}$ represent the normal and tangential relative velocities, respectively. ${m}_{{\rm{eff}}}=\tfrac{{m}_{i}{m}_{j}}{{m}_{i}+{m}_{j}}$ is the effective mass and ${\boldsymbol{\Delta }}{{\boldsymbol{s}}}_{{\boldsymbol{ij}}}$ 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 ${{\boldsymbol{F}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{t}}}$ is limited to the maximum of ${\mu }_{{\rm{c}}}| {{\boldsymbol{F}}}_{{\boldsymbol{ij}}}^{{\boldsymbol{n}}}|$ to satisfy the slipping criteria where μ_c is the frictional coefficient among particles. Table 1 describes the values of these simulation parameters.

Table 1.  The values of parameters used in the simulations.

Parameters	Value
Kn	$2.0\times {10}^{6}\rho {dg}$
Kt	$2.45\times {10}^{6}\rho {dg}$
${\gamma }_{{\rm{n}}}$	$3200g/{d}^{3/2}$
${\gamma }_{{\rm{t}}}$	$3200g/{d}^{3/2}$
${\mu }_{{\rm{c}}}$	0.2
dt	${10}^{-4}\sqrt{\tfrac{d}{g}}$

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 fixed as $51d$ where d = 10⁻³ m is the unit of length scale in the simulations. The cylinder is filled with particles of diameter D_g with a polydispersity of 10% by randomly pouring them from the top. The packing fraction ϕ in the simulations is varied roughly from 0.6 to 0.637. The number of particles poured in the vessel ranges from 15 850 to 53 500 depending on D_g. The height of the granular particles inside the cylinder is approximately $90d$ with the rod being inserted till a depth of 35d from the free surface as shown in figure 2. The density of the granular particles and the rod is set as ρ = 2500 kg m⁻³. The magnitude of gravitational acceleration is g = 9.8 m s⁻² in the negative z direction and the effective diameter of the rod is $D=1.125d$. The rod is constructed by using smaller particles of diameter $0.25d$ along the circumference of a cylinder of diameter 0.875d and thus, $D=(0.875+0.25)d=1.125d$.

Zoom In Zoom Out Reset image size

In the simulations, we pull the rod at a constant velocity of ${10}^{-3}\sqrt{{gd}}$ 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}_{{\rm{g}}}/D$ in the simulation study is limited in D_g/D ≥ 1.4.

We use Large-scale Atomic/Molecular Massively Parallel Simulator [26, 27] for the simulations. Ovito [28] and VMD [29] are used for the post-simulation visualization.

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⁻¹ 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.

Zoom In Zoom Out Reset image size

To characterize the steady granular friction force $\langle F\rangle$, 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 $\langle F\rangle$ are computed by the average in the range satisfying the condition $| F-{F}_{\mathrm{end}}| \lt 3{\sigma }_{\mathrm{end}}$ from the final measurement point.

The relation between steady resistance force $\langle F\rangle$ and initial packing fraction ϕ is shown in figure 4(a). For all ${D}_{{\rm{g}}}/D$ cases, one can confirm that $\langle F\rangle$ is an increasing function of ϕ. This positive correlation between $\langle F\rangle$ and ϕ is natural and has been observed in the similar previous experiments [8, 9]. The measured $\langle F\rangle$ can be transformed into the frictional coefficient by the relation $\mu =\langle F\rangle /N$. The simplest assumption for the normal force N is hydrostatic pressure times surface area of the rod. Then, the normal force acting on the rod surface is written as ${N}_{\mathrm{hydro}}=\pi {\rho }_{{\rm{g}}}{{gh}}_{0}^{2}D/2$, where ρ_g = ρϕ. Using these quantities, $\mu =\langle F\rangle /{N}_{\mathrm{hydro}}$ 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 ϕ 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.

Zoom In Zoom Out Reset image size

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 ϕ 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.

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.

Zoom In Zoom Out Reset image size

In this subsection, we present the results obtained by the numerical simulations. We report the average withdrawal force $\langle F\rangle$ as the average of the instantaneous withdrawal force once the steady state is achieved just like experiments. The average withdrawal force $\langle F\rangle$ 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^{\prime} ,{\rm{rod}}}$ where $h^{\prime}$ is the depth of the rod: the withdrawal force $\langle F\rangle =\int {S}{\rm{d}}{A}$ and the normal force ${N}_{{\rm{jam}}}=\int {P}_{h^{\prime} ,{\rm{rod}}}{\rm{d}}{A}=\int {P}_{h^{\prime} ,{\rm{rod}}}\pi {D}{\rm{d}}{h}^{\prime}$ where dA is a small area on the surface of the rod and ${\rm{d}}{h}^{\prime}$ is a small length of the rod along z direction. Although we calculate $\langle F\rangle$ by averaging the instantaneous force acting on the rod in the steady state regime, ${N}_{{\rm{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],

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{6}{\pi {D}_{i}^{3}}\displaystyle \sum _{i\ne j}^{{\rm{All}}\,{\rm{Contacts}}}{{\boldsymbol{f}}}_{{\boldsymbol{ij}}}\otimes {{\boldsymbol{r}}}_{{\boldsymbol{ij}}}.\end{eqnarray} \tag{ 3 }$$

Here, ${{\boldsymbol{f}}}_{{\boldsymbol{ij}}}$ is the contact force and ${{\boldsymbol{r}}}_{{\boldsymbol{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}_{i}=-\tfrac{1}{3}\mathrm{Tr}(\sigma )$. Based on this pressure, we calculate the ${N}_{{\rm{jam}}}$ on the rod as shown in figure 6(b). The increase of N_jam associated with ϕ can be observed. Then we show ${\mu }_{\mathrm{jam}}=\langle F\rangle /{N}_{{\rm{jam}}}$ in figure 6(c). In figure 6(c), μ_jam remains almost constant with respect to the packing fraction for a given ${D}_{{\rm{g}}}/D$. Although it could slightly depend on ${D}_{{\rm{g}}}/D$, the obtained values of μ_jam are close to μ_c = 0.2.

Zoom In Zoom Out Reset image size

Thus, the ϕ-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 $\langle F\rangle$ 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.

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 $\langle {P}_{\mathrm{hydro}}(z)\rangle$ is calculated by the average of pressure on particles that are located at a depth of $z\pm 2d$. Measured $\langle {P}_{\mathrm{hydro}}(z)\rangle /\rho {gd}$ with D_g/D = 1.422 and ϕ = 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. $\langle {P}_{\mathrm{hydro}}(z)\rangle$ 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.

Zoom In Zoom Out Reset image size

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^{\prime}$ 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^{\prime} \pm 2.55d$ and a depth $z\pm 2d;$ $\langle P(r^{\prime} ,z)\rangle$. 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 $\langle P(r^{\prime} ,z)\rangle$ with a hydrostatic pressure $\langle {P}_{\mathrm{hydro}}(z)\rangle$ by using a form,

$$\begin{eqnarray}&&{R}_{\mathrm{jam}}(z)=\displaystyle \frac{\langle P(r^{\prime} ,z)\rangle }{\langle {P}_{\mathrm{hydro}}(z)\rangle }-1,\end{eqnarray} \tag{ 4 }$$

which becomes zero when the depth-dependent pressure is simply identical to the hydrostatic one. In figure 7(b), we plot ${R}_{\mathrm{jam}}(z)$ with D_g/D = 1.422 and ϕ = 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^{\prime}$ regime, $\langle P(r^{\prime} ,z)\rangle$ approaches to $\langle {P}_{\mathrm{hydro}}(z)\rangle$. 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 $\langle P(r^{\prime} ,z)\rangle \simeq \langle {P}_{\mathrm{hydro}}(z)\rangle$ 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 ($\langle P(r^{\prime} ,z)\rangle \simeq \langle {P}_{\mathrm{hydro}}(z)\rangle$) point, due to the limitation of computational time for small ${D}_{{\rm{g}}}/D$ cases or the limited quality of statistics for large ${D}_{{\rm{g}}}/D$ case. In figure 7(c), the variation of the radius of shear cylinder nD_g against ϕ with D_g/D = 1.422 is displayed. It can be clearly seen that n increases with ϕ.

Now, we are sure about the physical reason for the large $\langle F\rangle$ in the large ϕ 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}_{{\rm{g}}}/D$ and ϕ 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.

As mentioned so far, the origin of large $\langle F\rangle$ 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 ϕ- and ${D}_{{\rm{g}}}/D$-dependent behaviors of $\langle F\rangle$, 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 ${\mu }^{* }={\mu }_{\mathrm{jam}}\kappa$ is used in the following calculation. Then, similarly to Coulomb–Amonton's frictional constitutive law, the force is expressed as $\langle F\rangle ={\mu }^{* }{M}_{\mathrm{jam}}g$, where M_jam is the mass of solidified region and ${\mu }^{* }$ can be regarded as a constant. Using this solidified cylinder model, the increase of $\langle F\rangle$ in large ϕ 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, ${\mu }^{* }$ is written as [7],

$$\begin{eqnarray}&&{\mu }^{* }=\displaystyle \frac{\langle F\rangle }{\pi {\rho }_{{\rm{g}}}{gh}\left[{\left({{nD}}_{{\rm{g}}}+D/2\right)}^{2}-{\left(D/2\right)}^{2}\right]}.\end{eqnarray} \tag{ 5 }$$

By solving this equation, n is written as

$$\begin{eqnarray}&&n=\displaystyle \frac{-{{AD}}_{{\rm{g}}}D+\sqrt{{\left({{AD}}_{{\rm{g}}}D\right)}^{2}+4{{AD}}_{{\rm{g}}}^{2}\langle F\rangle }}{2{{AD}}_{{\rm{g}}}^{2}},\end{eqnarray} \tag{ 6 }$$

where, $A={\mu }^{* }\pi {\rho }_{{\rm{g}}}{gh}$. The computed n values for all the experimental runs are shown in figure 8(a). To obtain the specific values of $n,{\mu }^{* }=0.1$ 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 ϕ. However, the data show systematic variation depending on the size ratio ${D}_{{\rm{g}}}/D$.

Zoom In Zoom Out Reset image size

To understand the underlying physics governing the behavior of n shown in figure 8(a), we analyze the ${D}_{{\rm{g}}}/D$ dependence of n. As a result, we find that the relation between n and ${D}_{{\rm{g}}}/D$ for various ϕ is scaled by the power-law relation $n\sim {({D}_{{\rm{g}}}/D)}^{-\alpha }$ 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}_{{\rm{g}}}/D$. Using the relation, the scaled size of shear-jammed cylinder, n(D_g/D)^α, is plotted as a function of ϕ in figure 8(b).

From the data trend observed in figure 8(b), we consider that n shows a critical divergence as ϕ approaches the critical (jamming) point ϕ_c. To characterize the critical divergence of the length scale n, the data are fitted to the form,

$$\begin{eqnarray}&&n={n}_{0}{\left(\displaystyle \frac{{D}_{{\rm{g}}}}{D}\right)}^{-\alpha }{\left({\phi }_{c}-\phi \right)}^{-\beta },\end{eqnarray} \tag{ 7 }$$

where β is a critical exponent characterizing the divergence and n₀ 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.35, ϕ_c = 0.62, and β = 1.3. In the inset of figure 8(b), the relation between $n{({D}_{{\rm{g}}}/D)}^{\alpha }$ and ϕ_c–ϕ is displayed in the log–log style. The data behavior is consistent with the critical divergence modeling.

Altogether, now we can conclude the effective constitutive law for the granular withdrawal frictional force as follows,

$$\begin{eqnarray}&&\langle F\rangle ={\mu }^{* }\pi {\rho }_{{\rm{g}}}{gh}\left[{\left({{nD}}_{{\rm{g}}}+\displaystyle \frac{D}{2}\right)}^{2}-{\left(\displaystyle \frac{D}{2}\right)}^{2}\right],\end{eqnarray} \tag{ 8 }$$

where n obeys the critical divergence form (equation (7)). This form is equivalent to equation (5). The parameter values for the glass beads and cylindrical stainless-steel rod case are ${\mu }^{* }(={\mu }_{\mathrm{jam}}\kappa )=0.1,{n}_{0}\,=\,0.35,{\phi }_{{\rm{c}}}=0.62,\alpha =0.56$, 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–6]), the current experiment actually focuses only on the very small inertial number regime to reveal the effects of ${D}_{{\rm{g}}}/D$ and ϕ. As discussed in [7], the rough estimate of I can be $I={v}_{\max }/\sqrt{{gh}}$. By substituting characteristic values (v_max ≃ 10⁻⁵ m s⁻¹ and h ≃ 5 mm), we can estimate I is less than 10⁻⁴. 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 $\langle F\rangle$ 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 $\langle F\rangle$ 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 $\langle F\rangle$. 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.