- Email: [email protected]
Discrete simulation and micromechanical analysis of two-dimensional saturated granular media
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Particuology xxx (2013) xxx–xxx
Contents lists available at ScienceDirect
Particuology journal homepage: www.elsevier.com/locate/partic
Discrete simulation and micromechanical analysis of two-dimensional saturated granular media Younes Khalili, Ahmad Mahboubi ∗ Department of Civil and Environmental Engineering, Power and Water University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 19 June 2012 Received in revised form 20 May 2013 Accepted 18 July 2013 Keywords: DEM Saturated granular media Numerical simulation Micromechanical investigation
a b s t r a c t In this study, a novel approach to incorporate the pore water pressure in the discrete element method (DEM) to comprehensively model saturated granular media was developed. A numerical model was constructed based on the DEM by implanting additional routines in the basic DEM code; pore water pressure calculations were used with a two-dimensional (2D) model to simulate the undrained behavior of saturated granular media. This model coupled the interaction of solid particles and the pore fluid in saturated granular media. Finally, several 2D undrained shear tests were simulated. The test results showed that the model could predict the response of the saturated granular soil to shear loading. The effect of initial compaction was investigated. Biaxial tests on dense and loose specimens were conducted, and the effect of the initial density on the change in shear strength and the volume change of the system was investigated. The overall behavior of loose and dense specimens was phenomenologically similar to the real granular material. Constant volume tests were simulated, and the results were compared to those from the coupled model. Induced anisotropy was micromechanically investigated by studying the contact force orientation. The change in anisotropy depended on the modeling scheme. However, the overall responses of the media obtained using the coupled and constant volume methods were similar. © 2013 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.
1. Introduction Granular materials are important to several industrial sectors, such as agriculture and food, pharmaceuticals, process engineering and construction. Soil is a major geological and granular material, which is used to construct roads, dams, and other earth structures because of its properties. The saturated response of a granular medium is of particular importance to practical applications, such as embankment dams. A saturated material is a two-phase mixture: the discrete phase consists of solid particles, and the continuum phase consists of the pore fluid. The interaction between these two phases constitutes the overall behavior of a saturated medium. Therefore, a saturated granular material shows complex behavior. The saturated soil behavior is generally modeled using analytical continuum models. These analytical models assume continuum behavior throughout the multiphase particulate medium. The continuum models are phenomenological and neglect the particulate behavior of granular materials, which certainly affects the deformation and failure modes of the material. Continuum models cannot
∗ Corresponding author. Tel.: +98 21 77312552; fax: +98 21 7700 6660. E-mail addresses: [email protected], ar [email protected] (A. Mahboubi).
explore the physical microscopic mechanisms that lead to liquefaction and can only simulate the macroscopic response if the material behavior is adequately emulated by the adopted constitutive law. An alternative approach to studying the behavior of a granular material accounts for the discrete nature of the particles. Early attempts to observe force distributions and particle behavior at microscopic levels were made by Dantu (1957) and Wakabayashi (1950) using optically sensitive materials. Later, De Josselin de Jong and Verruijt (1969) analyzed the force distribution in such assemblies by studying the individual particles. The results of these experiments provided a qualitative understanding of the mechanisms of load transfer and the deformation characteristics of granular materials. New experimental techniques using computer tomography have also provided important information on void ratio evolution in the fabric (Desrues, Chambon, Mokni, & Mazerolle, 1996). The discrete element method (DEM), which was originally developed by Cundall (1971) to numerically analyze rock blocks, assuming that each rock block satisfies an equation of motion, and later used to simulate an assembly of granular materials (Cundall & Strack, 1979), is a successful computational technique that treats a granular system as a collection of discrete granules. The method allows all of the required information to be extracted at the particle level at any stage of the deformation process. A comprehensive review of developments and applications of
1674-2001/$ – see front matter © 2013 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.partic.2013.07.005
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
2
DEM to simulate particulate systems has been published by Zhu, Zhou, Yang, and Yu (2007, 2008). Although this technique was originally developed for assemblies of idealized particles in a dry environment under shear, it can be extended to account for the effects of a pore fluid. The methods that incorporate fluid effects into DEM can generally be divided into (a) constant volume method and (b) coupled method. Numerical method that indirectly considers fluid effects is known as constant volume method. A constant volume method is based on the almost constant volume of saturated soil specimens under undrained loading because of the very low compressibility of the particles and the interstitial water. The alternative method, which considers fluid-particle interaction, is known as the coupled method. The inter-particle forces are caused by interactions between the solid particles, and the hydrostatic and hydrodynamic forces exerted on the solid particles result from changes in the pore volume and the fluid flow. The fluid flow in porous media is mainly controlled by the gradient of the pore pressure and the coefficient of permeability of the media. The coefficient of permeability of a granular medium depends on the porosity of the medium. The response of an individual particle is governed by the resultant of mechanical forces and hydromechanical forces acting on the particle. A complete review of applications of DEM to geomechanics has been presented by Jiang and Yu (2006). The main applications of DEM to simulate saturated granular materials are briefly explained in Table 1. In this study, a saturated soil model was formulated using the coupled method consisting of two interacting phases to model the mechanical behavior of the saturated granular materials. The elastic membrane in actual shear tests was simulated by a string of boundary particles. The boundary particles were connected to each other, but the boundary could deform freely like an elastic membrane. Imposing constant stress boundary conditions allowed the specimen to fail at the weakest plane. The bulk modulus and the permeability coefficient of each cell of the assembly of saturated particles were calculated at each time step from the particles situations. Constant volume tests were conducted, and the overall results and the local parameters were compared with those of the coupled method tests. A micromechanical investigation was also conducted. The differences between the constant volume method and the coupled method were examined.
The coupled model of the interactions between the solid phase and the fluid phase has three main components: • the fluid phase behavior and pressure generation in the fluid phase, • the solid phase dynamics, and • the interaction between the fluid phase and the solid phase. 2.1. Fluid phase behavior The governing equations for the liquid phase are the continuity and Darcy’s equations, which are given by:
vf = −k
∇ u w
,
2.1.1. Volume change Volume changes are calculated using the particle displacements in the cells. In this study, an indirect approach was used to calculate the volumetric strain from the porosity difference at the beginning and end of each time step: εv (i, j) =
V (i, j) n1 (i, j) − n2 (i, j) = , V (i, j) 1 − n2 (i, j)
(3)
where V(i, j) denotes the volume change during the time step t, V(i, j) denotes the volume of the cell (i, j), and n1 (i, j) and n2 (i, j) are the porosity of the assembly of particles in the cell (i, j) at the beginning and end of each time step, respectively. The porosity in a 2D medium can be defined as the void area divided by the total area. A row and column for a cell are denoted by i and j, respectively. 2.1.2. Pore pressure increment The volume change in cells that are saturated with a pore fluid in the inter-particle spaces can increase or decrease the pore pressure. The intensity of this effect depends on the bulk moduli of the particles and the fluid. Within a linear approximation, the bulk modulus of a saturated granular material can be calculated as: 1 nav (i, j) 1 − nav (i, j) + , = Bp Bw B(i, j)
(4)
where B(i, j) denotes the bulk modulus of the water-particle assembly in cell (i, j), and Bp and Bw are the bulk moduli of the particles and the water, respectively. The average porosity of the particles in cell (i, j) is denoted by nav (i, j) and can be calculated using Eq. (5): nav (i, j) =
n1 (i, j) + n2 (i, j) . 2
(5)
Thus, a change in the fluid volume can change the pore pressure. The increase in pore pressure in each cell is calculated by multiplying the volumetric strain by the bulk modulus of the saturated soil, as shown in Eq. (6).
2. Coupled fluid-particle model
∂n + ∇ · (nvf ) = 0, ∂t
where n is the porosity, vf is the fluid velocity vector, is the gradi is the pressure, and w is the ent operator, k is the permeability, u unit weight of the fluid. Implanting the fluid phase in DEM consists of three main computational steps involving three parameters, the volume change, pore pressure increment and fluid flow. These steps are implemented by dividing the fluid continuum into many cells using vertical and horizontal virtual lines. The fluid parameters and properties, such as the fluid velocity or the fluid pressure, are assumed to be constant in each cell.
(1)
(2)
u(i, j) = εv (i, j)B(i, j).
(6)
This pressure increase must be modified to account for the fluid flow between the cells during the time step. 2.1.3. Fluid flow The volume changes are not equal among the cells, and cells with different bulk moduli have different pressure increases and pore pressures. Thus, a pressure gradient is generated between the cells. The fluid flow depends on the coefficient of permeability of each cell and dissipates the pressure difference between the cells. In this study, Darcy’s law was assumed to be valid. Considering Fig. 1, which illustrates the pore fluid flow between adjacent cells, the continuity equation can be written as: u∗2 (i, j) − u1 (i, j) = B(i, j)
Vw (i, j) , V (i, j)
(7)
where u1 (i, j) denotes the pore pressure of cell (i, j) at the beginning of the time step (t). Considering the fluid flow between cells, the
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model
ARTICLE IN PRESS
PARTIC-603; No. of Pages 13
Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
3
Table 1 A brief list of DEM modeling of saturated granular materials. Model
Application
Remarks
References
Coupled (discrete–continuum model)
Sand liquefaction analysis
2D model; water flow between the pores are modeled by Darcy’s law and the buoyancy and hydrodynamic forces were neglected. Numerical modifications on Hakuno and Tarumi (1988) model
Hakuno and Tarumi (1988)
Liquefaction Liquefaction mechanism Liquefaction and quicksand simulation
3D model based on Hakuno’s model to simulate cyclic shear loading 3D model; using the continuity and Navier–Stokes equations for fluid modeling
Nakase, Ishikawa, and Fujitani (1996) Meguro and Ravichandran (2000) Zeghal and El Shamy (2004) and El Shamy and Zeghal (2005, 2007) Li, Chu, and Sheng (2007)
Biaxial test on saturated closely packed assemblage of particles Behavior of saturated soil mass Sand pile in water
DEM-SPH modeling; monodisperse hexagonal arrangement; the effect of Navier–Stokes and Darcy’s law on specimen behavior was studied. 2D model, laminar flow between the pores (Darcy’s law assumed to be valid) 3D model; CFD-DEM method, using Navier–Stokes equations for fluid modeling
Coupled (discrete–discrete)
Modeling the behavior of granular materials
2D model; fluid is modeled as overlapping cylinders, fluid flow was not considered
Lamei and Mirghasemi (2011)
Constant volume method
Monotonic and cyclic response of a granular soil
2D and 3D simulation of undrained cyclic simple shear test; pore pressure buildup to initial liquefaction and hysteresis loop formation are shown Effect of the particle shape on liquefaction behavior; significant effect of particle morphology on liquefaction susceptibility for sands with the same void ratio Studying the effect of various mechanical parameters on liquefaction
Dobry and Ng (1992) and Ng and Dobry (1994)
Liquefaction
Liquefaction
pore pressure of cell (i, j) at the end of the time step is denoted by u∗2 (i, j) in which the pore pressure increases caused by the cells’ volumetric strains are neglected. The volume (area) of cell (i, j) is denoted by V(i, j), and the total fluid volume that enters cell (i, j) from the adjacent cells is denoted by Vw (i, j), which is calculated using Eq. (8):
q1 =
Ashmawy, Sukumaran, and Hoang (2003) Sitharam (2003), Sitharam and Dinesh (2003), Dinesh, Sitharam, and Vinod (2004), and Sitharam and Vinod (2008)
2k (i, j − 1)k (i, j) u(i, ¯ j − 1) − u(i, ¯ j) x x kx (i, j − 1) + kx (i, j)
qk (i, j).
Zhao and Shan (2013)
j) are shown in Fig. 1. Note that an equivalent permeability coefficient was used for each pair of adjacent cells to account for the different permeability coefficients of adjacent cells. The first term in Eqs. (9)–(12) denotes the equivalent permeability coefficient.
4
Vw (i, j) = t
Shafipour and Soroush (2008)
(8)
q2 =
In Eq. (8), t is the time step, and q1 (i, j) to q4 (i, j) denote the discharges entering cell (i, j), which are calculated according to Darcy’s law as given in Eqs. (9)–(12). The finite difference method was used to calculate the fluid flow. The positive directions of q1 (i, j) to q4 (i,
q3 =
k=1
dy ,
¯ j) ¯ + 1, j) − u(i, u(i dx w dy
(9)
,
2k (i, j + 1)k (i, j) u(i, ¯ j) ¯ j + 1) − u(i, x x kx (i, j + 1) + kx (i, j)
q4 =
2ky (i + 1, j)ky (i, j) ky (i + 1, j) + ky (i, j)
w dx
2ky (i − 1, j)ky (i, j) ky (i − 1, j) + ky (i, j)
w dx
dy ,
¯ − 1, j) − u(i, ¯ j) u(i dx w dy
(10)
(11)
.
(12)
In the above equations, w denotes the unit weight of the pore fluid, dx and dy are the cell dimensions in the x and y directions, and kx (i, j) and ky (i, j) are the coefficients of permeability of cell (i, j) in the x and y directions, respectively. The kx and ky coefficients were assumed to be equal, which is realistic for granular materials. These coefficients are defined by the following equation (Das, 1985): 2
kx (i, j) = ky (i, j) = kinitial
e(i, j) , einitial 2
(13)
where e denotes the void ratio and k denotes the coefficient of permeability of each cell. The initial coefficient of permeability kinitial corresponds to einitial or the average void ratio at the beginning of the axial loading, i.e., axial strain = 0%. The equivalent pore pressure ¯ j), which is calculated as a linduring a time step is denoted by u(i, ear combination of the cell pressures at the beginning and at the end of the time step, as given in Eq. (14): ¯ j) = (1 − ˛)u1 (i, j) + ˛u∗2 (i, j), u(i, Fig. 1. Positive directions of discharges for cell (i, j) shown by q1 (i, j) to q4 (i, j).
(14)
where ˛ is a multiplier between 0.5 and 1. In this study, ˛ was taken to be equal to 0.6. The aforementioned equations were used
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
4
Table 2 Added subroutines. Subroutine No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Fig. 2. Contact model between solid particles.
to calculate u* for the cells; the calculation procedure needed to be repeated several times to reach an acceptable level of accuracy. Finally, the sum of u∗2 (i, j), as given by Eq. (7), and the change in pressure, u(i, j), as given by Eq. (6), yields the total pressure in cell (i, j) at the end of time step u(i, j)(t+t) . u(i, j)
(t+t)
= u(i, j) + u∗2 (i, j).
(15)
DEM models a soil mass as many interacting particles. The positions and velocities of all of the particles in DEM are determined by satisfying two main constraints: mv˙ = ˙ = Iω
f,
(16)
M,
(17)
where m denotes the particle mass, v denotes the particle linear velocity vector, f denotes the force vector acting on the particle, denotes the angular I denotes the particle moment of inertia, ω denotes the moment vector velocity vector of the particle, and M acting on particle. Time derivatives are denoted by dots. The particles were assumed to be rigid with soft contact areas; thus, particles in contact with each other may overlap. The contact law between the particles was implemented using two normal and tangential springs, a non-viscous damper in the normal direction and a slider in the tangential direction, as shown in Fig. 2. 2.3. Interaction between phases The pore fluid exerts interaction forces on particles because of the difference in the velocities of the fluid and the particles and the pressure gradient in the fluid phase. These forces are divided into two main types: hydrostatic and hydrodynamic forces. The hydrostatic force (i.e., the buoyancy force) can be calculated as follows: u(i, j))V p , fb = −(∇ p
(18)
where Vp is the volume of particle p. The hydraulic gradient vector u(i, j) and can be expressed in terms of for cell (i, j) is denoted by ∇ the pressure gradient between adjacent cells:
u(i, j) = ∇
Generate particle assemblage with appropriate characteristics Generate inter-connected boundary particles Divide specimen into cells Calculate cell porosity Calculate cell permeability Calculate volumetric strain of cells Calculate total volumetric strain of particles assemblage Calculate bulk modulus of cells Perform fluid flow calculations Calculate cell fluid velocity Calculate average particle velocity in cells Calculate hydrodynamic forces in X and Y directions Perform calculations of pore pressure using constant volume method Isotropic compaction of generated specimen Perform undrained simulation using coupled method Perform undrained simulation using constant volume method Perform drained simulation
Only drag forces were considered in this study, neglecting all other types of hydrodynamic forces. The drag force can be evaluated using Ergun’s equation as follows (Ergun, 1952):
fd (i, j) = (1−n(i, j))
2.2. Solid phase dynamics
17
Function
u1 (i, j + 1) − u1 (i, j − 1) u1 (i − 1, j) − u1 (i + 1, j) , . 2dx 2dy (19)
f |¯vf (i, j) − v¯ p (i, j)| [1−n(i, j)] 150 f + 1.75 2 ¯ d¯ p n(i, j)dp
[¯vf (i, j) − v¯ p (i, j)],
(20)
where n(i, j) denotes the porosity of cell (i, j), f the pore fluid viscosity, d¯ p the average particle diameter (d50 is used in this study), f the fluid density, and v¯ f (i, j) and v¯ p (i, j) are the average velocities of the fluid and the particles in cell (i, j), respectively. The hydrostatic and hydrodynamic forces acting on each particle were calculated using Eqs. (18)–(20). The coupled model describes the interaction between the fluid and solid phases by applying these forces to each particle at every time step. 3. Simulation with coupled method 3.1. Model implementation A 2D coupled continuum–discrete model of saturated granular soil was implemented using PFC2D (Itasca, USA), which can model the dry assemblies of particles. The subroutines, which were implanted in the code to model the pore fluid in the interstitial voids, calculated the volumetric strain, the velocities of the particles and the fluid, the hydrodynamic forces, the fluid flow and the pore pressure in the cells. Table 2 lists the most important subroutines that were added to the code. A general schematic of the calculation steps to model the saturated granular soil is shown in Fig. 3. The particle positions and particle velocities were specified at the start of the simulation. At each time step, the porosity in each cell was calculated. The bulk modulus and permeability coefficient of the soil in the cells were then updated because the bulk modulus and the permeability coefficient of the soil are porosity-dependent. The volumetric strain in the cells was calculated from the cell porosity. The calculated values for the volumetric strain and the bulk modulus were used to compute the pore pressure increase. This increment is updated following the fluid flow calculations using the previously described finite difference method. The permeability coefficient significantly affects the calculated pore water flow. The total pore pressure was then determined, from which the fluid
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
5
Fig. 3. Calculation procedure.
velocities were calculated. The calculated fluid velocities in the cells and the known values of the particle velocities were used in turn to calculate the drag forces from the Ergun’s equation; the hydrostatic forces were determined from the calculated pore pressure in the cells and the related pressure gradient. The contact forces were calculated from the particle positions using the contact law. The sum of the three types of forces was inserted into the equations of motion, and the particle positions and velocities were calculated as inputs for the next time step.
reduced to generate particles in random positions in a reasonable number of tries without overlapping between adjacent particles. In this study, the radius of any particle was divided by 1.6 in this step. In the second step, the particle radius was increased to its actual value. Particles were then allowed to overlap, producing large interparticle forces at the end of this step. The large inter-particle forces
3.2. Simulation scheme First, a dense specimen of particles was generated. Boundary particles were created at the edges of the generated specimen. The boundary conditions were applied to these boundary particles. The particle assembly was then divided into several cells. Increasing the number of cells increases the accuracy of the simulation but also increases the analysis run time. Therefore, the number of cells was maximized while ensuring a reasonable run time. Increasing the number of particles has the same advantages and limitations. After preparing the soil specimen, a confining pressure was applied to the boundary particles and maintained until the soil specimen was completely consolidated. The sample was then displacement loaded under undrained (or drained) conditions at a constant velocity. The lateral stress was maintained constant using a servo-mechanism to control the velocity of the boundary particles. After reaching a prescribed axial strain of the assembly, the loading was stopped. Note that a high loading velocity affects the simulation results and increases the shear strength parameters of the specimen. From studying the effect of the loading velocity on the overall results, a loading velocity was selected for the simulations that does not affect the global results. A saturated soil specimen with a height of 20 cm and a width of 10 cm was prepared to conduct the biaxial shear test described above. The specimen was divided into 18 pressure cells as shown in Fig. 4. The total number of boundary particles was 300. Table 3 lists the particle properties. Tables 4 and 5 present the properties of the top and bottom walls and the pore fluid (water), respectively. 3.2.1. Assembly generation and specimen preparation DEM uses a specific procedure to generate an assembly of particles. Different approaches can be used to generate the assembly. The generation method used in this study consists of three steps (see Fig. 5). In the first step, a primary assembly of particles with a reduced particle radius was generated. The particle radius was
Fig. 4. Assembly of particles.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model
ARTICLE IN PRESS
PARTIC-603; No. of Pages 13
Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
6
Fig. 5. Assembly generation and compaction: (a) Step I, generating an assembly of particles with reduced particle radii; (b) Step II, increasing particle radii to their actual values, allowing them to overlap; (c) Step III, allowing particles to move into an equilibrium state.
Table 3 Particle properties. Minimum diameter (mm) Maximum diameter (mm) Density (kg/m3 ) Normal stiffness (N/m) Shear stiffness (N/m) Bulk modulus (GPa) Initial permeability coefficient (mm/s)
3.6 5.58 2500 108 108 35 0.01
Table 4 Wall properties. 108 108 0.5
Normal stiffness (N/m) Shear stiffness (N/m) Friction coefficient
generated by particle overlap require that the velocity of particles be artificially reduced. Finally, the particles were allowed to move into an equilibrium state. At the end of the aforementioned steps, the inter-particle forces should be relatively small, and the particle assembly should be sufficiently compacted; otherwise, the number of generated particles must be adjusted. The inter-particle friction angle is an important parameter in the assembly generation. Selecting smaller friction angle values in the particle generation procedure results in the faster and more efficient compaction of the generated assembly. After attaining a desired compacted state, the inter-particle friction coefficient was changed to the actual value. Five different compacted specimens were prepared in this study. The specimens differed in their void ratios and consolidation pressures (see Table 6). The particle diameter distribution of the generated assembly is shown in Fig. 6.
Fig. 6. Particle diameter distribution.
triaxial test. In this study, the boundaries were modeled by a string of connected particles known as boundary particles (see Fig. 7). The connections between the boundary particles were defined such that they could only support tensile forces, i.e., the boundaries could be simply reshaped. This type of boundary has two main advantages compared to a rigid wall boundary. First, a shear band
3.2.2. Boundaries To realistically model a biaxial shear test, the boundaries of the test specimen should be deformable as an elastic membrane in a Table 5 Pore fluid (water) properties (Bardet & Sayed, 1993). Bulk modulus (GPa) Viscosity (Pa s)
2 0.001
Fig. 7. Elastic membrane formed by boundary particles.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model
ARTICLE IN PRESS
PARTIC-603; No. of Pages 13
Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
7
Table 6 Simulation conditions. Test No.
Simulation method
Initial void ratio, e
Confining pressure (kPa)
Inter-granular friction coefficient
Number of inner particles
1 2 3 4 5 6 7 8 9
U-C U-C U-C U-C U-C U-C U-C U-CV D
0.158 0.170 0.176 0.153 0.148 0.158 0.158 0.158 0.158
100 100 100 200 300 100 100 100 100
0.5 0.5 0.5 0.5 0.5 0.4 0.6 0.5 0.5
1228 1220 1209 1228 1228 1228 1228 1228 1228
D: drained and U: undrained test; C: coupled and CV: constant volume method.
can appear in the simulation results as has been observed in laboratory tests, and second, a confining pressure can be applied more uniformly. To facilitate specimen loading, two rigid walls were defined at the top and bottom of the specimen near the boundary particles. The sample was loaded by moving the top wall downwards at a constant velocity. The walls were also used to measure the axial force applied to the specimen. Boundary particles were used between the particle assembly and the rigid walls at the top and bottom of the specimen to transfer forces, i.e., the contact forces and the forces resulting from the pore fluid pressure, from the particle assembly to the rigid walls and to apply the predefined confining pressure to the assembly. 3.3. Simulation results As shown in Table 6, the numerical tests were performed under undrained and drained conditions. The objective of the tests was to numerically verify the effects of the initial compaction, the confining pressure, and the inter-particle friction coefficient on the global macroscopic behavior of the specimen. Simulations 1–7 were conducted using the coupled method to model the pore pressure. Test No. 9 was a drained test, which served as a reference for the drained behavior of an ideal granular material. Test No. 8 was a simulation of undrained loading using the alternative constant volume method. 3.3.1. Effect of pore pressure This section presents the study of the drained specimen behavior by simulation No. 9. The average void ratio and the total volume of the particle assembly in simulation No. 9 were not constant during the test. Fig. 8(a) shows the variation of the void ratio of the particle assembly in simulation Nos. 1–9. The void ratio of the drained specimen first decreased and then gradually increased to a maximum value. This behavior agreed with the experimental behavior of a dense sand specimen in a drained shear test. The void ratios of the undrained specimens (simulation Nos. 1–8) were almost constant, and the volume remained nearly constant during the shear test. Fig. 8(b) compares the variations of differential stress, ( 1 − 3 ) with the axial strain for simulation Nos. 1 and 9. Simulation No. 1 was a reference test to which the other simulations were compared. In the drained simulation (Test No. 9), the shear strength peaked before reaching a residual strength value, as expected from the dense sand behavior. At strains below 2%, a greater shear strength was observed for the drained specimen than that of the undrained test. This phenomenon can be attributed to the positive pore pressure of the undrained specimen (No. 1) at strains below 2%. The dilative behavior of the specimens at strains larger than 2% generated a negative pore pressure in the undrained specimen; thus, a higher differential stress was observed for the undrained specimens than for the drained specimen, as shown in Fig. 8(b).
Fig. 8(c) illustrates the change in the macroscopic effective internal friction angle ( 1 − 3 )/( 1 + 3 ) as a function of the axial strain, obtained from simulations 1 and 9. Both specimens showed similar effective behavior, which agreed well with the observed behavior of granular specimens in laboratory studies. Both diagrams show a peak at approximately 0.44 and a residual value of 0.36 of sin( ), which correspond to 26◦ and 21◦ , respectively, for the acceptable friction angles in dense assemblies of round (spherical) particles. 3.3.2. Effect of the initial void ratio on the behavior of undrained specimens Three identical specimens were prepared with different void ratios to study the effect of the initial void ratio on the specimens’ undrained behavior. The initial void ratios for Test Nos. 1, 2, and 3 were 0.158, 0.170, and 0.176, respectively. The results of the undrained biaxial tests for these specimens are shown in Fig. 9. The variations of the differential stress as a function of the axial strain for specimens 1, 2, and 3 are shown in Fig. 9(a). Increasing the void ratio at a given axial strain decreased the differential stress. This mechanical behavior is typical for granular materials. Fig. 9(b) shows the variations of the pore water pressure as a function of the axial strain for specimens 1, 2, and 3 during the biaxial test. Increasing the void ratio increased the maximum positive value of the pore pressure, and the negative pressure appeared at higher axial strains. The specimens with smaller void ratios showed higher negative pressures. These observations can be explained by the tendency of more compacted specimens toward higher dilation. The tendency of the specimens to increase in volume resulted in lower positive and higher negative pressures. Fig. 9(c) shows the variations of ( 1 − 3 )/( 1 + 3 ) as a function of the axial strain for specimens 1, 2, and 3 during the undrained biaxial test. The effective friction angles of the specimens in assemblies 1, 2, and 3 were 26.4◦ , 26.4◦ , and 22.3◦ , respectively, at a 1.6% axial strain, and 22.6◦ , 21.1◦ , and 22.6◦ , respectively, at a 10% axial strain. Comparing the three diagrams shows that the compression of specimens significantly affected the global friction angle of the particle assembly; however, the initial void ratio of the assembly did not significantly affect the global friction angles at high strains. 3.3.3. Effect of confining pressure Simulation Nos. 1, 4 and 5 investigated the effect of the confining pressure applied at the specimen boundaries on the global response of the particle assembly. The results of simulations 1, 4, and 5 were compared. The confining pressures for Test Nos. 1, 4, and 5 were 100, 200, and 300 kPa, respectively. Note that applying different confining pressures to homogenously compact the specimens in the drained state before axial loading changed the void ratios of specimens 4 and 5 relative to specimen 1. Fig. 10(a) shows the variations of the differential stress of specimen Nos. 1, 4, and 5 as a function of the axial strain of the assemblies
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
8
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Fig. 8. Variations of: (a) void ratio, (b) differential stress and (c) effective internal friction angle as a function of the axial strain for the tested specimens during undrained and drained tests.
during the undrained biaxial test. Applying larger confining pressure increased the differential stress at a given axial strain. This result can be explained by the increased inter-granular friction forces caused by larger confining pressures. Obviously, specimens with smaller void ratios were also affected by increasing the confining pressures, which increased the strength of the particle assembly and the differential stress. Fig. 10(b) shows the variations of the pore pressure of specimens 1, 4, and 5 as a function of the axial strain: increasing the confining pressure directly increased the positive pore pressure and
Fig. 9. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for specimens with different initial void ratios during undrained biaxial tests.
decreased the negative pressures. Increasing the confining pressure completely shifted the pore pressure diagram upwards. This result can be explained by the decrease in the tendency of the specimen to expand at higher confining pressures, which increased the pore pressure. Fig. 10(c) shows the evolution of ( 1 − 3 )/( 1 + 3 ) with increasing the axial strain for assemblies 1, 4, and 5. The effective internal friction angle was almost independent of the confining
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Fig. 10. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for specimens with different confining pressures during undrained biaxial tests.
pressure. Such behavior is typically observed in granular materials over this range of confining pressures. 3.3.4. Effect of inter-granular friction The effect of the inter-granular friction coefficient was studied for undrained simulation Nos. 1, 6 and 7. The results of simulations 1, 6, and 7 were compared. The inter-granular friction coefficient for test specimen Nos. 1, 6, and 7 were 0.5, 0.4, and 0.6, respectively.
9
Fig. 11. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for specimens with different friction coefficients during undrained biaxial tests.
Fig. 11(a) shows that the specimens with higher inter-granular friction coefficients exhibited higher differential stresses in the undrained biaxial tests. That is, increasing the inter-granular friction coefficient increased the shear strength of the particle assembly. Fig. 11(b) shows that the positive pore pressures of specimens 1, 6, and 7 do not significantly differ, while their negative pore
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
10
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
pressures show significant differences. This result can be explained by the significant effect of the inter-granular friction coefficient on the volume expansion of the particle assembly, which produced negative pore pressures, whereas contraction of the specimen produced positive pressures that were mainly affected by the elastic properties of the particles. Fig. 11(c) shows the variations of the ( 1 − 3 )/( 1 + 3 ) ratio as a function of the axial strain for assemblies 1, 6, and 7. The peak effective internal friction angles for assemblies 1, 6, and 7 (at 1.6% axial strain) were 26.4◦ , 24.1◦ , and 29.3◦ , respectively, and the residual friction angles (at 10% axial strain) were 22.6◦ , 21.0◦ , and 24.6◦ , respectively. The peak friction angle of the particle assembly clearly depends on the inter-granular friction coefficient, but at large strains (10% axial strain), a converging trend was observed and the difference between the friction angles decreased. Note that the effect of the inter-particle friction on the global friction angle was isolated by using identical parameters, such as the initial void ratios, for all of the specimens. 4. Comparison between coupled and constant volume methods Constant volume tests are often performed to simulate undrained behavior without accounting for the pore pressure. This method is an indirect approach that assumes that the sample volume remains unchanged during the undrained loading. This method has been used by Ng and Dobry (1994) and Sitharam (2003). Bonilla (2004) questioned the validity of using constant volume simulations to model undrained tests. Bonilla (2004) compared a DEM constant volume simulation, which did not account for the pore pressure, to an undrained DEM simulation with fluid coupling for 2D assemblies of elliptical particles. Bonilla found very similar stress–strain behaviors for both cases. 4.1. Macroscopic comparison of the simulation results by coupled and constant volume methods The effect of the simulation method on the overall and local behavior of the specimens was investigated for the results of simulations 1 and 8. Fig. 12(a) shows the differential stress versus the axial strain for simulations 1 and 8. A significant difference existed for axial strain around 6–8%, but nearly identical behavior was observed for axial strains below 4% and beyond 8%. However, the overall behavior of both specimens was similar. Fig. 12(b) also shows that similar trends were obtained for variations of the pore water pressures with the axial strain simulated with the two methods. Fig. 12(c) shows that the variations of ( 1 − 3 )/( 1 + 3 ) with the axial strain was reasonably similar for both simulation methods. The internal friction angle for both methods was almost identical at the peak and the residual states. 4.2. Micromechanical investigation of the simulation results by coupled and constant volume methods The contact orientation and magnitude of the contact force were micromechanically investigated. Fig. 13 illustrates the intergranular force distribution and the distribution of the contact orientation, which were used to investigate the microscopic mechanisms that occurred in specimens during simulations with the coupled and constant volume methods. The results show that the microscopic behavior of the specimen during the two simulations was very similar at small axial strains, while differed at large axial strains. The evolution of the anisotropy during both simulations was almost the same. At the beginning of the shear loading, the particle
Fig. 12. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for the identical undrained specimen, simulated with different (coupled and constant volume) methods.
assembly was isotropic, i.e., the inter-granular forces were almost uniformly distributed in all orientations, as shown in Fig. 13. A small increase in the axial strain dramatically changed the force distribution. The distribution of the contact orientation did not change significantly, and only the magnitude of the forces in the vertical direction increased. This behavior was observed because the particle positions did not change significantly. The specimen appeared to be completely anisotropic at 1% axial strain. Increasing the axial strain (to 6.1% and 10%) resulted in an anisotropy that persisted
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
11
Fig. 13. Micromechanical investigations of the development of anisotropy in the coupled method and constant volume method simulations: (a) polar distribution of the magnitude of the relative contact force (normalized to the maximum force), (b) polar distribution of the number of contacts ( = 5◦ ), (c) contact forces using the coupled method, and (d) contact forces using the constant volume method (the thickness of each line represents the magnitude of the contact force).
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
12
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Fig. 14. Specimens 1 and 8 at 10% axial strain: (a) particle assembly for specimen 1, (b) particle displacements in specimen 1, (c) particle assembly for specimen 8, and (d) particle displacements in specimen 8.
during the remainder of the test. The number of contacts also continuously decreased as the axial strain increased. Fig. 14 shows the assemblies of the particles and particle displacements for specimens 1 and 8, which were simulated using the coupled method and the constant volume method at 10% axial strain. The figure shows that a shear band formed in specimen 8, which was simulated using the constant volume method, and that the specimen failure was concentrated in the shear band. However, the specimen that was simulated using the coupled method failed mainly by bulging in the middle of the specimen. This behavior is also shown in Fig. 13, where the contact force distribution for specimen 8, which was simulated using the constant volume method at 6.1% and 10% axial strain, was more asymmetrical than in the other specimen. Finally, maintaining a constant specimen volume was unlikely to be adequate to comprehensively model the undrained behavior of saturated granular media. The two methods likely produced different results because the constant volume method neglected the pore fluid flow and its effects on the specimen behavior and assumed an average pore fluid pressure for the whole specimen.
5. Conclusions This paper focused on simulating the undrained behavior of saturated granular media, though the drained behavior was also simulated to verify the simulation method. Globally, the specimen behaved similarly to dense sand. At the initial stage of loading, the axial strength of the drained specimen was slightly higher than that of the undrained specimen; however, at large strains for which large negative pressures of pore fluid appeared, the axial strength of the undrained specimen was considerably higher than that of the drained specimen. As expected, similar effective friction angles were observed in both the drained and undrained simulations. The undrained behavior of saturated granular media under a biaxial shear test was simulated using two different methods: the coupled method and the constant volume method. Although the overall behavior of both specimens was similar, the behavior at high axial strains was different because of the simplified assumptions used in the constant volume method. The confining pressure significantly affected the material behavior. Applying larger confining pressures increased the axial strength and the effective friction angle of the specimens and shifted the pore pressure toward positive values. The simulation results showed that the void ratio significantly affected the observed friction angle of the particle assembly.
Increasing the void ratio increased the positive pore pressure and decreased the negative pressure values. The inter-granular friction substantially affected the specimen behavior. Increasing the inter-granular friction coefficient increased the values of the negative pressure. The inter-granular friction also increased the global friction angle of the particle assembly, but other parameters also affected the global friction angle. At the beginning of axial loading, the specimens were isotropic; however, after applying a small axial strain, an extreme increase in the anisotropy was observed. At high axial strains, the contact forces were more uniformly distributed and the anisotropy decreased to some extent. The anisotropy changes of specimens that were simulated using the coupled method and the constant volume method differed; however, the overall evolution of the anisotropy was the same for both methods. References Ashmawy, A. K., Sukumaran, B., & Hoang, V. V. (2003). Evaluating the influence of particle shape on liquefaction behavior using discrete element modeling. In Proceedings of the 13th international offshore and polar engineering conference (ISOPE 2003) (PCW-05) Honolulu, Hawaii, USA, (pp. 542–549). Bardet, J. P., & Sayed, H. (1993). Velocity and attenuation of compressional waves in nearly saturated soils. Soil Dynamics and Earthquake Engineering, 12, 391–401. Bonilla, R. R. O. (2004). Numerical simulations of undrained granular media. Ontario, Canada: Civil Engineering Department, University of Waterloo (Doctoral dissertation). Cundall, P. A. (1971). A computer model for simulating progressive, large-scale movement in blocky rock systems. In Proceedings of International Symposium on Rock Fractures (Vol. 2) Nancy, France, (pp. 129–136). Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Géotechnique, 29(1), 47–65. Dantu, P. (1957). Contribution à l’étude mécanique et géométrique des milieux pulvérulents. In Proceedings of the 4th international conference on soil mechanics and foundation engineering (Vol. 1) London, UK, (pp. 144–148). Das, B. M. (1985). Advanced soil mechanics. Singapore: McGraw Hill. De Josselin de Jong, G., & Verruijt, A. (1969). Etude photo-élastique d’un empilement de disques. Cahiers du Groupe Franc¸ais de Rheologie, 2, 73–86. Desrues, J., Chambon, R., Mokni, M., & Mazerolle, F. (1996). Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography. Géotechnique, 46(3), 529–546. Dinesh, S. V., Sitharam, T. G., & Vinod, J. S. (2004). Dynamic properties and liquefaction behaviour of granular materials using discrete element method. Current Science, 87(10), 1379–1387. Dobry, R., & Ng, T. (1992). Discrete modeling of stress–strain behavior of granular media at small and large strains. Engineering Computations, 9(2), 129–143. El Shamy, U., & Zeghal, M. (2005). Coupled continuum–discrete model for saturated granular soils. Journal of Engineering Mechanics, 131(4), 413–426. El Shamy, U., & Zeghal, M. (2007). A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils. Soil Dynamics and Earthquake Engineering, 27(8), 712–729. Ergun, S. (1952). Fluid flow through packed columns. Chemical Engineering Progress, 48(2), 89–94.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Hakuno, M., & Tarumi, Y. (1988). A granular assembly simulation for the seismic liquefaction of sands. Structural Engineering/Earthquake Engineering, 5(2), 333–342. Jiang, M. J., & Yu, H.-S. (2006). Application of discrete element method to geomechanics. In W. Wu, & H.-S. Yu (Eds.), Modern trends in geomechanics. Springer proceedings in physics (Vol. 106) (pp. 241–269). Berlin: Springer. Lamei, M., & Mirghasemi, A. A. (2011). A discrete element model for simulating saturated granular soil. Particuology, 9, 650–658. Li, X. K., Chu, X. H., & Sheng, D. C. (2007). A saturated discrete particle model and characteristic-based SPH method in granular materials. International Journal for Numerical Methods in Engineering, 72(7), 858–882. Meguro, K., & Ravichandran, N. (2000). 3-Dimensional distinct element simulation of liquefaction phenomena. Monthly Journal of Institute of Industrial Science, University of Tokyo, 52(12), 598–601. Nakase, H., Ishikawa, H., & Fujitani, M. (1996). A simulation study on liquefaction using distinct element method. In Proceedings of the Sixth Japan-U.S. Workshop on Earthquake Resistant Design of Lifeline Facilities and Countermeasures Against Soil Liquefaction MCEER, New York, (pp. 309–318). Ng, T.-T., & Dobry, R. (1994). Numerical simulations of monotonic and cyclic loading of granular soil. Journal of Geotechnical Engineering, 120(2), 388–403. Shafipour, R., & Soroush, A. (2008). Fluid coupled-DEM modeling of undrained behavior of granular media. Computers and Geotechnics, 35(5), 673–685.
13
Sitharam, T. G. (2003). Discrete element modeling of cyclic behavior of granular materials. Geotechnical and Geological Engineering, 21(4), 297–329. Sitharam, T. G., & Dinesh, S. V. (2003). Numerical simulation of liquefaction behaviour of granular materials using discrete element method. Journal of Earth System Science, 112(3), 479–484. Sitharam, T. G., & Vinod, J. S. (2008). Numerical simulation of liquefaction and pore pressure generation in granular materials using DEM. International Journal of Geotechnical Engineering, 2, 103–113. Wakabayashi, T. (1950). Photo-elastic method for determination of stress in powdered mass. Journal of the Physical Society of Japan, 5(5), 383–385. Zeghal, M., & El Shamy, U. (2004). A continuum–discrete hydromechanical analysis of granular deposit liquefaction. International Journal for Numerical and Analytical Methods in Geomechanics, 28(14), 1361–1383. Zhao, J., & Shan, I. (2013). Coupled CFD-DEM Simulation of fluid-particle interaction in geomechanics. Powder Technology, 239, 248–258. Zhu, H. P., Zhou, Z. Y., Yang, R. Y., & Yu, A. B. (2007). Discrete particle simulation of particulate systems: Theoretical developments. Chemical Engineering Science, 62, 3378–3396. Zhu, H. P., Zhou, Z. Y., Yang, R. Y., & Yu, A. B. (2008). Discrete particle simulation of particulate systems: A review of major applications and findings. Chemical Engineering Science, 63, 5728–5770.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
ARTICLE IN PRESS Particuology xxx (2013) xxx–xxx
Contents lists available at ScienceDirect
Particuology journal homepage: www.elsevier.com/locate/partic
Discrete simulation and micromechanical analysis of two-dimensional saturated granular media Younes Khalili, Ahmad Mahboubi ∗ Department of Civil and Environmental Engineering, Power and Water University of Technology, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 19 June 2012 Received in revised form 20 May 2013 Accepted 18 July 2013 Keywords: DEM Saturated granular media Numerical simulation Micromechanical investigation
a b s t r a c t In this study, a novel approach to incorporate the pore water pressure in the discrete element method (DEM) to comprehensively model saturated granular media was developed. A numerical model was constructed based on the DEM by implanting additional routines in the basic DEM code; pore water pressure calculations were used with a two-dimensional (2D) model to simulate the undrained behavior of saturated granular media. This model coupled the interaction of solid particles and the pore fluid in saturated granular media. Finally, several 2D undrained shear tests were simulated. The test results showed that the model could predict the response of the saturated granular soil to shear loading. The effect of initial compaction was investigated. Biaxial tests on dense and loose specimens were conducted, and the effect of the initial density on the change in shear strength and the volume change of the system was investigated. The overall behavior of loose and dense specimens was phenomenologically similar to the real granular material. Constant volume tests were simulated, and the results were compared to those from the coupled model. Induced anisotropy was micromechanically investigated by studying the contact force orientation. The change in anisotropy depended on the modeling scheme. However, the overall responses of the media obtained using the coupled and constant volume methods were similar. © 2013 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.
1. Introduction Granular materials are important to several industrial sectors, such as agriculture and food, pharmaceuticals, process engineering and construction. Soil is a major geological and granular material, which is used to construct roads, dams, and other earth structures because of its properties. The saturated response of a granular medium is of particular importance to practical applications, such as embankment dams. A saturated material is a two-phase mixture: the discrete phase consists of solid particles, and the continuum phase consists of the pore fluid. The interaction between these two phases constitutes the overall behavior of a saturated medium. Therefore, a saturated granular material shows complex behavior. The saturated soil behavior is generally modeled using analytical continuum models. These analytical models assume continuum behavior throughout the multiphase particulate medium. The continuum models are phenomenological and neglect the particulate behavior of granular materials, which certainly affects the deformation and failure modes of the material. Continuum models cannot
∗ Corresponding author. Tel.: +98 21 77312552; fax: +98 21 7700 6660. E-mail addresses: [email protected], ar [email protected] (A. Mahboubi).
explore the physical microscopic mechanisms that lead to liquefaction and can only simulate the macroscopic response if the material behavior is adequately emulated by the adopted constitutive law. An alternative approach to studying the behavior of a granular material accounts for the discrete nature of the particles. Early attempts to observe force distributions and particle behavior at microscopic levels were made by Dantu (1957) and Wakabayashi (1950) using optically sensitive materials. Later, De Josselin de Jong and Verruijt (1969) analyzed the force distribution in such assemblies by studying the individual particles. The results of these experiments provided a qualitative understanding of the mechanisms of load transfer and the deformation characteristics of granular materials. New experimental techniques using computer tomography have also provided important information on void ratio evolution in the fabric (Desrues, Chambon, Mokni, & Mazerolle, 1996). The discrete element method (DEM), which was originally developed by Cundall (1971) to numerically analyze rock blocks, assuming that each rock block satisfies an equation of motion, and later used to simulate an assembly of granular materials (Cundall & Strack, 1979), is a successful computational technique that treats a granular system as a collection of discrete granules. The method allows all of the required information to be extracted at the particle level at any stage of the deformation process. A comprehensive review of developments and applications of
1674-2001/$ – see front matter © 2013 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.
http://dx.doi.org/10.1016/j.partic.2013.07.005
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
2
DEM to simulate particulate systems has been published by Zhu, Zhou, Yang, and Yu (2007, 2008). Although this technique was originally developed for assemblies of idealized particles in a dry environment under shear, it can be extended to account for the effects of a pore fluid. The methods that incorporate fluid effects into DEM can generally be divided into (a) constant volume method and (b) coupled method. Numerical method that indirectly considers fluid effects is known as constant volume method. A constant volume method is based on the almost constant volume of saturated soil specimens under undrained loading because of the very low compressibility of the particles and the interstitial water. The alternative method, which considers fluid-particle interaction, is known as the coupled method. The inter-particle forces are caused by interactions between the solid particles, and the hydrostatic and hydrodynamic forces exerted on the solid particles result from changes in the pore volume and the fluid flow. The fluid flow in porous media is mainly controlled by the gradient of the pore pressure and the coefficient of permeability of the media. The coefficient of permeability of a granular medium depends on the porosity of the medium. The response of an individual particle is governed by the resultant of mechanical forces and hydromechanical forces acting on the particle. A complete review of applications of DEM to geomechanics has been presented by Jiang and Yu (2006). The main applications of DEM to simulate saturated granular materials are briefly explained in Table 1. In this study, a saturated soil model was formulated using the coupled method consisting of two interacting phases to model the mechanical behavior of the saturated granular materials. The elastic membrane in actual shear tests was simulated by a string of boundary particles. The boundary particles were connected to each other, but the boundary could deform freely like an elastic membrane. Imposing constant stress boundary conditions allowed the specimen to fail at the weakest plane. The bulk modulus and the permeability coefficient of each cell of the assembly of saturated particles were calculated at each time step from the particles situations. Constant volume tests were conducted, and the overall results and the local parameters were compared with those of the coupled method tests. A micromechanical investigation was also conducted. The differences between the constant volume method and the coupled method were examined.
The coupled model of the interactions between the solid phase and the fluid phase has three main components: • the fluid phase behavior and pressure generation in the fluid phase, • the solid phase dynamics, and • the interaction between the fluid phase and the solid phase. 2.1. Fluid phase behavior The governing equations for the liquid phase are the continuity and Darcy’s equations, which are given by:
vf = −k
∇ u w
,
2.1.1. Volume change Volume changes are calculated using the particle displacements in the cells. In this study, an indirect approach was used to calculate the volumetric strain from the porosity difference at the beginning and end of each time step: εv (i, j) =
V (i, j) n1 (i, j) − n2 (i, j) = , V (i, j) 1 − n2 (i, j)
(3)
where V(i, j) denotes the volume change during the time step t, V(i, j) denotes the volume of the cell (i, j), and n1 (i, j) and n2 (i, j) are the porosity of the assembly of particles in the cell (i, j) at the beginning and end of each time step, respectively. The porosity in a 2D medium can be defined as the void area divided by the total area. A row and column for a cell are denoted by i and j, respectively. 2.1.2. Pore pressure increment The volume change in cells that are saturated with a pore fluid in the inter-particle spaces can increase or decrease the pore pressure. The intensity of this effect depends on the bulk moduli of the particles and the fluid. Within a linear approximation, the bulk modulus of a saturated granular material can be calculated as: 1 nav (i, j) 1 − nav (i, j) + , = Bp Bw B(i, j)
(4)
where B(i, j) denotes the bulk modulus of the water-particle assembly in cell (i, j), and Bp and Bw are the bulk moduli of the particles and the water, respectively. The average porosity of the particles in cell (i, j) is denoted by nav (i, j) and can be calculated using Eq. (5): nav (i, j) =
n1 (i, j) + n2 (i, j) . 2
(5)
Thus, a change in the fluid volume can change the pore pressure. The increase in pore pressure in each cell is calculated by multiplying the volumetric strain by the bulk modulus of the saturated soil, as shown in Eq. (6).
2. Coupled fluid-particle model
∂n + ∇ · (nvf ) = 0, ∂t
where n is the porosity, vf is the fluid velocity vector, is the gradi is the pressure, and w is the ent operator, k is the permeability, u unit weight of the fluid. Implanting the fluid phase in DEM consists of three main computational steps involving three parameters, the volume change, pore pressure increment and fluid flow. These steps are implemented by dividing the fluid continuum into many cells using vertical and horizontal virtual lines. The fluid parameters and properties, such as the fluid velocity or the fluid pressure, are assumed to be constant in each cell.
(1)
(2)
u(i, j) = εv (i, j)B(i, j).
(6)
This pressure increase must be modified to account for the fluid flow between the cells during the time step. 2.1.3. Fluid flow The volume changes are not equal among the cells, and cells with different bulk moduli have different pressure increases and pore pressures. Thus, a pressure gradient is generated between the cells. The fluid flow depends on the coefficient of permeability of each cell and dissipates the pressure difference between the cells. In this study, Darcy’s law was assumed to be valid. Considering Fig. 1, which illustrates the pore fluid flow between adjacent cells, the continuity equation can be written as: u∗2 (i, j) − u1 (i, j) = B(i, j)
Vw (i, j) , V (i, j)
(7)
where u1 (i, j) denotes the pore pressure of cell (i, j) at the beginning of the time step (t). Considering the fluid flow between cells, the
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model
ARTICLE IN PRESS
PARTIC-603; No. of Pages 13
Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
3
Table 1 A brief list of DEM modeling of saturated granular materials. Model
Application
Remarks
References
Coupled (discrete–continuum model)
Sand liquefaction analysis
2D model; water flow between the pores are modeled by Darcy’s law and the buoyancy and hydrodynamic forces were neglected. Numerical modifications on Hakuno and Tarumi (1988) model
Hakuno and Tarumi (1988)
Liquefaction Liquefaction mechanism Liquefaction and quicksand simulation
3D model based on Hakuno’s model to simulate cyclic shear loading 3D model; using the continuity and Navier–Stokes equations for fluid modeling
Nakase, Ishikawa, and Fujitani (1996) Meguro and Ravichandran (2000) Zeghal and El Shamy (2004) and El Shamy and Zeghal (2005, 2007) Li, Chu, and Sheng (2007)
Biaxial test on saturated closely packed assemblage of particles Behavior of saturated soil mass Sand pile in water
DEM-SPH modeling; monodisperse hexagonal arrangement; the effect of Navier–Stokes and Darcy’s law on specimen behavior was studied. 2D model, laminar flow between the pores (Darcy’s law assumed to be valid) 3D model; CFD-DEM method, using Navier–Stokes equations for fluid modeling
Coupled (discrete–discrete)
Modeling the behavior of granular materials
2D model; fluid is modeled as overlapping cylinders, fluid flow was not considered
Lamei and Mirghasemi (2011)
Constant volume method
Monotonic and cyclic response of a granular soil
2D and 3D simulation of undrained cyclic simple shear test; pore pressure buildup to initial liquefaction and hysteresis loop formation are shown Effect of the particle shape on liquefaction behavior; significant effect of particle morphology on liquefaction susceptibility for sands with the same void ratio Studying the effect of various mechanical parameters on liquefaction
Dobry and Ng (1992) and Ng and Dobry (1994)
Liquefaction
Liquefaction
pore pressure of cell (i, j) at the end of the time step is denoted by u∗2 (i, j) in which the pore pressure increases caused by the cells’ volumetric strains are neglected. The volume (area) of cell (i, j) is denoted by V(i, j), and the total fluid volume that enters cell (i, j) from the adjacent cells is denoted by Vw (i, j), which is calculated using Eq. (8):
q1 =
Ashmawy, Sukumaran, and Hoang (2003) Sitharam (2003), Sitharam and Dinesh (2003), Dinesh, Sitharam, and Vinod (2004), and Sitharam and Vinod (2008)
2k (i, j − 1)k (i, j) u(i, ¯ j − 1) − u(i, ¯ j) x x kx (i, j − 1) + kx (i, j)
qk (i, j).
Zhao and Shan (2013)
j) are shown in Fig. 1. Note that an equivalent permeability coefficient was used for each pair of adjacent cells to account for the different permeability coefficients of adjacent cells. The first term in Eqs. (9)–(12) denotes the equivalent permeability coefficient.
4
Vw (i, j) = t
Shafipour and Soroush (2008)
(8)
q2 =
In Eq. (8), t is the time step, and q1 (i, j) to q4 (i, j) denote the discharges entering cell (i, j), which are calculated according to Darcy’s law as given in Eqs. (9)–(12). The finite difference method was used to calculate the fluid flow. The positive directions of q1 (i, j) to q4 (i,
q3 =
k=1
dy ,
¯ j) ¯ + 1, j) − u(i, u(i dx w dy
(9)
,
2k (i, j + 1)k (i, j) u(i, ¯ j) ¯ j + 1) − u(i, x x kx (i, j + 1) + kx (i, j)
q4 =
2ky (i + 1, j)ky (i, j) ky (i + 1, j) + ky (i, j)
w dx
2ky (i − 1, j)ky (i, j) ky (i − 1, j) + ky (i, j)
w dx
dy ,
¯ − 1, j) − u(i, ¯ j) u(i dx w dy
(10)
(11)
.
(12)
In the above equations, w denotes the unit weight of the pore fluid, dx and dy are the cell dimensions in the x and y directions, and kx (i, j) and ky (i, j) are the coefficients of permeability of cell (i, j) in the x and y directions, respectively. The kx and ky coefficients were assumed to be equal, which is realistic for granular materials. These coefficients are defined by the following equation (Das, 1985): 2
kx (i, j) = ky (i, j) = kinitial
e(i, j) , einitial 2
(13)
where e denotes the void ratio and k denotes the coefficient of permeability of each cell. The initial coefficient of permeability kinitial corresponds to einitial or the average void ratio at the beginning of the axial loading, i.e., axial strain = 0%. The equivalent pore pressure ¯ j), which is calculated as a linduring a time step is denoted by u(i, ear combination of the cell pressures at the beginning and at the end of the time step, as given in Eq. (14): ¯ j) = (1 − ˛)u1 (i, j) + ˛u∗2 (i, j), u(i, Fig. 1. Positive directions of discharges for cell (i, j) shown by q1 (i, j) to q4 (i, j).
(14)
where ˛ is a multiplier between 0.5 and 1. In this study, ˛ was taken to be equal to 0.6. The aforementioned equations were used
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
4
Table 2 Added subroutines. Subroutine No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Fig. 2. Contact model between solid particles.
to calculate u* for the cells; the calculation procedure needed to be repeated several times to reach an acceptable level of accuracy. Finally, the sum of u∗2 (i, j), as given by Eq. (7), and the change in pressure, u(i, j), as given by Eq. (6), yields the total pressure in cell (i, j) at the end of time step u(i, j)(t+t) . u(i, j)
(t+t)
= u(i, j) + u∗2 (i, j).
(15)
DEM models a soil mass as many interacting particles. The positions and velocities of all of the particles in DEM are determined by satisfying two main constraints: mv˙ = ˙ = Iω
f,
(16)
M,
(17)
where m denotes the particle mass, v denotes the particle linear velocity vector, f denotes the force vector acting on the particle, denotes the angular I denotes the particle moment of inertia, ω denotes the moment vector velocity vector of the particle, and M acting on particle. Time derivatives are denoted by dots. The particles were assumed to be rigid with soft contact areas; thus, particles in contact with each other may overlap. The contact law between the particles was implemented using two normal and tangential springs, a non-viscous damper in the normal direction and a slider in the tangential direction, as shown in Fig. 2. 2.3. Interaction between phases The pore fluid exerts interaction forces on particles because of the difference in the velocities of the fluid and the particles and the pressure gradient in the fluid phase. These forces are divided into two main types: hydrostatic and hydrodynamic forces. The hydrostatic force (i.e., the buoyancy force) can be calculated as follows: u(i, j))V p , fb = −(∇ p
(18)
where Vp is the volume of particle p. The hydraulic gradient vector u(i, j) and can be expressed in terms of for cell (i, j) is denoted by ∇ the pressure gradient between adjacent cells:
u(i, j) = ∇
Generate particle assemblage with appropriate characteristics Generate inter-connected boundary particles Divide specimen into cells Calculate cell porosity Calculate cell permeability Calculate volumetric strain of cells Calculate total volumetric strain of particles assemblage Calculate bulk modulus of cells Perform fluid flow calculations Calculate cell fluid velocity Calculate average particle velocity in cells Calculate hydrodynamic forces in X and Y directions Perform calculations of pore pressure using constant volume method Isotropic compaction of generated specimen Perform undrained simulation using coupled method Perform undrained simulation using constant volume method Perform drained simulation
Only drag forces were considered in this study, neglecting all other types of hydrodynamic forces. The drag force can be evaluated using Ergun’s equation as follows (Ergun, 1952):
fd (i, j) = (1−n(i, j))
2.2. Solid phase dynamics
17
Function
u1 (i, j + 1) − u1 (i, j − 1) u1 (i − 1, j) − u1 (i + 1, j) , . 2dx 2dy (19)
f |¯vf (i, j) − v¯ p (i, j)| [1−n(i, j)] 150 f + 1.75 2 ¯ d¯ p n(i, j)dp
[¯vf (i, j) − v¯ p (i, j)],
(20)
where n(i, j) denotes the porosity of cell (i, j), f the pore fluid viscosity, d¯ p the average particle diameter (d50 is used in this study), f the fluid density, and v¯ f (i, j) and v¯ p (i, j) are the average velocities of the fluid and the particles in cell (i, j), respectively. The hydrostatic and hydrodynamic forces acting on each particle were calculated using Eqs. (18)–(20). The coupled model describes the interaction between the fluid and solid phases by applying these forces to each particle at every time step. 3. Simulation with coupled method 3.1. Model implementation A 2D coupled continuum–discrete model of saturated granular soil was implemented using PFC2D (Itasca, USA), which can model the dry assemblies of particles. The subroutines, which were implanted in the code to model the pore fluid in the interstitial voids, calculated the volumetric strain, the velocities of the particles and the fluid, the hydrodynamic forces, the fluid flow and the pore pressure in the cells. Table 2 lists the most important subroutines that were added to the code. A general schematic of the calculation steps to model the saturated granular soil is shown in Fig. 3. The particle positions and particle velocities were specified at the start of the simulation. At each time step, the porosity in each cell was calculated. The bulk modulus and permeability coefficient of the soil in the cells were then updated because the bulk modulus and the permeability coefficient of the soil are porosity-dependent. The volumetric strain in the cells was calculated from the cell porosity. The calculated values for the volumetric strain and the bulk modulus were used to compute the pore pressure increase. This increment is updated following the fluid flow calculations using the previously described finite difference method. The permeability coefficient significantly affects the calculated pore water flow. The total pore pressure was then determined, from which the fluid
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
5
Fig. 3. Calculation procedure.
velocities were calculated. The calculated fluid velocities in the cells and the known values of the particle velocities were used in turn to calculate the drag forces from the Ergun’s equation; the hydrostatic forces were determined from the calculated pore pressure in the cells and the related pressure gradient. The contact forces were calculated from the particle positions using the contact law. The sum of the three types of forces was inserted into the equations of motion, and the particle positions and velocities were calculated as inputs for the next time step.
reduced to generate particles in random positions in a reasonable number of tries without overlapping between adjacent particles. In this study, the radius of any particle was divided by 1.6 in this step. In the second step, the particle radius was increased to its actual value. Particles were then allowed to overlap, producing large interparticle forces at the end of this step. The large inter-particle forces
3.2. Simulation scheme First, a dense specimen of particles was generated. Boundary particles were created at the edges of the generated specimen. The boundary conditions were applied to these boundary particles. The particle assembly was then divided into several cells. Increasing the number of cells increases the accuracy of the simulation but also increases the analysis run time. Therefore, the number of cells was maximized while ensuring a reasonable run time. Increasing the number of particles has the same advantages and limitations. After preparing the soil specimen, a confining pressure was applied to the boundary particles and maintained until the soil specimen was completely consolidated. The sample was then displacement loaded under undrained (or drained) conditions at a constant velocity. The lateral stress was maintained constant using a servo-mechanism to control the velocity of the boundary particles. After reaching a prescribed axial strain of the assembly, the loading was stopped. Note that a high loading velocity affects the simulation results and increases the shear strength parameters of the specimen. From studying the effect of the loading velocity on the overall results, a loading velocity was selected for the simulations that does not affect the global results. A saturated soil specimen with a height of 20 cm and a width of 10 cm was prepared to conduct the biaxial shear test described above. The specimen was divided into 18 pressure cells as shown in Fig. 4. The total number of boundary particles was 300. Table 3 lists the particle properties. Tables 4 and 5 present the properties of the top and bottom walls and the pore fluid (water), respectively. 3.2.1. Assembly generation and specimen preparation DEM uses a specific procedure to generate an assembly of particles. Different approaches can be used to generate the assembly. The generation method used in this study consists of three steps (see Fig. 5). In the first step, a primary assembly of particles with a reduced particle radius was generated. The particle radius was
Fig. 4. Assembly of particles.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model
ARTICLE IN PRESS
PARTIC-603; No. of Pages 13
Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
6
Fig. 5. Assembly generation and compaction: (a) Step I, generating an assembly of particles with reduced particle radii; (b) Step II, increasing particle radii to their actual values, allowing them to overlap; (c) Step III, allowing particles to move into an equilibrium state.
Table 3 Particle properties. Minimum diameter (mm) Maximum diameter (mm) Density (kg/m3 ) Normal stiffness (N/m) Shear stiffness (N/m) Bulk modulus (GPa) Initial permeability coefficient (mm/s)
3.6 5.58 2500 108 108 35 0.01
Table 4 Wall properties. 108 108 0.5
Normal stiffness (N/m) Shear stiffness (N/m) Friction coefficient
generated by particle overlap require that the velocity of particles be artificially reduced. Finally, the particles were allowed to move into an equilibrium state. At the end of the aforementioned steps, the inter-particle forces should be relatively small, and the particle assembly should be sufficiently compacted; otherwise, the number of generated particles must be adjusted. The inter-particle friction angle is an important parameter in the assembly generation. Selecting smaller friction angle values in the particle generation procedure results in the faster and more efficient compaction of the generated assembly. After attaining a desired compacted state, the inter-particle friction coefficient was changed to the actual value. Five different compacted specimens were prepared in this study. The specimens differed in their void ratios and consolidation pressures (see Table 6). The particle diameter distribution of the generated assembly is shown in Fig. 6.
Fig. 6. Particle diameter distribution.
triaxial test. In this study, the boundaries were modeled by a string of connected particles known as boundary particles (see Fig. 7). The connections between the boundary particles were defined such that they could only support tensile forces, i.e., the boundaries could be simply reshaped. This type of boundary has two main advantages compared to a rigid wall boundary. First, a shear band
3.2.2. Boundaries To realistically model a biaxial shear test, the boundaries of the test specimen should be deformable as an elastic membrane in a Table 5 Pore fluid (water) properties (Bardet & Sayed, 1993). Bulk modulus (GPa) Viscosity (Pa s)
2 0.001
Fig. 7. Elastic membrane formed by boundary particles.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model
ARTICLE IN PRESS
PARTIC-603; No. of Pages 13
Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
7
Table 6 Simulation conditions. Test No.
Simulation method
Initial void ratio, e
Confining pressure (kPa)
Inter-granular friction coefficient
Number of inner particles
1 2 3 4 5 6 7 8 9
U-C U-C U-C U-C U-C U-C U-C U-CV D
0.158 0.170 0.176 0.153 0.148 0.158 0.158 0.158 0.158
100 100 100 200 300 100 100 100 100
0.5 0.5 0.5 0.5 0.5 0.4 0.6 0.5 0.5
1228 1220 1209 1228 1228 1228 1228 1228 1228
D: drained and U: undrained test; C: coupled and CV: constant volume method.
can appear in the simulation results as has been observed in laboratory tests, and second, a confining pressure can be applied more uniformly. To facilitate specimen loading, two rigid walls were defined at the top and bottom of the specimen near the boundary particles. The sample was loaded by moving the top wall downwards at a constant velocity. The walls were also used to measure the axial force applied to the specimen. Boundary particles were used between the particle assembly and the rigid walls at the top and bottom of the specimen to transfer forces, i.e., the contact forces and the forces resulting from the pore fluid pressure, from the particle assembly to the rigid walls and to apply the predefined confining pressure to the assembly. 3.3. Simulation results As shown in Table 6, the numerical tests were performed under undrained and drained conditions. The objective of the tests was to numerically verify the effects of the initial compaction, the confining pressure, and the inter-particle friction coefficient on the global macroscopic behavior of the specimen. Simulations 1–7 were conducted using the coupled method to model the pore pressure. Test No. 9 was a drained test, which served as a reference for the drained behavior of an ideal granular material. Test No. 8 was a simulation of undrained loading using the alternative constant volume method. 3.3.1. Effect of pore pressure This section presents the study of the drained specimen behavior by simulation No. 9. The average void ratio and the total volume of the particle assembly in simulation No. 9 were not constant during the test. Fig. 8(a) shows the variation of the void ratio of the particle assembly in simulation Nos. 1–9. The void ratio of the drained specimen first decreased and then gradually increased to a maximum value. This behavior agreed with the experimental behavior of a dense sand specimen in a drained shear test. The void ratios of the undrained specimens (simulation Nos. 1–8) were almost constant, and the volume remained nearly constant during the shear test. Fig. 8(b) compares the variations of differential stress, ( 1 − 3 ) with the axial strain for simulation Nos. 1 and 9. Simulation No. 1 was a reference test to which the other simulations were compared. In the drained simulation (Test No. 9), the shear strength peaked before reaching a residual strength value, as expected from the dense sand behavior. At strains below 2%, a greater shear strength was observed for the drained specimen than that of the undrained test. This phenomenon can be attributed to the positive pore pressure of the undrained specimen (No. 1) at strains below 2%. The dilative behavior of the specimens at strains larger than 2% generated a negative pore pressure in the undrained specimen; thus, a higher differential stress was observed for the undrained specimens than for the drained specimen, as shown in Fig. 8(b).
Fig. 8(c) illustrates the change in the macroscopic effective internal friction angle ( 1 − 3 )/( 1 + 3 ) as a function of the axial strain, obtained from simulations 1 and 9. Both specimens showed similar effective behavior, which agreed well with the observed behavior of granular specimens in laboratory studies. Both diagrams show a peak at approximately 0.44 and a residual value of 0.36 of sin( ), which correspond to 26◦ and 21◦ , respectively, for the acceptable friction angles in dense assemblies of round (spherical) particles. 3.3.2. Effect of the initial void ratio on the behavior of undrained specimens Three identical specimens were prepared with different void ratios to study the effect of the initial void ratio on the specimens’ undrained behavior. The initial void ratios for Test Nos. 1, 2, and 3 were 0.158, 0.170, and 0.176, respectively. The results of the undrained biaxial tests for these specimens are shown in Fig. 9. The variations of the differential stress as a function of the axial strain for specimens 1, 2, and 3 are shown in Fig. 9(a). Increasing the void ratio at a given axial strain decreased the differential stress. This mechanical behavior is typical for granular materials. Fig. 9(b) shows the variations of the pore water pressure as a function of the axial strain for specimens 1, 2, and 3 during the biaxial test. Increasing the void ratio increased the maximum positive value of the pore pressure, and the negative pressure appeared at higher axial strains. The specimens with smaller void ratios showed higher negative pressures. These observations can be explained by the tendency of more compacted specimens toward higher dilation. The tendency of the specimens to increase in volume resulted in lower positive and higher negative pressures. Fig. 9(c) shows the variations of ( 1 − 3 )/( 1 + 3 ) as a function of the axial strain for specimens 1, 2, and 3 during the undrained biaxial test. The effective friction angles of the specimens in assemblies 1, 2, and 3 were 26.4◦ , 26.4◦ , and 22.3◦ , respectively, at a 1.6% axial strain, and 22.6◦ , 21.1◦ , and 22.6◦ , respectively, at a 10% axial strain. Comparing the three diagrams shows that the compression of specimens significantly affected the global friction angle of the particle assembly; however, the initial void ratio of the assembly did not significantly affect the global friction angles at high strains. 3.3.3. Effect of confining pressure Simulation Nos. 1, 4 and 5 investigated the effect of the confining pressure applied at the specimen boundaries on the global response of the particle assembly. The results of simulations 1, 4, and 5 were compared. The confining pressures for Test Nos. 1, 4, and 5 were 100, 200, and 300 kPa, respectively. Note that applying different confining pressures to homogenously compact the specimens in the drained state before axial loading changed the void ratios of specimens 4 and 5 relative to specimen 1. Fig. 10(a) shows the variations of the differential stress of specimen Nos. 1, 4, and 5 as a function of the axial strain of the assemblies
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
8
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Fig. 8. Variations of: (a) void ratio, (b) differential stress and (c) effective internal friction angle as a function of the axial strain for the tested specimens during undrained and drained tests.
during the undrained biaxial test. Applying larger confining pressure increased the differential stress at a given axial strain. This result can be explained by the increased inter-granular friction forces caused by larger confining pressures. Obviously, specimens with smaller void ratios were also affected by increasing the confining pressures, which increased the strength of the particle assembly and the differential stress. Fig. 10(b) shows the variations of the pore pressure of specimens 1, 4, and 5 as a function of the axial strain: increasing the confining pressure directly increased the positive pore pressure and
Fig. 9. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for specimens with different initial void ratios during undrained biaxial tests.
decreased the negative pressures. Increasing the confining pressure completely shifted the pore pressure diagram upwards. This result can be explained by the decrease in the tendency of the specimen to expand at higher confining pressures, which increased the pore pressure. Fig. 10(c) shows the evolution of ( 1 − 3 )/( 1 + 3 ) with increasing the axial strain for assemblies 1, 4, and 5. The effective internal friction angle was almost independent of the confining
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Fig. 10. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for specimens with different confining pressures during undrained biaxial tests.
pressure. Such behavior is typically observed in granular materials over this range of confining pressures. 3.3.4. Effect of inter-granular friction The effect of the inter-granular friction coefficient was studied for undrained simulation Nos. 1, 6 and 7. The results of simulations 1, 6, and 7 were compared. The inter-granular friction coefficient for test specimen Nos. 1, 6, and 7 were 0.5, 0.4, and 0.6, respectively.
9
Fig. 11. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for specimens with different friction coefficients during undrained biaxial tests.
Fig. 11(a) shows that the specimens with higher inter-granular friction coefficients exhibited higher differential stresses in the undrained biaxial tests. That is, increasing the inter-granular friction coefficient increased the shear strength of the particle assembly. Fig. 11(b) shows that the positive pore pressures of specimens 1, 6, and 7 do not significantly differ, while their negative pore
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
10
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
pressures show significant differences. This result can be explained by the significant effect of the inter-granular friction coefficient on the volume expansion of the particle assembly, which produced negative pore pressures, whereas contraction of the specimen produced positive pressures that were mainly affected by the elastic properties of the particles. Fig. 11(c) shows the variations of the ( 1 − 3 )/( 1 + 3 ) ratio as a function of the axial strain for assemblies 1, 6, and 7. The peak effective internal friction angles for assemblies 1, 6, and 7 (at 1.6% axial strain) were 26.4◦ , 24.1◦ , and 29.3◦ , respectively, and the residual friction angles (at 10% axial strain) were 22.6◦ , 21.0◦ , and 24.6◦ , respectively. The peak friction angle of the particle assembly clearly depends on the inter-granular friction coefficient, but at large strains (10% axial strain), a converging trend was observed and the difference between the friction angles decreased. Note that the effect of the inter-particle friction on the global friction angle was isolated by using identical parameters, such as the initial void ratios, for all of the specimens. 4. Comparison between coupled and constant volume methods Constant volume tests are often performed to simulate undrained behavior without accounting for the pore pressure. This method is an indirect approach that assumes that the sample volume remains unchanged during the undrained loading. This method has been used by Ng and Dobry (1994) and Sitharam (2003). Bonilla (2004) questioned the validity of using constant volume simulations to model undrained tests. Bonilla (2004) compared a DEM constant volume simulation, which did not account for the pore pressure, to an undrained DEM simulation with fluid coupling for 2D assemblies of elliptical particles. Bonilla found very similar stress–strain behaviors for both cases. 4.1. Macroscopic comparison of the simulation results by coupled and constant volume methods The effect of the simulation method on the overall and local behavior of the specimens was investigated for the results of simulations 1 and 8. Fig. 12(a) shows the differential stress versus the axial strain for simulations 1 and 8. A significant difference existed for axial strain around 6–8%, but nearly identical behavior was observed for axial strains below 4% and beyond 8%. However, the overall behavior of both specimens was similar. Fig. 12(b) also shows that similar trends were obtained for variations of the pore water pressures with the axial strain simulated with the two methods. Fig. 12(c) shows that the variations of ( 1 − 3 )/( 1 + 3 ) with the axial strain was reasonably similar for both simulation methods. The internal friction angle for both methods was almost identical at the peak and the residual states. 4.2. Micromechanical investigation of the simulation results by coupled and constant volume methods The contact orientation and magnitude of the contact force were micromechanically investigated. Fig. 13 illustrates the intergranular force distribution and the distribution of the contact orientation, which were used to investigate the microscopic mechanisms that occurred in specimens during simulations with the coupled and constant volume methods. The results show that the microscopic behavior of the specimen during the two simulations was very similar at small axial strains, while differed at large axial strains. The evolution of the anisotropy during both simulations was almost the same. At the beginning of the shear loading, the particle
Fig. 12. Variations of: (a) differential stress, (b) pore pressure and (c) effective internal friction angle as a function of the axial strain for the identical undrained specimen, simulated with different (coupled and constant volume) methods.
assembly was isotropic, i.e., the inter-granular forces were almost uniformly distributed in all orientations, as shown in Fig. 13. A small increase in the axial strain dramatically changed the force distribution. The distribution of the contact orientation did not change significantly, and only the magnitude of the forces in the vertical direction increased. This behavior was observed because the particle positions did not change significantly. The specimen appeared to be completely anisotropic at 1% axial strain. Increasing the axial strain (to 6.1% and 10%) resulted in an anisotropy that persisted
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
11
Fig. 13. Micromechanical investigations of the development of anisotropy in the coupled method and constant volume method simulations: (a) polar distribution of the magnitude of the relative contact force (normalized to the maximum force), (b) polar distribution of the number of contacts ( = 5◦ ), (c) contact forces using the coupled method, and (d) contact forces using the constant volume method (the thickness of each line represents the magnitude of the contact force).
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
12
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Fig. 14. Specimens 1 and 8 at 10% axial strain: (a) particle assembly for specimen 1, (b) particle displacements in specimen 1, (c) particle assembly for specimen 8, and (d) particle displacements in specimen 8.
during the remainder of the test. The number of contacts also continuously decreased as the axial strain increased. Fig. 14 shows the assemblies of the particles and particle displacements for specimens 1 and 8, which were simulated using the coupled method and the constant volume method at 10% axial strain. The figure shows that a shear band formed in specimen 8, which was simulated using the constant volume method, and that the specimen failure was concentrated in the shear band. However, the specimen that was simulated using the coupled method failed mainly by bulging in the middle of the specimen. This behavior is also shown in Fig. 13, where the contact force distribution for specimen 8, which was simulated using the constant volume method at 6.1% and 10% axial strain, was more asymmetrical than in the other specimen. Finally, maintaining a constant specimen volume was unlikely to be adequate to comprehensively model the undrained behavior of saturated granular media. The two methods likely produced different results because the constant volume method neglected the pore fluid flow and its effects on the specimen behavior and assumed an average pore fluid pressure for the whole specimen.
5. Conclusions This paper focused on simulating the undrained behavior of saturated granular media, though the drained behavior was also simulated to verify the simulation method. Globally, the specimen behaved similarly to dense sand. At the initial stage of loading, the axial strength of the drained specimen was slightly higher than that of the undrained specimen; however, at large strains for which large negative pressures of pore fluid appeared, the axial strength of the undrained specimen was considerably higher than that of the drained specimen. As expected, similar effective friction angles were observed in both the drained and undrained simulations. The undrained behavior of saturated granular media under a biaxial shear test was simulated using two different methods: the coupled method and the constant volume method. Although the overall behavior of both specimens was similar, the behavior at high axial strains was different because of the simplified assumptions used in the constant volume method. The confining pressure significantly affected the material behavior. Applying larger confining pressures increased the axial strength and the effective friction angle of the specimens and shifted the pore pressure toward positive values. The simulation results showed that the void ratio significantly affected the observed friction angle of the particle assembly.
Increasing the void ratio increased the positive pore pressure and decreased the negative pressure values. The inter-granular friction substantially affected the specimen behavior. Increasing the inter-granular friction coefficient increased the values of the negative pressure. The inter-granular friction also increased the global friction angle of the particle assembly, but other parameters also affected the global friction angle. At the beginning of axial loading, the specimens were isotropic; however, after applying a small axial strain, an extreme increase in the anisotropy was observed. At high axial strains, the contact forces were more uniformly distributed and the anisotropy decreased to some extent. The anisotropy changes of specimens that were simulated using the coupled method and the constant volume method differed; however, the overall evolution of the anisotropy was the same for both methods. References Ashmawy, A. K., Sukumaran, B., & Hoang, V. V. (2003). Evaluating the influence of particle shape on liquefaction behavior using discrete element modeling. In Proceedings of the 13th international offshore and polar engineering conference (ISOPE 2003) (PCW-05) Honolulu, Hawaii, USA, (pp. 542–549). Bardet, J. P., & Sayed, H. (1993). Velocity and attenuation of compressional waves in nearly saturated soils. Soil Dynamics and Earthquake Engineering, 12, 391–401. Bonilla, R. R. O. (2004). Numerical simulations of undrained granular media. Ontario, Canada: Civil Engineering Department, University of Waterloo (Doctoral dissertation). Cundall, P. A. (1971). A computer model for simulating progressive, large-scale movement in blocky rock systems. In Proceedings of International Symposium on Rock Fractures (Vol. 2) Nancy, France, (pp. 129–136). Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Géotechnique, 29(1), 47–65. Dantu, P. (1957). Contribution à l’étude mécanique et géométrique des milieux pulvérulents. In Proceedings of the 4th international conference on soil mechanics and foundation engineering (Vol. 1) London, UK, (pp. 144–148). Das, B. M. (1985). Advanced soil mechanics. Singapore: McGraw Hill. De Josselin de Jong, G., & Verruijt, A. (1969). Etude photo-élastique d’un empilement de disques. Cahiers du Groupe Franc¸ais de Rheologie, 2, 73–86. Desrues, J., Chambon, R., Mokni, M., & Mazerolle, F. (1996). Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography. Géotechnique, 46(3), 529–546. Dinesh, S. V., Sitharam, T. G., & Vinod, J. S. (2004). Dynamic properties and liquefaction behaviour of granular materials using discrete element method. Current Science, 87(10), 1379–1387. Dobry, R., & Ng, T. (1992). Discrete modeling of stress–strain behavior of granular media at small and large strains. Engineering Computations, 9(2), 129–143. El Shamy, U., & Zeghal, M. (2005). Coupled continuum–discrete model for saturated granular soils. Journal of Engineering Mechanics, 131(4), 413–426. El Shamy, U., & Zeghal, M. (2007). A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils. Soil Dynamics and Earthquake Engineering, 27(8), 712–729. Ergun, S. (1952). Fluid flow through packed columns. Chemical Engineering Progress, 48(2), 89–94.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005
G Model PARTIC-603; No. of Pages 13
ARTICLE IN PRESS Y. Khalili, A. Mahboubi / Particuology xxx (2013) xxx–xxx
Hakuno, M., & Tarumi, Y. (1988). A granular assembly simulation for the seismic liquefaction of sands. Structural Engineering/Earthquake Engineering, 5(2), 333–342. Jiang, M. J., & Yu, H.-S. (2006). Application of discrete element method to geomechanics. In W. Wu, & H.-S. Yu (Eds.), Modern trends in geomechanics. Springer proceedings in physics (Vol. 106) (pp. 241–269). Berlin: Springer. Lamei, M., & Mirghasemi, A. A. (2011). A discrete element model for simulating saturated granular soil. Particuology, 9, 650–658. Li, X. K., Chu, X. H., & Sheng, D. C. (2007). A saturated discrete particle model and characteristic-based SPH method in granular materials. International Journal for Numerical Methods in Engineering, 72(7), 858–882. Meguro, K., & Ravichandran, N. (2000). 3-Dimensional distinct element simulation of liquefaction phenomena. Monthly Journal of Institute of Industrial Science, University of Tokyo, 52(12), 598–601. Nakase, H., Ishikawa, H., & Fujitani, M. (1996). A simulation study on liquefaction using distinct element method. In Proceedings of the Sixth Japan-U.S. Workshop on Earthquake Resistant Design of Lifeline Facilities and Countermeasures Against Soil Liquefaction MCEER, New York, (pp. 309–318). Ng, T.-T., & Dobry, R. (1994). Numerical simulations of monotonic and cyclic loading of granular soil. Journal of Geotechnical Engineering, 120(2), 388–403. Shafipour, R., & Soroush, A. (2008). Fluid coupled-DEM modeling of undrained behavior of granular media. Computers and Geotechnics, 35(5), 673–685.
13
Sitharam, T. G. (2003). Discrete element modeling of cyclic behavior of granular materials. Geotechnical and Geological Engineering, 21(4), 297–329. Sitharam, T. G., & Dinesh, S. V. (2003). Numerical simulation of liquefaction behaviour of granular materials using discrete element method. Journal of Earth System Science, 112(3), 479–484. Sitharam, T. G., & Vinod, J. S. (2008). Numerical simulation of liquefaction and pore pressure generation in granular materials using DEM. International Journal of Geotechnical Engineering, 2, 103–113. Wakabayashi, T. (1950). Photo-elastic method for determination of stress in powdered mass. Journal of the Physical Society of Japan, 5(5), 383–385. Zeghal, M., & El Shamy, U. (2004). A continuum–discrete hydromechanical analysis of granular deposit liquefaction. International Journal for Numerical and Analytical Methods in Geomechanics, 28(14), 1361–1383. Zhao, J., & Shan, I. (2013). Coupled CFD-DEM Simulation of fluid-particle interaction in geomechanics. Powder Technology, 239, 248–258. Zhu, H. P., Zhou, Z. Y., Yang, R. Y., & Yu, A. B. (2007). Discrete particle simulation of particulate systems: Theoretical developments. Chemical Engineering Science, 62, 3378–3396. Zhu, H. P., Zhou, Z. Y., Yang, R. Y., & Yu, A. B. (2008). Discrete particle simulation of particulate systems: A review of major applications and findings. Chemical Engineering Science, 63, 5728–5770.
Please cite this article in press as: Khalili, Y., & Mahboubi, A. Discrete simulation and micromechanical analysis of two-dimensional saturated granular media. Particuology (2013), http://dx.doi.org/10.1016/j.partic.2013.07.005