Discrete Element Simulation of the Shear Behavior of Binary Mixtures Composed of Spherical and Cubic Particles

: This research paper presents an investigation into the shear behavior of binary mixtures composed of cubic and spherical particles, employing the discrete element method (DEM) through triaxial tests simulations. A range of binary particle samples with varying volume fractions of cubic and spherical particles is generated for analysis. The study primarily focuses on examining the contracting-dilatancy relationship of binary granular material samples by scrutinizing deviatoric stress and volumetric strain curves, while considering the inﬂuence of conﬁning pressure, initial porosity, and particle size ratio. Furthermore, the paper sheds light on the evolution of microstructures during the shearing process by presenting coordination numbers and rotational velocity ﬁelds for different particle types (overall particles, cubic particles, spherical particles), as well as between cubic-spherical particles. The ﬁndings demonstrate the substantial impact of both the volume fraction of cubic particles and the particle size ratio on the shear behavior of binary particles at both macroscopic and microscopic scales. Additionally, a comprehensive investigation reveals the dependence of anisotropy in normal contact forces, tangent contact forces, and contact orientations on the volume fraction of cubic particles.


Introduction
Granular materials have garnered significant attention as a distinct state of matter, finding wide-ranging applications in nature, everyday life, manufacturing, and technology. Examples include sand, soil, sugar, salt, and various industrial chemicals. The mechanical behavior of granular materials is influenced by multiple factors, such as particle size [1], particle shape [2], particle arrangement [3], and inter-particle friction [4]. However, due to the intricate interplay of these factors, isolating the specific impact of each one presents a challenge for researchers. To overcome this challenge, binary mixtures have emerged as a valuable approach for investigating the mechanical behavior of granular materials, wherein a specific influential factor is varied between two distinct states. By employing binary mixtures, researchers can effectively examine the isolated effects of individual factors, enabling a deeper understanding of granular material behavior.
Researchers have dedicated considerable efforts to investigating binary mixtures of granular materials using both experimental and numerical simulation approaches. Previous studies have highlighted the significant influence of volume fraction on the shear behavior of granular materials, wherein particle size heterogeneity affects packing density, shear strength, and energy dissipation. Moreover, the anisotropic fabric and stress tensors of granular media are directly related to integrated packing properties, which in turn, depend on particle size distribution. Additionally, it has been observed that fine and coarse particles play distinct roles in the formation of contact force chains. In the presence of coarse particles, fine particles dominate the contact force chains when in a floating state, whereas coarse shape. Binary samples comprising both spherical and cubic particles with varying volume fractions are prepared, and their shear strength and volumetric strain are analyzed from a macroscopic perspective. In addition to the macroscopic analysis, the microscopic metrics such as mechanical coordination number, contact force chains, and the anisotropy of normal and tangential forces are applied. These comprehensive analyses shed light on the underlying mechanisms governing the behavior of binary mixtures and provide valuable insights into the role of particle shape in shaping their macro and micromechanical properties.

Contact Law
In this study, the PFC 3D software [30] (https://www.itascacg.com/software/PFC, accessed on 20 June 2023) is applied to perform discrete element method (DEM) simulations. The Hertz contact model is utilized as the contact law, which consists of a nonlinear Hertzian component based on the theory of Mindlin and Deresiewicz to calculate the normal forces [31]. The interaction force between two contacting particles is determined by their overlap, and a dissipative component is also incorporated.
The normal force can be calculated as follows The coefficient h n depends on the geometrical and mechanical properties of the contacting pieces as where G, ν, andR are the effective shear modulus, Poisson's ratio, and contact radius, respectively. The effective shear modulus and Poisson's ratio can be input directly, or inherited from the contacting pieces as detailed below. The effective radius of the contact is computed via the radii of the contacting piece The tangential component F s is based on Coulomb's friction law, which is given by Equation (4) as follows: where (F h s ) 0 is the hertz shear force at the beginning of the timestep and ∆δ s is the relative shear increment. k s is the initial tangent shear stiffness, which depends on the current normal force as The value of (F h s ) 0 depends on the scaling mode M s : where (F h n ) 0 , (F h s ) 0 , and (k s ) 0 are the normal force, the shear force, and the tangent shear stiffness computed at the beginning of the timestep, respectively. M s = 1 means that the shear-force scaling is active.
Compute the shear strength: The Hertz shear force is calculated as Update the slip state If the slip state is true, then the contact is sliding. Whenever the slip states change, the slip change callback event occurs.

Generation of Cubic Particles
The cubic shape is in fact a variant of superellipsoid. The surface function of a super ellipsoidal particle in local Cartesian coordinates is defined in Equation (9): where x, y, z are the coordinates of a point on the surface of the superellipsoidal particle. The variables, a, b, c denote the half-length of the particle along the principal axis of the local coordinate system of the superellipsoidal particle. The shape indexes m and n are two shape factors to determine the particle shape and the surface blockness. The shape index m affects the upper and bottom margins, while the shape index n controls the four side edges [32,33]. Simultaneous variation of m and n transforms the sphere into a cube. Keeping m constant, increasing n leads to sphere-to-cylinder conversion. Constant n and increased m yield walnut-like tetrahedron. m influences upper/lower surfaces, and n defines lateral shape. Cross-sectionally, augmented m shifts the upper surface from circular to square; increased n shifts the side from circular to square. Thus, simultaneous m and n increase results in sphere-to-cube transformation. Figure 1 shows the evolution path of superellipsoids from sphere to cube. Figure 1. Evolution of superellipsoids from sphere to cube by changing two shape factors m and n to determine the particle shape and the surface blockness.

DEM Simulation Procedure
The aim of this study is to investigate the shear strength of binary mixture samples using three-dimensional discrete element method (DEM) simulations of triaxial tests under the boundary condition of a rigid cylindrical wall in dry condition. The influence of water has not been considered in the simulations. The simulation procedure consists of three steps: sample generation, isotropic consolidation, and shearing process. Six binary mixture samples with different volume fractions of cubes and spheres, denoted as CxSy (where x and y represent the volume fraction of cubes and spheres, respectively) are examined. The samples under investigation are C100S0, C80S20, C60S40, C40S60, C20S80, and C0S100 (see Figure 2 for sample visualization). Each sample has a diameter of 50 mm and a height of 125 mm, with approximately 5000 particles of the equivalent diameter of 4 mm for both spherical and cubic particles. The L/D ratio of the sample is 2.5. The initial porosity of the samples is approximately 0.45 monitored by measuring a circle in the middle of the sample. Default parameters for DEM simulations are employed (refer to Table 1 for specific DEM parameters). Following the sample generation, a stress servo-controlled rigid cylindrical wall boundary condition is applied to regulate the confining pressure throughout the isotropic consolidation and shearing stages. The servo-controlled boundary adjusts the radius of the rigid cylindrical wall in response to changes in lateral pressure, utilizing the Lame formula of Equation (10) [34] (see Figure 3) .
where ∆r is the radius variation, E is Young's modulus of the cylinder, r 1 and r 2 are, respectively, the internal and external radius of cylindrical boundary condition, and r is the radius between r 1 and r 2 . The confining pressure P 1 is determined by dividing the total force exerted by the particles in contact with the cylinder boundary by the lateral surface area of the cylinder. To achieve the desired pressure, the servo-controlled rigid cylindrical boundary adjusts P 1 at each time step to match the target pressure P 2 . The isotropic consolidation process is considered complete when P 1 reaches the target pressure P 2 . During the shearing process, the confining pressure P 1 is maintained at a constant value while the servo-controlled cylindrical wall is adjusted. Simultaneously, the sample is subjected to compression by the upper and bottom platens. The shearing rate is controlled by the inertia number I (Equation (11)), which determines the quasi-static state of the granular system [9].
whereε is the shear strain rate, M and d are the mass and diameter of the particle, and P 1 is the confining pressure. In our simulations, the shear strain rate is fixed at 1.0 /s to maintain an inertia number I below 10 −3 . This choice ensures that the granular system remains in a quasi-static state throughout the shearing process. The axial strain is continuously monitored to assess the progress of shearing. The simulation is concluded once the axial strain reaches a predetermined threshold of 30%. The effective mean stress p and the deviatoric stress q are defined as follows: where σ 1 denotes the axial stress obtained by dividing the applied force on the upper platen by the surface area of the upper platen. The axial strain ε z and volumetric strain ε v are derived from the movements of rigid walls as follows: where h 0 is the initial height of the sample, h is the current height of the sample, v 0 is the initial volume of the sample, and v is the current volume of the sample.

Influence of Volume Fraction of Cubic Particles
Figures 4a and 5a illustrate the evolution of deviatoric stress as a function of axial strain for binary mixtures under different confinement conditions of 1.0 MPa and 2.0 MPa, respectively. The data exhibit a strain-softening behavior with increasing axial strain, which aligns with previous investigations on dense granular materials [36][37][38]. The porosity values of the sample vary from 0.44 to 0.51 during the simulation. Notably, when the volume fraction of cubes (FC) decreases from 100% to 80%, the deviatoric stress at the peak rises from 3.8 MPa to 4.2 MPa. However, as FC continues to decrease, the deviatoric stress at the peak decreases from 4.2 MPa to 3.9 MPa. For FC values below or equal to 60%, the deviatoric stress at an axial strain of 30% remains relatively constant. In contrast, the deviatoric stress at the residual state generally exhibits an increasing trend as FC increases. Consequently, the volume fraction of cubes significantly influences both the peak and critical-state shear strength of the mixtures. Furthermore, as the confining pressure increases from 1.0 MPa to 2.0 MPa, the shear strength of granular mixtures nearly doubles. This finding underscores the substantial impact of confining pressure on the shear behavior of granular mixtures.   Figures 4b and 5b present the evolution of volumetric strain as a function of axial strain for the binary mixtures subjected to confining pressures of 1.0 MPa and 2.0 MPa, respectively. Initially, the mixture samples undergo slight contraction followed by dilation. As the sample approaches the residual state, the slope of the volume strain curve gradually decreases, consistent with previous findings. This behavior can be attributed to the easier compaction of spherical particles, which necessitates more void space for shear, leading to enhanced dilation in granular systems with lower cubic particle content. Consequently, the volumetric strain of the sample increases as the volume fraction of spheres during shearing increases. Moreover, the confining pressure influences the volumetric strain response. Samples exhibit greater resistance to deformation and dilation under higher confining pressures, as shown by the consistently higher volumetric strain under 1.0 MPa compared to that under 2.0 MPa.
The deviatoric stress and volumetric strain curves obtained from the triaxial shear test provide valuable insights into the intrinsic properties of granular materials, as demonstrated in Figure 6. These mechanical parameters, such as the friction angle (ϕ), dilatancy angle (ψ), Poisson's ratio (ν), and stiffness modulus (E) offer further understanding of the shear behavior of granular media under specific test conditions. The friction angle ϕ can be obtained from the deviatoric stress curve as a function of the effective stress ratio M , Equations (16) and (17): where q denotes the deviatoric stress, p indicates the effective average stress and M denotes the effective stress ratio. The dilatancy angle ψ is determined by Equation (18) at an inflection point of the axial strain-volumetric strain curve as illustrated in Figure 6.
The relationship between the peak internal friction angle (σ peak ) and critical internal friction angle (σ residual ) with FC is depicted in Figure 7. As FC increases, the friction angles at peak initially increase until reaching 80%, after which they start to decrease. However, the friction angle at the residual stage consistently increases with FC, indicating that particle shape has a greater influence on the strength of the granular material in the residual stage.
Additionally, the dilatancy angle shows a decreasing trend as FC increases.   Figure 9 illustrates the comparison of simulation results for different initial porosities, highlighting the significant impact of initial porosity on the shear behavior of granular materials. A clear inverse correlation between shear strength and initial porosity is observed, indicating that an increase in initial porosity results in a decrease in the deviatoric stress at the peak. However, the effect of initial porosity on shear strength at the residual state is not evident, as the shear strength eventually converges as shearing progresses. Moreover, the volumetric strain shows an inverse relationship with initial porosity, where the dilatation of the sample volume weakens with increasing initial porosity.

Influence of Particle Size Ratio
The influence of particle size ratio on the shear behavior of binary mixtures has received considerable attention in previous studies. However, there has been a noticeable gap in investigating the impact of particle shape on such behavior. To address this research gap, the effect of particle shape on the shear behavior of binary samples has been focused on, which consist of cubes and spheres with varying sizes. In this study, five mixtures are generated with a volume ratio of 1:1 for cubes and spheres, as outlined in Table 2. Sample SC represents a case where the particle sizes of cubes and spheres are equivalent. In Samples SS1 and SS2, the size of cubes is maintained while reducing the size of spheres. Conversely, for Samples SC1 and SC2, the size of spheres is kept constant while decreasing the size of cubes.  Figure 10a presents the deviatoric stress profiles of the five prepared mixture samples. It is observed that the deviatoric stress decreases as the particle size decreases, irrespective of the particle shape. For example, when comparing samples SS1/SC1 or SS2/SC2 with the same particle sizes, their peak strengths are nearly identical. However, the stress at the peak for samples SS1/SC1 is higher. Additionally, samples SC1 and SC2 exhibit a delayed attainment of peak strength compared to SS1 and SS2, indicating that smaller cube particle sizes result in a lower Young's modulus. Although the shear strengths of the five samples are almost identical at the residual state, suggesting that particle size has minimal effect on the residual state, the deviatoric stress curve exhibits the characteristic stick-slip phenomenon. This phenomenon can be attributed to variations in force networks under shear stress [39][40][41][42]. Previous studies have indicated that reducing the size of cube particles can amplify the amplitude of stick-slips. Figure 10b presents the volumetric strain profiles of the five mixture samples with equivalent volume fractions of cubic and spherical particles but varying particle sizes under a confining pressure of 1.0 MPa. The results demonstrate that reducing the particle size has a substantial influence on sample expansion, as samples with smaller particle sizes display smaller expansions. This effect can be attributed to the ability of small particles to efficiently fill the void spaces between larger particles without significantly displacing the surrounding particles. Furthermore, it is observed that cubic particles exhibit less dilation compared to spherical particles of the same size. This behavior arises from the inherent characteristics of the particle shapes: spherical particles have a higher tendency to form dense packing structures, while cubic particles of the same volume exhibit a comparatively looser arrangement.

Coordination Number
The coordination number, defined by Equation (20), serves as a metric for quantifying the contact number of individual particles and provides insights into the packing density at the particle level [8,9,43].
where Z is the coordination number, C is the total contact number, and N is the total particle number. However, it is worth noting that certain particles may have either no contact or only one contact with adjacent particles during the shearing process, rendering them less influential in maintaining a stable stress state. To address this issue, the mechanical average coordination number, denoted as Zm, is introduced (as shown in Equation (21)), which excludes the contribution of these floating particles.
where Z m is the number of total particles, Z mc is the number of cubic particles, and Z ms is the number of spherical particles. N 1 and N 0 are the numbers of particles with only one or no contacts, and the subscripts c and s represent cubic and spherical particles, respectively. Figure 11a presents the temporal evolution of the mechanical average coordination number, denoted as Z m , for mixture samples with different volume fractions of cubic and spherical particles under a confining pressure of 1.0 MPa. The initial coordination number, Z m , is relatively high and gradually decreases during the shearing process. Furthermore, it is observed that the coordination number, Z m , decreases as the volume fraction of cubic particles increases. To further investigate the coordination number, Figure 11b showcases the evolution of the mechanical average coordination number, Z mc , specifically between cubic particles, for samples with varying volume fractions of cubic particles under the confining pressure of 1.0 MPa. Interestingly, it is evident that Z mc decreases as the volume fraction of cubic particles decreases. In contrast, Figure 11c demonstrates the evolution of the mechanical average coordination number, Z ms , specifically between spherical particles, under the confining pressure of 1.0 MPa. The results indicate that Z ms increases with a decrease in the volume fraction of cubic particles. Moreover, Figure 11d illustrates the evolution of the mechanical average coordination number, Z msc , between cubic and spherical particles for samples with varying volume fractions of cubic particles under the confining pressure of 1.0 MPa. It is noteworthy that Z msc exhibits different initial values and initially increases, followed by a subsequent decrease as the volume fraction of cubic particles increases. This behavior arises due to the fact that the contact number between cubic and spherical particles is highest when the number of particles of both shapes is equal.

Rotational Velocity Field
The average rotational velocity is an important micromechanical metric for analyzing the mechanical behavior of binary mixture samples, and is defined as the average angular velocity of particles in a sample [43][44][45]. During shearing, the rotational velocity field of particles with different shapes is expected to differ. Figure 12 shows the variation of the average particle rotational velocity field for different mixture samples. The results indicate that the average particle rotational velocity increases as the FC decreases. Furthermore, Figure 13 depicts the rotational velocity of cubic and spherical particles during shearing. The rotational velocity of particles increases with the volume fraction of the corresponding particle shape. Previous studies [46,47] have suggested that particle rolling can offset the strengthening effect of inter-granular friction, while smaller particles tend to undermine the inter-particle locking effect due to their lubrication effect. In our study, cubic particles are demonstrated to tend to inhibit particle rotation.
The average rotational velocity serves as a crucial micromechanical parameter for analyzing the mechanical behavior of binary mixture samples, representing the average angular velocity of particles within a sample [43][44][45]. During shearing, it is expected that particles with different shapes will exhibit varying rotational velocity fields. Figure 12 presents the variation of the average particle rotational velocity field for different mixture samples. The results demonstrate that the average particle rotational velocity increases as the volume fraction of cubic particles decreases. Furthermore, Figure 13 displays the rotational velocity of cubic and spherical particles during shearing. It is observed that the rotational velocity of particles escalates with an increase in the volume fraction of the corresponding particle shape. Notably, previous studies [46,47] have indicated that particle rolling can counterbalance the strengthening effect of inter-granular friction, whereas smaller particles tend to weaken the inter-particle locking effect due to their lubrication effect. In our study, cubic particles have been found to tend to impede particle rotation, supporting these prior findings.    Figure 14 illustrates the contact force chains in cubic-spherical mixture samples with varying volume fractions of each particle shape. Each force is represented by a color line connecting the particles in contact, with the color indicating the magnitude of the contact force. Strong forces are depicted in red, while weak forces are shown in blue. The results demonstrate that as the volume fraction of cubic particles increases, the strength of a single force chain intensifies, but the force chain network becomes more dispersed. Conversely, as the volume fraction of spherical particles increases, the strength of a single force chain diminishes, but the force chain network becomes more compact. These findings indicate that a decrease in the volume fraction of cubic particles leads to a decrease in the strength of individual force chains, while simultaneously increasing the density of the force chain network.

Fabric Anisotropy
The utilization of DEM simulation provides a significant advantage in examining the micro-mechanical behavior at the particle scale by capturing crucial microscopic parameters, including the anisotropy of normal and tangential forces, as well as contact orientations. Following the shearing process, the anisotropy of normal and tangential contact forces, along with contact orientations, was analyzed for mixture samples comprising cube-cube, cube-sphere, and sphere-sphere contacts. To visualize the spatial distribution of contacts, a 3D rose diagram was constructed based on the collected data, as shown in  Each bar in the diagram represents the average contact number, average contact normal force, and average contact tangential force in the respective 3D direction. Figure 15 presents the anisotropy of normal stress among different contact types in the mixture samples. Although the anisotropy of normal stress gradually decreases with decreasing FC, the reduction is not significant. In contrast, the anisotropy of tangential stress increases as FC decreases, as depicted in Figure 16. The shape of the 3D rose contact diagram transitions from symmetrical cones to "X" shapes. Figure 17 showcases the anisotropy of contact orientations for mixture samples involving cube-cube, cube-sphere, and sphere-sphere contacts after shearing. The data indicates that contact orientations become more anisotropic when one particle type dominates the system, while they become more isotropic when the volume fractions of cubes and spheres are similar. The rose diagrams representing the anisotropy of these parameters offer valuable insights into the microscopic behavior of mixture samples.

Conclusions
The primary aim of this study is to explore the shear behavior of binary mixture samples using DEM simulations of triaxial tests. These samples are composed of cubic and spherical particles with different volume fractions and particle size ratios. The investigation analyzes the mechanical characteristics of the binary mixture samples from both macroscopic and microscopic viewpoints. The obtained results yield several significant conclusions, which can be summarized as follows: • The macroscopic shear behavior of binary particle samples is profoundly influenced by the volume fraction of cubic and spherical particles, as well as the particle size ratios. Variations in the volume fraction of these particles result in noticeable changes in deviatoric stress, volumetric deformation, as well as internal friction angles, dilatancy angles, stiffness modulus, and Poisson's ratio of the binary samples. • From a microscopic standpoint, the analysis reveals the influence of particle characteristics on the coordination number of the individual particles. As the volume fraction of coarse particles increases, the overall coordination number of the system decreases. More specifically, the coordination number Z mc of cubic particles decreases as the volume fraction of spheres decreases. Conversely, the coordination number Z ms of spherical particles increases as the volume fraction of spheres decreases. These findings demonstrate the intricate relationship between particle characteristics and the coordination number at a microscopic level. • The anisotropy of normal and tangential contact forces, as well as contact orientations, is influenced by the volume fraction of particles with different shapes. The variation in the volume fraction of particles with distinct shapes leads to changes in the directional distribution of contact forces and orientations.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.