In this paper, the software PFC 5.0 [26] is used for numerical simulations. The rolling resistance linear model is adopted [27, 28] for the contact law, which is commonly used to investigate the rolling effect on irregular-shaped particles (non-sphericity in 3D or non-circularity in 2D) [29] in pseudo-static system.

The force-displacement law for the rolling resistance linear model calculates the contact force and moment as in Eq. 1.

where F _l is the linear force, and F _d is the dashpot force.

The linear force F _l includes the normal force F _n and the shear force F _s . The normal force F _n is calculated by Eq. 2.

where k _n and k _s are respectively the normal and shear stiffness, g _s is the surface gap between contact particles at the beginning of the timestep.

The shear force F _s is calculated by Eq. 3.

The rolling resistance moment M _r is incremented as:

where Δθ _b is the relative bend-rotation increment.

where k _r is the rolling resistance stiffness, and k _s is the shear stiffness.

The normal and tangential stiffness can be calculated by the effective modulus E* and normal to shear stiffness ratio κ* at the contact as shown in Eqs. 6, 7, respectively.

where A is the surface area of the particle, L is the distance between the center of two contacting particles.

The effective radius R ̄ is calculated in Eq. 8.

The rolling resistance moment M _r is updated, but it cannot exceed the limiting torque M _limit calculated in Eq. 9.

where the rolling resistance coefficient μ _r corresponds to the tangent of the maximum angle of a slope on which the rolling resistance torque counterbalances the torque produced by gravity acting on the particle.

Binary mixture samples are prepared in different particle size ratios between coarse gravels and fines sands(SR) with SR = 1:2, 1:5, and 1:10. The coarse particles are distributed randomly in the samples with different volume fractions (VF). The simulation process involves three stages: sample preparation, isotropic consolidation, and shearing stage. In the beginning, the binary particle samples are generated in a rigid boundary condition with a diameter of 1,000 mm and a height of 2,000 mm. The particle sizes with different particle size ratios include 10, 20, 50, and 100 mm. The numbers of particles in the samples with different particle size ratio and volume fractions are different. The maximum number of particles in the sample with SR = 1:10 and VF = 0% is about 23,000, while the minimum number of particles in the sample with SR = 1:10 and VF = 100% is about 3,000. The particles are generated by using falling rain method, that is firstly created at the top of the boundary, and then deposited at the bottom of the sample under the influence of gravity. Each particle is created with irregular shapes by randomly selecting one of particle templates, shown in Figure 1. After reaching the equilibrium of the particle assembly, the confining pressures are exerted on the samples for isotropic consolidation. Then, samples are sheared under the same confining pressure. The main parameters of binary particle can be found in Table. 1 [30]. In this table, the damping is the ratio of the energy dissipated after collision to the energy before collision. The critical damping constant is calculated by C c = 2 k m . Figure 2 presents the numerical binary particle samples of SR = 1:10 with different volume fractions of coarse particles 0, 20, 40, 60, 80, and 100%. According to the different packing states, Figures 2A–C are in the coarse particle floating state, Figure 2D is in the transitional state, and Figures 2E,F are in the coarse particle non-floating state.

DEM simulations of the biaxial test are carried out in the membrane boundary condition, which allows to observe the irregular deformation of the samples during the shearing process. The membrane boundary condition consists of two rows of membrane particles, as shown in Figure 3A. The diameter of the sticky membrane ball is 5 × 10–3 m by default. The external force on the membrane particles generated by confining pressure and the composition of membrane boundaries are exhibited in Figure 3B. The external force exerted on the membrane particles by confining pressure is calculated by Eq. 10:

where p is the confining pressure, D is the diameter of the membrane ball, n A O and n O B are unit vectors normal to the central lines that connect membrane ball O and its adjacent membrane balls A and B, respectively.

The deviatoric stress q is calculated by Eq. 11.

where F is the external force applied on the upper platens, S is the surface of the upper platens, and p is the confining pressure.

The axial strain is defined by Eq. 12.

where h is the current difference between the upper and bottom platens, and h ₀ is the initial difference between the upper and bottom platens.

The evolution of the volumetric strain of samples is investigated in the simulations. In order to accurately trace the volume deformation, the area geometric division method is proposed in this study. This method calculates the volume of the sample by areas of four parts, as shown in Figure 3C. The upper part S1 and the bottom part S2 maintain triangles during the whole simulations, while the left part S3 and the right part S4 are two zones composed of a set of small triangles. Each small triangle consists of two adjacent membrane ball center points and the center point of the sample. During the shearing process, S3 and S4 deform irregularly, so the area of these two parts is calculated by accumulating the areas of small triangles. By summing up the four parts of the sample, the accurate volume of the deformed sample can be obtained. The volumetric strain is defined by Eq. 13.

where V is the current volume of the sample, and V ₀ is initial volume of the sample.

The simulations are carried out under different confining pressures of 100, 200, and 300 kPa, respectively. The compression rate is fixed at 0.005 m/s, which meets the criterion of the inertia I = ε ̇ m p d smaller than 10^–3, where ε ̇ is the shear strain rate, m and d are the average mass and average diameter of particle, and p is the confining pressure [31]. The simulation stops when the shearing strain reaches 15%.

Figures 4A–C show the deviatoric stress curves of binary particle samples of SR = 1:10 with different volume fractions of coarse particles under confining pressures of 100, 200, and 300 kPa. For each curve, the deviatoric stress first increases to the peak value then decreases to the residual value. The deviatoric stress curve is gentle at the residual stage. It can be seen that the volume fraction of coarse particles has an important impact on the deviatoric stress. For the granular system in the coarse particle floating state, the peak value of deviatoric stress decreases with the increase of volume fraction of coarse particle, as the mechanical strength is generally dominated by small particles. Fine particles play a lubrication role between coarse particles, which offsets the increase of the shear strength. The turning point happens in the transitional state. The lowest peak value is obtained at about VF = 60%. Then the peak value of deviatoric stress begins to increase with the volume fraction of coarse particles. In the coarse particle non-floating state, the peak value keeps increasing with the volume fractions of coarse particle. It is because that fine particles can not fill the voids between coarse particles, which weakens the lubrication effect of fine particles, and that the strength is dominated by the coarse particles with the increasing of its volume fraction.

Figures 4D–F present the volumetric strain curves of binary particle samples of SR = 1:10 with different volume fractions of coarse particles under confining pressures of 100, 200, and 300 kPa. If samples are in the coarse particle non-floating state or the coarse particle floating state, the samples shrink first and then dilate. On the other hand, if the samples are in or close to the transitional state, the initial porosity of samples is relatively small. The samples shrinkage while hardly dilate. The larger the particle size ratio is, the more obvious this phenomenon is. The finding is consistent with the evolution of the packing efficiency of the three states [25].

The intrinsic parameters of binary particle samples can be calculated from stress-strain curves. The friction angle φ is calculated by Eq. 14.

where M is the ratio between q and p′ by M = q/p′.

Figure 5 shows the friction angle of samples of SR = 1:10 at peak and at residual state under different confining pressures. The friction angles initially reduce with the increasing volume fraction of coarse particles and reach a minimum corresponding to the volume fraction of coarse particles of 60% for mixtures. A further increase of volume fraction of coarse particles causes an increase in friction angle up to the values obtained for GF = 100%. The variations of friction angle versus volume fraction of coarse particles with SR = 1:2 and SR = 1:5 are not shown here because they nearly follow the same trends as the SR = 1:10.

Figure 6 shows the deviatoric stress and volumetric strain of binary particle samples with different particle size ratios SR = 1:2, SR = 1:5, and SR = 1:10 under the confining pressure of 100 kPa. It can be seen that the difference of deviatoric stress of samples with the same particle size ratio under the same confining pressure is not apparent. However, the influence of the volume fraction of coarse particles on the deviatoric stress still follows the same trend of first decreasing and then increasing. Figure 7 shows the friction angle of binary particle samples with various particle size ratios SR = 1:2, SR = 1:5, and SR = 1:10 at peak and at residual state under the confining pressure of 100 kPa. The friction angle versus the volume fraction of coarse particles generally follows the same trend as the deviatoric stress, but this trend is not obvious when SR = 1:2, then the concave extent of curve increases obviously with the increase of particle size difference. The lowest point of friction angle at peak and at residual stage both appear at VF = 60%, but the change of internal friction angle at peak stage is slightly larger than that in residual stress stage.

To better understand the evolution law of contact force distribution between particles under shearing, the statistical distribution maps of normal contact force and tangential contact force are plotted. The distribution of orientations of normal and tangential contact forces are visualized. Ref. 32 proposed Eq. 15 of particle-to-particle contact force distribution to fit the statistical results of particle-to-particle contact force of anisotropic assemblies.

where f _n (θ) and f _t (θ) are angular distribution functions of normal and tangential forces between particles, respectively; f ₀ is the measure of average normal contact force over all contacts calculated by the formula of f 0 = ∑ 0 2 π f n ( θ ) d θ ; a _n and a _t are second order Fourier coefficients which define the magnitude of average normal forces anisotropy and average tangential forces anisotropy, respectively; θ _n and θ _t are second-order principal direction of average normal contact force anisotropy and average tangential contact force anisotropy.

Figure 8 exhibits the evolution of orientations of normal contact force distribution between particles in samples with different volume fractions of coarse particles at the axial strain of 15%. The distribution of the average normal forces changes from a peanut shape to an ellipse shape then returns to the peanut shape. The evolution of the tangential contact force distribution for samples of SR = 1:10 with different volume fractions of coarse particles at the axial strain of 15% is shown in Figure 9. The variation of the average tangential force distributions is obvious. The non-uniformity is more evident when volume fraction of either coarse or fine particles is low. As the content of the two tends to the average, the uniformity of the distribution of the average tangential force becomes larger.

The shear band can be regarded as a sign of failure for granular samples in the shearing process. DEM simulations of the shear band have been reported in the literature by particle rotation field [33, 34], particle displacement field [35, 36], porosity densification [35], as well as incremental deviatoric strain [37, 38]. As we know, particles in the shearing zone have a more frequently changed contact status than those in other zones, so the variations of contact status can be applied for the observation of shear bands. In DEM simulations, the contact force is represented by the segment of the line that connects the centroids of two contact particles. The existence of contact force chains has been demonstrated by photoelastic tests [39]. The thickness of the force chain represents the intensity of contact force. Generally, large force chains often exist between coarse particles, while small force chains often exist between fine particles.

During the shearing process, the arrangement and the contact status of particles are constantly changing, and the force chain between particles is incessantly generating and disappearing. Therefore, the change of the shear band can be determined by observing the newly generated contact force chain, which has the merit of presenting the shear band more accurately and sensitively. Figure 10 exhibits the evolution of the shear band in the monodispersed particle sample with VF = 0% under the confining pressure of 100 kPa. All the force chains in the samples have been counted. The evolution of the formation of shear bands can be clearly observed. As the axial strain increases, the non-uniform dilatation of the sample starts from the axial strain of 5%. When the axial strain increases to 7%, the first shear band starts to develop. When the axial strain grows to 9%, the first shear band is fully developes, and the top part of the second shear band appears on the top left corner of the sample. When, the axial strain grows to 13%, the bottom part of the second shear band starts to develop. Finally at 15%, two shear bands form an X-shape. The observation is consistent with the finding of Ref. 40 who indicated that the failure in all specimens can be characterized by two distinct and opposite shear bands.

Figure 11 shows the contact force chains of samples with different volume fractions of coarse particles at the axial strain of 15%. It can be seen that the shear band for the homogeneous sample without coarse particles is more clear. The contribution of the fine particles cannot be negligible even in small content. When the average spacing of coarse particles is around two times of the small particle size, the contribution of coarse particles disappears [4]. With the increase of the volume fraction of coarse particles, it becomes more and more difficult to form the shear bands. This could be due to the very small sample-particle size ratio. The shear band can be generated when the width of shear band is very small compared to the sample’s size. The number of force chains decreases, while the strength of force chains increases significantly. The coarse particles carry relatively large contact forces, which can be observed in Figure 11F. It is expected that the sample with VF = 100% has a higher shear strength than samples with other volume fractions of coarse particles.