Dilatancy of frictional granular materials under oscillatory shear with constant pressure

We perform numerical simulations of a two-dimensional frictional granular system under oscillatory shear confined by constant pressure. We found that the system undergoes dilatancy as the strain increases. We confirmed that compaction also takes place at an intermediate strain amplitude for a small mutual friction coefficient between particles. We also found that compaction depends on the confinement pressure while dilatancy little depends on the pressure.


Introduction
Granular materials consist of a collection of distinct macroscopic particles such as sands and glass beads. The evolution of the particles obeys Newton's equation with repulsive force between particles. In reality, due to the roughness of the particles, mutual friction between particles is unavoidable in granular systems. We call such particles frictional particles. However, we sometimes ignore the mutual friction as an idealistic model. We call such particles frictionless particles. So far many theoretical studies have only been conducted for frictionless particles. Previous studies, however, pointed out that the mutual friction causes drastic changes in rheology such as the emergence of discontinuous shear thickening (DST) [1][2][3][4]. Bi et al. found that the shear jammed state is only observable in frictional particles [5]. Various researchers studied shear jammed state and DST in systems taking into account the mutual friction [6][7][8][9][10][11][12][13].
It is well known that dilatancy which is the volume expansion of a collection of particles takes place in granular systems [14][15][16][17], when shear is applied to them. More dilatant structure is realized for the high shear rate [18][19][20][21]. Previous studies also reported the existence of compaction in which repeated oscillations to the system make particles denser and more ordered configurations [22][23][24][25][26][27]. Although some studies investigated the dilatancy and compaction at a fixed friction coefficient [23,26], it is not clear how the dilatancy and compaction depend on the friction coefficient. Thus, this paper, which is complementary to the previous paper [28], focuses on appearances of dilatancy and compaction of frictional particles under oscillatory shear.
The contents of this paper are as follows. In the next section, we explain the setup of our numerical simulation. * e-mail: daisuke.ishima@yukawa.kyoto-u.ac.jp In Section 3, we present our numerical results to clarify the conditions of occurrence of dilatancy and compaction. In the last section, we summarize our results.

Set up of our simulation
We consider a two-dimensional system containing N granular particles with the mutual friction characterized by Coulomb's friction coefficient µ. To prevent crystallization, the system contains four-dispersion particles whose diameters are d 0 , 0.9d 0 , 0.8d 0 , and 0.7d 0 , respectively. There are N/4 particles for each particle size [2,29,30]. Since we assume that the density of each particle is identical, the mass of each particle is proportional to the square of its diameter. In our study, m 0 stands for the mass of a particle with its diameter d 0 [28].
We ignore the effect of gravity throughout this paper. We adopt the periodic boundary condition in the shear (x−)direction and introduce mobile walls in the vertical (y−)direction. The top and bottom walls consist of particles of diameter d 0 , in which the number of particles of each wall is N w and their centers of gravity are located at x ± G (t) = ±A sin(Ωt) and y ± G (t) = ±L y (t)/2, with the amplitude A and angular frequency Ω of oscillation (see Fig.1).
Since the system size L y (t) under the influence of the confinement pressure P depends on time, we introduce an effective shear strain: where L 0 is the linear size of the shear direction. We adopt the equations of motion of the walls at y ± G (t): where ±, m w , v ± w,y , P ± w , and ξ d are indices referring to the wall at y ± G (t) , mass of the wall m w = N w m 0 , velocity in the y−direction of the wall at y ± G (t), internal pressures acting on the walls at y ± G (t), and damping coefficient, respectively 1 .
Let m i , r i , I i = m i d 2 i /8, and ω i be the mass, position, moment of inertia, and z−component of rotational velocity of particle i, respectively. The equations of motion for translation and rotation are, respectively, given by where we have introduced the contact force f i j and torque T i j exerted by particle j on particle i. The torque T i j in Eq. (4) satisfies where t i j is the tangential vector between particles i and j.
We adopt the Cundall-Strack model for the contact force between particles [31,32], where the contact force f i j is expressed as with d i j = (d i + d j )/2 and r i j = |r i − r j |. Here we have introduced the normal force f i j,n , the tangential force f i j,t , and Heaviside's step function Θ(x) satisfying Θ(x) = 1 for 1 The damping coefficient ξ d is neccesary to stabilize the motion of walls.
x ≥ 0 and Θ(x) = 0 otherwise. The normal repulsive force is modeled as an elastic force with a linear spring constant k n and the dissipative force with a damping constant ξ n . Then, we model the normal force as f i j,n = k n δ i j,n n i j − ξ n u i j,n , where we have introduced n i j = r i j /r i j and u i j,n = (u i j · n i j )n i j with r i j = r i − r j , u i j = dr i j /dt and δ i j,n = d i j − r i j 2 . On the other hand, the tangential force f i j,t takes two different expressions depending on the magnitude of the tangential force f i j,t = | f i j,t |: where δ i j,t is the tangential displacement between the particles. We have introduced the spring constant in the tangent direction k t , damping coefficient in the tangent direction ξ t , magnitude of the normal force f i j,n = | f i j,n |, and tangential velocity at the contact point c i j := u i j − u i j,n + t i j (d i ω i + d j ω j )/2. In our simulation, we adopt k t = k n /2, ξ n = (m 0 k n ) −1/2 , ξ t = ξ n , and ξ d = ξ n [13,19]. The control parameters of our system areP := P/k n , γ 0,eff , Ω √ m 0 /k n , and µ. We deal withP from 2.0 × 10 −5 to 2.0 × 10 −3 , γ 0,eff from 1.0 × 10 −6 to 1.0, and µ from 0 to 1.0. Furthermore, we mainly investigate that N = 4, 000, and Ω/(2π) = 1.0 × 10 −4 √ k n /m 0 in this paper. To know Ω dependences see Ref. [28].
To prepare the initial configuration, we place particles with diameters of 0.6d 0 , 0.5d 0 , 0.4d 0 , and 0.3d 0 at random and then increase the diameter of the particles to reach the area fraction φ ini = 0.82. Next, the confinement pressure is applied to both walls to achieve a steady state. Then, oscillatory shear is applied through the walls [28]. We use the symplectic Euler method with the time step ∆t = 0.05 √ m 0 /k n . We have averaged the data over 10 cycles after N ini = 10 cycles to remove the effect of specific initial configurations. We set t = 0 at the moment of the end of N ini cycles.
Let us introduce the area fraction φ(Ωt,P, γ 0,eff ) To characterize the density change we also introduce δφ: δφ(P, γ 0,eff ) :=φ(Ωt = 2nπ,P, γ 0,eff ) − φ 0 (P) where n and φ 0 (P) are a non-negative integer and the area fraction without shear atP, respectively. Since our system is confined by the constant pressure, the density depends on time. Then, we focus on the excess fraction δφ at the instance of zero strain and the maximum strain rate.

Results of our simulation
Let us present the results of our simulation for various µ. We found that dilatancy is always observable for large strain amplitude (see Fig. 2). We also found that dilatancy is enhanced as the friction coefficient increases. This result may be related to the fact that the region of shear jamming is more extensive for large friction coefficients [13]. Note that the excess fraction δφ defined by Eq. (10) does not linearly decrease with the shear rate, as in dilatancy in steady sheared systems under constant pressures [19,20]. It is noteworthy that compaction appears at an intermediate strain amplitude for small friction coefficients.
To clarify the dependence of compaction onP and µ, we plot δφ max , the maximum value of δφ at fixedP, against µ for two differentP in Fig. 3. We confirm that δφ max is almost independent ofP for µ ≤ 0.01 but it strongly depends on the pressure for µ ≥ 0.01. Moreover, compaction disappears for µ > 0.1 atP = 2.0 × 10 −3 while it survives even for µ = 0.4 atP = 2.0 × 10 −4 .
From now on, let us focus on the results for µ = 1.0 to discuss the pressure dependence. From Fig. 4, δφ is zero for γ 0,eff < 10 −4 , while it takes finite value for γ 0,eff > 10 −4 . It is remarkable that δφ seems to depend on the pressure only a little. In the previous paper [28],  we reported that the softening point characterized by the critical strain amplitude γ S corresponding to the bending point of shear modulus can be scaled by pressure, while the critical strain for dilatancy little depends on the pressure. Therefore, our result suggests that there is no direct connection between the density change and the behavior of shear modulus. Note that the trajectories of particles are reversible for the strain amplitude a little larger than γ S as in Refs. [34,35]. Our results suggest that the critical condition of dilatancy is not related to the softening point but the yielding point.

Concluding remarks
In conclusion, we conducted a numerical study of frictional granular systems under oscillatory shear confined by constant pressure. We found that compaction can be observed only for small µ at an intermediate strain amplitude, while dilatancy always takes place at large strain amplitude. We confirmed that dilatancy is a continuous nonequilibrium phase transition and little depends on the confinement pressure.
One of the most important findings in this paper is that the density changes such as dilatancy and compaction are not directly connected with the softening. It will be important to confirm the validity of this picture in the near future. As another important point, our preliminary study suggests that the results strongly depend on the initial packing fraction φ ini . Detailed studies on the φ ini -dependence will be reported elsewhere.