Study of Shear Behavior of Binary Mixtures by DEM Simulation of Biaxial Test in the Membrane Boundary Condition

This paper aims to study the shear behavior of binary particles with irregular shapes by discrete element method simulations of the biaxial test in the membrane boundary condition. Binary particle samples are generated according to different volume fractions of coarse and fine particles. The deviatoric stress and volumetric strain curves are plotted to describe the contracting-dilatancy relationship of binary samples under shearing conditions. The anisotropy of the normal and tangential contact forces are explored by visualization of the orientation of contact forces to describe the evolution of micro structures of samples during the shearing process. Besides, the formations of the shear band are observed by the visualization of the newly generated contact force chains between particles. The research shows that the volume fraction of coarse particles and particle size ratio have significant influences on the shear behavior of binary particles both in macroscopic and microscopic points of view. Moreover, the increased volume fraction of coarse particles leads to a more difficult formation of a shear band.


INTRODUCTION
Binary mixtures are special granular materials which consist of coarse particles and fine particles [1][2][3]. Due to the different attributes of coarse and fine particles (such as sizes, shapes, and volume fractions, etc.), the macroscopic behavior of binary mixtures distinguish with homogeneous materials [4][5][6]. Previous researches show that the strength of binary mixtures are significantly affected by the particle size ratio and the volume fraction of coarse particles [7,8]. However, the interactions between coarse and fine particles at microscopic scale has not been studied in depth. As experimental test is difficult to reveal the mechanism of granular materials in particle scale, and it often requires a lot of manpower and financial resources, the numerical simulations have become an alternative to study this problem.
In 1979, Cundall and Stack [9] first proposed the discrete element method (DEM), which has become an efficient tool to investigate the mechanical behavior of granular materials from both macroscopic and microscopic perspectives [10][11][12][13]. Many researchers have investigated the effect of particle size for various types of particles, including monodispersed [14,15], binary mixtures [16][17][18][19], polydispersed [20,21], as well as graded particles [22][23][24]. It is reported that there exist three packing states of coarse particles in the binary mixtures, including the coarse particle floating state, the transitional state, and the coarse particle non-floating state [2]. For the coarse particle floating state, most coarse particles are separated by fine particles, and the force chain networks are dominated by fine particles. For the coarse particle non-floating state, fine particles can hardly fill the voids between coarse particles, and the force chain networks are dominated by coarse particles. The transitional state is between the above two states. Ref. 25 pointed out that the maximum packing efficiency was obtained when there are enough fine particles to separate the coarse particles from each other. However, previous studies mainly studied the binary mixtures from macroscopic point of view, but the evolution of micro-structure of binary mixtures during the shearing process has not been fully discussed, and the formation mechanism of shear band in binary particle samples is still blank. Therefore, this issue is worthy of further study.
In this study, the shear behavior of binary particle samples is studied by DEM simulation of the biaxial test in a flexible membrane boundary condition in 2D. The coarse and fine particles are generated with irregular shapes. Binary samples with different volume fractions of coarse particles and different particle size ratios are investigated on the shear behavior. The rose diagrams are plotted to explore the anisotropy of orientations of normal and tangential contact forces, which helps us to better understand the macroscopic behavior of binary mixtures at particle scale. The comparisons of the shear band of binary samples in different states are visualized by the contact force chains as well.

Contact Law
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.
F n k n g s (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.

Preparation of Binary Particle Samples
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 km √ . 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.

Shearing Process of Biaxial Test
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:  Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 745721 7 where p is the confining pressure, D is the diameter of the membrane ball, n AO and n OB 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 0 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 0 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 pd 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%. 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. 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.

Anisotropy of Force Chain Orientation
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 0 is the measure of average normal contact force over all contacts calculated by the formula of f 0 2π 0 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.

Shear Band
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 nonuniform 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.

CONCLUSION
In this paper, numerical analysis on the shear behavior of binary particle samples is carried out using the DEM simulation of the biaxial triaxial test in the membrane boundary condition. The samples are prepared with different volume fractions of coarse particles and particle size ratios. The formation of the shear band has been systematically investigated through monitoring the disappearance and subsequent formation of contacts. Some essential conclusions can be drawn as follows: • The biaxial model in the membrane boundary condition can successfully describe the stress-strain relation of binary samples. Especially, the volumetric strain of irregular deformation of sample during the shearing process can be exactly recorded by the area geometric division method proposed in this research. • In the macroscopic point of view, the volume fraction of coarse particles and particle size ratios have been demonstrated to have important impacts on the macroscopic shear behavior of binary particle samples. For the granular system in the coarse particle floating state, the peak value decreases with the increase of volume fraction of coarse particles. For particles in the coarse particle non-floating state, the peak value increases monotonously with the volume fraction of coarse particles, the turning point happens in the transitional state. The friction angle at peak and at residual state are also affected by the volume fraction of coarse particles and particle size ratios. • In the microscopic point of view, The orientations of normal contact forces and tangential contact forces are applied to study the anisotropy of binary particle samples. An increased mixing degree between binary mixtures enhances the non-uniformity of the samples. • The evolution of the shear band can be followed via the vanishing and newly generated contacts between particles. The shear band in the sample with lower volume fraction of coarse particles is more obvious than samples with higher volume fraction of coarse particles. The binary particle samples with low coarse content are more likely to form shear bands.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
KW: Conceptualization and Writing-Original Draft WS: Methodology JD: Writing-Review and Editing XZ: Formal analysis SL: Supervision.