Drag of a cylindrical object in a two-dimensional granular environment

The drag of a cylindrical object in a two-dimensional granular environment is numerically studied. It is found that the drag law is fitted by the sum of the yield force and the dynamic force, the latter of which is reproduced by a simple collision model. The angular dependence of the radial stress on the surface of the object is given by the Gaussian below the yield force. The probability of the velocity drops of the object is investigated above the yield force, where this probability is independent of the packing fraction and the drag force.


Introduction
To understand the resistance force from granular materials is important for application [1]. The force acting on an object in granular materials is known to be different from that in fluids characterized by the drag coefficient. The drag in granular materials is studied in many setups and situations [2][3][4][5][6][7][8]. The previous studies [4] reported that the drag is proportional to the square of the moving speed in the high-speed regime. On the other hand, the low-speed regime shows various velocity dependence depending on the density of the system [5,6]. It is also noted that the granular jet experiments [7] and simulations [8] suggest the perfect fluidity of the granular jet flows.
We know that the object cannot move when the drag force is sufficiently smaller than a threshold, which is determined from interactions with the surrounding particles. This threshold is related to the Coulombic friction among particles or between particles and the bottom plate [5]. Although there are not so many papers studying this regime as far as we know, it must be important to understand this regime because there exist cracking processes or other phenomena. To this end, we study the stress field around the object by performing the discrete element method (DEM).
It is also noted that the high-speed regime is not simple. The movement of the object fluctuates, and there exist avalanches as time goes on. The information on the avalanche of a certain quantity might be related to the sheared granular systems [9]. We will also check this property in this paper.

Model and setup
In this section, we briefly explain our model and setup of the simulation. We prepare one spherical object with mass Figure 1. Setup of our simulation. The object is set at the origin when we start the simulation. The wall particles are aligned at y = ±L y /2. The arrow indicates the direction of the drag. For later discussion, we also introduce the polar coordinates (r, θ), whose origin is set at the center of the object. convenience, we regard i = 0 as the object, i.e., d 0 = D and m 0 = M. To avoid the crystalization, the surrounding particles are bidisperse with the ratio d : 1.4d. Correspondingly, the mass ratio is given by m : 1.4 2 m. Here, we choose the mass density of the object as the same as that of the surrounding particles, i.e., M = (D/d) 2 m. The equations of motion for the object and the surrounding particles are described by where F drag is the drag force acting on the object, µ b is the Coulombic friction between the surrounding particles and the bottom plate, g is the gravitational acceleration, andê x andû i are the unit vectors parallel to the x-direction and u i , respectively. Here, the interparticle force F i j between i-th and j-th particles is described by F i j = F n i j + F t i j . In this paper, we assume that the normal (tangential) interaction F n i j (F t i j ) is described by the sum of the linear spring with the spring constant k n (k t ) and the dashpot with the dissipation rate η n (η t ) [10]. We also consider the sliding friction in the tangential direction with the friction coefficient [10].
In the simulation, we basically choose D = 10d, N = 2000, k n = k t , η n = η t , µ = µ b = 0.3, and g = 1.0 × 10 −4 k n d/m. We apply the periodic boundary condition in the x-direction, and we put smaller particles at y = ±L y /2 as the bumpy boundary condition in the y-direction [5]. In the following, we perform ten simulations for each set of parameters to evaluate the ensemble averages.

Results
In this section, we present our results obtained from the simulations in the subsequent three subsections. First, we measure the steady speed V of the object when we control the drag force F drag . As shown in Fig.  2, there exists a threshold F Y , and the object cannot move for F drag F Y because the drag force is balanced with the interparticle forces from the surrounding particles. Here, we can regard F Y as the yield force [5,6]. From the simulations, F Y 1.1 × 10 −2 k n d and 1.8 × 10 −2 k n d for ϕ = 0.75 and 0.80, respectively. For F drag F Y , on the other hand, the object can move. From Fig. 2, the relation is fitted by with F dyn ∝ V 2 , which is also observed in high speed regime in many experiments and simulations [4][5][6]11]. Let us describe the second term F dyn (V) in Eq. (2) by a simple collision model [6,11,12]. We consider the coordinate that the object is set at the origin. In addition, we assume that the trajectories of the surrounding particles are ballistic before colliding with the object. In this situation, the impulse from the particle to the object in the xdirection is given by ∆p x = {[(2/3)+e n −e t /3] cos 2 θ +(1+ e t )/3}m r V (e.g., see Ref. [13]), where m r = Mm/(M + m) is the reduced mass with respect to the object and the particle, and θ is the angle from the x-axis (see Fig. 1). Here, the surrounding particles are regarded as monodisperse because the dispersity is much smaller than the object size. Summing up the impulse over the semicircle of the object (−π/2 ≤ θ ≤ π/2), we obtain Here, we have used that the (two-dimensional) volume of the collision cylinder is proportional to (D+d)V. The inset of Fig. 2 shows that this model (3) well reproduces the simulation results for ϕ = 0.75, but the deviation increases for ϕ = 0.80. This might be because our assumption of the trajectory is invalid for denser situations.

Stress field for
In this subsection, let us investigate the stress field around the object for F drag F Y . Let us define the stress obtained from the simulations as with the area S for calculating the stress. We note that the time and ensemble-averaged field is coarse-grained with the scale 2d for visibility [14].   Figure 3 shows the profile of σ rr on the surface of the object (r = D/2), where we have introduced the polar coordiantes (r, θ), whose origin is set at the center of the object (see Fig. 1). Here, we have found that this profile is well fitted by the Gaussian where the drag force dependencies of θ 0 and C are given in Fig. 4. The coefficient θ 0 ( 1.1) is insensitive to ϕ and F drag as shown in Fig. 4(a). On the other hand, the coefficient C is proportional to the drag force F drag (see Fig. 4(b)). This is understood as follows: by integrating σ rr cos θ over the surface of the object as 1 then, we get This shows that the coefficient C is proportional to the drag force F drag while it is independent of ϕ as shown in Fig.  4

(b).
Let us approximately rewrite Eq. (5) in Fourier series as Because the higher terms become small, we only consider the first two terms: At the surface of the object (r = D/2), this is consistent with the Airy stress function known in elasticity [15]: 1 We have used the following integral and the following approximation (θ 0 = 1.1) Here, the error caused from this approximation is sufficiently smaller than unity.
Using this function χ, we can obtain the stress components in the Cartesian coordinates as (a) (b) Figures 5 show the spatial distributions of the stress components σ xx and σ xy for ϕ = 0.80 and F drag = 1.6 × 10 −2 k n d. Here, σ xy obtained from the Airy function (11) well reproduces that from the simulation, while σ xx from Eq. (11) is approximately twice smaller than that from the simulation. We note that the magnitude of the discrepancy for the latter is almost independent of the drag force. This deviation might come from the fact that we cannot capture the behaviors of σ rθ and σ θθ from the Airy function (10).

Probability of the velocity drops for
In this subsection, we show the statistical properties of the velocity of the object for F drag F Y . As shown in Fig.  2, the object moves with the steady speed V (> 0) for F drag F Y while its instantaneous value fluctuates. Figure  6 shows a typical evolution of the velocity of the object in the x-direction, where the velocity drops appear intermittently.
Now, let us measure the drop probability. The velocity drop is defined by the magnitude of the decrease of the velocity, where the slope of the velocity is negative. This is analogous to the stress drop reported for sheared systems [9]. Here, because we control the drag force and measure the instantaneous velocity, we consider the velocity drop of the object. Figure 7 shows the probability distribution of the velocity drops for various object sizes, packing fraction, and drag force. Here, we have found that the distribution is independent of both ϕ and F drag . It is also found that the distribution has two peaks around {(D/d) 2 +1}∆V/V 10 0 and 10 −3 . The larger peak is understood as follows: When we consider an inelastic head-on collision between the object and one surrounding particle, the velocity difference of the object is written as where we have regarded 1 + e n ∼ 1 and m i m for simplicity. This description is true because the probability is scaled in terms of this as shown in Fig. 7.

Discussion and conclusion
In this paper, we have performed the two-dimensional drag simulations in terms of the DEM. First, we have clarified the threshold, below which the object cannot move due to the interaction with the surrounding particles. Above the threshold, the dynamic part of the drag force is found to be reproduced by a simple collision model when the system is not too dense. Below the threshold, the stress field around the object is analyzed. We have shown that the Airy stress function shows good agreements with the simulation results for σ xy , while the theoretical estimation of σ xx is twice smaller than that from the simulation. Above the threshold, on the other hand, the velocity drop statistics have been studied, where the probability shows two-peaks around {(D/d) 2 + 1}∆V/V 10 0 and 10 −3 , the former of which is understood from the simple model. In this paper, we have assumed that the Airy stress function is given by Eq. (10). However, various terms such as r −n cos(nθ) or r −n+2 cos(nθ) (n is an integer) are known to be included in the Michell solution [15,16]. Numerical estimation of these components will clarify the system more precisely, but this is our future work.