# Investigating the Combined Effects of Inherent and Stress-Induced Anisotropy on the Mechanical Behavior of Granular Materials Using Three-Dimensional Discrete Element Method

Authors: Xinran Chen; Jinsong Qian; Lei Zhang; Jianming Ling
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The three-dimensional discrete element method (DEM) was employed to investigate the combined effects of inherent and stress-induced anisotropy of granular materials. The particles were modeled following real particle shapes. Both isotropic and inherently anisotropic specimens were prepared, and then true triaxial numerical tests were conducted using different intermediate principle stress ratios (b). The results indicate that the oriented particles in the anisotropic specimens form strong contacts in their long axis direction in the early stages of shearing, which restrains the contraction of the specimens. As the strain increases, the oriented particles start to rotate and slide, which results in shorter contraction stages and fewer number of interparticle contacts with peak values compared to the isotropic specimens. In addition, the increase in b values aggravates the rotating and sliding of particles in the inherently anisotropic specimens and restrains the contraction of the granular and the increase of contact forces. As a result, the inherent anisotropy reduces the effects of stress-induced anisotropy on the mechanical behavior of granular materials.


---

## Body

## 1. Introduction

Granular materials are commonly used in the base layers of typical pavement structures. The directional distribution of particles in granular materials significantly affects their mechanical properties [1], which leads to an anisotropic phenomenon. The anisotropy of granular material can be divided into inherent anisotropy and stress-induced anisotropy [2]. Inherent anisotropy is induced by the directional arrangement of particles during compaction or deposition. Stress-induced anisotropy is caused by the nonuniform distribution of external loading, which leads to the rearrangement of the particles [3]. To investigate these two kinds of anisotropy, the long axis orientation (or “bedding angle”) of granular particles is widely used for inherent anisotropy analysis, and the intermediate principle stress ratio, denoted by b, which represents the relative magnitude of the intermediate principle stress, is used for stress-induced anisotropy analysis [4, 5].Several researchers have conducted laboratory tests to study the effect of anisotropy on the mechanical responses of particle assemblies. For example, Guo [6] conducted direct shear tests using specimens with different bedding angles. The results showed that the friction angle of granular materials varies with the orientation of the shear plane relative to the bedding plane. Besides, the degree of anisotropy is affected by particle shapes. Shi et al. [7] investigated the effect of intermediate principle stress on the macroscopic mechanical responses of a coarse-grained soil using true triaxial tests. The results showed that shear strength and intermediate principle strain increase as the b value increases, whereas the minor principle strain decreases. However, some experimental studies have shown uncertain relationships among the angle of shear resistance, b values, and other controversial behavior [8, 9], as laboratory test results can be affected by the initial fabric difference [10, 11], different sample preparation methods [12, 13], different boundary conditions, etc. [14]. In recent years, researchers used the discrete element method (DEM) [15] to investigate the behavior of granular materials from both macroscopic and microscopic perspectives [16, 17]. The DEM can simulate several test conditions using one specimen and thus avoid the initial fabric difference. Furthermore, microscopic characteristics at the particulate level are easy to obtain and can help to explain the stress and strain evolution during loading [18–20].Several studies have been undertaken to investigate the anisotropy of soil and granular materials using the DEM. For example, Hosseininia [21] studied the inherent anisotropy of granular soils using a two-dimensional DEM and found that the initial distribution of elongated particles and associated voids vary during shear deformation. Besides, the shear strength and deformability of granular materials are highly dependent on the initial fabric condition. Barreto and O’Sullivan [22] studied the stress-induced anisotropy of soil under a constant mean load using a three-dimensional (3D) DEM. Results showed that the friction coefficient affects the inherent stability of the strong force chains, whereas the intermediate stress ratio affects the lateral support provided to these force chains. Kuhn et al. [23] conducted numerical simulations using a 3D DEM and described the anisotropy of the voids through image analysis and Minkowski tensors. Their results showed that stress-induced anisotropy affects not only the mechanical stress-strain relationship but also the hydraulic properties. Some studies that employed DEM also considered the complex stress path [5, 24] and particle breakage [25, 26]. However, related studies have focused mainly on the effect of different b values on isotropic samples; research into the behavior of inherently anisotropic granular materials under different b values is limited.In real practice, the granular materials may be inherently anisotropic due to compaction or deposition, and may also work under complex stress states that cause stress-induced anisotropy. Therefore, the combined effects of the two types of anisotropy may affect the mechanical behavior of granular materials simultaneously and thus needs further investigations.For this purpose, 3D DEM simulations were conducted for this study. Isotropic and inherently anisotropic granular specimens were prepared, and then true triaxial numerical tests were conducted using differentb values. Finally, the macroscopic and microscopic responses, including stress-strain characteristics, interparticle contacts, and the anisotropies of contact normal and contact force, were analyzed in detail to address the mechanisms of anisotropy.


## 2. Numerical Model

### 2.1. Particle Modeling and Grading

In DEM simulations, the particle shape can obviously affect the anisotropic properties [3, 6, 27, 28]. Although the use of circular discs (in 2D) or spheres (in 3D) can significantly reduce the calculation time, they can also lead to excessive freedom of the particles and thus a higher degree of rotation and dilation compared with real granular particles [21, 29–31]. Therefore, noncircular/nonspherical particles are widely used and typically are oval or spherical particles made by clusters of bonded circles/spheres or overlapping rigid clusters [29, 32, 33]. Some researchers have tried another approach that uses convex polygon-shaped particles [34, 35]. The simulation accuracy and computational efficiency of these approaches are of vital importance and are a primary focus of related research.In this study, the numerical simulations utilized 3D DEM software, Itasca PFC3D (version 5.0). Ten typical shapes of granular particles were considered in the DEM simulations and were modeled in 3D following real particle shapes, as shown in Figure 1(a). Then, based on the bubble pack algorithm proposed by Taghavi [36], the 3D models were filled with different sized balls which partially overlapped each other without generating force to create complex clumps in PFC3D. During the filling process, the largest ball was placed first, followed by decreasingly smaller balls. The distance parameter, which controls the overlap of adjacent balls, and the ratio parameter, which is the ratio of the smallest to the largest ball, both affected the simulation accuracy and calculation efficiency. After several pretests, the distance and ratio parameters were set at 100 and 0.3, respectively. Figure 1(b) presents a comparison of a real particle, the original 3D model, and the PFC3D particle.Figure 1
(a) 3D granular models of different particle shapes and (b) comparison of the real particle, 3D model, and PFC3D particle.
(a)(b)In order to simulate an actual granular base that consists of particles of different sizes, a widely used grade was selected for the DEM simulations [31, 37]. In accordance with previous research studies, fine particles were deleted and the particle sizes were adjusted to consider computational efficiency [38, 39], as shown in Figure 2. It should be noted that, the DEM specimen without fine-graded particles may be less dense and more compressible than actual granular material [40]. This simplification would affect the values of some macro-microindexes, but the comparison among different anisotropic conditions would be less affected.Figure 2
Particle size distribution curve used in the DEM simulation.The total number of particles used for the numerical specimens were calculated using the target porosity of 0.33. To avoid the size effect, previous studies suggest that the diameter of the specimen should be at least 10 times larger than the maximum particle size [38, 41]. Therefore, the dimensions of each granular specimen were set to 15L × 15W × 15H cm3, and 6,884 particles were used. The ten different particle shapes were randomly distributed. Table 1 lists both the particle sizes and numbers of particles used in this study.Table 1
Particle size and number of particles for different specimens.
Particle size (mm)4.75–6.56.5–88–9.59.5–13.2513.25–15Number of particles39251346878531204


### 2.2. Specimen Preparation

Two specimen types were used in DEM simulation, which are the isotropic specimens and inherently anisotropic specimens. Both types of specimens were prepared using the gravitational deposition method [42, 43]. First, a deposition space of 15L × 15W × 30H cm3 was created and enclosed by six frictionless rigid walls, and the particles were then generated in the space, as shown in Figure 3. For the isotropic specimens, the long axis of the particles was randomly distributed, whereas for the inherently anisotropic specimens, the long axis of the particles was arranged in the intermediate principle stress (σ2) direction. During specimen preparation, the particle rotation was fixed to obtain the anticipated anisotropic conditions. Then, the gravity field was introduced using the gravitational forces applied to all particles in the vertical direction; thus, the particles would be deposited on the bottom boundary of the deposition space. To create a dense specimen and reduce the simulation time during specimen preparation, the particle friction coefficient was temporarily set as zero [5, 29]. At the end of the deposition, a space of 15L × 15W × 15H cm3 was chosen as the final specimen size. Finally, the specimen was isotropically compressed by giving the same loading speed to all six rigid walls until the target confining pressure of 100 kPa was reached. Figure 3(e) presents the isotropic and anisotropic specimens after preparation.Figure 3
The procedure of specimen preparation: (a) particles’ generation; (b) gravitational deposition; (c) resizing the walls; (d) isotropic compression; (e) specimens after preparation.
(a)(b)(c)(d)(e)The linear contact model [44] was adopted to simulate the interparticle contacts for computational efficiency. In accordance with related studies which also studied the mechanical behavior of granular materials [30, 39, 45, 46], the material parameters were selected for triaxial testing. The stiffness ratio κ∗, effective modulus E∗ friction coefficient μ, damping constant, and particle density were set at 1.333, 1 × 108 Pa, 0.5, 0.7, and 2.6 g/cm3, respectively.


### 2.3. True Triaxial Shearing

The true triaxial tests were conducted using a strain-controlled loading process that was performed by moving the rigid walls of the specimens. As suggested by Andrade et al. [47], increased loading speed increased the calculation efficiency but would affect the mechanical behavior of granular materials due to the dynamic effect. To select a proper loading speed, a sensitivity analysis was conducted using isotropic specimens under five loading speeds ranged from 0.03 m/s to 0.50 m/s, as depicted in Figure 4. It can be seen that the porosity, which reflects the volume change of the specimens, changes obviously under high loading speeds. By contrast, the porosity is similar for selected small loading speeds. This indicates that the loading speed had a threshold between quasi-static and dynamic. Thus, the loading speed for this study was set at 0.12 m/s to ensure the loading process was quasi-static.Figure 4
The effect of the loading speed sensitivity on the porosity.To simulate stress-induced anisotropy, the intermediate principle stress ratio,b, was induced and is defined as(1)b=σ2−σ3σ1−σ3,where σ1 (i = 1, 2, 3) are the principle stresses. In this paper, the stress and strain indexes take positive values when the granular is compressed.Fiveb values (b = 0.00, 0.25, 0.50, 0.75, 1.00) were selected for the numerical tests. The test under b = 0.00 (σ2 = σ3) represents the triaxial compression test, and the test under b = 1.00 (σ1 = σ2) is the triaxial extension test. The major principle stress σ1 changes when the top wall moves during testing, so the intermediate principle stress can be calculated using the following equation:(2)σ2=bσ1−σ3−σ3.The minor principle stressσ3 was set at 100 kPa for all conditions, and b values were also kept constant during the numerical simulation, so the relationship between σ1 and σ2 could be expressed by dσ2 = bdσ1. During triaxial shearing, the σ2 was dynamically adjusted by stress servo-control method to keep equation (2) satisfied.


### 2.4. Verification of the DEM Simulation

The DEM model was verified by the true triaxial test results conducted by Shi [48]. For the laboratory test, the fine gravel was used with the particle diameter ranged from 2–5 mm. The diameter of the laboratory test specimen is 70 mm × 35 mm × 70 mm, and the loading conditions are the same as the DEM simulation. The long axis of the particles in the laboratory test was randomly distributed, which consists of the isotropic DEM specimen.The comparison between the laboratory test and the numerical simulation results is depicted in Figure5. It can be seen that the stress-strain relationship and the peak stress ratio obtained from DEM simulation and laboratory test share the same trend, and the values are similar to each other. Therefore, it can be indicated that the DEM model is reasonable.Figure 5
Comparison between the DEM model and laboratory tests: (a) the deviator stress; (b) the peak stress ratio.
(a)(b)


## 2.1. Particle Modeling and Grading

In DEM simulations, the particle shape can obviously affect the anisotropic properties [3, 6, 27, 28]. Although the use of circular discs (in 2D) or spheres (in 3D) can significantly reduce the calculation time, they can also lead to excessive freedom of the particles and thus a higher degree of rotation and dilation compared with real granular particles [21, 29–31]. Therefore, noncircular/nonspherical particles are widely used and typically are oval or spherical particles made by clusters of bonded circles/spheres or overlapping rigid clusters [29, 32, 33]. Some researchers have tried another approach that uses convex polygon-shaped particles [34, 35]. The simulation accuracy and computational efficiency of these approaches are of vital importance and are a primary focus of related research.In this study, the numerical simulations utilized 3D DEM software, Itasca PFC3D (version 5.0). Ten typical shapes of granular particles were considered in the DEM simulations and were modeled in 3D following real particle shapes, as shown in Figure 1(a). Then, based on the bubble pack algorithm proposed by Taghavi [36], the 3D models were filled with different sized balls which partially overlapped each other without generating force to create complex clumps in PFC3D. During the filling process, the largest ball was placed first, followed by decreasingly smaller balls. The distance parameter, which controls the overlap of adjacent balls, and the ratio parameter, which is the ratio of the smallest to the largest ball, both affected the simulation accuracy and calculation efficiency. After several pretests, the distance and ratio parameters were set at 100 and 0.3, respectively. Figure 1(b) presents a comparison of a real particle, the original 3D model, and the PFC3D particle.Figure 1
(a) 3D granular models of different particle shapes and (b) comparison of the real particle, 3D model, and PFC3D particle.
(a)(b)In order to simulate an actual granular base that consists of particles of different sizes, a widely used grade was selected for the DEM simulations [31, 37]. In accordance with previous research studies, fine particles were deleted and the particle sizes were adjusted to consider computational efficiency [38, 39], as shown in Figure 2. It should be noted that, the DEM specimen without fine-graded particles may be less dense and more compressible than actual granular material [40]. This simplification would affect the values of some macro-microindexes, but the comparison among different anisotropic conditions would be less affected.Figure 2
Particle size distribution curve used in the DEM simulation.The total number of particles used for the numerical specimens were calculated using the target porosity of 0.33. To avoid the size effect, previous studies suggest that the diameter of the specimen should be at least 10 times larger than the maximum particle size [38, 41]. Therefore, the dimensions of each granular specimen were set to 15L × 15W × 15H cm3, and 6,884 particles were used. The ten different particle shapes were randomly distributed. Table 1 lists both the particle sizes and numbers of particles used in this study.Table 1
Particle size and number of particles for different specimens.
Particle size (mm)4.75–6.56.5–88–9.59.5–13.2513.25–15Number of particles39251346878531204


## 2.2. Specimen Preparation

Two specimen types were used in DEM simulation, which are the isotropic specimens and inherently anisotropic specimens. Both types of specimens were prepared using the gravitational deposition method [42, 43]. First, a deposition space of 15L × 15W × 30H cm3 was created and enclosed by six frictionless rigid walls, and the particles were then generated in the space, as shown in Figure 3. For the isotropic specimens, the long axis of the particles was randomly distributed, whereas for the inherently anisotropic specimens, the long axis of the particles was arranged in the intermediate principle stress (σ2) direction. During specimen preparation, the particle rotation was fixed to obtain the anticipated anisotropic conditions. Then, the gravity field was introduced using the gravitational forces applied to all particles in the vertical direction; thus, the particles would be deposited on the bottom boundary of the deposition space. To create a dense specimen and reduce the simulation time during specimen preparation, the particle friction coefficient was temporarily set as zero [5, 29]. At the end of the deposition, a space of 15L × 15W × 15H cm3 was chosen as the final specimen size. Finally, the specimen was isotropically compressed by giving the same loading speed to all six rigid walls until the target confining pressure of 100 kPa was reached. Figure 3(e) presents the isotropic and anisotropic specimens after preparation.Figure 3
The procedure of specimen preparation: (a) particles’ generation; (b) gravitational deposition; (c) resizing the walls; (d) isotropic compression; (e) specimens after preparation.
(a)(b)(c)(d)(e)The linear contact model [44] was adopted to simulate the interparticle contacts for computational efficiency. In accordance with related studies which also studied the mechanical behavior of granular materials [30, 39, 45, 46], the material parameters were selected for triaxial testing. The stiffness ratio κ∗, effective modulus E∗ friction coefficient μ, damping constant, and particle density were set at 1.333, 1 × 108 Pa, 0.5, 0.7, and 2.6 g/cm3, respectively.


## 2.3. True Triaxial Shearing

The true triaxial tests were conducted using a strain-controlled loading process that was performed by moving the rigid walls of the specimens. As suggested by Andrade et al. [47], increased loading speed increased the calculation efficiency but would affect the mechanical behavior of granular materials due to the dynamic effect. To select a proper loading speed, a sensitivity analysis was conducted using isotropic specimens under five loading speeds ranged from 0.03 m/s to 0.50 m/s, as depicted in Figure 4. It can be seen that the porosity, which reflects the volume change of the specimens, changes obviously under high loading speeds. By contrast, the porosity is similar for selected small loading speeds. This indicates that the loading speed had a threshold between quasi-static and dynamic. Thus, the loading speed for this study was set at 0.12 m/s to ensure the loading process was quasi-static.Figure 4
The effect of the loading speed sensitivity on the porosity.To simulate stress-induced anisotropy, the intermediate principle stress ratio,b, was induced and is defined as(1)b=σ2−σ3σ1−σ3,where σ1 (i = 1, 2, 3) are the principle stresses. In this paper, the stress and strain indexes take positive values when the granular is compressed.Fiveb values (b = 0.00, 0.25, 0.50, 0.75, 1.00) were selected for the numerical tests. The test under b = 0.00 (σ2 = σ3) represents the triaxial compression test, and the test under b = 1.00 (σ1 = σ2) is the triaxial extension test. The major principle stress σ1 changes when the top wall moves during testing, so the intermediate principle stress can be calculated using the following equation:(2)σ2=bσ1−σ3−σ3.The minor principle stressσ3 was set at 100 kPa for all conditions, and b values were also kept constant during the numerical simulation, so the relationship between σ1 and σ2 could be expressed by dσ2 = bdσ1. During triaxial shearing, the σ2 was dynamically adjusted by stress servo-control method to keep equation (2) satisfied.


## 2.4. Verification of the DEM Simulation

The DEM model was verified by the true triaxial test results conducted by Shi [48]. For the laboratory test, the fine gravel was used with the particle diameter ranged from 2–5 mm. The diameter of the laboratory test specimen is 70 mm × 35 mm × 70 mm, and the loading conditions are the same as the DEM simulation. The long axis of the particles in the laboratory test was randomly distributed, which consists of the isotropic DEM specimen.The comparison between the laboratory test and the numerical simulation results is depicted in Figure5. It can be seen that the stress-strain relationship and the peak stress ratio obtained from DEM simulation and laboratory test share the same trend, and the values are similar to each other. Therefore, it can be indicated that the DEM model is reasonable.Figure 5
Comparison between the DEM model and laboratory tests: (a) the deviator stress; (b) the peak stress ratio.
(a)(b)


## 3. Numerical Simulation Results

### 3.1. Macroscopic Analysis: Stress-Strain Characteristics

To evaluate the stress-strain characteristics of granular materials with different initial fabrics and under differentb values, several commonly used macroscopic indexes are selected based on the stress and strain in three principle stress directions. The peak internal friction angle (φmax), mean stress p, shear stress q, major principle strain ε1, and volumetric strain εv are defined, respectively, as equation (3)–(7) [49]:(3)sinφmax=σ1/σ3max−1σ1/σ3max+1,(4)p=σ1+σ2+σ33,(5)q=12σ1−σ22+σ2−σ32+σ1−σ32,(6)ε1=lnH0H,(7)εv=lnVV0,where H0 and V0 are the initial height and volume of the specimen, respectively, and H and V represent the height and volume of the specimen at time t, respectively.The relationship betweenφmax and b value is commonly used to investigate of peak strength in true triaxial testing. Figure 6(a) shows that, for both specimen types, φmax increases as the b value increases from 0.00 to 0.75 and then decreases as the b value increases from 0.75 to 1.00. Similar trend was also found by previous studies using laboratory tests [47, 49] and DEM simulations [4, 5, 11]. The peak value of the mean stress pmax that represents the overall stress level increases as the b value increases for both specimen types, as shown in Figure 6(b). Because all the numerical tests were conducted using the same σ3 value, the σ2 value increases as the b value increases and leads to the increase in the mean stress.Figure 6
Comparison of (a) peak internal friction angle (φmax) and (b) peak mean stress (pmax).
(a)(b)Compared to the isotropic specimens, theφmax and pmax values for anisotropic specimens are lower at the same b values and change less as the b value increases. This finding indicates that the horizontal alignment of the long axis of the particles reduces the effect of nonuniform stress distribution on the stress behavior of granular materials.To highlight the effects of the nonuniform distribution of external loading on the stress behavior of the granular materials, the stress ratio, defined asq/p, is introduced. Figure 7 shows that, for the isotropic specimens, the stress ratio deceases as the b value increases, which is consistent with related laboratory test results [50] and numerical simulation results [4, 14]. Also, the stress ratios for both specimen types are similar at the same b values. This finding indicates that the increased proportion of σ2 results in the decrease of the shear strength of granular materials, and this effect is similar for both isotropic and inherently anisotropic granular specimens.Figure 7
Evolution of stress ratio with the increase ofε1: (a) isotropic specimens; (b) anisotropic specimens.
(a)(b)Figure8 shows the variation of εv with the major principle strain ε1. All specimens exhibit the process from contraction to dilation as ε1 value increases. For the isotropic specimens, the amount of contraction tends to increase as the b value increases. By contrast, the anisotropic specimens are less contracted than the isotropic specimens for all the selected b values, and their contraction stages end with smaller ε1 values. This finding indicates that the inherently anisotropic granular materials are less compressible and may dilate earlier than isotropic granular materials and that the effect of b values on strain behavior is reduced.Figure 8
Variation of volumetric strain with the increase ofε1: (a) isotropic specimens; (b) anisotropic specimens.
(a)(b)


### 3.2. Microscopic Analysis: Contact Characteristics

In order to illustrate the macroscopic behavior of granular materials at the particle scale, the contact characteristics are analyzed. The coordination number (CN) is a commonly used index that is directly related to the stability of granular material [4, 51] and is defined as the average number of contacts per particle:(8)CN=2Nc−NcwNp,where CN is the coordination number, Nc is the total number of contacts, Ncw is the number of contacts between the particle and the wall, and Np is the total number of particles.Figure9 shows the evolution of the CN with the major principle strain. For both specimen types, the CN increases with ε1 to the peak value and then decreases to a residual value. The peak CNs increase obviously as the b value increases. Compared to the isotropic specimens, the inherently anisotropic specimens have higher initial CNs, lower peak CNs, and smaller ε1 values that correspond to the peak CNs. Besides, the increase of their peak CNs as b value increases is much smaller than isotropic specimens. This finding indicates that the oriented arrangement of particles in inherently anisotropic specimens form more contacts in the early stage of shearing, and the change of CNs during shearing is less affected by nonuniform stress distribution.Figure 9
Evolution of coordination number with the increase ofε1: (a) isotropic specimens; (b) anisotropic specimens.
(a)(b)During triaxial testing, the contact sliding typically occurs during the rearrangement of the particles, thus resulting in the permanent deformation of the granular materials. As defined by Gu et al. [39], contact slides when the inequality, shown as equation (9), is satisfied:(9)ftμfn>0.9999,where fn and ft are the normal contact force and tangential contact force, respectively, and μ is the contact friction coefficient.Figure10 presents a comparison of the initial and peak sliding contact percentages during triaxial shearing for all specimens. For both specimen types, the initial sliding contact percentages increase as b value increases from 0.00 to 0.25, and decrease as b value increases from 0.25 to 1.00. The peak sliding contact percentages for the isotropic specimens decrease as the b value increases, which is consistent with the DEM simulation results of Mahmud Sazzad et al. [4]. By contrast, the peak sliding contact percentages for the anisotropic specimens decrease as the b value increases from 0.00 to 0.75, then tend to increase as the b value increases from 0.75 to 1.00. Under the same b values, the anisotropic specimens have lower initial sliding contact percentages and higher peak sliding contact percentages compared to the isotropic specimens. This finding indicates that the sliding of contacts is restrained as the b values increase for isotropic specimens. By contrast, the anisotropic specimens restrain contact sliding in the early stages of testing compared with isotropic specimens, but contact sliding occurs more during shearing, especially with high b values.Figure 10
Comparison of (a) initial sliding contact percentages and (b) peak sliding contact percentages.
(a)(b)Both the CNs and the sliding contact percentages can be used to explain the change of volumetric strain (see Figure8). During the contraction stage, the particles are squeezed together in close contact with each other, which leads to an increase in the CN value. During the dilation stage, the particles start to rotate and slide, which results in the decrease of CN [51]. The specimen which is contracted mostly also has the largest peak CN and the lowest peak sliding contact percentage.


### 3.3. Microscopic Analysis: Fabric Tensor and Anisotropy

The directional distribution of particles, defined as the fabric, plays an important role in the shear behavior of a particle assembly [52–55]. To further investigate the combined effects of inherent and stress-induced anisotropy of granular materials, the contact normal, normal contact force, and tangential contact force are induced.The directional distribution of the contact normal is usually described by the fabric tensorRij [56, 57]:(10)Rij=∫ΩEΩninjdΩ=1N∑c∈Nninj,where ni is the unit contact normal in the i-direction, N is the total number of contacts, and EΩ is the distribution function on the unit sphere Ω and can be calculated by the second-order Fourier expansion as follows [58]:(11)EΩ=14π1+αijrninj,where αijr is the second-order anisotropic tensor and characterizes the fabric anisotropy. This tensor is determined by the deviatoric part of the fabric tensor Rij′ as(12)αijr=152Rij′.For normal and tangential contact force, the similar definition is induced, as expressed in equation (13)–(18):(13)Fijn=14π∫Ωf¯nΩninjdΩ=∑c∈NfnninjN1+αklrnknl,(14)f¯nΩ=f¯0n1+αijnninj,(15)αijn=152Fijn′f¯0n,(16)Fijt=14π∫Ωf¯itΩninjdΩ=∑c∈NftninjN1+αklrnknl,(17)f¯itΩ=f¯0nαijtnj−αkltnknlni,(18)αijt=5Fijt′f¯0n,where fn is the normal contact force, ft is the tangential contact force, and f¯0n is the average normal contact force over all the groups on Ω.Note thatαijr, αijn, and αijt are symmetric and deviatoric, so their second invariants αr, αn, and αt are used to quantify the degree of anisotropy as(19)α∗=32αij∗αij∗,where the super and subscript ∗ represents the contact normal r, tangential contact force t, and normal contact force n, respectively.Figure11 shows that the αr, αn, and αt values increase at the beginning and then decrease to a residual value as ε1 value increases. As the b value increases, the initial slope of αr, shows a small difference, whereas the peak values and residual values decrease obviously. This finding indicates that the contact distribution tends to be more uniform with a higher b value. Besides, the αn values are about 3 to 4 times higher than the αt values for all conditions, which suggests that contact force anisotropy is mainly caused by normal contact forces [51].Figure 11
Effects ofb values on the anisotropy parameter αr, αn, and αt. (a) αr Iso, (b) αr Ani, (c) αn Iso, (d) αn Ani, (e) αt Iso, and (f) αt Ani.
(a)(b)(c)(d)(e)(f)The initialαr values for the two specimen types are significantly different. The isotropic specimens have an initial αr at zero, which means the initial contact distribution is uniform. By contrast, the anisotropic specimens have an initial αr, at about 1.0, which means the oriented arrangement of particles make the initial contact distribution nonuniform. However, the peak and residual αr values for two specimen types are similar, which suggests that inherent anisotropy changes the initial contact distribution, but will not obviously affect the extent of contact normal anisotropy during shearing. For αn and αt, the values change with different b values for anisotropic specimens is smaller than the isotropic specimens, which suggests that inherent anisotropy significantly reduces the effect of b values on the distribution of contact force of granular materials.To address the reason for the difference in anisotropy parametersαr and αn, between the two specimen types, the principle values corresponding to them are induced. The αt is not analyzed here because it is much smaller than αn.The principle values of the fabric tensorR are used to quantify the cluster extent of contacts in the principle stress directions [14]. R1, R2, and R3 represent the contact clustering extent in the σ1, σ2, and σ3 direction, respectively. As shown in Figure 12, the R principle values show similar trends with the increase of ε1 for both specimen types. Before the peak of R principle values, R1 increases, whereas R3 decreases as ε1 value increases for all the selected b values. R2 decreases with the increase of ε1 when b < 0.5 and increases as ε1 value increases when b > 0.5. Such trends are similar to the precious DEM simulation results [11, 14, 22], which also indicates that contacts tend to cluster in the σ1 direction to carry most of the external load. As the b value increases, the contacts also cluster in the σ2 direction.Figure 12
Evolution of theR principle values with an increase in b values: (a) b = 0.00; (b) b = 0.25; (c) b = 0.50; (d) b = 0.75; (e) b = 1.00.
(a)(b)(c)(d)(e)The anisotropic specimens have higherR3 values and lower R2 values compared with the isotropic specimens for all selected b values; thus, the contact normal anisotropy αr for the two specimen types is similar. In addition, for the anisotropic specimens, the R2 value is lower than the R3 value under selected small b values, whereas for the isotropic specimens, the relationship is opposite. This finding indicates that the inherently anisotropic specimens form less contact in the direction of the long axis orientation of particles during shearing, compared with isotropic specimens.Figure13 shows the principle values of the normal contact force tensor Fn (Fn1, Fn2, Fn3). The trends are similar for the three Fn principle changes as the ε1 increases. As b value increases, all three Fn principle values increase obviously, which leads to the increase of peak internal friction angle and peak mean stress (Figure 6).Figure 13
Evolution of theFn principle values with the increase in b values: (a) b = 0.00; (b) b = 0.25; (c) b = 0.50; (d) b = 0.75; (e) b = 1.00; (f) comparison of peak Fn1 values.
(a)(b)(c)(d)(e)(f)Compared to the isotropic specimens, the anisotropic specimens have lower principleFn values (except for Fn2 when b = 0), and the growth rate of contact force for anisotropic specimens is lower as the b value increases, as shown in Figure 13(f). Therefore, the difference in αn for the anisotropic specimens among different b values is smaller (Figure 11), and their peak internal friction angle and peak mean stress are lower than isotropic specimens (Figure 6).For the anisotropic specimens, the clustering of the contact normal in theσ2 direction (R2) is lower than in the σ3 direction (R3) for the selected small b values (Figure 12), whereas the clustering of the contact force in the σ2 direction (Fn2) is higher than the σ3 direction (Fn3), which suggests that the oriented particles form strong contacts in their long axis direction and thus comprise of a strong force subnetwork [51, 56, 59, 60].


### 3.4. Discussion of the Combined Effects of Anisotropy

In summary, the combined effects of inherent anisotropy and stress-induced anisotropy can be illustrated in Figure14. For selected small b values, the long axis of the particles in the anisotropic specimens are initially perpendicular to the σ1 direction to carry most of the external load [21] and resist the sliding tendency in σ1 direction. Also, these particles form strong contacts in the σ2 direction, as explained in Section 3.3. Therefore, the anisotropic specimens form a support frame in the early stages of testing, which leads to higher initial CN values and lower initial sliding contact percentages and makes the specimens less compressible compared with isotropic specimens. As the strain value increases, the isotropic specimens are obviously compressed, whereas for the anisotropic specimens, particle rotating and sliding occur due to the increased deformation in the σ2 direction, which leads to higher peak sliding contact percentages, restrains the increase in the CNs, and shortens the contraction stage.Figure 14
Schematics of the combined effects of inherent and stress-induced anisotropy.The increase inb values increases the proportion of σ2, thus leads to a greater sliding tendency in σ2 direction. For isotropic specimens, the increase in b values leads to higher peak CNs, restrains the sliding of the contacts, and thus makes the specimens more compressible. For inherently anisotropic specimens, the increased b value significantly aggravates the rotating and sliding of the oriented particles and offsets its compressive effect. Therefore, the anisotropic specimens have higher peak sliding contact percentages and lower growth rate for the coordination number and contact force and are less compressible than the isotropic specimens. As a result, the effect of b values on the mechanical behavior of anisotropic specimens is significantly reduced when combined with the inherently anisotropic condition, compared with isotropic specimens.


## 3.1. Macroscopic Analysis: Stress-Strain Characteristics

To evaluate the stress-strain characteristics of granular materials with different initial fabrics and under differentb values, several commonly used macroscopic indexes are selected based on the stress and strain in three principle stress directions. The peak internal friction angle (φmax), mean stress p, shear stress q, major principle strain ε1, and volumetric strain εv are defined, respectively, as equation (3)–(7) [49]:(3)sinφmax=σ1/σ3max−1σ1/σ3max+1,(4)p=σ1+σ2+σ33,(5)q=12σ1−σ22+σ2−σ32+σ1−σ32,(6)ε1=lnH0H,(7)εv=lnVV0,where H0 and V0 are the initial height and volume of the specimen, respectively, and H and V represent the height and volume of the specimen at time t, respectively.The relationship betweenφmax and b value is commonly used to investigate of peak strength in true triaxial testing. Figure 6(a) shows that, for both specimen types, φmax increases as the b value increases from 0.00 to 0.75 and then decreases as the b value increases from 0.75 to 1.00. Similar trend was also found by previous studies using laboratory tests [47, 49] and DEM simulations [4, 5, 11]. The peak value of the mean stress pmax that represents the overall stress level increases as the b value increases for both specimen types, as shown in Figure 6(b). Because all the numerical tests were conducted using the same σ3 value, the σ2 value increases as the b value increases and leads to the increase in the mean stress.Figure 6
Comparison of (a) peak internal friction angle (φmax) and (b) peak mean stress (pmax).
(a)(b)Compared to the isotropic specimens, theφmax and pmax values for anisotropic specimens are lower at the same b values and change less as the b value increases. This finding indicates that the horizontal alignment of the long axis of the particles reduces the effect of nonuniform stress distribution on the stress behavior of granular materials.To highlight the effects of the nonuniform distribution of external loading on the stress behavior of the granular materials, the stress ratio, defined asq/p, is introduced. Figure 7 shows that, for the isotropic specimens, the stress ratio deceases as the b value increases, which is consistent with related laboratory test results [50] and numerical simulation results [4, 14]. Also, the stress ratios for both specimen types are similar at the same b values. This finding indicates that the increased proportion of σ2 results in the decrease of the shear strength of granular materials, and this effect is similar for both isotropic and inherently anisotropic granular specimens.Figure 7
Evolution of stress ratio with the increase ofε1: (a) isotropic specimens; (b) anisotropic specimens.
(a)(b)Figure8 shows the variation of εv with the major principle strain ε1. All specimens exhibit the process from contraction to dilation as ε1 value increases. For the isotropic specimens, the amount of contraction tends to increase as the b value increases. By contrast, the anisotropic specimens are less contracted than the isotropic specimens for all the selected b values, and their contraction stages end with smaller ε1 values. This finding indicates that the inherently anisotropic granular materials are less compressible and may dilate earlier than isotropic granular materials and that the effect of b values on strain behavior is reduced.Figure 8
Variation of volumetric strain with the increase ofε1: (a) isotropic specimens; (b) anisotropic specimens.
(a)(b)


## 3.2. Microscopic Analysis: Contact Characteristics

In order to illustrate the macroscopic behavior of granular materials at the particle scale, the contact characteristics are analyzed. The coordination number (CN) is a commonly used index that is directly related to the stability of granular material [4, 51] and is defined as the average number of contacts per particle:(8)CN=2Nc−NcwNp,where CN is the coordination number, Nc is the total number of contacts, Ncw is the number of contacts between the particle and the wall, and Np is the total number of particles.Figure9 shows the evolution of the CN with the major principle strain. For both specimen types, the CN increases with ε1 to the peak value and then decreases to a residual value. The peak CNs increase obviously as the b value increases. Compared to the isotropic specimens, the inherently anisotropic specimens have higher initial CNs, lower peak CNs, and smaller ε1 values that correspond to the peak CNs. Besides, the increase of their peak CNs as b value increases is much smaller than isotropic specimens. This finding indicates that the oriented arrangement of particles in inherently anisotropic specimens form more contacts in the early stage of shearing, and the change of CNs during shearing is less affected by nonuniform stress distribution.Figure 9
Evolution of coordination number with the increase ofε1: (a) isotropic specimens; (b) anisotropic specimens.
(a)(b)During triaxial testing, the contact sliding typically occurs during the rearrangement of the particles, thus resulting in the permanent deformation of the granular materials. As defined by Gu et al. [39], contact slides when the inequality, shown as equation (9), is satisfied:(9)ftμfn>0.9999,where fn and ft are the normal contact force and tangential contact force, respectively, and μ is the contact friction coefficient.Figure10 presents a comparison of the initial and peak sliding contact percentages during triaxial shearing for all specimens. For both specimen types, the initial sliding contact percentages increase as b value increases from 0.00 to 0.25, and decrease as b value increases from 0.25 to 1.00. The peak sliding contact percentages for the isotropic specimens decrease as the b value increases, which is consistent with the DEM simulation results of Mahmud Sazzad et al. [4]. By contrast, the peak sliding contact percentages for the anisotropic specimens decrease as the b value increases from 0.00 to 0.75, then tend to increase as the b value increases from 0.75 to 1.00. Under the same b values, the anisotropic specimens have lower initial sliding contact percentages and higher peak sliding contact percentages compared to the isotropic specimens. This finding indicates that the sliding of contacts is restrained as the b values increase for isotropic specimens. By contrast, the anisotropic specimens restrain contact sliding in the early stages of testing compared with isotropic specimens, but contact sliding occurs more during shearing, especially with high b values.Figure 10
Comparison of (a) initial sliding contact percentages and (b) peak sliding contact percentages.
(a)(b)Both the CNs and the sliding contact percentages can be used to explain the change of volumetric strain (see Figure8). During the contraction stage, the particles are squeezed together in close contact with each other, which leads to an increase in the CN value. During the dilation stage, the particles start to rotate and slide, which results in the decrease of CN [51]. The specimen which is contracted mostly also has the largest peak CN and the lowest peak sliding contact percentage.


## 3.3. Microscopic Analysis: Fabric Tensor and Anisotropy

The directional distribution of particles, defined as the fabric, plays an important role in the shear behavior of a particle assembly [52–55]. To further investigate the combined effects of inherent and stress-induced anisotropy of granular materials, the contact normal, normal contact force, and tangential contact force are induced.The directional distribution of the contact normal is usually described by the fabric tensorRij [56, 57]:(10)Rij=∫ΩEΩninjdΩ=1N∑c∈Nninj,where ni is the unit contact normal in the i-direction, N is the total number of contacts, and EΩ is the distribution function on the unit sphere Ω and can be calculated by the second-order Fourier expansion as follows [58]:(11)EΩ=14π1+αijrninj,where αijr is the second-order anisotropic tensor and characterizes the fabric anisotropy. This tensor is determined by the deviatoric part of the fabric tensor Rij′ as(12)αijr=152Rij′.For normal and tangential contact force, the similar definition is induced, as expressed in equation (13)–(18):(13)Fijn=14π∫Ωf¯nΩninjdΩ=∑c∈NfnninjN1+αklrnknl,(14)f¯nΩ=f¯0n1+αijnninj,(15)αijn=152Fijn′f¯0n,(16)Fijt=14π∫Ωf¯itΩninjdΩ=∑c∈NftninjN1+αklrnknl,(17)f¯itΩ=f¯0nαijtnj−αkltnknlni,(18)αijt=5Fijt′f¯0n,where fn is the normal contact force, ft is the tangential contact force, and f¯0n is the average normal contact force over all the groups on Ω.Note thatαijr, αijn, and αijt are symmetric and deviatoric, so their second invariants αr, αn, and αt are used to quantify the degree of anisotropy as(19)α∗=32αij∗αij∗,where the super and subscript ∗ represents the contact normal r, tangential contact force t, and normal contact force n, respectively.Figure11 shows that the αr, αn, and αt values increase at the beginning and then decrease to a residual value as ε1 value increases. As the b value increases, the initial slope of αr, shows a small difference, whereas the peak values and residual values decrease obviously. This finding indicates that the contact distribution tends to be more uniform with a higher b value. Besides, the αn values are about 3 to 4 times higher than the αt values for all conditions, which suggests that contact force anisotropy is mainly caused by normal contact forces [51].Figure 11
Effects ofb values on the anisotropy parameter αr, αn, and αt. (a) αr Iso, (b) αr Ani, (c) αn Iso, (d) αn Ani, (e) αt Iso, and (f) αt Ani.
(a)(b)(c)(d)(e)(f)The initialαr values for the two specimen types are significantly different. The isotropic specimens have an initial αr at zero, which means the initial contact distribution is uniform. By contrast, the anisotropic specimens have an initial αr, at about 1.0, which means the oriented arrangement of particles make the initial contact distribution nonuniform. However, the peak and residual αr values for two specimen types are similar, which suggests that inherent anisotropy changes the initial contact distribution, but will not obviously affect the extent of contact normal anisotropy during shearing. For αn and αt, the values change with different b values for anisotropic specimens is smaller than the isotropic specimens, which suggests that inherent anisotropy significantly reduces the effect of b values on the distribution of contact force of granular materials.To address the reason for the difference in anisotropy parametersαr and αn, between the two specimen types, the principle values corresponding to them are induced. The αt is not analyzed here because it is much smaller than αn.The principle values of the fabric tensorR are used to quantify the cluster extent of contacts in the principle stress directions [14]. R1, R2, and R3 represent the contact clustering extent in the σ1, σ2, and σ3 direction, respectively. As shown in Figure 12, the R principle values show similar trends with the increase of ε1 for both specimen types. Before the peak of R principle values, R1 increases, whereas R3 decreases as ε1 value increases for all the selected b values. R2 decreases with the increase of ε1 when b < 0.5 and increases as ε1 value increases when b > 0.5. Such trends are similar to the precious DEM simulation results [11, 14, 22], which also indicates that contacts tend to cluster in the σ1 direction to carry most of the external load. As the b value increases, the contacts also cluster in the σ2 direction.Figure 12
Evolution of theR principle values with an increase in b values: (a) b = 0.00; (b) b = 0.25; (c) b = 0.50; (d) b = 0.75; (e) b = 1.00.
(a)(b)(c)(d)(e)The anisotropic specimens have higherR3 values and lower R2 values compared with the isotropic specimens for all selected b values; thus, the contact normal anisotropy αr for the two specimen types is similar. In addition, for the anisotropic specimens, the R2 value is lower than the R3 value under selected small b values, whereas for the isotropic specimens, the relationship is opposite. This finding indicates that the inherently anisotropic specimens form less contact in the direction of the long axis orientation of particles during shearing, compared with isotropic specimens.Figure13 shows the principle values of the normal contact force tensor Fn (Fn1, Fn2, Fn3). The trends are similar for the three Fn principle changes as the ε1 increases. As b value increases, all three Fn principle values increase obviously, which leads to the increase of peak internal friction angle and peak mean stress (Figure 6).Figure 13
Evolution of theFn principle values with the increase in b values: (a) b = 0.00; (b) b = 0.25; (c) b = 0.50; (d) b = 0.75; (e) b = 1.00; (f) comparison of peak Fn1 values.
(a)(b)(c)(d)(e)(f)Compared to the isotropic specimens, the anisotropic specimens have lower principleFn values (except for Fn2 when b = 0), and the growth rate of contact force for anisotropic specimens is lower as the b value increases, as shown in Figure 13(f). Therefore, the difference in αn for the anisotropic specimens among different b values is smaller (Figure 11), and their peak internal friction angle and peak mean stress are lower than isotropic specimens (Figure 6).For the anisotropic specimens, the clustering of the contact normal in theσ2 direction (R2) is lower than in the σ3 direction (R3) for the selected small b values (Figure 12), whereas the clustering of the contact force in the σ2 direction (Fn2) is higher than the σ3 direction (Fn3), which suggests that the oriented particles form strong contacts in their long axis direction and thus comprise of a strong force subnetwork [51, 56, 59, 60].


## 3.4. Discussion of the Combined Effects of Anisotropy

In summary, the combined effects of inherent anisotropy and stress-induced anisotropy can be illustrated in Figure14. For selected small b values, the long axis of the particles in the anisotropic specimens are initially perpendicular to the σ1 direction to carry most of the external load [21] and resist the sliding tendency in σ1 direction. Also, these particles form strong contacts in the σ2 direction, as explained in Section 3.3. Therefore, the anisotropic specimens form a support frame in the early stages of testing, which leads to higher initial CN values and lower initial sliding contact percentages and makes the specimens less compressible compared with isotropic specimens. As the strain value increases, the isotropic specimens are obviously compressed, whereas for the anisotropic specimens, particle rotating and sliding occur due to the increased deformation in the σ2 direction, which leads to higher peak sliding contact percentages, restrains the increase in the CNs, and shortens the contraction stage.Figure 14
Schematics of the combined effects of inherent and stress-induced anisotropy.The increase inb values increases the proportion of σ2, thus leads to a greater sliding tendency in σ2 direction. For isotropic specimens, the increase in b values leads to higher peak CNs, restrains the sliding of the contacts, and thus makes the specimens more compressible. For inherently anisotropic specimens, the increased b value significantly aggravates the rotating and sliding of the oriented particles and offsets its compressive effect. Therefore, the anisotropic specimens have higher peak sliding contact percentages and lower growth rate for the coordination number and contact force and are less compressible than the isotropic specimens. As a result, the effect of b values on the mechanical behavior of anisotropic specimens is significantly reduced when combined with the inherently anisotropic condition, compared with isotropic specimens.


## 4. Conclusion

This paper discussed the results of numerical simulations to investigate the combined effects of inherent and stress-induced anisotropy on the mechanical behavior of granular materials using 3D DEM. The macroscopic stress-strain relationships, microscopic contact characteristics, and anisotropic parameters were analyzed in detail and could be used to provide insights into the mechanism of anisotropy. The following conclusions can be drawn from this study:(i)
The macroscopic and microscopic responses of granular materials are significantly affected byb values. The increase in b values makes isotropic specimens more compressible and leads to lower shear strength at the same stress level, and the directional distribution of contact normal and contact force becomes more uniform.(ii)
Compared with isotropic specimens, the oriented granular particles in inherently anisotropic specimens form a support frame in the early stages of triaxial shearing, and thus make the granular material less compressible. As the strain increases, the oriented particles start to rotate and slide, which result in the redistribution of interparticle contacts and make the granular tend to dilate.(iii)
When applying on the inherently anisotropic specimens, increasedb values break the strong contacts parallel to the long axis of particles, aggravate the rotating and sliding of particles, and offset the compressive effect of b values. Thus, the contraction of granular materials and the increase of shear strength as the b value increases are significantly restrained. As a result, the inherent anisotropy reduces the effects of stress-induced anisotropy on the mechanical behavior of granular materials.



---

*Source: 7841824-2020-09-17.xml*