CFD-DEM model to assess stress-induced anisotropy in undrained granular material

Computers and Geotechnics xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

CFD-DEM model to assess stress-induced anisotropy in undrained granular material Talat Foroutan, Ali Asghar Mirghasemi

⁎

School of Civil Engineering, College of Engineering, University of Tehran, Iran

A R T I C LE I N FO

A B S T R A C T

Keywords: True triaxial test Granular soils Undrained shear strength CFD-DEM Stress-force-fabric

This paper presents CFD-DEM simulations of undrained true triaxial tests on a granular soil. A CFDEM solver that considers compressible fluid and moving meshes has been developed and calibrated against the results of laboratory testing on Nevada sand. The effect of the intermediate stress ratio (b) on the macro- and micro-mechanical features of undrained granular assemblies was studied. The tests were repeated at different confining pressures, and the trends were similar to each other. The lowest effective frictional resistance was observed for specimens for which b = 0 and the highest value was obtained for specimens with b = 0.5. An increase in b caused samples to exhibit the greatest dilative behavior and increased the absolute maximum negative excess pore pressure. The stress-force-fabric equation that represents the stress state of the assembly in terms of the fabric anisotropy and contact force anisotropy tensors was used to study the micromechanics of the undrained granular samples. The tests were varied from compression to extension in order to qualitatively and quantitatively analyze the evolution of anisotropy under anisotropic loading conditions. Observations on the microstructural evolution of undrained samples were in agreement with those already reported for drained samples.

1. Introduction The current research aimed to develop CFD-DEM simulations of true triaxial testing on a granular soil under undrained conditions. In this type of test, three different principal stresses are applied to a soil sample under conditions that are most likely to occur in reality. The results of true triaxial testing revealed that the stress path could influence the shear strength and deformation characteristics of granular soil [1–4]. Study of the effects of the stress path in drained soil using the discrete element method (DEM) has a strong background in the literature. Ng [5,6] investigated the behavior of ellipsoid granular material under different stress paths. Thornton [7,8] examined granular assemblies under different loading paths and investigated the contributions of normal and tangential contact forces to the evolution of deviator stress. Barreto and O'Sullivan [9] investigated the combined effect of the intermediate stress ratio of b = (σ2 − σ3)/(σ1 − σ3) and the coefficient of friction between particles on strong force chain formation and its resistance to collapse. O'Sullivan et al. [10] used DEM data and an analytical model to study the sensitivity of shear strength versus b by simulating the buckling of strong force chains. Huang et al. [11] studied the critical state response of granular assemblies under a generalized stress condition and observed that the critical state line depends on b. Zhou et al. [12,13] investigated the influence of the intermediate stress ⁎

ratio of granular samples under different anisotropic consolidation conditions. They studied the influence of the intermediate stress ratio and initial stress on the critical state and stress-dilatancy relationship of granular material. Xiao et al. [14] carried out a series of true triaxial tests to investigate the critical state of coarse-grained soil. Sazzad and Suzuki [15] studied the behavior of loose and dense granular material at different intermediate stress ratios. Dorostkar and Mirghasemi [16,17] investigated the micromechanics of a granular sample using stress-force-fabric (SFF) relationships at different intermediate stress ratios. On the other hand, little research has focused on true triaxial testing under undrained conditions. Huang et al. [11] used the constant volume method to simulate undrained triaxial compression (b = 0) and extension (b = 1) tests. Yimsiri and Soga [18] investigated the effects of the initial fabric of the soil, including specimen reconstitution methods, large pre-shearing, and anisotropic characteristics in undrained triaxial compression and extension tests. In their studies, information about the midrange of b was unavailable. Studies that have been done to investigate the undrained shear strength of granular soil using DEM fall into two main groups. In the first group, the constant volume approach is applied, in which the fluid phase is not explicitly considered. Sitharam et al. [19–21] used this method to investigate the micromechanical behavior of loose and dense

Corresponding author. E-mail addresses: [email protected] (T. Foroutan), [email protected] (A.A. Mirghasemi).

https://doi.org/10.1016/j.compgeo.2019.103318 Received 16 May 2019; Received in revised form 15 October 2019; Accepted 16 October 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Talat Foroutan and Ali Asghar Mirghasemi, Computers and Geotechnics, https://doi.org/10.1016/j.compgeo.2019.103318

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

and is termed cfdemSolverSonicLiquidDyM. The model was calibrated using experimental results from drained and undrained triaxial tests on Nevada sand [44]. Several specimens were used, with intermediate stress ratios of 0 to 1 at different confining pressures. The shear strength of granular assemblies in terms of the effective mobilized friction angle and pore pressure evolution with respect to b has been depicted for undrained samples. The SFF equation presented by Chantawarangul [45] was used to monitor the evolution of the anisotropy coefficients and the contribution of each anisotropy coefficient to the mobilized shear strength of the samples in the undrained condition.

undrained 3D granular material. They also investigated soil liquefaction under cyclic load and the critical state behavior of undrained granular assemblies. Gong et al. [22,23] studied simulations of axisymmetric undrained tests on loose assemblies. Zhao and Guo [24] investigated the rotational resistance effect among spherical particles in undrained triaxial tests. Liu et al. [25] presented the influence of grading on the undrained behavior of granular material. Nguyen et al. [26] studied the micromechanical and macro-mechanical features of 3D granular assemblies of ellipsoid particles at characteristic and phase transformation states. In the second group, the fluid phase is computationally modeled using Darcy's equation. Among researches that have extended DEM and computationally modeled particle fluid interaction, the pioneer work is by Hakuno and Tarumi [27], who for the first time modified and developed DEM using Darcy's law to include pore fluid effects; pore scale meshes were used to solve fluid flow equations. Bonilla [28] used Darcy's equation to simulate undrained biaxial tests of ellipsoid particle assemblies, and fluid flow was investigated in the pore scale. Shafipour and Soroush [29] and Khalili and Mahboubi [30] used a coarse grid technique in studying 2D saturated granular materials by Darcy's equation. Okada and Ochiai [31] simulated the liquefaction phenomenon in 3D granular soil assemblies. In their study, a measurement sphere surrounded each particle and fluid flow occurred among these spheres. Liu et al. [32] improved Okada and Ochiai's approach and used it in studying undrained triaxial tests. An overview of this literature indicates that great efforts have been made to study the undrained shear strength of granular soil. However, investigation into the effect of intermediate principal stress on undrained tests has been neglected. The focus of the current study is to develop an undrained true triaxial test to study the intermediate stress ratio effect. As mentioned earlier, the undrained shear strength of the granular soils is mostly studied via the constant volume approach. In studies that modeled the fluid phase explicitly, Darcy's empirical equation has been used because of its simplicity. Because of shortcomings in the development of an ideal model, the present study developed a coupled method to simulate undrained true triaxial testing. This coupled method uses Newton's second law of motion to solve the particulate phase and the Navier-Stokes momentum equation and the continuity equation for solving the fluid phase. To solve the fluid phase using the Navier-Stokes equations in the presence of particles, the most popular method is the coarse grid technique first used by Tsuji [33] to model fluidized beds. This approach was employed in the current study. In the recent past, the interest in using the CFD-DEM approach to solve geomechanical problems has increased. Zeghal and El Shamy [34,35] have implemented this approach in geomechanics and investigated the liquefaction of loose and cemented granular deposits. Climent et al. [36] studied sand production in oil wells using a coupled method. Jiang et al. [37] investigated the micromechanical properties of saturated sand under cyclic load. Zhang et al. [38] modeled seepage erosion around shield tunnels for loose, medium dense and dense silty sands. Guo and Yu [39] used different mathematical CFD-DEM models to simulate the surface erosion of granular soil inside a pipe flow. To the author’s knowledge, no work has been carried out to use this approach to investigate undrained triaxial tests. In this study, the CFD-DEM approach is employed to study the undrained shear strength of a granular soil in triaxial tests. With this aim, LIGGGHTS [40], an open source DEM code on the LAMMMPS [41] platform, and open-source software package OpenFOAM [42] were used. The coupling framework of CFD-DEM was realized by CFDEMcoupling [43]. Both codes are executed in parallel and the communication among processors and data exchange between CFD and DEM are done by Message Passing Interface functionality. The basic CFD-DEM solver cfdemSolverPiso does not support compressible fluid and dynamic meshes. For this reason, a new CFD-DEM solver that considers compressible fluid and dynamic meshes has been developed

2. Methodology and formulation 2.1. DEM Particle interaction is solved using DEM as proposed by Cundall and Strack [46] in the following governing equations for particle i motion:

mi

Ii

dvi = fpf , i + dt

dωi = dt

nc

∑ j=1 fc,ij

+ mi g

(1)

nc

∑ j=1 (Mij)

(2)

where vi and ωi are, respectively, the translational and angular velocities of particle i , nc is the number of particles in contact with particle i , fpf , i is the particle-fluid interaction forces, fc, ij is the particle-particle contact forces, including force induced by contact law and damping force, mi g is the gravitational force and is ignored in this study, and Mij is the torque acting on particle i through other particles. According to Eq. (2), it is clear that fluid does not contribute torque to particles. Parameters mi and Ii are the mass and moment of inertia of particle i , respectively. Different contact laws are implemented in LIGGGHTS. In this study, the nonlinear Hertz contact law was used. Normal and tangential contact forces can be calculated using Eqs. (3) and (4), respectively. Eqs. (5) and (6) are the normal and tangential damping coefficients, respectively.

Fn =

4 ∗ ∗ 3/2 E R δn − γn vn 3

(3)

Ft = 8G∗ R∗δn δt − γt vt γn = −2

5 6

ln2 (e ) + π 2

ln (e )

γt = −2

5 6

ln2 (e ) + π 2

ln (e )

(4)

2m∗E ∗ R∗δn

8m∗G∗ R∗δn

(5)

(6)

where Fn is the normal contact force, Ft is the tangential contact force, δn is the overlap distance between two particles, δt is the tangential displacement vector between two spherical particles which is truncated to satisfy a frictional yield criterion, vn and vt are the normal and the tangential component of the relative velocity of the two particles, respectively, γn and γt are the normal and tangential damping coefficients, respectively and R∗, E ∗, G∗ and m∗ are the equivalent radius, equivalent Young's modulus, equivalent shear modulus and equivalent mass of two particles in contact, respectively, and e is the coefficient of restitution. The Coulomb friction criterion is implemented for the tangential force:

|Ft | ≤ μs |Fn|

(7)

where μs is the coefficient of sliding friction. To take into account the particle shape, a rolling contact model was implemented that transfers a moment through contact [47]. This consists of a spring moment (Eq. (9)) and a damping moment (Eq. (13)).

Mr = Mrk + Mrd 2

(8)

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

Initial configuration

Deformed configuration

Fig. 1. Schematic diagram of deforming CFD cells, including DEM particles.

ΔMrk = −kr Δθr

kr = βkn R∗2 kn =

(10)

and support mesh deformation. The new solver is denoted cfdemSolverSonicLiquidDym. The governing equations for this solver based on the volume averaged Navier-Stokes equations are: Continuity equation:

(11)

∂ (ρf αf )

(9)

4 ∗ ∗ E R δn 3

∂t

ΔMrk

+ ∇ . (ρf α u f ) = 0 f

(16)

is the increase in spring moment during Δt , Δθr is the rewhere lative rotation increment, kr is the rolling stiffness, β is the rolling stiffness coefficient and kn is the elastic constant for normal contact. The spring moment is limited as:

where αf is the fluid volume fraction, u f is the fluid velocity and ρf is the fluid density. Momentum equation:

|Mrk| ≤ μr R∗ |Fn|

∂ (αf ρf u f )

(12)

∂t

where μr is the rolling friction coefficient. The damping moment Mrd is realized as:

Mrd =

k ∗ ⎧ − Cr θṙ if Mr < μr R Fn ⎨− fCr θṙ if Mrk = μr R∗Fn ⎩

Cr = 2ηr Ir kr

(17)

where p is the fluid pressure, T is the fluid stress tensor and rpf is the fluid-particle interaction force. Here, the density of the fluid is a function of pressure and is defined as: (13)

∂ρf

(14)

∂t

= Cf

∂p ∂t

(18)

where Cf is the compressibility. Because the continuity and momentum equations are coupled, there is a need for special treatment to solve these equations. In this case, the PISO algorithm, a predictor-corrector method, has been applied. In the predictor step, the momentum equation is solved implicitly, pressure and density are used from the previous time step and a first approximation for the velocity field is obtained. By solving the pressure equation, a new pressure field is estimated. Because of the new pressure distribution, the velocity field should be updated. The velocity is corrected in an explicit manner; this is the corrector step [49]. In order to simulate the undrained triaxial test, it is necessary for the solver to support the movement of boundaries. The ability to squeeze and stretch mesh during the run is possible in OpenFOAM. The existing dynamic mesh library should be added to the new solver. To provide a realistic pressure distribution during shear or to capture the decrease and increase in pressure during the soil softening and hardening phases accurately, the CFD boundaries should be synchronized with the DEM boundaries. A mapping technique was employed to update the new node position in the CFD mesh with the current position of the particles. The CFD boundaries move in such a way that they are always tangential

−1

1 1 ⎞ Ir = ⎜⎛ + 2 2⎟ I + m R I + m R i i j j i j ⎝ ⎠

+ ∇ . (αf ρf u f u f ) = −αf ∇p + ∇ . (αf T ) − rpf

(15)

where f is a non-dimensional coefficient that represents the effects of rolling damping after full mobilization, θṙ is the relative rotational velocity of two particles in contact, Cr is the rolling viscous damping coefficient, ηr is the rolling damping ratio, and Ir is the equivalent moment of inertia. In this study f = 0 , which implies that in the case of full mobilization (Mrk = μr R∗Fn ) , rolling damping will be disabled; in other words, when the spring moment reaches its maximum value according to Eq. (12), the damping moment vanishes. 2.2. CFD Goren [48] found that pore pressure evolution in granular material is dependent on system drainage. When drainage conditions are poor, pore pressure variation is highly influenced by compressibility and volumetric strain. In the present study, the basic solver cfdemSolverPiso from the CFDEM package has been modified to consider compressibility 3

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

3. Calibration of model

to the outer marginal particles in the assembly. Fig. 1 shows a schematic 2D diagram of the CFD cells deforming according to particle motion. D represents a domain configuration at the given time of t , and during a time interval of dt , D changes shape to Dn . The x component of the CFD node positions for each time step is updated as:

Pxn = xcn +

l xn (Px − x c ) lx

Considering the resemblance of the results obtained from conventional triaxial compression tests to the results of true triaxial tests, where the minor principal stress and intermediate principal stress are equal (b = 0), the model was calibrated using the results of a series of drained and undrained triaxial tests done on Nevada sand [44]. To keep the computational effort to a minimum, the particle size distribution of Nevada sand was upscaled ten-fold [54–56]. A total of 18,000 particles with diameters ranging from 0.75 to 2.5 mm were generated. Fig. 3 shows the particle size distribution curve of the assembly and the Nevada sand. The ratio of the side length of the cubical box to mean particle diameter d50 was set to 25, which is a reasonable value to eliminate the size effect [57]. The particle density was set to 2670 kg/m3 , the same as for Nevada sand. The time step has been identified in accordance with the critical time step [58]; the relevant values for the DEM time step are shown in Table 1. Six rigid walls are used to model the cubical test cell. When particles are generated randomly in a loose state, consolidation is done by moving opposing walls at the same rate in opposing directions until the intended confining stress and the void ratio are achieved. A view of the compacted assembly at a confining pressure of 160 kPa is shown in Fig. 4. Under shear, a constant stress value was assigned to direction 3 to provide minor principal stress σ3 and a vertical loading system implemented at a constant strain rate to prepare major principal stress σ1, a horizontal loading system specified to produce intermediate principal stress σ2 , in such a way as to keep b constant during shear. This mechanism of controlling stresses in the true triaxial test is known as the mixed boundary condition [59]. For the quasi-static simulation, the unbalanced forces index should be kept under 0.01. The results for this value from the current study are shown in Fig. 5 for specimens at a confining pressure of 160 kPa and b = 0, 0.5 and 1. It is clear that, for confining pressures above this value, the samples are more stable. The unbalanced forces index is calculated as [60]:

(19)

n n − x min )/2 , l x = (x max − x min ) where x c = (x max − x min )/2, xcn = (x max n n n and l x = (x max − x min ) , x min and x max are the minimum and maximum x component of the particle’s position in the particle assembly, Px is the x component of a node position and n denotes the next time steps. The y and z components of the node position can be calculated in the same way.

2.3. Fluid-particle interaction The particle phase is fully saturated and surrounded by fluid; the drag and pressure gradient forces are considered to be exerted on the particles. The other components of hydrodynamic forces, such as lift force and virtual mass force, are ignored as is common in geomechanical problems [50]. For numerical reasons, the momentum exchange term has an implicit and an explicit term [40,43]. The last term in Eq. (17) considers the particle force on the flow field and, in accordance with the formulation of momentum, is defined as: n

rpf =

∑i = 1 (fd, i ) (20)

ΔV

where

fd, i =

ΔVβpf αp

(u − v )

(21)

where αp = 1 − αf , fd, i is the drag force on particle i , n is the number of particles in the CFD cell, ΔV is the volume of the computational cell, u is the fluid velocity and v is the average particle velocity in the CFD cell. βpf is the drag force correlation. For purposes of evaluation of the drag force, there are a number of correlations based on experimental and numerical studies which have been summarized by Zhu et al. [51]. Gidaspow [52] presented a general correlation for drag force estimation in fluid-particle systems. His correlation for dense systems is the same as that in Ergun’s equation [53] and was used in this study as:

βpf =

150(1 − αf )2μf αf dp2

+

1. 75(1 − αf ) ρf dp

|u − v|

αf ≤ 0. 8

Iuf =

1 Nc

N

p 2 ) ∑p =p 1 (fres N

∑c =c 1 (f c )2

≤ 0.01 (24)

where Nc and Np are the total number of contacts and particles, rep spectively, f c is the contact force and fres is the resultant force on each particle. The required DEM micro-parameters were E (Young's modulus), μs (coefficient of sliding friction), μr (coefficient of rolling friction), β (rolling stiffness coefficient), ηr (rolling damping ratio), ν (Poisson's ratio) and e (coefficient of restitution). First, the experimental results of the drained triaxial test were used to identify the DEM micro-parameters. Next, the experimental results of the undrained triaxial tests were used to calibrate the characteristics of the coupled method. Fig. 6 shows diagrams for deviator stress and volumetric strain versus axial strain for an isotropically consolidated drained compression test. During shear, the mean confining stress (p' ) is held constant at 160kPa and the initial void ratio is set at 0.65. As seen, the results of DEM generally agree with the experimental results. Table 1 presents the finalized micro-parameters. The coefficient of friction between the walls and particles is assumed to be zero, thus, the stresses applied to each wall will remain normal. In accordance with the averaging process used to calculate the particle properties for the coupled unresolved CFD–DEM method (flow field among particles is not resolved and the impact of particles is averaged in one cell), the CFD cell size should be large enough to contain a reasonable number of particles [40,43]. Sensitivity analysis was carried out to determine if there were enough particles in each cell. The fluid phase was divided to 16 × 16 × 16 cells and each cell contains 4–5 particles. The CFD time step was set at five times the DEM time step

(22)

where μf is fluid viscosity, dp is particle diameter. The particle-fluid interaction force implemented in Eq. (1) is defined as:

fpf , i = fd, i + f∇p, i

1 Np

(23)

where f∇p, i is the pressure gradient force on particle i . In the coupling scheme, the CFD and DEM solvers run simultaneously and data exchange between them occurs at specific time intervals. In this iterative loop, the particle positions and velocities are calculated by LIGGGHTS and are transferred to the CFD solver. Then the corresponding CFD cell for every particle is determined and the mean particle velocity and void fraction for the cells are assigned. Using this information, the momentum exchange term is calculated and fluid flow is evaluated, then the fluid forces on the particles are determined and returned to the DEM solver. Fig. 2 is a flowchart showing the coupling procedure. 4

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

Force displacement law

Synchronize CFD mesh boundaries

Contact forces

Void fraction and mean particle velocity in each cell is determined

Equations of motion Fluid pressure and velocity in each cell is calculated Particle positions and velocities

Fluid forces on particles is calculated

Fig. 2. The algorithm of adopted continuum-discrete model.

100 90 Passing percent (%)

80 70 60 50

Experiment

40

DEM

30 20 10 0 0.001

0.01

0.1 1 10 Particle diameter (mm)

100

1000

Fig. 4. Simulation cell for true triaxial test at confining pressure of 160 kPa.

0.004

Fig. 3. Particle size distribution.

b=0 b=0.5 b=1

0.003

Parameter

Value

Young's modulus (GPa) Poisson's ratio Coefficient of restitution Coefficient of friction between particles Coefficient of rolling friction between particles Coefficient of rolling stiffness Rolling damping ratio Fluid compressibility (sec2/m2) Fluid density (kg/m3) Fluid viscosity (Pa sec) DEM time step (sec) CFD time step (sec)

40 0.22 0.87 0.6 0.35 0.3 0.13 4.4e−07 1000 1.5e−3 1e−7 5e−7

Iuf

Table 1 Parameters used in simulations.

0.002

0.001

0 0

5

10 Axial strain %

15

20

Fig. 5. Unbalanced forces index versus axial strain at confining pressure of 160 kPa and b = 0, 0.5 and 1.

initial void ratio of 0.65 and confining pressure of 80 kPa are presented in Fig. 8. There is good agreement between the experimental results and those obtained from the CFD-DEM approach. The proposed method is used to study the evolution of anisotropy during undrained true triaxial testing in the generalized principal stress condition. To that end, granular assemblies are generated and sheared at initial confining pressures of 160, 300 and 600 kPa with the same initial void ratio of 0.65, and an intermediate stress ratio (b) of 0, 0.25, 0.5, 0.75 and 1.

because of the need to sync the CFD boundaries to the DEM, and the coupling time interval was selected to be the lowest value that is equal to the CFD time step. The isotropically consolidated undrained triaxial test results on Nevada sand are presented numerically and experimentally in Fig. 7. The effective consolidation pressure is 160kPa and the initial void ratio is 0.65. The required fluid parameters including Cf (compressibility), ρf (density) and μf (viscosity) are presented in Table 1. To verify the adopted parameters, the deviatoric stress-strain response and excess pore pressure evolution of another isotropically consolidated undrained triaxial test on the same Nevada sand with an 5

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

300

0

b)

-1

250 Volumetric strain %

Deviatoric stress (kPa)

a)

200 150 100 Experiment

50

0

5

10

15

5

10

20

-3 -4 -5 -6 Experiment DEM

-8

20

15

-2

-7

DEM

0

0

Axial strain %

Axial strain%

Fig. 6. Comparison of experimental and numerical results of drained triaxial test: (a) deviator stress-strain response; (b) volumetric strain versus axial strain.

4. Discussion of results

higher confining pressure had a lower effective frictional resistance. These observations are qualitatively similar to what has already been reported for drained samples [3,17,61,62]. In Fig. 10(a), 10(b), and 10(c), with an increase in the value of b, the undrained samples tended to show the most dilative behavior. The maximum primary positive development of pore pressure and its rate of development increased as b increased. The samples with higher b values decreased more quickly to the negative regime. The absolute value of maximum negative excess pore pressure increased when transferred from the compression test (b = 0) to the extension test (b = 1). Variation of the maximum negative excess pore pressure with the b value at different confining pressures is shown in Fig. 12. Similar trends were obtained for different confining pressures. At a given b value, the sample with the higher confining pressure had a higher absolute value of maximum negative excess pore pressure. When b = 0, there is a small difference among the evolved maximum negative excess pore pressures for different confining pressures. As the b value increased, the difference among the values of maximum negative excess pore pressures increased. These observations highlight the stress path effect on the strength characteristics of the undrained soil samples.

4.1. Macro parameters

3000 2500

4.2. Micro parameters In granular material, deformation of the system leads to considerable change in the magnitude and orientation of the contact forces. The contact force magnitudes and contact orientations vary randomly from one contact to another, but, by averaging contacts over a class having the same orientation, regular behavior can be seen [63]. Such averaging

b) Excess pore water pressure (kPa)

a) Deviator stress (kPa)

The sine of the effective mobilized friction angle and the excess pore pressure evolution versus axial strain with different b values are shown in Figs. 9 and 10, respectively. These results are presented for simulated tests under confining pressures of 160, 300, and 600 kPa. As can be seen in Fig. 9(a), 9(b), and 9(c), by implementing the axial strain in the major principal direction, the effective mobilized frictional resistance began to increase and attained its maximum magnitude around a strain of 6% for compression tests (b = 0). At b values other than 0, the sine of the effective mobilized friction angle reached its maximum value at lower strains. Compression tests show a smooth peak, but for b values greater than 0 after the peak point, the strength dropped somewhat and then reached the steady state. At residual state, the sine of the effective mobilized friction angle for all tests with different b values is almost the same. Loading in the intermediate direction increases the frictional resistance of the assembly, although for small intermediate stress ratios such as 0.25 there is a significant increase. The maximum values of the sine of effective mobilized friction angle for b values of 0.25, 0.5, 0.75 and 1 are comparable. To improve the comparison, the maximum ef' fective mobilized friction angle φmax versus b is presented in Fig. 11. At all confining pressures, the maximum effective friction angle increased from b = 0 to b = 0.5 and then started to decrease at a smaller rate to b = 1; thus, the lowest friction angle occurred at b = 0 and the highest value was obtained at b = 0.5. At a given b value, the sample with a

2000 1500 1000 Experiment CFD-DEM

500 0 0

5

10 15 Axial strain %

20

100 0 -100 0

5

10

15

20

-200 -300 -400 -500

Experiment CFD-DEM

-600 -700 -800 Axial strain %

Fig. 7. Comparison of experimental and numerical results of undrained triaxial testing at confining pressure of 160 kPa: (a) deviator stress-strain response; (b) excess pore pressure versus axial strain. 6

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

b) Excess pore water pressure (kPa)

Deviatoric stress (kPa)

a) 3000 2500 2000 1500 1000 Experiment CFD-DEM

500 0 0

5

10 15 Axial strain %

20

100 0 -100

0

5

10

15

20

-200 -300

Experiment CFD-DEM

-400 -500 -600 -700 -800

Axial strain %

Fig. 8. Comparison of experimental and numerical results of undrained triaxial test at confining pressure of 80 kPa: (a) deviator stress-strain response; (b) excess pore pressure versus axial strain.

a)

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

b=0 b=0.25 b=0.5 b=0.75 b=1

0

5 10 15 Axial strain (%)

20

b)

b=0 b=0.25 b=0.5 b=0.75 b=1

0.4 0.3 0.2 0.1 0

Sine of effective mobilized friction angle

0

5 10 15 Axial strain (%)

10

15

20

-700 -1200

b=0 b=0.25 b=0.5 b=0.75 b=1

-1700

-200 0

5

20

10

15

20 b=0 b=0.25 b=0.5 b=0.75 b=1

-700 -1200 -1700 -2200 Axial strain%

c)

b=0 b=0.25 b=0.5 b=0.75 b=1

5 10 15 Axial strain (%)

5

20

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

-200 0

300

Excess pore pressure (kPa)

0.6

Excess pore pressure (kPa)

Sine of effective mobilized friction angle

0.7 0.5

300

Axial strain%

0.8

b)

c)

Excess pore pressure (kPa)

Sine of effective mobilized friction angle

a)

300 -200 0

5

10

15

-700

20 b=0 b=0.25 b=0.5 b=0.75 b=1

-1200 -1700 -2200 -2700 Axial strain%

Fig. 10. Evolution of excess pore pressure versus axial strain in undrained true triaxial test for initial confining pressures of: (a) 160 kPa; (b) 300 kPa and (c) 600 kPa.

Fig. 9. Sine of effective mobilized friction angle versus axial strain in undrained true triaxial test for initial confining pressures of: (a) 160 kPa; (b) 300 kPa; (c) 600 kPa.

7

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

Maximum effective mobilized friction angle

48

1

r ⎧ E (Ω) = 4π (1 + aij ni nj ) ⎨ aijr = ajir aiir = 0 ⎩

46 44

ni and nj are components of contact normal and are defined based on the azimuthal φ and polar θ angles of the spherical coordinate system as Eq. (26).

42 40

n1 = cosθ n2 = cosφsinθ n3 = sinφsinθ

160 kPa 300kPa 600kPa

38 36 0

0.25

0.5

is the fabric anisotropy tensor which can be derived from the deviatoric fabric tensor as Eq. (27). Rij is the fabric tensor of the assembly and has the information about the geometrical arrangements of particles.

0.75

1

b

aijr = 7.5Rij'

Fig. 11. Maximum effective mobilized frictional angle vs. b in undrained true triaxial test under confining pressures of 160, 300 and 600 kPa.

1 Nc

Rij =

0

(27)

∑ nic njc

(28)

c∈V

Rij'

where = Rij − δij /3, δij is the unit matrix, and Nc is the total number of contacts. Similar to contact normal orientation, the distribution of average normal contact forces f¯ n (Ω) could be approximated by a second-order Fourier series (Eq. (29)). This distribution function specifies the average magnitude of normal contact forces among contacts with the same orientation.

160 kPa 300kPa 600kPa

-500

(26)

aijr

34

Maximum negative excess pore pressure

(25)

-1000 -1500 -2000

n n n ⎧ f¯ (Ω) = f¯o (1 + aij ni nj ) n n n ⎨ aij = aji , aii = 0 ⎩

-2500

n where f¯o is the average normal contact force in the assembly, which is defined as Eq. (30).

0

0.25

0.5 b

0.75

1

1 n f¯o = 4π

Fig. 12. Maximum generated negative pore pressure vs. b in undrained true triaxial test under confining pressures of 160, 300 and 600 kPa.

∑Ng

(29)

f¯ n (Ω g )ΔΩ g

(30)

¯n

f (Ω g ) is the average normal force component acting at contacts with an orientation of Ω g . Ω g is a partitioned group of the unit sphere, Ng is for the discretized form. aijn is the normal contact force anisotropy tensor which is related to the average normal contact force tensor Fijn . aijn = 7.5

Fij'n n f¯

(31)

o

Fijn =

1 4π

∑Ng

Fij'n

f¯ n (Ω g ) ni nj ΔΩ g

Fijn

(32)

n f¯o /3.

where = − The distribution function of average tangential contact forces can also be represented by a second-order Fourier series: t n t t ⎧ f¯i (Ω) = f¯o (aij nj − (akl nk nl ) ni )

⎨ ⎩

f¯ t (Ω) =

t t (f¯i × f¯i )

(33)

aijt

is defined in Eq. The tangential contact force anisotropy tensor (34), which is related to the average tangential contact force tensor Fijt .

Fig. 13. Spherical coordinate system.

aijt = 5

can be accomplished on a unit spherical coordinate system that is divided into fragments at intervals of ΔΩ [45]. The unit spherical coordinate system and the intervals of ΔΩ are shown in Fig. 13. The angular distribution of E (Ω) is the distribution function of contact normal orientation and specifies the portion of contacts having the orientation of ΔΩ. An acceptable approximation of E (Ω) may be provided by a second-order Fourier series, which has been described by Eq. (25) [45].

Fijt n f¯

(34)

o

Fijt =

1 4π

∑Ng

t f¯i (Ω g ) nj ΔΩ g

(35)

t f¯i (Ω g ) is the average tangential force component acting at contacts with orientation of Ω g . Figs. 14, 15 and 16, respectively, show the distributions of contact normals (E (Ω) ), average normal contact forces ( f¯ n (Ω) ), and tangential contact forces ( f¯ t (Ω) ). As mentioned earlier, E (Ω) specifies the portion

8

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

b)

c)

b)

c)

b)

c)

b=1

b=0.5

b=0

a)

Fig. 14. Evolution of 3D histograms of contact normal at b = 0, 0.5 and 1 for a confining pressure of 160 kPa: (a) initial state; (b) peak state; (c) residual state.

of contacts in the assembly having the orientation of ΔΩ. f¯ n (Ω) specifies the average magnitude of normal contact forces among contacts having the orientation of ΔΩ. And f¯ t (Ω) specifies the average magnitude of tangential contact forces among contacts having the orientation of ΔΩ. In this study, the unit sphere at the spherical coordinate system was divided into 36(azimuthal) × 18(polar) fractions, and every group of contacts is averaged over one fraction. Histograms of Figs. 14, 15 and 16 are presented at b = 0, 0.5 and 1 for the initial, peak and residual states and a confining pressure of 160 kPa. In each figure, the image in panel (a) is just after isotropic compaction and before implementation of deviatoric stress, which is the same for all b values. The images in panel (b) are at peak shear strength and the images in panel (c) are at the residual states. It is clear that under other confining pressures, these histograms are qualitatively similar. At the initial hydrostatic stress condition, for all b values, the distributions of contact normals (Fig. 14(a)) and average normal contact forces (Fig. 15(a)) are isotropic and spherical. The average tangential contact force (Fig. 16(a)) has a magnitude of zero under the hydrostatic stress condition. Implementing deviatoric stress increased the contact normal and contact force anisotropies. At b = 0, the loading condition is axisymmetric, i.e., there is one axis of loading and the sample expands in the two lateral directions

(σ1 > σ2 = σ3) . At b = 0.5, the loading condition is general (σ1 > σ2 > σ3) , and the histograms will no longer be axisymmetric. At b = 1, the loading condition is also axisymmetric, there are two major principal axes of loading (σ1 = σ2 > σ3), and the sample expands in one minor direction. When a granular assembly is under loading, the contact orientations evolve parallel to the axis of loading and loss of contacts occurs at minor stress directions. Therefore, for different stress paths, the histograms of contact orientations and the average contact forces may be different. Fig. 14(b) shows histograms of contact normals at peak state. For b = 0, the histogram evolves into an axisymmetric peanut shape with a principal axis that is coincident with the axis of loading. For b = 0.5, the histogram evolves from a perfect peanut shape to a bi-directional peanut shape. For this bi-directional peanut, a major principal axis is coincident with the major axis of loading and a second principal axis is in the direction of the intermediate axis of loading. For b = 1, the histogram is axisymmetric and has become donut-shaped without the central hole. The axis of symmetry for this donut shape is the minor principal stress direction. Fig. 15(b) shows the average normal contact force distribution at peak state. The shape of the evolution of the average normal contact forces histograms is similar to the histograms of the contact normals. For b = 0 the average normal contact force distribution evolves to a 9

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

b)

c)

b)

c)

b=0.5

b=0

a)

b)

b=1

c)

Fig. 15. Evolution of 3D histograms of average normal contact forces at b = 0, 0.5 and 1 for a confining pressure of 160 kPa: (a) initial state; (b) peak state; (c) residual state.

evolves in the 45° direction with respect to the major principal stress plane for all stress paths. These observations are qualitatively similar to what has been reported for drained true triaxial tests [45]. For quantitative analysis of the microstructural anisotropy of the samples, scalar measures of the fabric anisotropy tensor and contact force anisotropy tensors are used. The deviator coefficients of the anisotropy tensors are defined as:

perfect peanut. For b = 0.5 it evolves to a bi-directional peanut and for b = 1 the histogram becomes like a donut without the central hole. Fig. 16 (b) shows the tangential contact force distribution at peak state. This histogram, for b = 0, evolves to an axisymmetric dumbbell shape with an axis of symmetry that is coincident with the major principal stress axis. The tangential contact force distribution for b = 0.5 is not axisymmetric, but has evolved differently in different directions. For b = 1, the tangential contact force distribution is an axisymmetric dumbbell shape. The axis of symmetry for this dumbbell shape is the minor principal stress direction. In the residual states in panel (c), the histograms maintain the same shapes as in the peak state, but the size of the histograms decreases somewhat. The difference in histogram size between the peak and residual states is more pronounced for the average tangential contact forces. For all states and all values of b, the magnitude of normal contact force anisotropy becomes greater than the magnitude of the fabric and tangential contact force anisotropies. The lower magnitude of anisotropy relates to tangential contact force anisotropy. For all stress paths, the coincidence of the direction of anisotropy with the major principal stress direction can be seen. The maximum tangential contact force

adr =

3 r r aij aij 2

(36)

adn =

3 n n aij aij 2

(37)

adt =

3 t t aij aij 2

adr ,

adn

(38)

adt

and are the deviatoric coefficients of anisotropy tensors where of the fabric, normal and tangential contact forces, respectively. The change in these coefficients with axial strain at different values of b and a confining pressure of 160 kPa is shown in Fig. 17. These coefficients are representative of the shearing capacity of the system and their evolutionary trend is very similar to the assembly deviatoric stress10

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

b)

c)

b)

c)

b)

c)

b=1

b=0.5

b=0

a)

Fig. 16. Evolution of 3D histograms of average tangential contact forces with b = 0, 0.5 and 1 for a confining pressure of 160 kPa: (a) initial state; (b) peak state; (c) residual state.

magnitudes of the deviatoric coefficients of normal and tangential contact force anisotropies are the same at large strains for different b values with the exception of b = 0. Despite the deviatoric coefficient of fabric anisotropy, the deviatoric coefficient of (normal and tangential) contact force anisotropy decreased with an increase in b. It can be deduced that the coefficients of fabric anisotropy and contact force anisotropy are complementary and preserve the stability of the system. When the contact force anisotropy increases, the magnitude of the contact forces decreases in the minor principal direction. This indicates a need for more contacts in the minor direction to compensate for this. Quantitative comparison of Fig. 17 suggests that normal contact force anisotropy contributes most to the evolution of shear strength, and tangential contact force anisotropy contributes the least. These observations are qualitatively similar to previous reports for drained true triaxial tests [17,45]. Several expressions can be used to formulate the stress tensor of a particulate system subject to external loading at the boundaries. For an assembly in static equilibrium, the basic representation of the stress tensor has been presented by Rothenburg [64]. This form of stress tensor is the volume average of the micromechanical parameters of contact force fic and position vector l jc of the particle-to-particle contact point relative to the particle center in volume V . This form of stress

strain curve evolution. Fig. 17(a) shows that the deviatoric coefficient of fabric anisotropy increased gradually to a strain of about 5% and then decreased to the ultimate value. This coefficient is relevant to the number of contacts and represents anisotropy in the contact orientation. For compression testing at b = 0, loss of contact can occur in two lateral directions perpendicular to the major principal stress direction. For extension testing at b = 1, there are two major principal directions of loading and loss of contact occurs in one minor direction. It can be seen that an increase in b will increase the maximum deviatoric coefficient of fabric anisotropy. At large strains, the magnitude of this coefficient for different intermediate stress ratios differs slightly. Fig. 17(b) shows the more rapid evolution of the deviatoric coefficient of normal force anisotropy. It evolves to a strain that is similar to the strain of development of positive pore pressure, and then begins to reduce to its ultimate value. Unlike the deviatoric coefficient of fabric anisotropy, an increase in b decreased the deviatoric coefficient of normal force anisotropy. Fig. 17(c) shows the rapid growth of the deviatoric coefficient of tangential contact force anisotropy, which then declines slowly to the ultimate value. Similar to the deviatoric coefficient of normal force anisotropy, the deviatoric coefficient of the tangential contact force anisotropy decreased as b increased. The 11

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

b) 0.6

0.4

0.2

b=0 b=0.5 b=1

0 0

Dev. coeff. of tangential force anisotropy

c)

Dev. coeff. of normal force anisotropy

Dev. coeff. of fabric anisotropy

a)

5

b=0.25 b=0.75

10 15 Axial strain %

1 0.8 0.6 0.4 b=0 b=0.5 b=1

0.2

b=0.25 b=0.75

0

20

0

5

10 15 Axial strain %

20

0.4

0.3

0.2

0.1

b=0 b=0.5 b=1

0 0

5

b=0.25 b=0.75

10 15 Axial strain %

20

Fig. 17. Deviatoric coefficient of: (a) fabric anisotropy; (b) normal force anisotropy; (c) tangential force anisotropy versus axial strain at a confining pressure of 160 kPa and b values of 0, 0.25, 0.5, 0.75 and 1.

the results occurred at b = 0.75 and 1 (Fig. 18(d) and 18(e)) at strains of 1% to 7%. The maximum difference was under 10%. When the deviatoric stress increases in granular media, strong force chains evolve parallel to the major principal stress directions. Fig. 19 shows the force chain network at peak strength at b = 0 and 1 for a confining pressure of 160 kPa. For the compression test with one axis of loading in the direction Z, the strong force chains in the ZX and ZY planes show evolution that is almost parallel to the Z-axis. In the X and Y directions, the assembly expanded, the particles did not transmit high contact forces, and no strong force chain evolved. For the extension test, with two axes of loading (directions Z and X) in the ZX plane, strong force chain evolution created an orthogonal network parallel to the Z and X directions. In the ZY plane, strong force chains formed only in the direction parallel to the Z-axis because the Y-axis is the tensile axis of assembly and interparticle contact forces in this direction are slight. Barreto and O'Sullivan [9] divided the fabric tensor of assembly into strong and weak fabrics and concluded that a strong force chain network is relevant to changes in the principal stress and becomes anisotropic when increasing the anisotropic load; however, the weak force network that supports the strong force chains remains almost isotropic during loading. They studied drained true triaxial tests. In the current study, strong and weak fabric tensors were used to quantify the force chain evolution during shear in the undrained true triaxial test. The fabric tensor has already been defined by Eq. (28). The strong contacts are those that are equal to or greater than the average contact force and weak contacts are those that are lower than the average contact force according to Thornton and Antony [66]. Fig. 20(a) and 20(b), respectively, illustrate the changes in the major and minor principal components of the strong and weak fabric tensors for compression testing at a confining pressure of 160 kPa. Fig. 20(a)

representation is called the general stress tensor and is represented as:

σij =

1 V

3

∑

fic l jc

i, j = 1

(39)

The SFF equation relates the stress tensor of the granular assembly to micromechanical features. This equation identifies the contributions of fabric anisotropy and normal and tangential contact force anisotropies to the shear strength of granular material. It was proposed by Rothenburg and Bathurst [63] and generalized to ellipsoid-shaped particles by Ouadfel [65]. The SFF relationship proposed by Chantawarangul is used in the current study as [45]:

σij =

n m v f¯o l¯o

3

⎡δij + 2 ⎛aijr + aijn + 3 aijt ⎞ + 2 {(akln − aklt ) aklr δij 5⎝ 2 ⎠ 35 ⎣

+ (4ailn + 3ailt ) aljr } ⎤ ⎦

(40)

where m v is the contact density, or the ratio of twice the number of n physical contacts to the assembly volume, f¯o is the average normal contact force, l¯o is the contact vector length averaged over all contacts, aijr is the fabric anisotropy tensor (Eq. (27)), aijn is the normal contact force anisotropy tensor (Eq. (31)) and aijt is the tangential contact force anisotropy tensor (Eq. (34)). Fig. 18 compares the effective frictional strengths as calculated by these two approaches for a confining pressure of 160 kPa and b values of 0, 0.25, 0.5, 0.75 and 1. The results of the general stress tensor equation and SFF equation show good consistency. The evolutionary trend of the curves is similar; however, the SFF method underestimated the frictional resistance to some extent. The greatest difference between 12

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

General stress tensor SFF equation 0

Sine of effective mobilized friction angle

c)

Sine of effective mobilized friction angle

5

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

10 15 Axial strain %

5

10

15

General stress tensor SFF equation 0

d)

20

5

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

10 15 Axial strain %

20

General stress tensor SFF equation 0

Axial strain % 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

20

General stress tensor SFF equation

0

e)

Sine of effective mobilized friction angle

b)

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

Sine of effective mobilized friction anggle

Sine of effective mobilized friction angle

a)

5

10

15

20

Axial strain %

General stress tensor SFF equation 0

5

10 15 Axial strain %

20

Fig. 18. Sine of effective mobilized friction angle vs. axial strain at a confining pressure of 160 kPa at: (a) b = 0; (b) b = 0.25; (c) b = 0.5; (d) b = 0.75, (e) b = 1.

and micro-mechanical features of undrained true triaxial tests by considering the effect of the b parameter. Newton's second law of motion has been used to solve the particulate phase and the Navier-Stokes momentum equation and continuity equation were employed for solving the fluid phase. A new CFDEM solver with the ability to consider compressible fluid and dynamic meshes has been developed and named cfdemSolverSonicLiquidDyM. The results of drained and undrained triaxial tests on Nevada sand have been used to calibrate the model. Variation in the intermediate stress ratio affected the macro- and micro-mechanical behavior of the samples. When implementing deviatoric stress at small strains, the sample tended to compact, increasing the pore pressure. As the value of b varied from 0 toward 1, positive pore pressure growth occurred more rapidly and at greater values; however, the positive pore pressure dropped faster at higher b values. In other words, as the intermediate stress ratio increased from the compression test toward the extension test, the undrained granular samples tended to show the greatest dilative behavior and the absolute maximum negative excess pore pressure increased. The frictional resistance of the undrained samples was also dependent on the value of b. The maximum effective frictional angle increased from compression testing (b = 0) toward b = 0.5 and, when

shows that after implementing deviatoric loading, the minor principal component of strong fabric tensor R3s decreased and the major principal component of strong fabric tensor R1s increased. In Fig. 20(b), R1w and R3w , respectively, the major and minor principal components of the weak fabric tensor, remain constant during shear and the difference between them (R1w − R3w ) is nearly zero. It could be concluded that the strong fabric shows anisotropic evolution and weak fabric stays isotropic. For tests with b values greater than 0, this observation is qualitatively the same. This conclusion can also be qualitatively seen in Fig. 19. For compression testing at b = 0, where one axis of major principal stress exists, the ZY and ZX planes show a weak network of force chains with isotropic evolution that form a lateral support for the strong force chains. In the extension test at b = 1, the weak force chains evolved isotropically in the ZX and ZY planes. As noted, the Y-axis is the direction of tensile strain and the weak force network in this direction forms a lateral support for the strong force chains.

5. Conclusion This paper has made a contribution to understanding the macro13

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

ZY plane

b=1

b=0

ZX plane

Fig. 19. Contact force networks at peak strength for b = 0 and 1 at a confining pressure of 160 kPa.

b)

0.6 0.5

Weak principal fabric

Strong principal fabric

a)

0.4 0.3 0.2 major principal of strong fabric minor principal of strong fabric

0.1 0

0.6 0.5 0.4 0.3 0.2 major principal of weak fabric minor principal of weak fabric

0.1 0

0

5

10 15 Axial strain %

20

0

5

10 15 Axial strain %

20

Fig. 20. Major and minor principal values of: (a) strong; (b) weak fabric in true triaxial compression testing at a confining pressure of 160 kPa.

evolution of the micro features in the granular assembly. For quantitative comparison, the evolution of the scalar measures of the fabric, normal and tangential contact force anisotropies for different values of b were evaluated. The deviatoric coefficient of fabric anisotropy increased with an increase in the value of b and a reversal in behavior was observed for the deviatoric coefficient of (normal and tangential) contact force anisotropy. It can be inferred that these coefficients are complementary and counteract each other to keep the system stable. An additional microstructure, which was focused on, is the formation of strong and weak force chains. It has been qualitatively and quantitatively confirmed that strong force chains become anisotropic in correlation with principal stress anisotropy, but weak force chains show isotropic behavior during shear.

moving toward extension testing (b = 1), it decreased. The highest effective frictional angle was recorded for tests at b = 0.5. These trends were qualitatively similar at different confining pressures. For a given b value, the sample with a higher confining pressure had lower effective frictional resistance and a higher absolute maximum negative excess pore pressure. The results obtained using the SFF equation, which relates the stress tensor of the assembly to the anisotropy tensors of the fabric, normal and tangential contact forces, were compared to the results obtained using the general stress tensor and good agreement was observed. The distribution of contact normals, average normal contact forces, and average tangential contact forces are presented in a spherical coordinate system for b values of 0, 0.5 and 1. As the b increased from the compression test toward the extension test, the histograms of contact normals and average normal contact forces evolved from a peanutshaped distribution to a donut shape without the central hole. The histogram of tangential contact forces evolved from a dumbbell with an axis of symmetry of the major principal stress axis to a dumbbell with an axis of symmetry of the minor principal stress axis. These illustrations help to visually comprehend the effect of the stress path on the

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 14

Computers and Geotechnics xxx (xxxx) xxxx

T. Foroutan and A.A. Mirghasemi

of saturated soil by a fluid coupled-DEM model. Eng Geol 2015;193:256–66. [33] Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol 1993;77(1):79–87. [34] Zeghal M, El Shamy U. A continuum-discrete hydromechanical analysis of granular deposit liquefaction. Int J Numer Anal Meth Geomech 2004;28(14):1361–83. [35] Zeghal M, El Shamy U. Liquefaction of saturated loose and cemented granular soils. Powder Technol 2008;184(2):254–65. [36] Climent N, Arroyo M, O’Sullivan C, Gens A. Sand production simulation coupling DEM with CFD. Eur J Environ Civ Eng 2014;18(9):983–1008. [37] Jiang M, Liu J, Sun C. An upgraded CFD-DEM investigation on the liquefaction mechanism of sand under cyclic loads. International, conference on discrete element methods. Springer; 2016. p. 609–17. [38] Zhang D-M, Gao C-P, Yin Z-Y. CFD-DEM modeling of seepage erosion around shield tunnels. Tunn Undergr Space Technol 2019;83:60–72. [39] Guo Y, Yu XB. Comparison of the implementation of three common types of coupled CFD-DEM model for simulating soil surface erosion. Int J Multiph Flow 2017;91:89–100. [40] Kloss C, Goniva C, Hager A, Amberger S, Pirker S. Models, algorithms and validation for opensource DEM and CFD–DEM. Prog Comput Fluid Dyn Int J 2012;12(2–3):140–52. [41] Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117(1):1–19. [42] Available from: www.openfoam.org. [43] Goniva C, Kloss C, Deen NG, Kuipers JA, Pirker S. Influence of rolling friction on single spout fluidized bed simulation. Particuology 2012;10(5):582–91. [44] Arulmoli K, Muraleetharan K, Hossain M, Fruth L. VELACS verification of liquefaction analyses by centrifuge studies laboratory testing program soil data report. Irvine, CA: Earth Technology Corporation; 1992. [45] Chantawarangul K. Numerical simulation of three dimensional granular assemblies. university of Waterloo; 1993. [46] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique 1979;29(1):47–65. [47] Ai J, Chen J-F, Rotter JM, Ooi JY. Assessment of rolling resistance models in discrete element simulations. Powder Technol 2011;206(3):269–82. [48] Goren L, Aharonov E, Sparks D, Toussaint R. Pore pressure evolution in deforming granular material: a general formulation and the infinitely stiff approximation. J Geophys Res: Solid Earth 2010;115(B9). [49] Issa RI. Solution of the implicitly discretised fluid flow equations by operatorsplitting. J Comput Phys 1986;62(1):40–65. [50] O'Sullivan C. Particulate discrete element modelling: a geomechanics perspective. CRC Press; 2011. [51] Zhu HP, Zhou ZY, Yang RY, Yu AB. Discrete particle simulation of particulate systems: theoretical developments. Chem Eng Sci 2007;62(13):3378–96. [52] Gidaspow D. Multiphase flow and fluidization: continuum and kinetic theory descriptions. Academic press; 1994. [53] Ergun S, Orning AA. Fluid flow through randomly packed columns and fluidized beds. Ind Eng Chem 1949;41(6):1179–84. [54] Calvetti F, Nova R. Micromechanical approach to slope stability analysis. Degradations and instabilities in geomaterials. Springer 2004:235–54. [55] Utili S, Nova R. DEM analysis of bonded granular geomaterials. Int J Numer Anal Meth Geomech 2008;32(17):1997–2031. [56] Gabrieli F, Cola S, Calvetti F. Use of an up-scaled DEM model for analysing the behaviour of a shallow foundation on a model slope. Geomech Geoeng 2009;4(2):109–22. [57] Cheng K, Wang Y, Yang Q, Mo Y, Guo Y. Determination of microscopic parameters of quartz sand through tri-axial test using the discrete element method. Comput Geotech 2017;92:22–40. [58] Li Y, Xu Y, Thornton C. A comparison of discrete element simulations and experiments for ‘sandpiles’ composed of spherical particles. Powder Technol 2005;160(3):219–28. [59] Alshibli KA, Williams HS. A true triaxial apparatus for soil testing with mixed boundary conditions. Geotech Test J 2005;28(6):534–43. [60] Ng T-T. Input parameters of discrete element methods. J Eng Mech 2006;132(7):723–9. [61] Wang Q, Lade PV. Shear banding in true triaxial tests and its effect on failure in sand. J Eng Mech 2001;127(8):754–61. [62] Xiao Y, Liu H, Chen Y, Zhang W. Particle size effects in granular soils under true triaxial conditions. Geotechnique 2014;64(8):667–72. [63] Rothenburg L, Bathurst RJ. Analytical study of induced anisotropy in idealized granular materials. Geotechnique 1989;39(4):601–14. [64] Rothenburg L. Micromechanics of idealized granular systems. Carleton University; 1981. [65] Ouafdel H. Numerical simulations of granular assemblies with three-dimensional ellipsoid-shaped particles. University of Waterloo; 1998. [66] Thornton C, Antony S. Quasi–static deformation of particulate media. Philos Trans R Soc London Ser A 1998;356(1747):2763–82.

References [1] Choi C, Arduino P, Harney M. Development of a true triaxial apparatus for sands and gravels. Geotech Test J 2008;31(1):32–44. [2] Chu J, Lo S-CR, Lee IK. Strain softening and shear band formation of sand in multiaxial testing. Géotechnique 1996;46(1):63–82. [3] Lade PV, Duncan JM. Cubical triaxial tests on cohesionless soil. J Soil Mech Found Div 1973;99(SM10):793–812. [4] Lade PV. Analysis and prediction of shear banding under 3d conditions in granular materials. Soils Found 2003;43(4):161–72. [5] Ng T-T. Macro-and micro-behaviors of granular materials under different sample preparation methods and stress paths. Int J Solids Struct 2004;41(21):5871–84. [6] Ng T-T. Behavior of gravity deposited granular material under different stress paths. Can Geotech J 2005;42(6):1644–55. [7] Thornton C. Numerical simulations of deviatoric shear deformation of granular media. Géotechnique 2000;50(1):43–53. [8] Thornton C, Zhang L. On the evolution of stress and microstructure during general 3D deviatoric straining of granular media. Géotechnique 2010;60(5):333–41. [9] Barreto D, O’Sullivan C. The influence of inter-particle friction and the intermediate stress ratio on soil response under generalised stress conditions. Granular Matter 2012;14(4):505–21. [10] O’sullivan C, Wadee M, Hanley K, Barreto D. Use of DEM and elastic stability analysis to explain the influence of the intermediate principal stress on shear strength. Géotechnique 2013;63(15):1298–309. [11] Huang X, Hanley KJ, O’Sullivan C, Kwok CY, Wadee MA. DEM analysis of the influence of the intermediate stress ratio on the critical-state behaviour of granular materials. Granular Matter 2014;16(5):641–55. [12] Zhou W, Wu W, Ma G, Huang Y, Chang X. Study of the effects of anisotropic consolidation on granular materials under complex stress paths using the DEM. Granular Matter 2017;19(4):76. [13] Zhou W, Liu J, Ma G, Chang X. Three-dimensional DEM investigation of critical state and dilatancy behaviors of granular materials. Acta Geotech 2017;12(3):527–40. [14] Xiao Y, Sun Y, Liu H, Yin F. Critical state behaviors of a coarse granular soil under generalized stress conditions. Granular Matter 2016;18(2):17. [15] Sazzad MM, Suzuki K. Density dependent macro-micro behavior of granular materials in general triaxial loading for varying intermediate principal stress using DEM. Granular Matter 2013;15(5):583–93. [16] Dorostkar O, Mirghasemi AA. Micro-mechanical study of stress path and initial conditions in granular materials using DEM. Comput Particle Mech 2016;3(1):15–27. [17] Dorostkar O, Mirghasemi AA. On the micromechanics of true triaxial test, Insights from 3D DEM study. Iran J Sci Technol Trans Civ Eng 2018:1–15. [18] Yimsiri S, Soga K. DEM analysis of soil fabric effects on behaviour of sand. Géotechnique 2010;60(6):483–95. [19] Sitharam TG, Dinesh SV, Shimizu N. Micromechanical modelling of monotonic drained and undrained shear behaviour of granular media using three-dimensional DEM. Int J Numer Anal Meth Geomech 2002;26(12):1167–89. [20] Sitharam TG, Dinesh SV. Numerical simulation of liquefaction behaviour of granular materials using Discrete Element Method. J Earth Syst Sci 2003;112(3):479. [21] Sitharam TG, Vinod JS. Critical state behaviour of granular materials from isotropic and rebounded paths: DEM simulations. Granular Matter 2009;11(1):33–42. [22] Gong G. DEM [Discrete Element Method] simulations of drained and undrained behaviour. University of Birmingham; 2008. [23] Gong G, Thornton C, Chan AHC. DEM simulations of undrained triaxial behavior of granular material. J Eng Mech 2012;138(6):560–6. [24] Zhao J, Guo N. Rotational resistance and shear-induced anisotropy in granular media. Acta Mech Solida Sin 2014;27(1):1–14. [25] Liu Y-J, Li G, Yin Z-Y, Dano C, Hicher P-Y, Xia X-H, et al. Influence of grading on the undrained behavior of granular materials. Comptes Rendus Mécanique 2014;342(2):85–95. [26] Nguyen HBK, Rahman MM, Fourie AB. Characteristic behavior of drained and undrained triaxial compression tests: DEM study. J Geotech Geoenviron Eng 2018;144(9):04018060. [27] Hakuno M, Tarumi Y. Sand liquefaction analysis by granular assembly simulation. Ninth world conference on earthquake engineering. 1988. p. 231–6. [28] Olivera Bonilla RR. Numerical simulations of undrained granular media. University of Waterloo; 2004. [29] Shafipour R, Soroush A. Fluid coupled-DEM modelling of undrained behavior of granular media. Comput Geotech 2008;35(5):673–85. [30] Khalili Y, Mahboubi A. Discrete simulation and micromechanical analysis of twodimensional saturated granular media. Particuology 2014;15:138–50. [31] Okada Y, Ochiai H. Coupling pore-water pressure with distinct element method and steady state strengths in numerical triaxial compression tests under undrained conditions. Landslides 2007;4(4):357–69. [32] Liu G, Rong G, Peng J, Zhou C. Numerical simulation on undrained triaxial behavior

15