# 3D Mechanical Characters and Their Fabric Evolutions of Granular Materials by DEM Simulation

Authors: Xuefeng Li; Weinan Lu; Yuqi He
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

To quantitatively describe the influence of the contact characteristics of granular materials on their mechanical response, the true-triaxial tests with different particle shapes are simulated by the discrete element method (DEM), and the connection between the evolutions of particle contact fabric and the anisotropic mechanical responses is studied. The contact normal vector of the particle in 3D space is described by two independent angles, by which the contact fabric tensor is defined. The amplitude parameters in three orthogonal directions are defined by the invariants of the plane fabric tensor, which are scalars and describe the degree of anisotropy of the contact fabric in each plane. The expression of orthotropic fabric tensor is derived from the amplitude parameters, with the change of geometric space of contact normal vector, which is reduced to the different tensor of transverse isotropic naturally. The fabric tensor can be directly applied to the constitutive equation to describe the effects of the particle contact on the mechanical response. For verifying the rationality of contact characteristics described by fabric tensor, four particle shapes are clumped by PFC3D. The mechanical properties of specimens with different particles are simulated under the true-triaxial loading path, and the data of contact normal vector is extracted in real time. The simulation results showed that the particle shapes have a significant effect on the 3D stress-strain relationship and strength, which showed apparent anisotropy, and the invariants of fabric tensor can be used to describe the evolution of particle contact in the loading process.


---

## Body

## 1. Introduction

Due to the unique properties of granular materials, such as dispersion, for a long time, their macroscopic mechanical properties and microscopic mechanical properties are often studied separately. Taking geotechnical engineering as an example, the strength and stability of granular materials in slope, foundation, and underground engineering studies from a macroscopic perspective can often meet engineering needs. However, under complex stress conditions such as cyclic loading, the failure’s forms are quite different macroscopically, due to anisotropy, strain localization, strain softening, etc. The fundamental reason for these phenomena is the microscopic properties of the material, such as particle morphology, motion, fabric, and evolution [1–4]. Therefore, it is of theoretical and practical significance to study the relationship between macroscopic and microscopic mechanical properties of granular materials [5–9].Fabric tensor is a great connection tool for granular materials, and quantitative description is helpful in revealing the influence of microscopic fabric on macroscopic mechanical properties [10]. Since Oda et al. [11–13] put forward the definition of fabric tensors, the study on the mechanical properties of granular materials has gone from macroscopic to microscopic, and several constitutive models and strength criteria based on fabric tensors have been established. Guo [14] and Yao and Kong [15] established the strength criterion of granular materials using fabric tensors. Ouadfel and Rothenburg [16] researched the relationship of the stress-force-fabric. Liu et al. [17] presented the relationship between macroscopic stress, microscopic contact force, and anisotropic fabric. Petalas et al. [18] introduced fabric anisotropy parameters into the constitutive model and studied the influence of fabric anisotropy evolution on mechanical properties. Pouragha et al. [19] combined the contact force, contact direction, and strength criterion between particles. Hu et al. [20] presented a boundary surface model containing fabric tensors, in which the anisotropy was represented by introducing back stress opposite to the contact normal direction and contrasted the predicted results of the model between simulation results by DEM. The author of this paper also established the potential theory [21] and the strength criterion [22] using the fabric tensor [23].The above research results enrich the description performance and scope of the existing constitutive models and strength criteria based on fabric tensors. However, the characteristics of macroscopic and microscopic connections and their evolutions under complex stress conditions need to be further studied, such as the relationship between the evolution of contact fabric and macroscopic stress-strain change [24], and the relationship with the evolution of stress-induced anisotropy [25]. In the loading process, the evolution of the particle contact normal determines the evolution direction of macroscopic deformation to a certain extent. Liu et al. [26] studied the influence of stress path on fabric evolution by analyzing the evolution of partial fabric. Zheng et al. [27] studied the relationship between strength and fabric evolution. Hu et al. [28] studied the relationship between strain rate and fabric evolution. Nie et al. [29] believed that particle shape was also one of the factors affecting the evolution of fabric. Ng [30] and Oda et al. [31, 32] found that the contact normal direction between particles tended to be consistent with the maximum principal stress direction. Gu et al. [33] analyzed the evolutions of mechanical quantities such as contact number, contact force, and anisotropic parameters between particles. Vijayan et al. [34] studied the evolution of fabric and the average coordination number in the shearing process of granular materials. Wang et al. [35] and Zhou et al. [36] believed that the initial fabric influenced the evolutions of fabric.The true-triaxial test can determine the macroscopic mechanical properties of granular materials in the 3D space and combine them with the measurement of the materials’ microscopic particle information. We can explain the mechanism of the macro-microscopic relations. Under the true-triaxial condition, the granular materials in the stress state of nonhydrostatic pressure lead to stress-induced anisotropy [37–40], which is helpful in studying the loading influences on contact fabric anisotropy [41–44]. In practice, granular materials are composed of many irregular particles and voids. It is difficult to quantify their contact properties, especially during loading. DEM is a widely used discontinuous method [45–47], which can better simulate the macroscopic mechanical response of different particle morphology materials under complex loading conditions and, at the same time, can extract the microscopic characteristics of the contact in real time. Therefore, it is more suitable for studying the mechanical properties of granular materials from the macroscopic and microscopic perspectives [48–50]. Suhr and Six [51] studied the relationship between particle shapes and fabric evolution by DEM. Dorostkar et al. [52], Yuan et al. [53], and Sazzad and Suzuki [54] used DEM to simulate a triaxial test and believed that loading with different ratios in three orthogonal directions had a significant influence on the evolution of anisotropic fabric. He and Jiang [55] conducted a true-triaxial simulation and analyzed the influence of intermediate principal stress coefficient (b) for fabric tensor. The above research has laid an excellent foundation for granular materials’ macro-microscopic relations and evolutions.In this paper, the influence of the four particle shapes on the macroscopic and microscopic properties of granular materials and the evolution trend of anisotropic fabric will be studied under the condition of true three-dimensional. Based on the author’s research on fabric tensor, the contact fabric tensor will be defined by the normal vector of the contact point, and three amplitude parameters of the orthogonal direction can be defined with the invariant of the fabric tensor. Then, the expression of the orthotropic contact fabric with scalar parameters can be derived. The influence of particle shapes on the strength of granular materials will be analyzed by extracting the stress-strain relationship, peak internal friction angle, and other mechanical parameters. The evolution of anisotropic fabric and the relationship between microscopic fabric and macroscopic stress can be explored by extracting contact fabric, contact number of points, contact force, and other statistics of the contact points.


## 2. Contacts Fabric

As one of the three granular material fabrics, the contact fabric describes the statistical characteristics of the normal direction at the contact points. The force inside the material is transmitted by the contact point on the microscopic scale, for the complexity of the spatial distribution of the contact points caused different forces in a different direction, which shows anisotropy on a macroscopic scale. The degree of anisotropy is related to the contact direction between particles.


### 2.1. The Definition of Fabric

To describe the pattern of particle-to-particle contact direction, Oda et al. [11–13] and Tobita and Kuhn [56] have defined the expression of the second-order fabric tensor as(1)Fij=12N∑k=12Nniknjk,i,j=1,2,3,where N is the number of the particle, and nik and njk are the components of the unit contact normal vector on the coordinate axis, respectively.


### 2.2. Definition of Amplitude Parameters of Orthogonal Planes

Two independent angles can represent the components of the contact normal vector on three orthogonal axes. Take a certain contact point as an example, at the contact point, one angle between the contact normal vector and thex1 axis is θ1k, the other between the projection on the horizontal plane and the x3 axis is α(k), and then the direction vector of equation (1) is(2)n=cosθ1k,sinθ1ksinαk,sinθ1kcosαk.It is difficult to analyze the three-dimensional fabric of sand, so it is necessary to use two-dimensional graph analysis and then carry out reasonable three-dimensional. According to equations (1) and (2), the expression on the orthogonal plane is(3)Fij=12N∑k=12Nniknjki,j=1,2.In the interior of granular materials, the particles are in contact with each other, and the directions of the contact normal vector are normally positively distributed. The two angles can also represent the projection components of particles onto three orthogonal planes. For example, on the plane ofx1-x3, the 2D fabric tensor can be expressed as(4)F¯ij=F¯11F¯13F¯31F¯33=12N∑k=12Nsin2θk12N∑k=12Ncosθksinθk12N∑k=12Ncosθksinθk12N∑k=12Ncos2θk,where F¯ij is a second-order plane-symmetric tensor, and θ(k) is the angle between the x1-x3 plane and the x3 axis projection. As the properties of plane-symmetric stress, the plane fabric tensor can also be described by its invariant. The two principal values of the F¯1 and F¯3 can be expressed as(5)F¯1,3=12F¯11+F¯33±14F¯11−F¯332+F¯1321/2=12±a12,where, according to the fabric tensor definition, the trace of Fij is equal to 1, so the direction vector describes a scalar parameter a1in equation (5) as(6)a1=12N∑k=12Ncos2θ1k−sin2θ1kcos2αk2+∑k=12Nsin2θ1kcosαk2,where a1 is the anisotropic amplitude parameter of contact normal fabric, which describes the statistical probability distribution of the normal direction of the contact point. According to the theoretical definition, a1 value range is [0, 1].Similarly, the amplitude parameters ofx1-x2 plane and x2-x3 plane can be obtained as follows:(7)a2=12N∑k=12Ncos2θ1k−sin2θ1ksin2αk2+∑k=12Nsin2θ1ksinαk2,(8)a3=12N∑k=12Nsin2θ1kcos2αk2+∑k=12Nsin2θ1ksin2αk2.As shown in Figure1, a1, a2, and a3 are the anisotropic amplitude parameters on three orthogonal planes, such as Fx1- Fx3, Fx1- Fx2, and Fx2- Fx3. Three parameters can be determined by the normal vector of all particle contact points from equations (6)–(8), and the range of the three parameters is 0∼1, describing the anisotropy degree of materials on each surface.Figure 1
The schematic diagram for determining the amplitude parameters of orthogonal fabrics.According to the fabric definition, the trace of the fabric tensor is equal to 1. There are only two independent variables among the three amplitude parameters on the orthogonal planes. Hence, the fabric tensorFij (i, j = 1, 2, 3) can be derived from any two of three amplitude parameters defined by equations (6)–(8). This paper uses the orthotropic fabric tensor derived from reference [20]. For a detailed derivation process, refer to reference [57].(9)F=1+a1+a2+a1a23+a1+a2−a1a20001+a1−a2−a1a23+a1+a2−a1a20001−a1+a2−a1a23+a1+a2−a1a2,(10)F=1+a1−a3−a1a33−a1−a3−a1a30001−a1+a3−a1a33−a1−a3−a1a30001−a1−a3+a1a33−a1−a3−a1a3.In equation (9), when α(k) = π/4, sin α(k) = cos α(k), then a1 = a2. The amplitude parameters of vertical planes and in the horizontal direction are equal. Then, the materials are shown as transversely isotropic, and then equation (9) is depredated to transversely isotropic fabric:(11)Fij=13−a1′1+a1′0001−a1′0001−a1′,where(12)a1′=12N∑k=12Ncos2θ1k−sin2θ1k22+12∑k=12Nsin2θ1k2,where a1′ is the amplitude parameter after degradation of a1 and a2, and the a1′ value range is [0, 1].In equation (10), an amplitude parameter is obtained according to any one of the two planes that can be used to describe the transverse isotropy of natural soil, and the amplitude parameter value can be obtained by equation (12). In equation (10), if θ1k=π/2 in equation (6) and α(k) = 0 in equation (8), then a1 = a3. Equation (10) has degenerated into the transversely isotropic form.(13)Fij=13+a3′1+a3′0001+a3′0001−a3′,where(14)a3′=12N∑k=12Ncos2θk2+∑k=12Nsin2θk2,where a3′ is the amplitude parameter after degradation of a1 and a3.Equation (14) is the same as the formula proposed by Oda [13, 58] if the variable satisfies θk=π/2−θ1k. Equation (13) is the only expression that quantitatively describes the fabric tensor of particle orientation at present. It has been widely used in both quantitative tests and macro-microscopic mechanical models. The difference is that equation (13) is derived from the 3D definition of fabric, while the expression of Oda [13, 58] is derived directly from the assumption of transverse isotropy.


## 2.1. The Definition of Fabric

To describe the pattern of particle-to-particle contact direction, Oda et al. [11–13] and Tobita and Kuhn [56] have defined the expression of the second-order fabric tensor as(1)Fij=12N∑k=12Nniknjk,i,j=1,2,3,where N is the number of the particle, and nik and njk are the components of the unit contact normal vector on the coordinate axis, respectively.


## 2.2. Definition of Amplitude Parameters of Orthogonal Planes

Two independent angles can represent the components of the contact normal vector on three orthogonal axes. Take a certain contact point as an example, at the contact point, one angle between the contact normal vector and thex1 axis is θ1k, the other between the projection on the horizontal plane and the x3 axis is α(k), and then the direction vector of equation (1) is(2)n=cosθ1k,sinθ1ksinαk,sinθ1kcosαk.It is difficult to analyze the three-dimensional fabric of sand, so it is necessary to use two-dimensional graph analysis and then carry out reasonable three-dimensional. According to equations (1) and (2), the expression on the orthogonal plane is(3)Fij=12N∑k=12Nniknjki,j=1,2.In the interior of granular materials, the particles are in contact with each other, and the directions of the contact normal vector are normally positively distributed. The two angles can also represent the projection components of particles onto three orthogonal planes. For example, on the plane ofx1-x3, the 2D fabric tensor can be expressed as(4)F¯ij=F¯11F¯13F¯31F¯33=12N∑k=12Nsin2θk12N∑k=12Ncosθksinθk12N∑k=12Ncosθksinθk12N∑k=12Ncos2θk,where F¯ij is a second-order plane-symmetric tensor, and θ(k) is the angle between the x1-x3 plane and the x3 axis projection. As the properties of plane-symmetric stress, the plane fabric tensor can also be described by its invariant. The two principal values of the F¯1 and F¯3 can be expressed as(5)F¯1,3=12F¯11+F¯33±14F¯11−F¯332+F¯1321/2=12±a12,where, according to the fabric tensor definition, the trace of Fij is equal to 1, so the direction vector describes a scalar parameter a1in equation (5) as(6)a1=12N∑k=12Ncos2θ1k−sin2θ1kcos2αk2+∑k=12Nsin2θ1kcosαk2,where a1 is the anisotropic amplitude parameter of contact normal fabric, which describes the statistical probability distribution of the normal direction of the contact point. According to the theoretical definition, a1 value range is [0, 1].Similarly, the amplitude parameters ofx1-x2 plane and x2-x3 plane can be obtained as follows:(7)a2=12N∑k=12Ncos2θ1k−sin2θ1ksin2αk2+∑k=12Nsin2θ1ksinαk2,(8)a3=12N∑k=12Nsin2θ1kcos2αk2+∑k=12Nsin2θ1ksin2αk2.As shown in Figure1, a1, a2, and a3 are the anisotropic amplitude parameters on three orthogonal planes, such as Fx1- Fx3, Fx1- Fx2, and Fx2- Fx3. Three parameters can be determined by the normal vector of all particle contact points from equations (6)–(8), and the range of the three parameters is 0∼1, describing the anisotropy degree of materials on each surface.Figure 1
The schematic diagram for determining the amplitude parameters of orthogonal fabrics.According to the fabric definition, the trace of the fabric tensor is equal to 1. There are only two independent variables among the three amplitude parameters on the orthogonal planes. Hence, the fabric tensorFij (i, j = 1, 2, 3) can be derived from any two of three amplitude parameters defined by equations (6)–(8). This paper uses the orthotropic fabric tensor derived from reference [20]. For a detailed derivation process, refer to reference [57].(9)F=1+a1+a2+a1a23+a1+a2−a1a20001+a1−a2−a1a23+a1+a2−a1a20001−a1+a2−a1a23+a1+a2−a1a2,(10)F=1+a1−a3−a1a33−a1−a3−a1a30001−a1+a3−a1a33−a1−a3−a1a30001−a1−a3+a1a33−a1−a3−a1a3.In equation (9), when α(k) = π/4, sin α(k) = cos α(k), then a1 = a2. The amplitude parameters of vertical planes and in the horizontal direction are equal. Then, the materials are shown as transversely isotropic, and then equation (9) is depredated to transversely isotropic fabric:(11)Fij=13−a1′1+a1′0001−a1′0001−a1′,where(12)a1′=12N∑k=12Ncos2θ1k−sin2θ1k22+12∑k=12Nsin2θ1k2,where a1′ is the amplitude parameter after degradation of a1 and a2, and the a1′ value range is [0, 1].In equation (10), an amplitude parameter is obtained according to any one of the two planes that can be used to describe the transverse isotropy of natural soil, and the amplitude parameter value can be obtained by equation (12). In equation (10), if θ1k=π/2 in equation (6) and α(k) = 0 in equation (8), then a1 = a3. Equation (10) has degenerated into the transversely isotropic form.(13)Fij=13+a3′1+a3′0001+a3′0001−a3′,where(14)a3′=12N∑k=12Ncos2θk2+∑k=12Nsin2θk2,where a3′ is the amplitude parameter after degradation of a1 and a3.Equation (14) is the same as the formula proposed by Oda [13, 58] if the variable satisfies θk=π/2−θ1k. Equation (13) is the only expression that quantitatively describes the fabric tensor of particle orientation at present. It has been widely used in both quantitative tests and macro-microscopic mechanical models. The difference is that equation (13) is derived from the 3D definition of fabric, while the expression of Oda [13, 58] is derived directly from the assumption of transverse isotropy.


## 3. True-Triaxial Test Simulations by PFC3D

### 3.1. Generation of Nonspherical Particles

In DEM, spherical particle (3D) with different diameters is used, which is too ideal to describe the natural geometry of granular materials, especially when the microscopic particle characteristics are used to describe the influence of their fabric on macroscopic mechanism. Therefore, the four particle shapes, that is, spheres, elongate clump, pyramid clump, and cube clump (Figure2), are used to establish a true-triaxial specimen, respectively, and bonding 1∼4 spherical particles generate the clumps. In order to distinguish the differences in the geometric shapes of particles, Yang and Luo [59] defined the overall regularity value to describe the geometric characteristics of the particles. The smaller the value is, the more complex the particle morphology is. The overall regularity values for spherical, ellipsoid, pyramidal, and cubic particles are 1, 0.91, 0.88, and 0.89, respectively.Figure 2
The four particle shapes in the specimens. (a) Spherical. (b) Elongate clump. (c) Pyramid clump. (d) Cube clump.
(a)(b)(c)(d)Elongate, pyramid, and cube clumps are used to replace spherical particles with clamp function. Figure3 shows replacing a single spherical particle with an elongate clump. It should be pointed out that particle replacement is carried out one by one in which particles are arranged. In the process of particle replacement, only the particle shape is changed, but other material parameters of the particle are not changed, and the volume of the DEM model is constant during the process of generation irregular particles.Figure 3
Particle replacement process.


### 3.2. Simulation Processes

In this paper, the rectangular models are established with a ratio of length, width, and height of 1 : 1:2 (Figure4). Rigid walls are set on six faces of the models (the boundaries are the rigid boundaries) [60, 61]. The contact model is linear. The friction coefficient between particles is 0.4, and the friction coefficient between particles and the rigid walls is 0. According to the target porosity of 0.35, 12000 spherical particles are generated in the model with the random distribution method. The spherical particle parameters are adopted in Table 1.Figure 4
Schematic of the true-triaxial specimen.Table 1
The parameters used in DEM simulation.
Density/ρ (kg/m3)Modulus of elasticityE (MPa)Normal contact stiffnessKn (N/m)Normal to tangential stiffness ratioDamping coefficientParticle radiusd (mm)1600201 × 1082.000.71.8∼2.5In true-triaxial test simulations, intermediate principal stress coefficientb=σ2−σ3/σ1−σ3 and mean normal stress p=σ1+σ2+σ3/3 are required to be kept constant throughout the loading process, and then the loading path in the three principal directions is(15)dσ1=dσ1,dσ2=1−2bb−2dσ1,dσ3=1+bb−2dσ1,where dσ1, dσ2 and dσ3 are the increments of σ1, σ2, and σ3 respectively. For ease of understanding, Figure 5 is the schematic of the loading.Figure 5
Schematic illustration of the loading process.In the true-triaxial simulation test, isotropic consolidation is applied earlier, and then shearing is conducted. For isotropic consolidation, 500 kPa confining pressureσc (σc = p) is applied to the six rigid walls. The model gradually achieved stress equilibrium under the confining pressure, and the first figure in Figure 5 shows the isotropic consolidation. Then, the vertical downward velocity is applied to the rigid wall at the top of the model, and compressive stress is increased on the rigid walls at the four sides of the model. In other words, under the condition of σc = 500 kPa, dσ1 is loaded in the vertical direction, and then four lateral stress increments are obtained by the relationship between dσ2, dσ3 and dσ1 (equation (15)). The shearing process ended when the vertical strain ε1 reached 20%. The shearing process can be shown in the second picture of Figure 5. The test schemes are shown in Table 2. The stress paths of the tests are constant p and b.Table 2
Schemes of true-triaxial simulation.
Particle shapesIntermediate principal stress coefficientbMean normal stressp (kPa)Spherical0, 0.2, 0.4, 0.6, 0.8,1500Elongate clumpPyramid clumpCube clump


### 3.3. Influences of Particle Shapes on the Strength

For obtaining the influences of the particle shapes on macroscopic strength, a true-triaxial test is carried out by PFC3D. Figure 6 shows the stress-strain relationships of the specimens of the four particle shapes under the same b. The purpose is to get the variation of generalized shear stress q=σ1−σ22+σ2−σ32+σ3−σ12/2. The horizontal axis of the curves is the vertical strain ε1 of the specimens in Figure 6.Figure 6
Stress-strain relationships with the four particle shapes under different intermediate principal stress coefficients. (a)b = 0 (b) b = 0.2 (c) b = 0.4 (d) b = 0.6 (e) b = 0.8 (f) b = 1.0.
(a)(b)(c)(d)(e)(f)It can be seen from Figure6 that the strain softening and hardening characteristics of the same particle shape increase with the increase of b. The more complex the particle shape is under the same b, the greater the stress-strain relationship’s influence is. The spherical particle specimen has the most significant slope in the strain softening process, and the pyramids have the slightest slope. This law is consistent with the experimental and simulation results carried out by Yang and Luo [59], Jerves et al. [62], and Lashkari et al. [63].Figures7 and 8 show the relationships between the peak strength (qmax) and peak internal friction angle (φmax) of the four particle shapes with different b, respectively. The peak strength and the peak internal friction angles are the maximum values of the generalized shear stress and the internal friction angle, which can be expressed as sinφ=σ1−σ3/σ1+σ3. The results of all the specimens in Figure 7 show a monotonous decreasing trend of qmax with the increase of b, which is consistent with the true-triaxial test trend of granular materials. Moreover, with the complexity of particle shapes increasing, the qmax is greater at the same b. Similarly, the relationships between the φmax and b in Figure 8 are also consistent with the test results of granular materials. The φmax increases and then decreases. With the particle shapes becoming complex, this trend became more apparent. Some existing research has also reached similar conclusions. The specimens of spherical particles have a lower friction angle and shear strength [64]. Some experiments and simulations revealed that granular materials comprising irregular particles usually have higher macroscopic shear strength at a given stress path [65, 66].Figure 7
Relationships between peak strength and intermediate principal stress coefficient with the four particle shapes.Figure 8
Evolutions of peak friction angle with intermediate principal stress coefficient with the four particle shapes.


### 3.4. Stress Path and Fabric Variations

The evolution of the fabric in three principal directions is extracted during the test under the true-triaxial stress path. From shear to failure,F11 first increases, decreases, and then stabilizes. The rules are the same as the pattern of stress-strain relationship from hardening, peaking to softening, finally. Oda et al. [32] conclude that the contact normal tends to align in the maximum principal stress direction to resist greater external forces during the shearing. Ultimately, the main direction of the fabric will be consistent with the direction of principal stress. Due to space limitations, this article only gives the data at the critical state.Figure9(a) shows the true-triaxial stress path of pressure control, and Figure 9(b) ∼ 9(e) shows the evolutions of the principal fabric components of the four particle shapes at the peak state. With the variations of b, F11 and F33 gradually decrease, while F22 increases. The evolutions of contact fabric at the peak state are consistent with the true-triaxial test of stress path under the constant b and p. The varying degree of the principal fabric with spherical particles is far less than that of other nonspherical particles, and the contact fabric changes continuously as the stress. The more complex the particle shape is, the greater the discontinuity of the fabric is when b = 1. It can be found that the complexity of particle shapes directly affected the evolutions of the contact normal vector.Figure 9
Evolutions of principal stresses and fabric components with intermediate principal stress coefficient. (a) The stress of true triaxial test at the critical state. (b) The fabric of the true triaxial test at the critical state.
(a)(b)It can be seen from Figure6 that the strain softening and hardening characteristics of the same particle shape increase with the increase of b. The more complex the particle shape is under the same b, the greater the stress-strain relationship’s influence is. The spherical particle specimen has the most significant slope in the strain softening process, and the pyramids have the slightest slope. This law is consistent with the experimental and simulation results carried out by Yang and Luo [59], Jerves et al. [62], and Lashkari et al. [63].Figures7 and 8 show the relationships between the peak strength (qmax) and peak internal friction angle (φmax) of the four particle shapes with different b, respectively. The peak strength and the peak internal friction angles are the maximum values of the generalized shear stress and the internal friction angle, which can be expressed as sinφ=σ1−σ3/σ1+σ3. The results of all the specimens in Figure 7 show a monotonous decreasing trend of qmax with the increase of b, which is consistent with the true-triaxial test trend of granular materials. Moreover, with the complexity of particle shapes increasing, the qmax is greater at the same b. Similarly, the relationships between the φmax and b in Figure 8 are also consistent with the test results of granular materials. The φmax increases and then decreases. With the particle shapes becoming complex, this trend became more apparent. Some existing research has also reached similar conclusions. The specimens of spherical particles have a lower friction angle and shear strength [64]. Some experiments and simulations revealed that granular materials comprising irregular particles usually have higher macroscopic shear strength at a given stress path [65, 66].


### 3.5. Evolutions of Fabric

Both inherent anisotropy and stress-induced anisotropy are closely related to the microscopic properties of the granular materials. Casagrande and Carillo [67] first discriminate the difference between them. The sedimentation of a particle will trigger inherent anisotropy, while the induced anisotropy is mainly caused by nonelastic deformation [68]. Oda et al. [32] revealed that induced anisotropy is mainly produced by changing the distribution of internal contact normal of materials. The four particle shapes can make different inherent anisotropic specimens in our simulation. Changing the loading ratio of the three main directions of the specimen can simulate stress-induced anisotropy. The relationships between the two anisotropies can be explored by analyzing the particle contact evolution data.Based on the simulation results above, it is considered thatFij evolutions of granular materials are the function of principal stress. Hu et al. [69] and Lashkari and Norouzi [70] all established the relationship between stress ratio and fabric evolution in order to analyze the distribution of fabric evolution and anisotropy. In this paper, the relationship between the fabric tensor and the deviatoric stress ratio suggested by Guo et al. [14] is used to analyze the law of fabric evolution. Guo [14] supposed that Fij is proportional to the components of the deviatoric stress ratio tensor:(16)dFij=λdηij,where λ is a fabric evolution parameter that could describe the stress-fabric relationship, dηij denotes deviatoric stress ratio, dηij=dsij/p, sij is deviatoric stress tensor, and sij=σij−pδij, δij is Kronecker notation. Under constant p loading, p is a constant, so equation (16) can be turned into(17)dFij=λdσij−pp=λdσijp.We can get the fabric-stress evolution relationship by integrating equation (17):(18)Fij=λpσij+C,where C is an arbitrary constant. Assume that the sample has an initial fabric Fij|0 in the isotropic state, and then equation (18) can be written as(19)Fij=λpσij+Fij|0.For the initial isotropic sample, the initial fabric is equal to 1/3. The fabric evolution parameterλ of granular material can be calculated by substituting Fij|0=1/3 into equation (19). The experimental evidence by Oda et al. [32] revealed that induced anisotropy is mainly produced by changing the distribution of internal contact of materials. To describe the effect of fabric evolution on anisotropy, we used the invariant of fabric tensor as a reference quantity as follows:(20)Fq=22F11−F222+F22−F332+F11−F332.Figure10 shows the variation of λ of the four particle shapes with different b at the peak state. When the b is small, the change rate of the λ-b curve of spherical particles is small, while the change rate of cubic and pyramidal particles is larger. When the b is larger, the rate of change of λ-b curves of different shapes is approximately constant. The overall change of λ is smaller and around 1. For the nonspherical particle, λ decreased first as b increased, and it decreased slowly after b = 0.6; then, the rate of change tends to be constant. The more irregular the particle shape is, the more considerable the λ is, and the more significant the reduction when b < 0.6. It can be seen that the spherical particle is easy to rotate in the loading process, but due to the regular shape, the spatial probability and distribution probability of the particle contact point tend to be even. However, for the nonspherical particles, with the change of particle shapes and principal stress loading ratio, the spatial distribution of the contact points is different, and the variation rules of fabric evolution parameters are quite different.Figure 10
Evolutions of fabric evolution parameter with the intermediate principal stress coefficient.Figure11 shows the evolutions of Fq with b of the four particle shapes in the peak state. The expression of Fq is the same as the expression of q, which is the second invariant of fabric and stress partial tensor, respectively. Fq can describe the degree of anisotropy of contact fabric. The greater the Fq is, the higher the degree of anisotropy is. For the spherical particles, Fq increases monotonously with b, which is contrary to the evolution that q changed with b in Figure 7. When the particle morphology becomes complicated, Fq does not increase monotonously with b but shows a different anisotropy. Therefore, Fq can be used to quantify the evolution of fabrics caused by stress-induced anisotropy.Figure 11
Evolutions of deviatoric fabric with intermediate principal stress coefficient.Figure12 shows the rose diagram on the XZ plane of the distribution of contact number and contact force with the different b. The black curves in the figure are the distribution of contact characteristics after isotropic consolidation, and it can be seen that the curve is close to the isotropic distribution. Other curves show apparent anisotropy distribution, indicating that the contact fabric shows different degrees of anisotropy distribution under the conditions of different proportions, and that can directly reflect the stress-induced anisotropy. The distribution of contact numbers in Figure 12(a) and contact force in Figure 12(b) are similar for spherical particles, but the anisotropic degrees are different. There is a similar rule for other irregular particles, only to a different degree. This rule indicates that the magnitude of the force on the contact points aggravate the degree of anisotropy in the distribution direction of the same contact number. Therefore, a more reasonable method to describe the contact fabric is to consider the distribution of contact points’ numbers and the contact force’s magnitude.Figure 12
Rose diagrams of contact number and contact force. (a) Distribution of contact number of spherical. (b) Distribution of contact number of the elongate. (c) Distribution of contact number of the pyramid. (d) Distribution of contact number of the cube. (e) Distribution of contact force of the spherical. (f) Distribution of contact force of the elongate. (g) Distribution of contact force of the pyramid. (h) Distribution of contact force of the cube.
(a)(b)(c)(d)(e)(f)(g)(h)Figure13 shows the evolutions of Fq with stress ratio (η = q/p) of the four particle shapes. The results show that Fq and η are increased to peak point first and then reduced to a critical state under the stress path of constant p. Fq increases with the increase of b, and η decreases with the increase of b, which is consistent with the research results of Yuan and Yu [53]. As the irregularity of the particles increases, η and Fq will also gradually increase. The micro-macroscopic mechanical interpretation can be given as follows. Due to the uniform distribution of the fabric and the smoothness of spherical particle specimens, the particles are easy to rotate, so the shear stress is relatively lower, while for irregular particles, the spatial distribution of fabric is more complex, and particle movement requires greater stress. It is consistent with the observation that interlocking between spherical particles is more unstable than between angular particles in reference [71].Figure 13
Evolutions of deviatoric fabric with stress ration. (a) Evolutions ofFq with η of spherical specimens. (b) Evolutions of Fq with η of elongate specimens. (c) Evolutions of Fq with η of pyramid specimens. (d) Evolutions of Fq with η of cube specimens.
(a)(b)(c)(d)


## 3.1. Generation of Nonspherical Particles

In DEM, spherical particle (3D) with different diameters is used, which is too ideal to describe the natural geometry of granular materials, especially when the microscopic particle characteristics are used to describe the influence of their fabric on macroscopic mechanism. Therefore, the four particle shapes, that is, spheres, elongate clump, pyramid clump, and cube clump (Figure2), are used to establish a true-triaxial specimen, respectively, and bonding 1∼4 spherical particles generate the clumps. In order to distinguish the differences in the geometric shapes of particles, Yang and Luo [59] defined the overall regularity value to describe the geometric characteristics of the particles. The smaller the value is, the more complex the particle morphology is. The overall regularity values for spherical, ellipsoid, pyramidal, and cubic particles are 1, 0.91, 0.88, and 0.89, respectively.Figure 2
The four particle shapes in the specimens. (a) Spherical. (b) Elongate clump. (c) Pyramid clump. (d) Cube clump.
(a)(b)(c)(d)Elongate, pyramid, and cube clumps are used to replace spherical particles with clamp function. Figure3 shows replacing a single spherical particle with an elongate clump. It should be pointed out that particle replacement is carried out one by one in which particles are arranged. In the process of particle replacement, only the particle shape is changed, but other material parameters of the particle are not changed, and the volume of the DEM model is constant during the process of generation irregular particles.Figure 3
Particle replacement process.


## 3.2. Simulation Processes

In this paper, the rectangular models are established with a ratio of length, width, and height of 1 : 1:2 (Figure4). Rigid walls are set on six faces of the models (the boundaries are the rigid boundaries) [60, 61]. The contact model is linear. The friction coefficient between particles is 0.4, and the friction coefficient between particles and the rigid walls is 0. According to the target porosity of 0.35, 12000 spherical particles are generated in the model with the random distribution method. The spherical particle parameters are adopted in Table 1.Figure 4
Schematic of the true-triaxial specimen.Table 1
The parameters used in DEM simulation.
Density/ρ (kg/m3)Modulus of elasticityE (MPa)Normal contact stiffnessKn (N/m)Normal to tangential stiffness ratioDamping coefficientParticle radiusd (mm)1600201 × 1082.000.71.8∼2.5In true-triaxial test simulations, intermediate principal stress coefficientb=σ2−σ3/σ1−σ3 and mean normal stress p=σ1+σ2+σ3/3 are required to be kept constant throughout the loading process, and then the loading path in the three principal directions is(15)dσ1=dσ1,dσ2=1−2bb−2dσ1,dσ3=1+bb−2dσ1,where dσ1, dσ2 and dσ3 are the increments of σ1, σ2, and σ3 respectively. For ease of understanding, Figure 5 is the schematic of the loading.Figure 5
Schematic illustration of the loading process.In the true-triaxial simulation test, isotropic consolidation is applied earlier, and then shearing is conducted. For isotropic consolidation, 500 kPa confining pressureσc (σc = p) is applied to the six rigid walls. The model gradually achieved stress equilibrium under the confining pressure, and the first figure in Figure 5 shows the isotropic consolidation. Then, the vertical downward velocity is applied to the rigid wall at the top of the model, and compressive stress is increased on the rigid walls at the four sides of the model. In other words, under the condition of σc = 500 kPa, dσ1 is loaded in the vertical direction, and then four lateral stress increments are obtained by the relationship between dσ2, dσ3 and dσ1 (equation (15)). The shearing process ended when the vertical strain ε1 reached 20%. The shearing process can be shown in the second picture of Figure 5. The test schemes are shown in Table 2. The stress paths of the tests are constant p and b.Table 2
Schemes of true-triaxial simulation.
Particle shapesIntermediate principal stress coefficientbMean normal stressp (kPa)Spherical0, 0.2, 0.4, 0.6, 0.8,1500Elongate clumpPyramid clumpCube clump


## 3.3. Influences of Particle Shapes on the Strength

For obtaining the influences of the particle shapes on macroscopic strength, a true-triaxial test is carried out by PFC3D. Figure 6 shows the stress-strain relationships of the specimens of the four particle shapes under the same b. The purpose is to get the variation of generalized shear stress q=σ1−σ22+σ2−σ32+σ3−σ12/2. The horizontal axis of the curves is the vertical strain ε1 of the specimens in Figure 6.Figure 6
Stress-strain relationships with the four particle shapes under different intermediate principal stress coefficients. (a)b = 0 (b) b = 0.2 (c) b = 0.4 (d) b = 0.6 (e) b = 0.8 (f) b = 1.0.
(a)(b)(c)(d)(e)(f)It can be seen from Figure6 that the strain softening and hardening characteristics of the same particle shape increase with the increase of b. The more complex the particle shape is under the same b, the greater the stress-strain relationship’s influence is. The spherical particle specimen has the most significant slope in the strain softening process, and the pyramids have the slightest slope. This law is consistent with the experimental and simulation results carried out by Yang and Luo [59], Jerves et al. [62], and Lashkari et al. [63].Figures7 and 8 show the relationships between the peak strength (qmax) and peak internal friction angle (φmax) of the four particle shapes with different b, respectively. The peak strength and the peak internal friction angles are the maximum values of the generalized shear stress and the internal friction angle, which can be expressed as sinφ=σ1−σ3/σ1+σ3. The results of all the specimens in Figure 7 show a monotonous decreasing trend of qmax with the increase of b, which is consistent with the true-triaxial test trend of granular materials. Moreover, with the complexity of particle shapes increasing, the qmax is greater at the same b. Similarly, the relationships between the φmax and b in Figure 8 are also consistent with the test results of granular materials. The φmax increases and then decreases. With the particle shapes becoming complex, this trend became more apparent. Some existing research has also reached similar conclusions. The specimens of spherical particles have a lower friction angle and shear strength [64]. Some experiments and simulations revealed that granular materials comprising irregular particles usually have higher macroscopic shear strength at a given stress path [65, 66].Figure 7
Relationships between peak strength and intermediate principal stress coefficient with the four particle shapes.Figure 8
Evolutions of peak friction angle with intermediate principal stress coefficient with the four particle shapes.


## 3.4. Stress Path and Fabric Variations

The evolution of the fabric in three principal directions is extracted during the test under the true-triaxial stress path. From shear to failure,F11 first increases, decreases, and then stabilizes. The rules are the same as the pattern of stress-strain relationship from hardening, peaking to softening, finally. Oda et al. [32] conclude that the contact normal tends to align in the maximum principal stress direction to resist greater external forces during the shearing. Ultimately, the main direction of the fabric will be consistent with the direction of principal stress. Due to space limitations, this article only gives the data at the critical state.Figure9(a) shows the true-triaxial stress path of pressure control, and Figure 9(b) ∼ 9(e) shows the evolutions of the principal fabric components of the four particle shapes at the peak state. With the variations of b, F11 and F33 gradually decrease, while F22 increases. The evolutions of contact fabric at the peak state are consistent with the true-triaxial test of stress path under the constant b and p. The varying degree of the principal fabric with spherical particles is far less than that of other nonspherical particles, and the contact fabric changes continuously as the stress. The more complex the particle shape is, the greater the discontinuity of the fabric is when b = 1. It can be found that the complexity of particle shapes directly affected the evolutions of the contact normal vector.Figure 9
Evolutions of principal stresses and fabric components with intermediate principal stress coefficient. (a) The stress of true triaxial test at the critical state. (b) The fabric of the true triaxial test at the critical state.
(a)(b)It can be seen from Figure6 that the strain softening and hardening characteristics of the same particle shape increase with the increase of b. The more complex the particle shape is under the same b, the greater the stress-strain relationship’s influence is. The spherical particle specimen has the most significant slope in the strain softening process, and the pyramids have the slightest slope. This law is consistent with the experimental and simulation results carried out by Yang and Luo [59], Jerves et al. [62], and Lashkari et al. [63].Figures7 and 8 show the relationships between the peak strength (qmax) and peak internal friction angle (φmax) of the four particle shapes with different b, respectively. The peak strength and the peak internal friction angles are the maximum values of the generalized shear stress and the internal friction angle, which can be expressed as sinφ=σ1−σ3/σ1+σ3. The results of all the specimens in Figure 7 show a monotonous decreasing trend of qmax with the increase of b, which is consistent with the true-triaxial test trend of granular materials. Moreover, with the complexity of particle shapes increasing, the qmax is greater at the same b. Similarly, the relationships between the φmax and b in Figure 8 are also consistent with the test results of granular materials. The φmax increases and then decreases. With the particle shapes becoming complex, this trend became more apparent. Some existing research has also reached similar conclusions. The specimens of spherical particles have a lower friction angle and shear strength [64]. Some experiments and simulations revealed that granular materials comprising irregular particles usually have higher macroscopic shear strength at a given stress path [65, 66].


## 3.5. Evolutions of Fabric

Both inherent anisotropy and stress-induced anisotropy are closely related to the microscopic properties of the granular materials. Casagrande and Carillo [67] first discriminate the difference between them. The sedimentation of a particle will trigger inherent anisotropy, while the induced anisotropy is mainly caused by nonelastic deformation [68]. Oda et al. [32] revealed that induced anisotropy is mainly produced by changing the distribution of internal contact normal of materials. The four particle shapes can make different inherent anisotropic specimens in our simulation. Changing the loading ratio of the three main directions of the specimen can simulate stress-induced anisotropy. The relationships between the two anisotropies can be explored by analyzing the particle contact evolution data.Based on the simulation results above, it is considered thatFij evolutions of granular materials are the function of principal stress. Hu et al. [69] and Lashkari and Norouzi [70] all established the relationship between stress ratio and fabric evolution in order to analyze the distribution of fabric evolution and anisotropy. In this paper, the relationship between the fabric tensor and the deviatoric stress ratio suggested by Guo et al. [14] is used to analyze the law of fabric evolution. Guo [14] supposed that Fij is proportional to the components of the deviatoric stress ratio tensor:(16)dFij=λdηij,where λ is a fabric evolution parameter that could describe the stress-fabric relationship, dηij denotes deviatoric stress ratio, dηij=dsij/p, sij is deviatoric stress tensor, and sij=σij−pδij, δij is Kronecker notation. Under constant p loading, p is a constant, so equation (16) can be turned into(17)dFij=λdσij−pp=λdσijp.We can get the fabric-stress evolution relationship by integrating equation (17):(18)Fij=λpσij+C,where C is an arbitrary constant. Assume that the sample has an initial fabric Fij|0 in the isotropic state, and then equation (18) can be written as(19)Fij=λpσij+Fij|0.For the initial isotropic sample, the initial fabric is equal to 1/3. The fabric evolution parameterλ of granular material can be calculated by substituting Fij|0=1/3 into equation (19). The experimental evidence by Oda et al. [32] revealed that induced anisotropy is mainly produced by changing the distribution of internal contact of materials. To describe the effect of fabric evolution on anisotropy, we used the invariant of fabric tensor as a reference quantity as follows:(20)Fq=22F11−F222+F22−F332+F11−F332.Figure10 shows the variation of λ of the four particle shapes with different b at the peak state. When the b is small, the change rate of the λ-b curve of spherical particles is small, while the change rate of cubic and pyramidal particles is larger. When the b is larger, the rate of change of λ-b curves of different shapes is approximately constant. The overall change of λ is smaller and around 1. For the nonspherical particle, λ decreased first as b increased, and it decreased slowly after b = 0.6; then, the rate of change tends to be constant. The more irregular the particle shape is, the more considerable the λ is, and the more significant the reduction when b < 0.6. It can be seen that the spherical particle is easy to rotate in the loading process, but due to the regular shape, the spatial probability and distribution probability of the particle contact point tend to be even. However, for the nonspherical particles, with the change of particle shapes and principal stress loading ratio, the spatial distribution of the contact points is different, and the variation rules of fabric evolution parameters are quite different.Figure 10
Evolutions of fabric evolution parameter with the intermediate principal stress coefficient.Figure11 shows the evolutions of Fq with b of the four particle shapes in the peak state. The expression of Fq is the same as the expression of q, which is the second invariant of fabric and stress partial tensor, respectively. Fq can describe the degree of anisotropy of contact fabric. The greater the Fq is, the higher the degree of anisotropy is. For the spherical particles, Fq increases monotonously with b, which is contrary to the evolution that q changed with b in Figure 7. When the particle morphology becomes complicated, Fq does not increase monotonously with b but shows a different anisotropy. Therefore, Fq can be used to quantify the evolution of fabrics caused by stress-induced anisotropy.Figure 11
Evolutions of deviatoric fabric with intermediate principal stress coefficient.Figure12 shows the rose diagram on the XZ plane of the distribution of contact number and contact force with the different b. The black curves in the figure are the distribution of contact characteristics after isotropic consolidation, and it can be seen that the curve is close to the isotropic distribution. Other curves show apparent anisotropy distribution, indicating that the contact fabric shows different degrees of anisotropy distribution under the conditions of different proportions, and that can directly reflect the stress-induced anisotropy. The distribution of contact numbers in Figure 12(a) and contact force in Figure 12(b) are similar for spherical particles, but the anisotropic degrees are different. There is a similar rule for other irregular particles, only to a different degree. This rule indicates that the magnitude of the force on the contact points aggravate the degree of anisotropy in the distribution direction of the same contact number. Therefore, a more reasonable method to describe the contact fabric is to consider the distribution of contact points’ numbers and the contact force’s magnitude.Figure 12
Rose diagrams of contact number and contact force. (a) Distribution of contact number of spherical. (b) Distribution of contact number of the elongate. (c) Distribution of contact number of the pyramid. (d) Distribution of contact number of the cube. (e) Distribution of contact force of the spherical. (f) Distribution of contact force of the elongate. (g) Distribution of contact force of the pyramid. (h) Distribution of contact force of the cube.
(a)(b)(c)(d)(e)(f)(g)(h)Figure13 shows the evolutions of Fq with stress ratio (η = q/p) of the four particle shapes. The results show that Fq and η are increased to peak point first and then reduced to a critical state under the stress path of constant p. Fq increases with the increase of b, and η decreases with the increase of b, which is consistent with the research results of Yuan and Yu [53]. As the irregularity of the particles increases, η and Fq will also gradually increase. The micro-macroscopic mechanical interpretation can be given as follows. Due to the uniform distribution of the fabric and the smoothness of spherical particle specimens, the particles are easy to rotate, so the shear stress is relatively lower, while for irregular particles, the spatial distribution of fabric is more complex, and particle movement requires greater stress. It is consistent with the observation that interlocking between spherical particles is more unstable than between angular particles in reference [71].Figure 13
Evolutions of deviatoric fabric with stress ration. (a) Evolutions ofFq with η of spherical specimens. (b) Evolutions of Fq with η of elongate specimens. (c) Evolutions of Fq with η of pyramid specimens. (d) Evolutions of Fq with η of cube specimens.
(a)(b)(c)(d)


## 4. Conclusions

The normal direction of particle contact substantially influences the macroscopic mechanical properties of granular materials. In this paper, the true-triaxial tests with four particle shapes are simulated by DEM, the microscopic evolutions of a particle are studied using the novel defined fabric tensor, and the relationships between the evolution of contact fabric and the macroscopic mechanical behaviors are explored. The main conclusions are as follows:(1)
The contact fabric tensor is defined by the contact normal vector, which describes the probabilistic and statistical laws of the contact characteristics of microscopic particles. Three amplitude parameters in the orthogonal direction are defined using the invariant of the plane fabric tensor, by which the scalar expression of orthotropic fabric is derived. With the change of the geometric relationship of the contact points, the orthotropic fabric can naturally degenerate into different forms of transverse isotropy. The fabric tensor defined in this paper can be directly applied to the macroscopic constitutive equation.(2)
The four particle shapes are constructed using PFC3D software. Under the condition of the same parameters and loading paths of the true-triaxial test, simulation results show that particle morphology changes directly affect the anisotropy of the stress-strain relationship and the strength. The more complex the particle shape is, the more significant the influence on the anisotropy is. These showed the influence of particle morphology on the macro-mechanics of granular materials.(3)
The contact characteristics of particles directly affect their macroscopic mechanical response. The contact fabric can be used to describe the contact characteristics and the evolution of particles. Under the true-triaxial loading path, the spherical particles are easy to rotate, the distributions of contact points are uniform, and the variations of fabric evolution parameters are small. Even if the three principal stress directions are loaded in different proportions, the effect on anisotropy is small. Irregular particles greatly influence the spatial distributions of contact points, and the more pronounced the fabric evolution is, the more the anisotropy changes. The distribution of contact points, contact force, and the evolution of the four particle shapes show that the distribution of particle contact points and the magnitude of the contact forces should be considered in the fabric tensor.



---

*Source: 4765887-2022-08-09.xml*