Discrete Element Modelling of the Mechanical Behavior of Sand–Rubber Mixtures under True Triaxial Tests

Sand–rubber mixtures (SRMs) consisting of stiff sand particles and soft rubber particles are typical binary mixture materials that possess a variety of complicated properties. The complexity of the properties of sand–rubber mixtures is increased when complex stress path is involved. This study investigates the mechanical behavior of sand–rubber mixtures under generalized loading conditions using the discrete element method. A series of numerical true triaxial shear tests were conducted on pure sand and sand–rubber mixtures. The effect of rubber content and loading path on both of the macroscopic and microscopic performances of sand–rubber mixtures was investigated, and the associated microscale mechanism was also discussed. Numerical simulations show that the relationship between the peak friction angle ϕp and the intermediate principal stress ratio b is influenced by the addition of rubber particles, and a suggested explanation of this phenomenon is that the rubber particles mainly affect the inherent stability of the strong network. Particle-scale observations, including the coordinate number, the proportion of strong contacts, and the fabric anisotropy, are also presented in this study. Microscopic results confirm the explanation above, and explore the force transmission characteristics of sand–rubber mixtures under generalized loading conditions. This research can provide a reference for the constitutive model development of sand–rubber mixtures.


Introduction
Sand-rubber mixtures (SRM) recently have received considerable amount of interest from both the research community and practicing engineers, and have been widely used in geotechnical engineering, including lightweight backfill, retaining walls, highway embankments, and road construction, due to its light weight, high damping, and high permeability [1][2][3][4][5][6]. Sand-rubber mixtures are unconventional geo-material and possess a variety of unique properties. The mechanical behavior of sand-rubber mixtures depends not only on the properties of host sands but also on the characteristics of rubber particles. A full understanding of the mechanical properties of sand-rubber mixtures can provide a solid backup for the application of sand-rubber mixtures in geotechnical projects.
A lot of experimental and numerical research has been performed to investigate the mechanical behavior of sand-rubber mixtures. Previous results show that the mechanical behavior of sand-rubber where k n and k s are the normal and shear stiffness, respectively, g s is the overlap, (F s ) 0 is the linear shear force at the beginning of the timestep, and ∆δ s is the adjusted relative shear-displacement increment which is nonzero only when g s is negative. The shear force follows a slip model, which is defined by the friction coefficient µ at the contact. The contact is sliding when |F s | ≥ µ·F n is satisfied. The rolling resistance moment M r is added to resist the relative rotation between particles, and it is calculated in the incremental form as follows: where ∆θ b is the relative bend-rotation increment, and k r is the rolling resistance stiffness. There are many different definitions of k r , and in this study k r is defined as: and R is the contact effective radius, which is defined as: where R (1) and R (2) are the radii of the contact ends. If one of the ends is a wall, then R (2) is infinite. Similar to the sliding model, a rolling model is applied when the magnitude of rolling resistance moment exceeds a threshold limit, which is calculated as: where µ r is the rolling friction coefficient.

Model Generation and Loading Paths
To simulate the behavior of sand-rubber mixtures in generalized stress paths, a numerical true triaxial shear test sample is established. The numerical sample is enclosed by six rigid, frictionless walls with a size of 80 mm × 39.1 mm × 39.1 mm. The top and bottom walls are used as loading platens. The particle size distribution (PSD) of the simulated sand-rubber mixtures is similar to that of the materials reported by Deng et al. [39]. Deng et al. [39] and Gong et al. [31] can be referred to for a brief description of the sand-rubber mixtures and the background of the laboratory testing. The particle size distribution curves of sand and rubber are shown in Figure 1. The PSD of the simulated rubber matches that of Deng et al. [39] following a uniform distribution with minimum and maximum sizes of 4 mm and 5 mm. For numerical assemblies containing small, stiff particles, the computation times could be extremely long if the actual sizes were used. To increase the computation efficiency and reduce the calculation time, model scaling is necessary. In general, there are three primary methods of model scaling [40]: density scaling [41,42], gravity scaling [43,44], and mass scaling [45,46]. In the density scaling method, particle density increases by several orders of magnitude, while the particle sizes remain the same. The mass scaling method is just the opposite of the density scaling method. In the mass scaling method, particle sizes increases, while the particle density remain constant. The gravity scaling method increases gravity field, and keep the particle size and particle density constant. The mass scaling method has been proved to be suitable for sand-rubber mixtures [47], and has been successfully introduced into the DEM simulation of sand-rubber mixtures under different loading conditions [47,48]. In the current research, the mass scaling [40] was applied to scaling the numerical sand particles for the sake of the computation efficiency. There are two phases of the preparation of the particle assembly. In the first phase, the sand particles and rubber particles are initially generated with random locations in the box enclosed by six walls following the size distributions shown in Figure 1. Note that there may be large particle-particle overlaps in the particle assembly. To reduce the overlaps to an acceptable level, a very low value of inter-particle friction coefficient is set to the particles at first. The numerical particle assembly is solved to allow the particles to rearrange until an equilibrium and isotropic state is achieved.
In the second phase, the isotropically compressed sample is prepared by compressing the generated sample with a stress servo-control algorithm. The servo-controlled algorithm in PFC3D is called in each computational cycle to determine the current wall stresses and to adjust the wall velocities so as to reduce the difference between the current stress and the required stress. Sand particles are assigned an inter-particle friction coefficient ( ) of 0.01 while rubber particles are assigned a of 1.5 during the isotropic compression stage. The numerical samples with different rubber contents at the end of an isotropic compression of 100 kPa are shown in Figure 2.   There are two phases of the preparation of the particle assembly. In the first phase, the sand particles and rubber particles are initially generated with random locations in the box enclosed by six walls following the size distributions shown in Figure 1. Note that there may be large particle-particle overlaps in the particle assembly. To reduce the overlaps to an acceptable level, a very low value of inter-particle friction coefficient is set to the particles at first. The numerical particle assembly is solved to allow the particles to rearrange until an equilibrium and isotropic state is achieved.
In the second phase, the isotropically compressed sample is prepared by compressing the generated sample with a stress servo-control algorithm. The servo-controlled algorithm in PFC3D is called in each computational cycle to determine the current wall stresses and to adjust the wall velocities so as to reduce the difference between the current stress and the required stress. Sand particles are assigned an inter-particle friction coefficient (µ) of 0.01 while rubber particles are assigned a µ of 1.5 during the isotropic compression stage. The numerical samples with different rubber contents at the end of an isotropic compression of 100 kPa are shown in Figure 2. There are two phases of the preparation of the particle assembly. In the first phase, the sand particles and rubber particles are initially generated with random locations in the box enclosed by six walls following the size distributions shown in Figure 1. Note that there may be large particle-particle overlaps in the particle assembly. To reduce the overlaps to an acceptable level, a very low value of inter-particle friction coefficient is set to the particles at first. The numerical particle assembly is solved to allow the particles to rearrange until an equilibrium and isotropic state is achieved.
In the second phase, the isotropically compressed sample is prepared by compressing the generated sample with a stress servo-control algorithm. The servo-controlled algorithm in PFC3D is called in each computational cycle to determine the current wall stresses and to adjust the wall velocities so as to reduce the difference between the current stress and the required stress. Sand particles are assigned an inter-particle friction coefficient ( ) of 0.01 while rubber particles are assigned a of 1.5 during the isotropic compression stage. The numerical samples with different rubber contents at the end of an isotropic compression of 100 kPa are shown in Figure 2.   Table 1 summarizes the simulations carried out in this study. Samples are labelled as CS-XXX-b for Clean Sand, with XXX indicating the confining pressure σ 3 and with b indicating the intermediate principle stress ratio. Samples are labelled as RS-RC-XXX-b for sand-rubber mixtures, with RC indicating the content of rubber in percentage, and similarly, with XXX indicating the confining pressure σ 3 and with b indicating the intermediate principle stress ratio. The numerical tests are conducted on dense sand-rubber mixtures (D r = 70%) with different RC of 0%, 10%, and 30%. RC is the percentage of rubber particles by weight and is defined as RC = (m rubber /m total ) × 100%, where m rubber is the mass of rubber particles and m total is the total mass of sand-rubber sample. The number of particles for specimens with different rubber content is also listed in Table 1. Once the isotropic compression state is achieved, the numerical specimens are sheared under true triaxial shear conditions. In the true triaxial stress state, the major principal stress σ 1 , the intermediate principal stress σ 2 , and the minor principal stress σ 3 can be respectively manipulated to simulate the complex loading paths.
In this paper, the true triaxial shear tests under the constant σ 3 and constant intermediate principal stress ratio b westress state are presented. For the samples under the constant σ 3 and constant b stress state, σ 1 is coincident with the moving directions of the top and the bottom walls, and the value of σ 1 varies during the test. The intermediate principal stress σ 2 can be estimated from: The stress servo-control method is adopted to keep the minor principal stress and the intermediate principal stress ratio constant. To explore the effect of b, six drained true triaxial tests with different b (b = 0.0, 0.2, 0.4, 0.6, 0.8 and 1.0) are conducted on particle assemblies. The confining pressure levels simulated in this work are 100 and 200 kPa.

Calibration Procedure
The mechanical properties of DEM samples are controlled by the micro-parameters assigned to the sand particles and the rubber particles. To get the micro-parameters, a calibration procedure is given in this study. Since the rolling resistance model is applied to simulating the behavior of contacting particles, the micro-parameters calibrated in this study are the normal contact stiffness, the shear contact stiffness, friction coefficient, and the rolling friction coefficient. The micro-parameters of the numerical sand-rubber mixtures model, while involving complex interactions between two different materials, i.e., sand to sand, rubber to rubber or sand to rubber, are calibrated by using a strategy as follows. The microscopic parameters of sand-sand contact in pure sand samples are first calibrated by iteratively changing microscopic parameters until the numerical results match the stress-strain curves measured in the laboratory. Then the micro-parameters of sand-rubber contacts are calibrated in samples with rubber content of 10%. Note that this is possible only because the number of rubber-rubber contacts is very small in samples with rubber content of 10%, and rubber-rubber contacts have little influence on the mechanical behavior of samples with 10% rubber particles. The micro-parameters of rubber-rubber contacts are set as the same as those of sand-rubber contacts in this step. The micro-parameters of rubber-rubber contacts are finally calibrated Materials 2020, 13, 5716 6 of 24 in samples with rubber content of 30% based on the calibrated micro-parameters of sand-sand contacts and sand-rubber contacts. After that, the microscopic parameters are adjusted to be more acceptable. The micro-parameters calibrated in this study are list in Table 2. The density of rubber particle is 1330 kg/m 3 , and the density of sand particles is 2620 kg/m 3 .  Figure 3 presents the comparisons between the experimental data and the corresponding DEM simulation results. It is observed that the stress-strain behavior of the simulation results matches well with the experimental data, suggesting that the calibrated micro-parameters in Table 2 Figure 3 presents the comparisons between the experimental data and the corresponding DEM simulation results. It is observed that the stress-strain behavior of the simulation results matches well with the experimental data, suggesting that the calibrated micro-parameters in Table 2 are suitable.

Deviatoric Stress and Volumetric Strain Against the Axial Strain
A series of DEM simulations of true triaxial tests with different rubber contents (0%, 10% and 30%) are conducted with intermediate principal stress ratio b varying from 0.0 to 1.0. The curves of the macroscopic response of the numerical samples under a constant minimum principal stress σ 3 of 100 kPa are shown in Figure 4. The variation trend of b in each curve is marked with an arrow line.

Deviatoric Stress and Volumetric Strain Against the Axial Strain
A series of DEM simulations of true triaxial tests with different rubber contents (0%, 10% and 30%) are conducted with intermediate principal stress ratio varying from 0.0 to 1.0. The curves of the macroscopic response of the numerical samples under a constant minimum principal stress of 100 kPa are shown in Figure 4. The variation trend of in each curve is marked with an arrow line.
Data in Figure 4a show that the stress ratio / is initially independent from the intermediate principal stress ratio b until the axial strain reaches a specific value. After that point the stress ratio shows a clear dependency on intermediate principal stress ratio b. An increase in b value leads to a decrease in the value of stress ratio at the same axial strain . The same trend is also observed in Figure 4b,c, which indicate the evolution of versus axial strain for sand-rubber samples with rubber contents of 10% and of 30%, respectively. This observation is also in line with previous experimental and numerical results [20,22,26,28]. By comparing the curves of sand and sand-rubber mixtures with different rubber contents, it can be observed that the pure sand samples show the typical response of strain softening, while sand rubber mixture materials show the response of strain hardening. Figure 5 shows the comparison between the peak deviatoric stress in experimental tests and that in DEM simulations on specimens ( = 0) with various rubber contents under confining pressures of 100 kPa and of 200 kPa. It is observed that the experimental results and the DEM simulation results of the relationship between the peak friction angle and rubber contents share the same trend. The peak friction angle first decreases with the rubber content varying from 0% to 10%, then increases with the rubber content varying from 10% to 30%. This observation can also be found in other The evolution of stress ratio η = q/p along the axial strain ε 1 of pure sand samples with different b values is shown in Figure 4a, where q and p are the deviatoric stress and the mean stress, and can be calculated as:

Effect of Rubber Content and Intermediate Principal Stress Ratio on the Stress Ratio in the Peak State
Data in Figure 4a show that the stress ratio q/p is initially independent from the intermediate principal stress ratio b until the axial strain ε 1 reaches a specific value. After that point the stress ratio η shows a clear dependency on intermediate principal stress ratio b. An increase in b value leads to a decrease in the value of stress ratio η at the same axial strain ε 1 . The same trend is also observed in Figure 4b,c, which indicate the evolution of η versus axial strain ε 1 for sand-rubber samples with rubber contents of 10% and of 30%, respectively. This observation is also in line with previous experimental and numerical results [20,22,26,28]. By comparing the curves of sand and sand-rubber mixtures with different rubber contents, it can be observed that the pure sand samples show the typical response of strain softening, while sand rubber mixture materials show the response of strain hardening. Figure 5 shows the comparison between the peak deviatoric stress in experimental tests and that in DEM simulations on specimens (b = 0) with various rubber contents under confining pressures of 100 kPa and of 200 kPa. It is observed that the experimental results and the DEM simulation results of the relationship between the peak friction angle and rubber contents share the same trend. The peak friction angle first decreases with the rubber content varying from 0% to 10%, then increases with the rubber content varying from 10% to 30%. This observation can also be found in other experimental results of sand-rubber mixtures containing large rubber particles [49]. The influences of rubber content on the relationship between peak friction angle and intermediate principal stress ratio b are also explored and demonstrated in Figure 6a. The peak friction angle is the maximum value of the mobilized friction angle during shearing. The mobilized friction angle can be calculated by:

Effect of Rubber Content and Intermediate Principal Stress Ratio on the Stress Ratio in the Peak State
As shown in Figure 6a, the peak friction angle of each specimen first increases with intermediate principal stress ratio b until a special value is reached, and then decreases with b. The samples of = 1.0 have a greater value of peak friction angle than those of = 0.0. This trend is in line with the previous experimental and numerical results available, suggesting that the numerical true triaxial shearing test conducted on sand-rubber mixtures in this study is reasonable. It is noted that the special value of b for pure sand is = 0.6, while for sand-rubber mixtures the special value of b is = 0.8, and the peak friction angles for sand-rubber samples with = 0.2, 0.4, 0.6, and 0.8 nearly fall on a straight line, which are quite different from those for pure sand samples. This means that the failure pattern can be changed by the addition of rubber particles, as confirmed in Figure 6b where the responses of sand-rubber mixtures are quite different from those of pure sand when normalized by for = 0.0.
(a) (b) Figure 6. Influence of rubber content on (a) the peak friction angle and (b) the normalized peak friction angle in true triaxial compression. Normallized peak friction angle, The influences of rubber content on the relationship between peak friction angle φ p and intermediate principal stress ratio b are also explored and demonstrated in Figure 6a. The peak friction angle φ p is the maximum value of the mobilized friction angle φ m during shearing. The mobilized friction angle can be calculated by: The influences of rubber content on the relationship between peak friction angle and intermediate principal stress ratio b are also explored and demonstrated in Figure 6a. The peak friction angle is the maximum value of the mobilized friction angle during shearing. The mobilized friction angle can be calculated by: As shown in Figure 6a, the peak friction angle of each specimen first increases with intermediate principal stress ratio b until a special value is reached, and then decreases with b. The samples of = 1.0 have a greater value of peak friction angle than those of = 0.0. This trend is in line with the previous experimental and numerical results available, suggesting that the numerical true triaxial shearing test conducted on sand-rubber mixtures in this study is reasonable. It is noted that the special value of b for pure sand is = 0.6, while for sand-rubber mixtures the special value of b is = 0.8, and the peak friction angles for sand-rubber samples with = 0.2, 0.4, 0.6, and 0.8 nearly fall on a straight line, which are quite different from those for pure sand samples. This means that the failure pattern can be changed by the addition of rubber particles, as confirmed in Figure 6b where the responses of sand-rubber mixtures are quite different from those of pure sand when normalized by for = 0.0.
(a) (b) Figure 6. Influence of rubber content on (a) the peak friction angle and (b) the normalized peak friction angle in true triaxial compression. Normallized peak friction angle, As shown in Figure 6a, the peak friction angle of each specimen first increases with intermediate principal stress ratio b until a special value is reached, and then decreases with b. The samples of b = 1.0 have a greater value of peak friction angle than those of b = 0.0. This trend is in line with the previous experimental and numerical results available, suggesting that the numerical true triaxial shearing test conducted on sand-rubber mixtures in this study is reasonable. It is noted that the special value of b for pure sand is b = 0.6, while for sand-rubber mixtures the special value of b is b = 0.8, and the peak friction angles for sand-rubber samples with b = 0.2, 0.4, 0.6, and 0.8 nearly fall on a straight line, which are quite different from those for pure sand samples. This means that the failure pattern can be changed by the addition of rubber particles, as confirmed in Figure 6b where the responses of sand-rubber mixtures are quite different from those of pure sand when normalized by φ p for b = 0.0.

Coordination Number
An important parameter widely used to evaluate the internal contact response is the coordinate number. The coordinate number is the average number of the contacts of each particle. Various authors have noted that in granular materials there are particles with no contact or only one contact and they do not contribute to the force transmission. Thornton [41] defines the mechanical coordination number to neglect these particles as: where N c is the number of total contacts, N b is the number of total particles, while N 1 and N 0 are the number of particles with one and with zero contacts, respectively. The evolution of the mechanical coordinate number Z m against axial strain ε 1 for sand-rubber specimens with different b values under constant confining pressure 100 kPa is plotted in Figure 7a.

Coordination Number
An important parameter widely used to evaluate the internal contact response is the coordinate number. The coordinate number is the average number of the contacts of each particle. Various authors have noted that in granular materials there are particles with no contact or only one contact and they do not contribute to the force transmission. Thornton [41] defines the mechanical coordination number to neglect these particles as: where is the number of total contacts, is the number of total particles, while and are the number of particles with one and with zero contacts, respectively. The evolution of the mechanical coordinate number against axial strain for sand-rubber specimens with different b values under constant confining pressure 100 kPa is plotted in Figure 7a. As shown in Figure 7a Figure 4. It is also observed that of samples with 30% rubber particles is higher than of those with 10% rubber particles. This may be attributed to the fact that with the increase of rubber content, the total number of particles decreases, while the contact number of each particle increases. For binary mixture materials with two types of particles, Minh and Cheng [50] has defined three types of coordinate number relating to each type of contact. Following Minh and Cheng [50], the coordinate numbers of rubber-sand mixture for sand-sand contacts , rubber-sand contacts and rubber-rubber contacts are defined as:  Figure 7b illustrates the evolution of Z m for specimens with different rubber contents under triaxial compression and extension tests. The pure sand samples demonstrate a typical behavior of a dense material, while the sand-rubber samples show a behavior of a loose one. This trend is in agreement with the results shown in Figure 4. It is also observed that Z m of samples with 30% rubber particles is higher than Z m of those with 10% rubber particles. This may be attributed to the fact that with the increase of rubber content, the total number of particles decreases, while the contact number of each particle increases.
For binary mixture materials with two types of particles, Minh and Cheng [50] has defined three types of coordinate number relating to each type of contact. Following Minh and Cheng [50], the coordinate numbers of rubber-sand mixture for sand-sand contacts Z s−s , rubber-sand contacts Z r−s and rubber-rubber contacts Z r−r are defined as: where N c,s−s , N c,r−s , and N c,r−r are respectively the number of sand-sand, rubber-sand, and rubber-rubber contacts, while N b,s and N b,r are the number of sand and rubber particles. Figure 8 shows the evolution of Z s−s , Z r−s , and Z r−r for samples with rubber content of 10% in different loading paths with constant confining pressure of 100 kPa. Each contact type of coordinate number has the same initial value. The initial value of Z s−s is the largest, followed by Z r−r , and then Z r−s . Z s−s and Z r−s are initially independent from intermediate principal stress ratio b and quickly divided from each other at a very small axial strain. After that Z s−s and Z r−s generally increase as b increases at the same axial strain. The values of Z r−r are very volatile because of the limited number of rubber particles and rubber-rubber contacts. It is obvious that because of the limitation of the number of rubber particles, the micro structure and the mechanical behavior are dominated by sand-sand particles when rubber content is 10%. Figure 9 shows the evolutions of Z s−s , Z r−s , and Z r−r for samples with rubber content of 10% and 30% under triaxial compression and tension tests. As shown in Figure 9, Z s−s decreases when rubber content increases from 10% to 30%, while Z r−s and Z r−r increase with rubber content. The values of Z s−s and Z r−s increase with b for samples with rubber content of 30%, which agrees with the trend in samples with rubber content of 10%.

Proportion of Strong Contact in Different Types of Contacts
Gong et al. [31] has investigated the role of different contact type in a strong force network and has concluded that rubber-rubber contacts share a very low proportion in the strong network, and contribute less to stress transmission within the system. In this study, we investigate the role of different contact types in a different way. The strong contact ratio is proposed and defined as the ratio of the number of strong contacts to the number of all contacts (sand-sand, sand-rubber, rubber-rubber, and overall contacts). According to Radjai et al. [51] as well as Shi and Guo [52], the contacts can be divided into strong contacts and weak contacts by the average normal contact force, which is the average value of the normal contact force of all contacts.
The evolutions of strong contact ratios in different types of contacts with rubber content of 10% and of 30% are shown in Figure 10. For sand-sand contact, the strong contact ratios nearly equal to those of overall contacts at the same axial strain for all samples. The strong contact ratios of rubber-rubber contact present a different performance. They are much higher than those of overall contacts, and their initial values are nearly 1.0, which means nearly all rubber-rubber contacts are strong contacts. When the shearing begins, the strong contact ratios of rubber-rubber contacts present fluctuate declining with axial strain ε 1 , but their values keep no less than 0.8 for all the samples. This means rubber-rubber contacts mainly participate in strong force networks, and even with rubber content of 10%, the number of rubber particles is very small, rubber particles can significantly affect the micro structure and force transmission in the contact network. This finding also can be used to explain why the failure behavior is changed by adding rubber particles. A suggested explanation of this phenomenon is that the added large rubber particles mainly affect the inherent stability of the strong network. , and (c) for sand-rubber mixtures with 10% rubber particles under different loading paths Figure 9 shows the evolutions of , , and for samples with rubber content of 10% and 30% under triaxial compression and tension tests. As shown in Figure 9, decreases when rubber content increases from 10% to 30%, while and increase with rubber content. The  values of and increase with b for samples with rubber content of 30%, which agrees with the trend in samples with rubber content of 10%. Zs-s_RC30%, b0.0 Zr-r_RC30%,b0.0 Zr-s_RC30%, b0.0 Zs-s_RC30%, b1.0 Zr-r_RC30%, b1.0 Zr-s_RC30%, b1.0 Zs-s_RC10%, b0.0 Zr-r_RC10%, b0.0 Zr-s_RC10%, b0.0 Zs-s_RC10%, b1.0 Zr-r_RC10%, b1.0 Zr-s_RC10%, b1.0

Proportion of Strong Contact in Different Types of Contacts
Gong et al. [31] has investigated the role of different contact type in a strong force network and has concluded that rubber-rubber contacts share a very low proportion in the strong network, and contribute less to stress transmission within the system. In this study, we investigate the role of different contact types in a different way. The strong contact ratio is proposed and defined as the ratio of the number of strong contacts to the number of all contacts (sand-sand, sand-rubber, rubberrubber, and overall contacts). According to Radjai et al. [51] as well as Shi and Guo [52], the contacts can be divided into strong contacts and weak contacts by the average normal contact force, which is the average value of the normal contact force of all contacts.
The evolutions of strong contact ratios in different types of contacts with rubber content of 10% and of 30% are shown in Figure 10. For sand-sand contact, the strong contact ratios nearly equal to those of overall contacts at the same axial strain for all samples. The strong contact ratios of rubberrubber contact present a different performance. They are much higher than those of overall contacts, and their initial values are nearly 1.0, which means nearly all rubber-rubber contacts are strong contacts. When the shearing begins, the strong contact ratios of rubber-rubber contacts present fluctuate declining with axial strain , but their values keep no less than 0.8 for all the samples. This means rubber-rubber contacts mainly participate in strong force networks, and even with rubber content of 10%, the number of rubber particles is very small, rubber particles can significantly affect the micro structure and force transmission in the contact network. This finding also can be used to explain why the failure behavior is changed by adding rubber particles. A suggested explanation of this phenomenon is that the added large rubber particles mainly affect the inherent stability of the strong network. For the case of sand-rubber contacts, the strong contact ratios are greater than those of overall contacts for samples with 10% rubber particles, while equal to those of overall contacts for samples with 30% rubber particles. This implies that force chains consisting only of rubber particles may exist in the sand-rubber specimens when rubber content exceeds 30%. It should be noted that the coordinate number of rubber-sand contacts increase with as shown in Figures 8b and 9 in the For the case of sand-rubber contacts, the strong contact ratios are greater than those of overall contacts for samples with 10% rubber particles, while equal to those of overall contacts for samples with 30% rubber particles. This implies that force chains consisting only of rubber particles may exist in the sand-rubber specimens when rubber content exceeds 30%. It should be noted that the coordinate number of rubber-sand contacts increase with b as shown in Figures 8b and 9 in the previous section. Considering the effect of intermediate principal stress ratio on the strong contact ratios and coordinate numbers of rubber-sand contacts together, we can find that with an increase of b, the force transmitted through the rubber particles increases, and more rubber-sand contacts are needed to laterally support the force chain.

Fabric Tensor and Anisotropy
Stress-induced anisotropy is one of the most important properties associated with the key features of granular materials. The concept of fabric, which was first proposed by Oda [53], is a useful statistical method to describe the geometrical anisotropy induced by applied stress in granular materials. The fabric anisotropy can be quantified by a fabric tensor based on the contact normal distribution. The fabric tensor of the whole contact networks is defined by Satake [54] as: where n k i is the unit contact normal component in the ith direction. Similar to the fabric tensor, the fabric tensor of strong contacts and weak contacts can be defined as: where N s C and N w C are respectively the numbers of strong and weak contacts. φ s 1 , φ s 2 , and φ s 3 are the principal values of the fabric tensor of strong contacts associated with the directions of maximum, intermediate, and minimum principal stresses, respectively. Figure 11 shows the evolutions of φ s 1 , φ s 2 , and φ s 3 of sand-rubber samples with rubber content of 10% under constant confining stress of 100 kPa. φ s 1 first increases with axial strain and then nearly keeps constant, and an increase in the value of b leads to a decrease in the φ s 1 . The value of φ s 2 first decreases with the axial strain until b reaches 0.4, and after that the value of φ s 2 increases with the axial strain (Figure 11b). It is also seen that the value of φ s 2 rises as the intermediate principal stress ratio b increases. The evolution of φ s 3 has the reverse trend compared with that of φ s 1 : its value first decreases with the axial strain and then remains nearly constant. On the contrary, the intermediate principal stress ratio b has the same effect on the evolutions of φ s 3 and φ s 1 , and an increase in the value of b leads to a decrease in the φ s 3 . To clarify the effect of rubber content on the principal values of fabric tensor of strong contacts, Figure 12 presents the evolutions of φ s 1 , φ s 2 , and φ s 3 for samples with 0%, 10%, and 30% rubber particles under triaxial compression and tension conditions. As shown in Figure 12a, the value of φ s 1 of pure sand increases with the axial strain until a peak value is reached, and then slightly decreases to the end of shearing, while the behavior of sand-rubber mixtures has the same trend mentioned before. It is noted that the value of φ s 1 is first dependent on rubber content until the axial strain reaches a certain value, and after that the φ s 1 is nearly independent from the rubber content. Similar trends are observed in the evolutions of φ s 2 and φ s 3 , as shown in Figure 12b,c.
value of b leads to a decrease in the . The value of first decreases with the axial strain until b reaches 0.4, and after that the value of increases with the axial strain (Figure 11b). It is also seen that the value of rises as the intermediate principal stress ratio b increases. The evolution of has the reverse trend compared with that of : its value first decreases with the axial strain and then remains nearly constant. On the contrary, the intermediate principal stress ratio b has the same effect on the evolutions of and , and an increase in the value of b leads to a decrease in the . To clarify the effect of rubber content on the principal values of fabric tensor of strong contacts, Figure 12 presents the evolutions of , , and for samples with 0%, 10%, and 30% rubber particles under triaxial compression and tension conditions. As shown in Figure 12a, the value of of pure sand increases with the axial strain until a peak value is reached, and then slightly decreases to the end of shearing, while the behavior of sand-rubber mixtures has the same trend mentioned before. It is noted that the value of is first dependent on rubber content until the axial strain reaches a certain value, and after that the is nearly independent from the rubber content. Similar trends are observed in the evolutions of and , as shown in Figure 12b,c. particles under triaxial compression and tension conditions. As shown in Figure 12a, the value of of pure sand increases with the axial strain until a peak value is reached, and then slightly decreases to the end of shearing, while the behavior of sand-rubber mixtures has the same trend mentioned before. It is noted that the value of is first dependent on rubber content until the axial strain reaches a certain value, and after that the is nearly independent from the rubber content. Similar trends are observed in the evolutions of and , as shown in Figure 12b,c. The evolution of the strong deviatoric fabric for the samples with 10% rubber particles under different stress conditions are presented in Figure 13a. The deviatoric fabric is initially independent on the intermediate principal stress ratio b until the axial strain reaches nearly 5%. After that the deviatoric fabric becomes dependent on the b values, and an increase in b value leads to a decrease in the value of deviatoric fabric. The effects of rubber content on the evolutions of the deviatoric fabric of strong contacts is shown in Figure 13b. It can be found that a higher rubber content leads to a lower value of deviatoric fabric. This means the samples with a higher rubber content shows lower level of fabric anisotropy.
The evolution of the strong deviatoric fabric for the samples with 10% rubber particles under different stress conditions are presented in Figure 13a. The deviatoric fabric is initially independent on the intermediate principal stress ratio b until the axial strain reaches nearly 5%. After that the deviatoric fabric becomes dependent on the b values, and an increase in b value leads to a decrease in the value of deviatoric fabric. The effects of rubber content on the evolutions of the deviatoric fabric of strong contacts is shown in Figure13b. It can be found that a higher rubber content leads to a lower value of deviatoric fabric. This means the samples with a higher rubber content shows lower level of fabric anisotropy.    It can be observed that the peak deviatoric fabric decreases with intermediate principal stress ratio b for all the samples, and the peak deviatoric fabric of pure sand is higher than those ones of sandrubber samples with the same b values. The samples with rubber content of 10% and of 30% nearly have the same peak value of deviatoric fabric.

Normal Contact Force and the Probability Density Function of Normal Contact Force
The probability density function (PDF) of contact force can be used as a statistical tool to explore the characteristic of the force transmission behavior in granular materials. The contact forces are partitioned into strong contact forces and weak contact forces by the average contact force. The probability of the strong forces is presented as an exponential equation, while the probability of the weak forces is described by a power equation, as follows:

Normal Contact Force and the Probability Density Function of Normal Contact Force
The probability density function (PDF) of contact force can be used as a statistical tool to explore the characteristic of the force transmission behavior in granular materials. The contact forces are partitioned into strong contact forces and weak contact forces by the average contact force. The probability of the strong forces is presented as an exponential equation, while the probability of the weak forces is described by a power equation, as follows: where f x is the normal contact force or shear contact force; f x is the average value of f x ; A and B are the parameters that imply the inhomogeneity of forces in discrete element model. Figure 15 shows the probability density function (PDF) of the normal contact forces at the shear strain of 15% for pure sand and sand-rubber mixtures with different intermediate principal stress ratio in a log-linear space. It is clear that the normal contact forces show a high agreement with the probability density distribution. It can also been found, the probability density decreases with the contact forces, and the magnitude of the contact forces increases with b values.
The probability density distributions of pure sand and sand-rubber mixtures under the true triaxial loading paths of b = 0 and b = 1.0 are plotted in Figure 16. As shown in Figure 16, the rubber content has a significant effect on the probability density distributions of sand-rubber mixtures. A wider distribution of contact force in the sample is shown as the rubber content increases, and the discrepancies in the probability density distribution of normal contact force associated with b = 1.0 are larger than the ones associated with b = 0.0.
The average normal contact force is defined as the average value of the normal contact force of all contacts. The average normal contact forces of samples with different rubber contents under various loading paths at the axial strain of 15% are shown in Figure 17. It can be seen that for the samples with the same rubber content, the average normal contact forces increase with b values, and for the same b value, the average normal contact forces increase with rubber content.
where is the normal contact force or shear contact force; 〈| |〉 is the average value of | |; A and B are the parameters that imply the inhomogeneity of forces in discrete element model. Figure 15 shows the probability density function (PDF) of the normal contact forces at the shear strain of 15% for pure sand and sand-rubber mixtures with different intermediate principal stress ratio in a log-linear space. It is clear that the normal contact forces show a high agreement with the probability density distribution. It can also been found, the probability density decreases with the contact forces, and the magnitude of the contact forces increases with b values. The probability density distributions of pure sand and sand-rubber mixtures under the true triaxial loading paths of = 0 and = 1.0 are plotted in Figure 16. As shown in Figure 16, the rubber content has a significant effect on the probability density distributions of sand-rubber mixtures. A wider distribution of contact force in the sample is shown as the rubber content increases, and the discrepancies in the probability density distribution of normal contact force associated with b =1.0 are larger than the ones associated with = 0.0. The probability density distributions of pure sand and sand-rubber mixtures under the true triaxial loading paths of = 0 and = 1.0 are plotted in Figure 16. As shown in Figure 16, the rubber content has a significant effect on the probability density distributions of sand-rubber mixtures. A wider distribution of contact force in the sample is shown as the rubber content increases, and the discrepancies in the probability density distribution of normal contact force associated with b =1.0 are larger than the ones associated with = 0.0. The average normal contact force is defined as the average value of the normal contact force of all contacts. The average normal contact forces of samples with different rubber contents under various loading paths at the axial strain of 15% are shown in Figure 17. It can be seen that for the samples with the same rubber content, the average normal contact forces increase with b values, and for the same b value, the average normal contact forces increase with rubber content.

Conclusions
Several DEM simulated true triaxial shear tests have been conducted to investigate the effect of intermediate principal stress ratio on the mechanical behavior of sand-rubber mixtures. Both the macroscale and microscale performance has been discussed. The main conclusions drawn from the numerical study are the following: 1. The peak strengths of samples under conventional triaxial tests first decrease with 10% rubber particles added, and then increase when the proportion of rubber particles rises up to 30%, but the peak strengths sand-rubber mixtures with either 10% or 30% rubber particles are lower than those of pure sand. This trend is in agreement with previous experimental results and numerical simulations on sand-rubber mixtures with large rubber particles, which confirms the feasibility

Conclusions
Several DEM simulated true triaxial shear tests have been conducted to investigate the effect of intermediate principal stress ratio on the mechanical behavior of sand-rubber mixtures. Both the macroscale and microscale performance has been discussed. The main conclusions drawn from the numerical study are the following: 1.
The peak strengths of samples under conventional triaxial tests first decrease with 10% rubber particles added, and then increase when the proportion of rubber particles rises up to 30%, but the peak strengths sand-rubber mixtures with either 10% or 30% rubber particles are lower than those of pure sand. This trend is in agreement with previous experimental results and numerical simulations on sand-rubber mixtures with large rubber particles, which confirms the feasibility of the simulations conducted in this study. The same trend can also be observed in the peak friction angles for the samples at each intermediate principal stress ratio.

2.
For sand-rubber mixtures, the relationship between the peak friction angle and the intermediate principal stress ratio is quite different from that for pure sand, which means adding rubber particles can change the failure behavior of sand under complex loading conditions. A suggested explanation of this phenomenon is that the added large rubber particles mainly affect the inherent stability of the strong network. This study can provide a reference for the constitutive model development of sand-rubber mixtures. 3.
The investigation on the strong contact ratio of different types of contacts show that nearly all of the rubber-rubber contacts of sand-rubber mixtures are strong contacts, no matter what the rubber contents and the values of intermediate principal stress ratio are. While the strong contact ratio of rubber-sand contacts is higher than that of overall contacts for specimens with 10% of rubber particles, and becomes nearly equal to that of the overall contacts when rubber content rises up to 30%. It can be concluded that rubber content can significantly influence the micro structure and the force transmission in the contact network. This finding also confirms the explanation in the Conclusion 2.

4.
For samples with the same rubber content, the strong contact ratio of rubber-sand contacts decrease with the principal stress ratio b, while the coordinate number of rubber-sand contacts increase with b. It means that with an increase of b, the force transmitted through the rubber particles increases, and more rubber-sand contacts are needed to support the force chain. 5.
The analysis of the fabric anisotropy shows that the deviatoric fabric of strong contacts demonstrates a decline by adding large rubber particles, and the deviatoric fabric of strong contacts also decreases with b, which is in line with the previous numerical simulations.

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

RSM Sand-rubber mixtures DEM
Discrete element method PSD Particle size distribution RC Rubber content PDF Probability density function k n The normal stiffness k s The shear stiffness µ The friction coefficient M r The rolling resistance moment k r The rolling resistance stiffness µ r The rolling friction coefficient σ 1 The major principal stress σ 2 The intermediate principal stress σ 3 The minor principal stress The major principle strain q The deviatoric stress p The mean stress b The Intermediate principal stress ratio φ p The peak friction angle φ m The mobilized friction angle Z m The mechanical coordination number Z s−s The coordination number of sand-sand contacts Z r−s The coordination number of rubber-sand contacts Z r−r The coordination numbers of rubber-rubber contacts φ s