Rheology of two-dimensional crushable granular materials

The rheology of two-dimensional crushable granular materials under shear is numerically studied using the discrete element method. We find that the mean fragment size changes as the shear strain increases while the shear stress is almost independent of this mean size. The fragment size distribution is found to follow a power law. In particular, the exponent in the intermediate fragment size regime becomes approximately −11/6, which is almost independent of the shear rate.


Introduction
To understand the rheology of granular materials is important not only in industrial usages but also in physics because they exhibit various phenomena such as jamming transition [1,2]. The existence of the discontinuous change and the hysteresis of the shear stress for frictional granular materials are reported when the shear rate changes [3]. However, these previous studies focused on the systems where the component particles are never broken into pieces. For practical usage of granular materials such as soil or rock engineerings [4], we should treat the effect of the breakage of particles. To this end, we should investigate how the brittleness of particles affects the rheology of the system.
There are some papers studying the effect of the fragmentation of particles [5][6][7][8]. Ben-Nun et al. [5] considered a particle, which is replaced by a group of smaller particles when the compression pressure acting on the particle is larger than a threshold. They reported that the contact force distribution follows log-normal, which is different from that observed in uncrushable systems [9]. Tsoungui et al. [6] performed similar simulations and reported that the relation between the external pressure and the strain shows good agreement with experiments.
Cheng et al. [7] used another model. They collect many small particles and glue them together, and then they regard them as one "macroscopic" particle. The bonds between particles are removed depending on the probability determined by the stress. They performed uniaxial compression simulations using this model and reported that the distribution of the crushing stress is in good agreement with that in laboratory experiments. McDowell et al. [8] used a similar model to investigate how the fragmentation pattern of a single macroscopic particle is affected by the * Corresponding author, e-mail: takada@go.tuat.ac.jp defect fraction inside the particle when uniaxial compression is applied to the particle. They reported that the fragmentation mode changes after the defect fraction reaches 30% of the total number of constituent particles.
In this paper, we consider the hybrid model of the previous two approaches. We consider many small particles and glue them together. The bond between particles is broken when the force acting on them satisfies a certain threshold. If the bond is broken, this bond is never reproduced in the subsequent simulation. We apply shear to a collection of the above particles and investigate how the brittleness of the particles affects the rheology.
The organization of this paper is as follows: we explain our model and setup of the simulation in the next section. Section 3 is the main part of this paper, where we show various results obtained from the simulations. In Sec. 4, we discuss our results and draw some conclusions.

Model and simulation setup
In this section, let us explain our model and setup of the discrete element method (DEM) simulation [10][11][12]. First, we prepare macroscopic two-dimensional crushable particles composed of small particles, whose mass and diameter are m i and d i , respectively. To avoid crystallization, we prepare two sizes of component particles (d and 1.4d) with the same number ratio, where the corresponding masses are m and 1.4 2 m, respectively. When two particles interact with each other, the interparticle force in the normal direction is given by Fig. 1, where k n and η n are the spring constant and the magnitude of the dashpot in the normal direction, respectively, d i j = (d i + d j )/2 is the sum of the radii of the two colliding particles, r i j = |r i j | = |r i − r j | is the relative displacement between i-th and j-th particles,r i j = r i j /r i j , u i j = u i − u j , and Θ(x) is the step function. Here, ε represents the existence of the bond between particles, i.e., ε = 1 (0) means that there exists (does not exist) a bond. The parameter u relates to the strength of the bond, where u is assumed to be small. In the same manner, the tangential force is described by where k t and η t are the spring constant and the magnitude of the dashpot in the tangential direction, respectively, µ is the Coulombic friction coefficient, dt is the tangential displacement from the beginning of the contact [11,12].
At the beginning of the simulation, the bonds between particles are set if the relative distance between particles is less than (1 + u)d i j . During the simulation, a bond between i-th and j-th particles is broken when is satisfied with certain thresholds F n 0 and F t 0 for normal and tangential directions, respectively [13]. Once this occurs, this bond is removed subsequently [14].
Using these particles, we prepare macroscopic crushable particles [13]. Hereafter, we refer to these macroscopic particles as "clusters" for simplicity. Even in this case, we adopt two sizes to avoid crystallization, where each crushable particle is made of 500 or 1000 small component particles with the packing fraction ϕ = 0.80 with the number ratio of clusters 1 : 1. It is noted that the corresponding sizes of the clusters are 30.4d and 43.0d, respectively. When we put component particles to make clusters, the maximum overlap between component particles is smaller than 10 −6 d using the same procedure as reported in Ref. [15].
When the simulation starts, shear is applied to the system with the shear rateγ with the aid of the Lees-Edwards boundary condition [16]. Hereafter, we measure the quantities as a function of the shear strain γ =γt.
In this paper, we use k t = k n , η n = η t = 1.0 √ mk n , µ = 0.3, F n 0 = F t 0 = 5.0 × 10 −3 k n d, and u = 5.0 × 10 −3 . The corresponding normal and tangential restitution coefficients are e n = e t = 4.3 × 10 −2 . We note that we perform the simulations 10 times for each set of parameters to evaluate the ensemble averages of the quantities.

Results
In this section, the simulation results are presented. Figure 2 shows the typical snapshots of the simulations foṙ γ * = 1.0 × 10 −3 and 10 −4 , where we have introduced the dimensionless shear rateγ * ≡γ √ k n /m. Initially, almost all particles are glued together with the bonds (Fig. 2(a)). When the shear rate is small (γ * = 1.0 × 10 −4 ), the macroscopic particles are compressed, and the particles whose brittleness is relatively large are broken as shown in Fig.  2(b). Therefore, there exist both large and small particles after a sufficiently large shear strain. On the other hand, when the shear rate is large (γ * = 1.0 × 10 −3 ), almost all macroscopic particles are broken due to the compression as shown in Fig. 2(c) and (d). The above trend is also observed by the evolution of the mean fragment size k = k kN k / k N k with the shear strain as shown in Fig. 3, where N k is the number of fragments whose size is k. Here, two particles belong to the same fragment if there remains a bond between them. We note that the initial value is approximately given by 500, because there exist isolated particles even in the initial condition. As the shear is applied, the mean size starts decreasing by around γ = 0.5. In the transient regime (γ 0.5), the non-affine deformation becomes dominant, and it enhances collisions between clusters. Accordingly, the number of bond breakages also increases. Indeed, the mean fragment size decreases with k (γ) ∼ γ −β with 0 β 2 depending on the magnitude of the shear rate in the range 10 −5.0 ≤γ * ≤ 10 −2.5 . As the mean size becomes small, this exponent gradually decreases to zero. After large shear strain, the system behaves as non-crushable granular materials.  We also evaluate the fragment size distribution after applying the shear strain γ = 5 as shown in Fig. 4. The distribution for k 10 2 is well fitted by a power law distribution, although large fragments are observed for lower shear rate (seeγ * = 10 −4.5 in Fig. 4). For the larger shear rate (γ * 10 −2.5 ), the exponent of the distribution becomes larger because the particle breakage is dominant. On the other hand, the exponent in the wide shear rate regime is well fitted by −11/6. It should be noted that this value is the same as the steady-state distribution in the collision cascade [17,18] and observed in our solar system such as the distributions of asteroids [17,19] and ring particles [20]. This result might suggest that there is some commonality between two different systems although our system does not include the effect of coagulation processes. Figure 5. Shear stress (P xy ) field forγ * = 10 −3 and γ = 1.0, where P xy is evaluated from Eq. (4) and we have introduced P * xy ≡ P xy /k n . The color indicates the magnitude of the shear stress.
Next, let us check the shear rate dependence of the shear stress. The two-dimensional stress tensor is defined by with the area S for calculating the stress tensor. We note that the kinetic contribution of the stress is much smaller than its collisional contribution as far as we have checked.
The shear stress is localized as shown in Fig. 5. This yields the following: When the magnitude of the shear stress locally increases, the fragmentation process occurs at that place, and then the shear stress decreases. Eventually, this produces that the global shear stress is almost independent of the shear strain. Here, the global shear stress is evaluated by Eq. (4) with S = L 2 (L = 233.4d is the system size). This trend is observed in Fig. 6(a), where the evolution of the global shear stress is plotted as a function of the shear strain for several shear rates. It should be noted that the mean fragment size is decreasing in this regime as shown in Fig. 3. We also show the shear rate dependence on the shear stress when we fix the value of the shear strain. When the applied shear strain is small, the shear stress shows a weak dependence on the shear rate (see Fig. 6(b)), whose shear rate dependence is much weaker than the Bagnold scaling. We note that this trend is also observed in the solid regime for non-crushable granular materials [3]. On the other hand, as the applied shear strain increases, the shear stress slightly depends on the shear rate forγ * 10 −4 .

Discussion and Conclusion
In this paper, we have investigated the rheology of crushable granular materials using two-dimensional simulations. As the shear strain increases, the cluster breakage proceeds, where the breakage depends on the magnitude of the shear. It is also noted that a universality class of the fragmentation process is reported for thin glass beads [21] and glass tubes [22]. The relation between these systems and ours is important, but this is for future work.
In this system, we have performed two-dimensional DEM simulations, where the particles are bonded together when the initial distance is less than a threshold. As the shear is applied, the fragmentation proceeds, and the decrease of the mean size of fragments follows a power law in the transient regime. On the other hand, the shear stress is found to be almost independent of the fragmentation. We note that a similar evolution of the shear stress is observed in rock experiments [23]. We have also found that the fragment size distribution shows a power law, where the exponent for the intermediate shear rate is equal to −11/6. Here we have considered the rheology under constant volume conditions, where the packing fraction of the whole system is not changed. Indeed, the previous studies investigated the rheology of the system under a constant pressure condition [24], and the µ-I rheology is discussed [25,26]. To provide an insight into experimental setups, the comparison with the rheology under constant pressure is our important next work.