Scaling law of the drag force in dense granular media
Sonu Kumar K. Anki Reddy Satoshi Takada Hisao Hayakawa
cond-mat.soft




Department of Chemical Engineering, Indian Institute of Technology, Guwahati, Assam, India

anki.reddy@iitg.ernet.in
Department of Chemical Engineering, Indian Institute of Technology, Guwahati, Assam, India

Earthquake Research Institute, University of Tokyo, 1-1-1, Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan
Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan

hisao@yukawa.kyoto-u.ac.jp
Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan





Making use of the system of pulling a spherical intruder in static three-dimensional granular media, we numerically study the scaling law for the drag force F_ drag acting on the moving intruder under the influence of the gravity.
Suppose if the intruder of diameter D immersed in a granular medium consisting of grains of average diameter d is located at a depth h and moves with a speed V, we find that F_ drag can be scaled as (D+d)^ϕ_μ h^α_μ with two exponents ϕ_μ and α_μ, which depend on the friction coefficient μ and satisfy an approximate sum rule ϕ_μ+α_μ≈ 3.
This scaling law is valid for the arbitrary Froude number (defined by Fr=2 V√(2D/g)/(D+d)), if h is sufficiently deep.
We also identify the existence of three regimes (quasistatic, linear, and quadratic) at least for frictional grains in the velocity dependence of drag force.
The crossovers take place at Fr≈ 1 between the quasistatic to the linear regimes and at Fr≈ 5 between the linear to the quadratic regimes.
We also observe that Froude numbers at which these crossovers between the regimes happen are independent of the depth h and the diameter of the intruder D.
We also report the numerical results on the average coordination number of the intruder and average contact angle as functions of intruder velocity.



Scaling law of the drag force in dense granular media
    Hisao Hayakawa
    December 30, 2023
=====================================================






§ INTRODUCTION

Granular medium, an assembly of discrete particles behaves as a solid or a fluid depending on its density <cit.>.
Variety of interesting phenomena exhibited by granular materials <cit.> has attracted many physicists and engineers in the last few decades.
There are many studies on the motion of the objects in granular media (through analysis of drag and lift forces) to improve the understanding of rheology of granular flows.

There are both experimental and simulation studies related to the characterization of the drag force on a moving passive object.
In most of these studies the objects considered are either spheres <cit.>, cylinders <cit.>, or plate-like objects <cit.>.
Through the analysis of slow drag on a cylinder in aluminium oxide polydisperse granular mixtures under high pressure <cit.>, it is clarified that the drag resistance in granular media is related to the packing effects.
The studies of drag in a quasi-two-dimensional dilute supersonic granular flow <cit.> and in monolayered (two-dimensional) dense granular media <cit.> have provided more quantitative descriptions of drag force as the effects of flow fields around the intruder.
A quite recent study on the drag force of a spherical intruder moving through sedimented granular hydrogels <cit.> reports that drag force is a constant up to a critical velocity and starts increasing with velocity quadratically.
In a study of granular drag inspired by self-burrowing rotary seeds <cit.>, the drag reduction by rotation for a vertically penetrating intruder is observed.
Similarly the drag reduction by mechanical vibration for the penetration of an intruder into a dry granular medium has been reported <cit.> and this is attributed to the local fluidization which could rupture the force chains.
Such studies on the drag force in granular media are relevant to understand the animal locomotion through granular media <cit.>.
We expect that our study can stimulate the robotics or animal locomotion moving in sands.
  
  

Though many studies exists for characterization of the drag force, we need to establish a scaling law for the drag forces in granular media similar to the case of viscous fluid flows.
If there is an object in a viscous fluid and if there exists a relative speed V between the object and the fluid, it is well-known that F_ drag can be given as (1/2)ρ C_D A V^2, where A is a reference area, ρ is the fluid density, and C_D is the drag coefficient.
Such kind of universal law does not exist for granular fluids.
In a recent study <cit.> for the drag force acting on a disk at high velocity near the jamming point, it is found that the average drag force can be expressed as F_drag= F_0(Φ) + α(Φ)V^2, where Φ is the area fraction and V is the pulling speed of the intruder disk.
This scaling  has also been observed in the numerical study <cit.> of similar system.
Nevertheless, the situation is not simple for three-dimensional cases.
Indeed, the simulation in Ref. <cit.> suggests F_ drag for a pulling sphere in a three-dimensional granular media under the influence of the gravity is expressed as the summation of a yield force and the dynamical term which is proportional to V, at least, for V<3√(g D/2), where D and g are, respectively, the diameter of the intruder and the  acceleration due to gravity.
On the other hand, an experiment on the moving cylinder in a three-dimensional granular media using Taylor-Couette configuration  <cit.> suggests that the drag law obeys V∝exp[(F_ drag-F_ c)/F_ c G], where F_ c and G are the yield force and a function that depends on the location (or the shear).
Therefore, we should clarify what is the proper scaling law to characterize the moving object in granular media in three-dimensional situations.

The recent work to study the effects of vibration on the drag force of an intruder revealed that the drag force does not depend on the velocity for low vibration level and depends linearly for high vibrations <cit.>.
Similarly velocity dependence of F_ drag has been discussed in both experimental and numerical works <cit.> for various shaped intruders and various packing conditions either in two-dimensional or three-dimensional systems.
Based on the results of these studies, it seems that there exist three or four regimes such as velocity independent (quasistatic), weakly dependent (logarithmic), linear, and quadratic regimes, though the distinction between quasistatic and the logarithmic regimes is not easy.
In the present work with the help of extensive computer simulations, we numerically analyze the drag force for a wide range of intruder velocities characterized by the Froude number defined by Fr=2 V√(2D/g)/(D+d), where d is the average diameter of the granular bed particles.
Plotting the numerical data against the Froude number, we clarify the existence of crossovers between these regimes.
Moreover we also attempt to know the role of depth of the immersed object on the crossover.
We discuss how the drag forces are scaled by the intruder diameter and the depth.
We also discuss how the scaling exponents depend on the friction coefficient.

The present study makes use of the simulation of pulling a spherical intruder in dense granular medium to answer the questions discussed above.
In Sec. <ref>, we explain the simulation methodology employed for the present study.
We illustrate that the drag force can be represented by a product of the yield force and the dynamical part.
We also discuss how the yield force depends on the intruder size, the depth of the intruder, and the friction coefficient in Sec. <ref>.
Finally we summarize our results with the possible future work in Secs. <ref> and <ref>.





§ SIMULATION METHODOLOGY

We employ the discrete element method (DEM) as our simulation method <cit.>, which is widely used to simulate the granular materials.
The simulations are carried out for a three-dimensional system. Initial configurations for the simulation are generated by fixing the intruder sphere at a given position while granular particles of diameter uniformly distributed between 0.9d to 1.1d are poured from the top by the gravity.
The particles are allowed to settle and reach almost zero kinetic energy due to the presence of dissipating and damping forces.
In our simulations, physical quantities are scaled by the density ρ, the average diameter d, and the gravity acceleration g.
Therefore, for example, the time is scaled by √(d/g), the mass is scaled by ρ d^3, the force is scaled by ρ d^3 g, etc.



Periodic boundary conditions are applied across x and y-directions, where x coordinate is assigned to the moving direction of the intruder, and z is the vertical direction (anti-parallel to the direction of the gravity).
A rough base is used for the bottom of the simulation box.
The gravity acceleration acts along negative z-direction throughout the simulation.
We have considered two systems in our study.
The first system consists of 70,001 particles (including one intruder) in a simulation box of 40d×40d×38d in the x, y, and z-directions.
The intruder is initially located at x=y=20 d and z=13.5 d.
In this system, the depth of the intruder from the surface is kept constant while the diameter of the intruder is varied.
Six cases of the intruder diameter are considered with the intruder diameter varying from 1.0d to 6.0d at various moving velocities inside the granular bed.
It is ensured that the system is large enough for all the simulations carried out and that the periodic boundary does not affect the forces on the intruder due to its periodic image.
The initial configuration for the intruder of diameter 2.0d is shown in Figs. <ref>(a) and <ref>(b).
We have additionally considered the intruders of diameters D=8d and 10d in the first system to confirm the scaling law on the yield force.
The second system consists of 300,001 particles including the intruder in a system of 80d×40d×80d in the x, y and z-directions (see Fig. <ref>(c)).
The intruder diameter is fixed as 4.0d while the depth of the intruder is varied as 10d, 25d, 45d, 50d, 55d, 60d, and 65d from the surface in the negative z-direction.
The first system I is used to study the dependence of drag force on the diameter of the intruder at constant depth for both frictional and frictionless cases, the drag force in the quasistatic regime and the yield force of the intruder (here, we define the yield force as the drag force acting on the intruder in the zero velocity limit).
The second system II, on the other hand, is used to study the drag forces on the intruder of a constant diameter 4.0d at various depths and the dependence of the drag forces on the other factors such as non-linear spring stiffness and damping coefficient both of which are defined in the next paragraph.
The benefit of using two systems is as follows.
The system I can be used to determine the yield force on the intruder because we need huge computation time (the number of particles in the system II is four times or more than that in the system I), while the system II is adequate to know the depth dependence of the drag force because the location of the intruder should be away from the bottom rough base and the surface of the granular bed.

The normal and tangential contact forces <cit.> between the two interacting particles or between the particle and the wall are expressed as:

    F_n=√(R_effδ) (K_nδ-m_effγ_n v_n),
       F_t=-min(μ F_n,√(R_effδ) (K_tΔs_t+m_effγ_t v_t)),

Here, K_n and K_t are non-linear spring constants for the normal and the tangential contact forces, respectively, with the dimensions of the pressure, γ_n and γ_t are damping factors for the normal and the tangential forces, respectively, δ refers to the overlap distance of the colliding particles and is given by δ=(d_1+d_2)/2-|d_12| where d_1 and d_2 are the diameters of the colliding particles 1 and 2, respectively, and d_12 is the instantaneous distance between the two colliding particles.
Moreover, Δ s_t is the tangential displacement between two particles when in contact, v_n and v_t are the normal and the tangential components of the relative velocity, respectively.
R_eff is expressed as d_1d_2/(2(d_1+d_2)) and the reduced mass m_eff is expressed as m_eff=m_1m_2/(m_1+m_2), where m_1 and m_2 are the masses of the colliding particles 1 and 2, respectively.
Here μ is the coefficient of friction, for particle-particle interactions, which is varied from 0.0 to 0.9 in the present simulation study.
There is an upper limit μ F_n on the tangential component of contact force to express the slip motion between two contacting particles.




We have taken K_n=2×10^8 ρ dg and the value of K_t=2.45×10^8 ρ dg for our simulations.
These values of K_n and K_t correspond to Young's modulus, E=200  GPa and ν=0.3 if we choose ρ=10  g/cm^3, d=7.5  mm, and g=9.8  m/s^2.
The identical K_n and K_t can be used to represent a material of Young's modulus, E=10  GPa and ν=0.3 if we choose ρ=1  g/cm^3, d=3.75  mm, and g=9.8  m/s^2.
The value of γ_n is chosen as 50000 √(g/d^3) to obtain a coefficient of restitution curve that represents a realistic value.
It is seen that the coefficient of restitution depends on the colliding velocity <cit.>.
We have carried out the simulation of dropping a sphere on a flat wall to obtain coefficient of restitution for various V_n as shown in Fig. <ref> with the interaction parameters identical to that of the sphere-sphere interaction.
The theoretical value of the coefficient of restitution is also plotted (calculated taking diameter of wall as infinity) using the well-known expression obtained by Kuwabara and Kono <cit.> with the first four constants of the infinite series <cit.> and is seen to converge with the one obtained from our simulations at various velocities (see Fig. <ref>).
The coefficient of restitution in our simulation satisfies the well known curve of the form 1-k_1 (v_n)^1/5+k_2 (v_n)^2/5+..., where k_i are constants for a given K_n and γ_n, and v_n is the normal component of velocity along the line of impact.
Note that the coefficient of restitution curve remains the same if γ_n/K_n^0.6 ratio is kept constant <cit.>.
This has been used in some part of our study to maintain the same curve for the coefficient of restitution while varying the value of K_n.
Basically, we fix K_n as K_n=2×10^8 ρ d g except for the places where we change K_n explicitly.
The value of γ_t is identical to γ_n <cit.> in our simulation.




The value of a time step is chosen as 1.0× 10^-5√(d/g) <cit.>.
The intruder is moved with a constant speed in x-direction at a constant depth.
We have two time scales t_1 and t_2 in our system as shown in the Fig. <ref>.
The first time scale t_1 characterizes the forward motion of the intruder and is given by t_1=(D+d)/2V.
The second time-scale t_2 represents the characteristic time for the particle to fall from the top to the bottom, which is given by t_2=√(2D/g).
The Froude number Fr introduced in Sec. <ref> corresponds to:

    Fr=t_2/t_1=2V/(D+d)√(2D/g).


We are interested in this definition of Froude number because of the transition that happens at t_1=t_2, i.e., Fr=1.
The simulations are performed for Fr≤ 10.
The coefficient of the friction μ is assumed to be a constant for all contacts.
The standard velocity Verlet time  integration algorithm  is used to update the particles positions and velocities in our DEM simulation.
The simulation is carried out until the intruder to cover the whole simulation box length, i.e., to move one complete length of the simulation box in the x direction.
The simulations are carried out using Large-scale Atomic/Molecular Massively Parallel Simulator <cit.>. OVITO <cit.> and VMD <cit.> are used for post-simulation visualization and analysis.





§  RESULTS AND ANALYSIS





This section is the main part to present the results of our simulation study.
This section consists of 8 subsections.
In Sec. <ref>, we study the yield force on an intruder and find a scaling law for the yield force for several μ and D.
In Sec. <ref>, we study the yield force on an intruder for various μ at different intruder depth h and propose an exponential scaling for the depth dependence.
In Sec. <ref>, we show the results of the drag force at μ=0.5 for various intruder diameters D.
In Sec. <ref>, we present the results of our simulation on the drag force for frictionless grains (μ=0) for various intruder diameters.
In Sec. <ref>, we examine how the drag force depends on the friction coefficient μ in the range 0≤μ≤ 0.3 for Fr≤ 1.
In Secs. <ref>, we present the detailed results of the drag force for various depths for a constant intruder size D=4.0d.
In Sec. <ref>, we examine how the drag force depends on the stiffness constant K_n and the dissipation constant γ_n.
In Sec. <ref>, we demonstrate the existence of a scaling law for the drag force with respect to the yield force. We have plotted the instantaneous drag force acting on the intruder at a given time as shown in Fig. <ref> for D=4.0d and the pulling speed 1.0√(gd). Although the instantaneous force (thin line) acting on the intruder fluctuates with time, the cumulative average of this drag force is well defined (thick line).
Then, we call this average force  the drag force on the intruder in the subsequent subsections.




 §.§ Yield force and its dependence on μ and the intruder size D


In the present study, we have introduced the yield force as the drag force on an intruder moving through granular bed in the zero limit of Fr.
As can be seen in the velocity field around the intruder (Fig. <ref>), there is a bulk reorganization of the particles around the intruder.
To obtain the yield force, we perform simulations at V=0.01√(gd) and regard the average force acting on the intruder as the yield force.
Note that it is computationally very expensive to move an intruder below this velocity.
The simulations are carried out for nine cases for various coefficients of friction μ= 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.75 and 0.9 for the intruder diameters D/d=1, 2, 3, 4, 5, 6, 8 and 10 in this study using System I.




The dimensionless yield forces F_ Y^*=F_ Y/(ρ d^3 g) with respect to μ and D^*+1 with D^*=D/d are plotted in Figs. <ref>(a) and (b), respectively.
As can be seen in these figures, the yield force increases with μ. At any moment, the force acting on the intruder can be split into two components: the normal forces due to the K_n part and the tangential forces due to the K_t part.
Here, the velocity dependent damping forces whose components are proportional to γ_n and γ_t does not play important role for small Fr.
However, we notice that the contact force between the intruder and the particle in contact increases several times as μ increases from zero to a finite value.
It should be noticed that the yield force starts to saturate beyond a certain coefficient of friction μ.
This saturation might be related to the upper limit of tangential force in the form of μ F_n where F_n is the normal component of force, i.e., if the tangential force exceeds this particular value, then the tangential force becomes constant.
However, this switching of force does not take place for high μ because the tangential force does not become as large as μ F_n. Hence, the yield force saturates after a certain value of μ.

To understand the mechanism behind the saturation of this yield force with μ, we have considered three different parameters. First is <|F_t|/|μ F_n|> which represents the time-averaged ratio of tangential and μ times the normal collision force for each intruder contact. Second is <F_n.x̂/(F_n+F_t).x̂> which is the fraction of drag force that comes from the normal collision force (x̂ is the unit vector along the moving direction of the intruder). It can be seen in Fig. <ref> that the second parameter decreases while the first one increases with μ. Moreover, we have plotted the third parameter (P(|F_t|/|μ F_n|≈1) which represents the fraction of particles in contact with the intruder undergoing slipping. It can be seen that only around 5 % of the particles undergo slipping at a high μ. This is in accordance with the saturation of F_Y with μ we discussed earlier. Importantly, we also see that the contribution of normal force to the drag changes from 100 % at the frictionless case to approximately 80 % at high μ although the increase in drag occurs several times than that. This proves that the presence of drag not only increases the frictional forces, but also the normal forces. This is possible as the particle in collision with the intruder shall be held strongly by the particles supporting it if μ is high which will also lead to a stronger force chain network.














Next, we plot the yield force against D^*+1 in a log-log plot in Fig. <ref>(b).
The reason we have used D^*+1 instead of D^* is that we want to plot the yield force against the effective diameter instead of just the intruder diameter.
In our work, the effective diameter is twice as much as the maximum distance between the intruder and a particle in the granular bed (of the average diameter d) that are in contact.
We find that yield force is proportional to (D+d)^ϕ_μ where ϕ_μ is a constant for a given coefficient of friction estimated from Fig. <ref>(b), where the dimensionless yield force F_ Y^* satisfies the scaling F_ Y^*/(f(μ)(D^*+1)^ϕ_μ)= const.
We also plot μ-dependence to express the function f(μ) used in Fig. <ref>(b) in the inset of Fig. <ref>(c).
The exponent ϕ_μ decreases with μ as shown in the Fig. <ref>.
The exponent ϕ_μ for various K_n is examined while keeping K_t/K_n (which is a function of the Poisson's ratio only), γ_n/K_n^0.6 (the term that characterizes coefficient of restitution curve), and γ_t/γ_n (the ratio of tangential and normal dissipation constants) as constants.
We can indicate that the exponent ϕ_μ is almost independent of K_n if we fix the other simulation parameters except μ.

Note that ϕ_μ=0≈ 2 in the frictionless limit can be understood if the yield force is proportional to the cross section of the effective sphere surrounding the intruder.
In the absence of the frictional forces, the gravity plays the major role i.e., the intruder moves away from surrounding grains which are settled in a stable position by the gravity.
The work done by the intruder for changing the potential energy of the intruder is not retrieved back leading to a finite yield force for μ=0.
In fact, if the intruder moves by a distance of x, all the particles inside a volume of π (D+d)^2 x/4 in front of the intruder has to be rearranged.
Since, in the zero velocity limit there is no momentum based collision forces or the frictional forces (μ=0) thus leading to the square dependence for μ=0.

It is important to understand how a given intruder contact displaces to make a way for the intruder. In a dense system with gravity such as ours, any given contact collides several times with the intruder. As the intruder moves, a given particle remains in contact for some time after which it separates. Due to gravity and a push from the particles at the back, the same particle may again be in contact with the intruder after some time. Let us say that t_total is the total time a given particle spends in contact with the intruder summed for several small collisions, t_coll is the average collision time for each contact (averaged for each such small collision) and N_coll is the average number of collisions assuming the particle has collided at least once with the intruder. We see that although t_total decreases with μ, N_coll is seen to increase with μ several times (refer to Fig. <ref>). Moreover, t_coll is also seen to decrease with μ. The decrease in ϕ_μ can be a consequence of a change in these parameters.  In Fig. <ref>, it can be seen that although the total number of collisions increases with μ, the overall time of contact t_coll decreases, thereby implying that as μ increases, the intruder experiences an overall less time of contact with the other particles in the system. This leads to deviation from a perfect square dependence of the drag on the Froude number.




 §.§ Yield force dependence on μ and the intruder depth h


In this subsection, we have considered the depth dependence of the yield force.
As shown in Fig. <ref>(a), we have verified that the yield force is proportional to h^α_μ for h/d=10, 25, 45, 50, 55, and 60 using system II.
Indeed, F_ Y^*/h^*α_μ is independent of h/d as shown in Fig. <ref>(c), which can be scaled as a universal horizontal line with the introduction of g(μ).
We also plot μ-dependence of the g(μ) in the inset of Fig. <ref>(c).
It is noted that the function g(μ) is different from f(μ) used in Sec. <ref> (see inset of Fig. <ref>(c)).
Similar to the case with ϕ_μ it is possible to scale the yield force with respect to the depth of the intruder. Then we study the scaling factor (D+d)^ϕ_μ for various h and D for few cases to confirm that ϕ_μ is independent of depth. With this in mind, we find that there exists a depth scaling factor α_μ as F_Y∝ h^α_μ(D+d)^ϕ_μ.



We also examine the α_μ for various μ and the non-linear normal spring constant K_n for a constant value of K_t/K_n, γ_t/γ_n, and γ_n/K_n^0.6 at seven different depths in our study (h/d=10, 25, 45, 50, 55, 60, and 65).
Unlike the case with the scaling factor ϕ_μ, α_μ increases as μ increases.
However, α_μ weakly depends on the non-linear spring stiffness K_n as shown in Fig. <ref>(b).
The simplest contribution of the depth is the hydrostatic pressure expressed by ρ g h which is roughly true for frictionless case, i.e., α_μ=0≃ 1.
Changing the value of the non-linear spring constant K_n or the Young's modulus E causes the change of α_μ because the overlap between two granular particles δ increases as the K_n decreases. Hence, we propose

    F_Y/ρ (D+d)^ϕ_μ h^α_μ g=Ψ

where Ψ is a function dependent of μ. We shall call Ψ as the yield parameter.


 §.§ Drag force in frictional system

In this subsection, we present the results of the drag force acting on a moving intruder moving through the granular bed at various velocities (in the positive x-direction) for μ=0.5 and various D with a constant depth of the intruder.
All the results in this subsection are based on the simulation of the system I.





The drag force with respect to the Froude number for μ=0.5 has been plotted in the Fig. <ref>(a), where F_ drag^*=F_ drag/(ρ d^3 g).
We find three regimes, I, II and III in the drag force trend.
The drag force in the regime I (Fr<1) is almost independent of Fr as can be seen in Fig. <ref>(b).
The force in this regime strongly depends upon the coefficient of friction μ (as shown in Secs. <ref> and <ref>) and D.
We call this regime the quasistatic regime as in Ref. <cit.>, though we cannot exclude the weakly velocity dependent drag from our data.
In Ref. <cit.>, a similar quasistatic regime was observed for Fr<1 for various μ.
We will study the dependence of drag force in quasistatic regime on the coefficient of friction later.
In a different experimental setup in Ref. <cit.>, a quasistatic regime was observed for Fr<1 though their definition of Fr is different from ours.
We observe that the drag force in this regime is almost identical to the yield force which is the minimum force to make the intruder move inside the granular medium in the  limit of zero Fr.
This yield force is synonymous to the yield stress in Bingham fluid which is the minimum stress to begin to flow under the shear.
The yield force is originated from the local jamming and buckling of force chain networks inside the granular bed.

In the regime II or the intermediate regime, i.e., 1<Fr<5, it seems that the drag force increases linearly with Fr.
Least square linear fit has been applied in the regime II as shown in the graph of the form a+b Fr, where a and b are fitting parameters.
Note that the slight deviation from linearity observed for the intruder of diameter D=6.0d is mostly because the rough bottom wall is only a few particle diameters apart from the intruder.
The regime II which nicely follows linear trend as a function of Fr has been observed for all cases in the depth study (Sec. <ref>) where the  system size of our simulation is larger and the wall effect can be suppressed by keeping intruder far away from the bottom wall.
This regime was also observed in Ref. <cit.>.
This regime II can be considered as a transition zone from the quasistatic regime (the regime I) to the regime III as mentioned in the next paragraph.

In the regime III, i.e., Fr>5, a quadratic dependence of drag force on Fr can be observed due to the direct momentum transfer by collisions of grains to the intruder.
We call this regime  the inertial regime. A quadratic fit of the form a^'+b^'Fr^2 (where a^' and b^' are fitting parameters for the quadratic expression) has been applied in this regime.
In a similar work for a two-dimensional system, where an intruder disk is dragged at a constant velocity under gravity <cit.>, a quadratic dependence of drag force on velocity of intruder was observed.
There was absence of quasistatic and linear regime at least for the range of parameters (d=0.32  cm, velocity = 10.3–309  cm/sec, depth of immersion of intruder =3.75–37.5  cm) they investigated.
The reason they missed  the first two regimes could be due to the limitation of number of simulation runs at the low Froude number. We also plot the first two parameters, <|F_t|/|μ F_n|> and <F_n.x̂/(F_n+F_t).x̂>,  discussed in Sec. <ref> using intruder diameter D=2.0 d in Fig. <ref> (c) and we observed that the two parameters remain the same as that of yield limit in the range of Fr we studied. This is important as we shall see in scaling of the drag force in the last subsection.


Figure <ref> shows the velocity field of the surrounding particles around the intruder for various D and the moving speed V_0.
One can see that the surrounding particles for V_0=0.1 √(gd) have more particles moving at speeds comparable to that of the intruder than for V_0=1.0 √(gd).
The same trend can be observed between the intruder velocities 1.0√(gd) and 3.0√(gd).
Figure <ref> shows the larger number of rearrangements at lower velocity around the intruder than that for larger V_0.




 §.§ Drag force in frictionless system

In the presence of the friction, the drag force is determined by two forces, the normal contact force F_n and the tangential contact force F_t.
In the presence of the tangential contact forces, stick-slip motion in granular media plays an important role, i.e., when the intruder is in contact with a particle, the tangential contact force keeps building up (stick) until it reaches a value of μ F_n, after which it starts slipping.
This allows the strain to be built up and then release in the form of rearrangements.
This phenomenon is quite active in the quasistatic regime. It is necessary to see how the drag force changes in the absence of these tangential contact forces.






All the simulations are carried out for the system I with μ=0 for all the contacts.
We have obtained the result for the drag force on the intruder for frictionless cases as shown in Fig. <ref>.
Even in the frictionless case, three regimes exist.
However, contrary to the case with friction with a constant force, the drag force in the first regime can be fitted by the form a^''+b^''Fr^0.5 as shown in Fig. <ref>.
The reason for the absence of the quasistatic regime here is that the velocity dependent  dissipative forces play dominant roles throughout the frictionless case because of the absence of the tangential contact force.
In other words, the main contribution is from the frictional force in this regime in the presence of friction.
However, we should stress that there exists a finite yield force even for the frictionless case.
This indicates that tangential force alone is not responsible for the yield force.
Nevertheless, we should note that the presence of tangential forces does increase the yield force multiple times (discussed earlier).
In the regimes II and III, we observe similar trends to those observed in the case with friction, i.e., we observe a linear and a quadratic dependences in II and III, respectively.
This is because the dominant force in these regimes is the normal contact forces as has already been discussed in the previous subsections.
The fact that the drag force follows the quadratic trend even in the case without friction also shows that the dominant force in this regime indeed is the collision based force and would be present even in the absence of the friction.
In the frictionless case, we observe that the crossover between the regimes (the drag force dependence on velocity of intruder) takes place at the identical Froude numbers as observed in the case with the friction.




 §.§ Drag force dependence on μ in the quasistatic regime





In the previous subsections, we have reported the quasistatic regime only for the frictional case at μ=0.5, while the quasistatic regime disappears in the frictionless case μ=0.
Therefore, we should clarify whether a quasistatic regime exists for small μ.
In this subsection, we study the drag forces in the quasistatic regime keeping the diameter of the intruder as D=2.0 d for various μ.

We have confirmed the existence of a constant force regime or a quasistatic regime for Fr≤ 1 except for μ=0 as seen in Fig. <ref>, where we have examined the cases with μ=0, 0.05, 0.1, 0.2, and 0.3.
Figure <ref> indicates the existence of the quasistatic regime even for the small coefficient of friction such as μ=0.05.
We suggest that the drag force on an intruder in the quasistatic regime is constant unless μ is zero.
Another important quantity to be addressed is the yield force (the force required to initiate the motion of an intruder), which is  non-zero even for frictionless systems.
From Fig. <ref>, the drag force in this quasistatic regime is equal to yield force for the cases of μ≥ 0.05.





 §.§ Drag force at various depth h




In a recent two-dimensional simulation study<cit.> on a depth dependence of drag force on an intruder disk reports that the onset of inertial regime (the quadratic dependence of the drag force on V) clearly depends on the depth of the intruder.
Even in three-dimensional system, the onset of the quasistatic, the intermediate, and the quadratic regimes might depend on the amount of grains located above the intruder in granular media under gravity, which acts as a hydrostatic pressure proportional to ρ g h, where h is the depth of the intruder.
Nevertheless, this simple picture might be only applicable to the frictionless case, because the exponent ϕ_μ for μ>0 is smaller than 2, and thus, we expect that the effective pressure is proportional to h^α_μ with α_μ>1 for μ>0.
Note that the sum rule ϕ_μ+α_μ≈3 should be satisfied because of the argument on the yield force in the end of the last subsection.
Therefore, it is important to study the scaling of the drag force on the depth of the intruder.

All the results presented in this subsection are based on the simulations for the system II.
We have performed the simulations at depths h/d=10, 25, 45, 50, 55, 60, and 65 from the free surface in negative z-direction with the intruder diameter D=4.0 d in a much larger system of dimensions 80d× 40d× 80d.
The value of μ is fixed to be 0.5, and all the other parameters such as K_n, K_t, γ_n, γ_t are the same as mentioned in the simulation methodology section for the most of arguments in this subsection.
Based on Fig. <ref>, we have confirmed that the onset of the different drag regimes starts at roughly the identical Froude numbers as before.
It is reasonable that if there are insufficient numbers of granular layers above the intruder, the drag may undergo a regime change differently and the quadratic regime may start quicker as has been seen in Ref. <cit.> for a two-dimensional system.
However, if the intruder is located sufficiently deep from the free surface, Froude number at which crossover between the regimes happens, seems to be almost independent of the depth.
Therefore, we expect the effect of the depth appears as the yield force.











 §.§ Drag force dependence on non-linear spring stiffness K_n, the dissipation coefficient γ_n, and the coefficient of restitution curve





In this subsection, we study how the values of the non-linear spring stiffness constant K_n changes the onset of regimes.
The simulations are performed at constant depth using the system II with h=45d.
We control K_n with keeping K_n/K_t, γ_n/K_n^0.6, and γ_n/γ_t constant.
Because K_n/K_t is purely a function of the Poisson's ratio, we fix the Poisson's ratio for the material.
Similarly fixing γ_n/K_n^0.6 corresponds to the fixed curve for the coefficient of restitution e (Fig. <ref>) for the three cases of K_n as shown in Fig. <ref>(a).
Finally, γ_n/γ_t is kept equal to unity to maintain uniformity in the normal and tangential dissipation factor in all the three cases.

We observe that, all the graphs for F_drag versus Fr seem to converge in the inertial regime, i.e., for Fr>5 if we choose an identical restitution coefficient curve (by fixing the value of γ_n/K^0.6_n).
We, however, notice that the drag force depends on the spring stiffness in the first two regimes.
The drag force becomes smaller for larger K_n in these regimes.
This is reasonable because the softer materials may form clusters (force chains) under the compression.
We know that the coefficient of restitution depends on the Young's Modulus E or the normal spring stiffness K_n <cit.>, but it is possible to vary the dissipation term γ_n under unchanged curve of the coefficient of restitution.
This is what we have carried out by keeping γ_n/K_n^0.6 constant for various K_n in Fig. <ref>(a) while we have varied γ_n/K_n^0.6 term in Fig. <ref>(b) for a fixed K_n=2×10^8 ρ d g and thus changing the coefficient of restitution curve.
As can be seen in our results, the drag force in the quadratic regime for various K_n is similar as long as the γ_n/K^0.6_n is kept constant irrespective of the value of K_n (see Fig. <ref>(a)) and it is different for various γ_n/K^0.6_n (Fig. <ref>(b)).
Hence, we can say that the forces in the quadratic regime depend on γ_n/K^0.6_n but not on the normal non-linear spring constant K_n or the normal dissipation γ_n term separately.
We also observe that the forces in the quasistatic regime depend on K_n but not on γ_n as shown in Figs. <ref>(a) and <ref>(b).
This is to be expected as the damping terms containing γ_n are associated with the velocity part and hence the drag independence on γ_n in the quasistatic regime. This was also discussed in the yield force study (which lies in the quasistatic regime).




From both Figs. <ref>(a) and <ref>(b), we also conclude that the crossover in the drag regime that we have observed throughout our study is independent of the material properties such as K_n and γ_n.
This is important as it sets a general rule for the regime changes in the drag force with respect to the Froude number used in our study, i.e., the drag force regime crossovers occur at Fr≈ 1 between the quasistatic to the linear regimes and at Fr≈ 5 between the linear to the quadratic regimes.



 §.§ Scaling of the drag force with respect to the yield force


In this subsection, we try to scale the drag force with respect to the intruder diameter and the depth. In the Sec. <ref>, we saw that the contribution of the normal collision force and the ratio of tangential and normal collision force remains constant against velocity of the intruder. Hence, indicating the presence of yield force in the scaling of the drag. Note that the data for the dynamical simulation with finite Fr are only for D^*≤ 6.
In Sec. <ref>, we have found that the yield force can be scaled as ρ (D+d)^ϕ_μ h^α_μ, where ϕ_μ and α_μ obey Fig. <ref> and Fig. <ref>(b), respectively.
We see that the yield force can be used as a scaling for the drag forces up to Fr=10, i.e., the scaled drag force (F_drag/F_Y-1) becomes independent of the size of the intruder (see Figs. <ref>(a) and <ref>(b)) as well as the depth (Fig. <ref>(c)).
If we are interested in the case for D≫ d, F_drag∝ D^ϕ_μ for any given Fr.
The existence of the scaling via the the yield force suggests that the drag force can be expressed as a product of the yield force and the dynamical part depending on Fr.
Indeed, three regimes and their crossover values do not depend on D and μ.
We, therefore, propose that the drag force can be expressed as (for μ≥ 0.05):

    F_drag/F_Y-1=
    
    0,     Fr<1 
    Linear ,     1<Fr<5 
    Quadratic    Fr>5,



We verify this by fitting a function of type F_drag/F_Y-1=a_i×Fr^b_i on a log-log plot for the frictional cases in all the regimes where a_i and b_i are fitting parameters for the regime i. We see that b_III is 1.7 and 1.9 for Fig. <ref>(a) and (c), respectively. In frictionless case, although the scalings are valid with respect to F_Y, these are applicable in regime II and III only when the fitting function is modified to the form a_i×Fr^b_i+c_i where c_i is additional fitting parameter.













§  DISCUSSION

The scaling law that we address in our work can play an important role in miniaturization of the drag forces on intruders to move through granular media.
The miniaturization of the drag in the granular media keeping the dimensionless number, would be important in industrial applications.
We, however, want to scale the diameter of the intruder D by the diameter of the grains d, because the numerical simulation can treat only dimensionless quantities such as D/d.
This can be successfully used to do a complete miniaturization of the intruders on a laboratory scale while keeping the diameter of the granular bed particles same.
In the present work, we have found that the drag force can be represented by a product of the dynamical part and the yield force, if the intruder is located enough deep.
Thus, the scaling law for the moving intruder in frictional grains can be expressed by

    F_drag/ρ h^α_μ(D+d)^ϕ_μg Ψ=
    
    1,     Fr<1 
    25/4λ/Ψ (Fr-1) + 1,     1<Fr<5 
    λ/ΨFr^2 + 1,     Fr>5

where the exponents α_μ and ϕ_μ satisfy the approximate sum rule α_μ+ϕ_μ≈ 3.
Unfortunately, we cannot explain μ-dependences of the nontrivial exponents α_μ and ϕ_μ.
The explanation of the exponents would be an important subject of our future work.

To express Eq. (<ref>), we make the following assumptions:
(i) in the quasistatic regime, i.e., for Fr<1, the force remains constant throughout and is equal to its yield force.
We have seen that this is not true for the case without friction.
However, for μ≥0.05, this is found to be true.
(ii) The graphs are continuous at both Fr=1 and Fr=5. This is what we have used to obtain the function in regime II using parameters from regime I and III (Ψ and λ, respectively).
(iii) The intruder is not located in shallow region.
In other word, if the location of the intruder is shallow, we cannot use Eq. (<ref>) because the drag force cannot be expressed as a product of the yield force and the dynamical part.

Moreover, based on our work in Sec. <ref>, <ref> and <ref>, we can say that 
    Ψ≡Ψ(K_n^*,μ)


    λ≡λ(γ_n^*/K_n^*0.6, μ)


We already saw in Fig. <ref>(a) that Ψ is heavily dependent on μ, while Fig. <ref>(a) and (b) showed that it is weakly dependent on K_n while is independent of γ_n. The weak dependence on K_n was seen because the yield force could be affected by the ability of the particles to form force chains which is affected by the softness or hardness of the particles. The independence from γ_n is a consequence of the fact that the damping term is coupled with the velocity component and thus gets nullified in the zero velocity limit. λ is seen to be heavily dependent on all the parameters. However, since it is associated with the inertial regime, we can couple two of the parameters γ_n^*/K_n^*0.6 together since the forces remain same in inertial regime if this ratio is kept constant irrespective of the K_n or γ_n (see Fig. <ref> (a) and (b)).





§ CONCLUSION

We have computed the drag force acting on an intruder moving with a constant velocity through the granular bed made up of slightly polydispersed particles. We have found that the drag force in the zero velocity limit could be scaled by the size and depth of the particle in a manner similar to that in fluids. We have seen how the parameters associated with these scalings change with frictional coefficient and other material properties. We have also investigated the saturation of yield force with frictional coefficient and related it to very low fraction of particles undergoing slipping at high μ. Alternatively, we have also studied the time spent by an average particle as an intruder contact and its successive collisions with the intruder for various μ. Also, we have related the depth scaling factor α_μ with μ for various material properties. It was observed that ϕ_μ+α_μ≈ 3. Moreover, ϕ_μ=2 and α_μ=1 for frictionless system suggesting a projected area dependence like in fluids.

We have found the drag force can be expressed as a product of the yield force and the dynamical part, if the location of the intruder is sufficiently deep.
In the dynamical part of the drag force, there are three regimes: the regime I, the intermediate and, the inertial regimes, at least, for the frictional cases, where the crossovers take place at Fr≈1 and Fr≈5.
For μ≥0.05, the regime I for Fr≤ 1 is a quasistatic regime, while the drag force in the regime I for the frictionless case depends on Fr.
Moreover, we have seen that the drag force in the inertial regime remains unchanged if the term that characterizes the coefficient of restitution curve, i.e., γ_n/K^0.6_n  is kept constant. A drag force in the velocity range we have studied can always be scaled with respect to the yield force.




§ ACKNOWLEDGMENT

We would like to thank IITG for providing excellent computational facilities through Param-Ishan.
K. Anki Reddy would like to thank YITP, Kyoto University for funding one month visit.
This work is partially supported by the Grant-in-Aid of MEXT for Scientific Research (Grant No. 16H04025).


plain

99
jaeger96
	H. M. Jaeger, S. R. Nagel, and R. P. Behringer,
	Rev. Mod. Phys. 68, 1259 (1996).
aranson06
	I. S. Aranson and L. S. Tsimring,
	Rev. Mod. Phys. 78, 641 (2006).
hilton
	J. E. Hilton and A. Tordesillas,
	Phys. Rev. E 88, 062203 (2013).
albert99
	R. Albert, M. A. Pfeifer, A. L. Barabasi, and P. Schiffer,
	Phys. Rev. Lett. 82, 205 (1999).
chehata03
	D. Chehata, R. Zenit, and C. R. Wassgren,
	Phys. Fluids 15, 1622 (2003).
wassgren03
	C. R. Wassgren, J. A. Cordova, R. Zenit, and A. Karion,
	Phys. Fluids 15, 3318 (2003).
bharadwaj06
	R. Bharadwaj and C. Wassgren,
	Phys. Fluids 18, 043301 (2006).
reddy11
	K. A. Reddy, Y. Forterre, and O. Pouliquen,
	Phys. Rev. Lett. 106, 108301 (2011).
guillard13
	F. Guillard, Y. Forterre, and O. Pouliquen,
	Phys. Rev. Lett. 110, 138303 (2013).
guillard15
	F. Guillard, Y. Forterre, and O. Pouliquen,
	Phys. Rev. E 91, 022201 (2015).
zhang15
	X. Zhang, D. Sheng, G. P. Kouretzis, K. Krabbenhoft, and S. W. Sloan,
	Phys. Rev. E 91, 022204 (2015).
gravish10
	N. Gravish, P. B. Umbanhowar, and D. I. Goldman,
	Phys. Rev. Lett. 105, 128301 (2010).
guo12
	H. Guo, J. Goldsmith, I. Delacruz, M. Tao, Y. Luo, and S. A. Koehler,
	J. Stat. Mech. Theory Exp. 7, P07013 (2012).
zhou05
	F. Zhou and S. G. Advani,
	Phys. Rev. E 71, 061304 (2005).
boudet10
	J. F. Boudet and H. Kellay,
	Phys. Rev. Lett. 105, 104501 (2010).
takehara10
	Y. Takehara, S. Fujimoto, and K. Okumura,
	EPL 92, 44003 (2010).
panaitescu17
	A. Panaitescu, X. Clotet, and A. Kudrolli,
	Phys. Rev. E 95, 032901 (2017).
jung17
	W. Jung, S. M. Choi, W. Kim, and H. Y. Kim,
	Phys. Fluids 29, 041702 (2017).
texier17
	B. D. Texier, A. Ibarra, and F. Melo,
	PLoS ONE 12, e0175412 (2017).
maladen09
	R. D. Maladen, Y. Ding, and C. L. D. I. Goldman,
	Science 325, 314 (2009).
Li09
	C. Li, P. B. Umbanhowar, H. Komsuoglu, D. E. Koditschek, and D. I. Goldman,
	PNAS 106, 3029 (2009).
mazouchova10
	N. Mazouchova, N. Gravish, A. Savu, and D. I. Goldman,
	Biol. Lett. 6, 398 (2010).
maladen11
	R. D. Maladen, Y. Ding, P. B. Umbanhowar, and D. I. Goldman,
	Int. J. Robotics Res. 30, 793 (2011).
pak11
	O. S. Pak and E. Lauga,
	Phys. Fluids 23, 081702 (2011).
marvi14
	H. Marvi, C. Gong, N. Gravish, H. Astley, M. Travers, R. L. Hatton, J. R. Mendelson, H. Choset, D. L. Hu, and D. I. Goldman,
	Science 346, 224 (2014).
sharpe15
	S. S. Sharpe, R. Kuckuk, and D. I. Goldman,
	Phys. Biol. 12, 012904 (2015).
hosoi15
	A. E. Hosoi and D. I. Goldman,
	Annu. Rev. Fluid Mech. 47, 431 (2015).
aguilar16
	J. Aguilar and D. I. Goldman,
	Nat. Phys. 12, 278 (2016).
takehara14
	Y. Takehara and K. Okumura,
	Phys. Rev. Lett. 112, 148001 (2014).
takada17
	S. Takada and H. Hayakawa,
	J. Eng. Mech. 143, C4016004 (2017).
seguin17
	A. Seguin and P. Gondret,
	Phys. Rev. E 96, 032905 (2017).
geng05
	J. Geng and R. P. Behringer,
	Phys. Rev. E 71, 011302 (2005).
baldassarri06
	A. Baldassarri, F. Dalton, A. Petri, S. Zapperi, G. Pontuale, and L. Pietronero,
	Phys. Rev. Lett. 96, 118002 (2006).
potiguar13
	F. Q. Potiguar and Y. Ding,
	Phys. Rev. E 88, 012204 (2013).
faug15
	T. Faug,
	Eur. Phys. J. E 38, 34 (2015).
seguin16
	A. Seguin, C. Coulais, F. Martinez, Y. Bertho, and P. Gondret,
	Phys. Rev. E 93, 012904 (2016).
candelier09
	R. Candelier and O. Dauchot,
	Phys. Rev. Lett. 103, 128001 (2009).
cundall79
	P. A. Cundall and O. D. L. Strack,
	Géotechnique 29, 47 (1979).
brilliantov96
	N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, and T. Pöschel,
	Phys. Rev. E 53, 5382 (1996).
zhang05
	H. P. Zhang and H. A. Makse,
	Phys. Rev. E 72, 011301 (2005).
kuwabara87
	G. Kuwabara and K. Kono,
	Jpn. J. Appl. Phys. 26, 1230 (1987).
ramirez99
	R. Ramírez, T. Pöschel, N. V. Brilliantov, and T. Schwager,
	Phys. Rev. E 60, 4465 (1999).
schwager98
	T. Schwager and T. Pöschel,
	Phys. Rev. E 57, 650 (1998).
silbert01
	L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton,
	Phys. Rev. E 64, 051302 (2001).
schwager08
	T. Schwager and T. Pöschel,
	Phys. Rev. E 78, 051304 (2008).
plimpton95
	S. Plimpton,
	J. Comput. Phys. 117, 1 (1995).
plimpton07
	S. Plimpton, P. Crozier, and A. Thompson,
	Sandia National Laboratories 18 (2007).
stukowski10
	A. Stukowski,
	Model. Simul. Mater. Sci. Eng. 18, 015012 (2010).
humphrey96
	W. Humphrey, A. Dalke, and K. Schulten,
	J. Mol. Graph. 14, 33 (1996).
kohlmeyer17
	A. Kohlmeyer,
	TopoTools (2017).
guillard14
	F. Guillard, Y. Forterre, and O. Pouliquen,
	Phys. Fluids 26, 043301 (2014).


