Wei Zhang, Leo Rothenburg

Department of Civil and Environmental Engineering, University of Waterloo, Waterloo, N2L 3G1, Canada

Keywords:Micromechanics Constant volume method Fluid-coupled discrete element method(DEM)Undrained behavior

A B S T R A C T The fluid-coupled discrete element method (DEM) and the constant volume method as two types of discrete modeling methods for fundamental study of undrained responses of granular materials, have been discussed by many researchers. The fluid-coupled DEM, which couples the motions of discrete particles with pore fluid movements, is theoretically robust although it requires a large amount of computation time. As a substitution for the complex fluid-coupled DEM, the constant volume method simulates an undrained condition for a saturated granular material by simply preserving the total volume of a granular assembly without considering interactions between fluids and particles;hence,the validity of its results is questionable. In this paper, the undrained behaviors of granular assemblies simulated using the aforementioned two methods are compared. Based on a comparison of both macroscopic and microscopic responses given by the two methods, it is demonstrated that the constant volume method may reasonably simulate the responses of a loose saturated granular material with very coarse grains,which has a high permeability, and thus a rapid pore pressure equalization. However, it is ineffective in simulating the responses of a loose material with fine components due to its failure to capture the process of a slow dissipation of the excess pore pressure among the individual pores.With regard to the dense material adopted,similar behaviors at the early and intermediate shearing stages given by the two methods are displayed.

1. Introduction

Numerical simulations using discrete element method (DEM)have been used in the past to explore the behavior of granular materials on microscale (Cundall and Strack, 1979; Thornton and Randall, 1988; Ng, 1989; Rothenburg and Bathurst, 1989, 1992;Ting et al., 1989; Bathurst and Rothenburg, 1992; Dobry and Ng,1992; Iwashita and Oda, 2000; Thornton, 2000; McDowell and Harireche, 2002; Sitharam et al., 2002; Yimsiri and Soga, 2010). In an undrained condition, a saturated granular medium is usually simulated by coupling the effect of fluid on its particles, the socalled fluid-coupled DEM (Hakuno and Tarumi, 1988; Nakase et al., 1999; Olivera, 2004; Zeghal and El Shamy, 2004; Boutt et al., 2007; El Shamy and Zeghal, 2007; Okada and Ochiai, 2007;Shafipour and Soroush, 2008; Shimizu, 2011; Han and Cundall,2013; Lomine et al., 2013; Catalano et al., 2014; Liu et al., 2015),or by preserving the total volume of the granular material assuming that the pore fluid is incompressible, which is denoted “the constant volume method”(Thornton and Barnes,1986;Ng and Dobry,1994; Dubujet and Dedecker, 1998; El-Metskawy, 1998; Cheng et al., 2003; Sitharam, 2003; Sitharam and Dinesh, 2003; Gong et al., 2011; Gong, 2015). Due to the advantage of saving computation time by eliminating complex fluid coupling calculations,the constant volume method was adopted in many studies and results from these studies reasonably match the corresponding laboratory results(Ng and Dobry,1994;Sitharam,2003;Sitharam and Dinesh,2003;Gong et al.,2011;Gong,2015).However,because it does not evaluate the excess pore pressure and its dissipation, this method lacks a sound theoretical basis. In comparison, the fluid-coupled DEM mechanically couples the motions of granular particles with the movements of pore fluids and thus can depict the undrained response of a granular material more realistically.

Fig.1. Flow chart of DEM.

The fluid-coupled DEM used in the past studies falls into three categories: the continuum-based method, the Lattice Boltzmann method, and the pore fluid network method. Representative studies using the continuum-based method are those of Shimizu(2004), Zeghal and El Shamy (2004), El Shamy and Zeghal (2005),El Shamy and Zeghal (2007), and Suzuki et al. (2007). In their studies,computation of the fluid motion is based on the continuity and Navier-Stokes equations and follows the idea of spatial average quantities over the fixed cells proposed by Anderson and Jackson (1967). This method has the advantage of saving computational cost. However, since it cannot specifically depict the local pore pressure variations and interactions between particles and fluid on a pore scale, it is not applicable to the micromechanical study of the responses of saturated granular materials. Different from the continuum-based method,the Lattice Boltzmann method is based on a very fine discretization of the fluid phase. Representative studies are those of Boutt et al. (2007), Han and Cundall(2011), Han and Cundall (2013), Wang et al. (2016), and Wang et al. (2018). This method is effective in simulating the response of a granular material at the pore scale;however,due to its very fine discretization, this method is only applicable to a granular sample with a limited number of particles. As opposed to the former two methods, the pore fluid network method identifies the pore pressure changes of each pore, and it considers the fluid transfer between each pair of interconnected pores. Representative studies using this method are those of Hakuno and Tarumi(1988),Hakuno(1995), Olivera (2004), Shimizu (2011), Chareyre et al. (2012),Catalano et al. (2014) and Zhang (2018).

In this paper, the individual pore-based fluid-coupled DEM (a pore fluid network method) proposed by Zhang (2018) is used to conduct undrained simulations. Zhang’s model is built on a preceding model proposed by Olivera (2004) and is verified by comparing it with Biot’s linear poroelastic continuum (Zhang,2018). In Olivera’s model, the range of application is limited to cases in which contact creation and disintegration do not induce drastic modifications of the pore space, i.e. undrained shear deformation in a relatively dense granular sample.This becomes a practical limitation when considering liquefaction problems. During the occurrence of static liquefaction in a loose sample, very complex pore space modifications occur,and usually affect several pores simultaneously. Zhang’s model includes a “pore group”scheme, which removes this limitation and is capable of overcoming the difficulty of preserving a fluid mass balance during dramatic pore structure changes in the strain-softening phase;therefore, it is employed in this study. Both the conventional stress-strain relationships and the micromechanical responses computed using Zhang’s model are compared to those given by the constant volume method.The similarities and differences between the undrained behaviors given by the two methods are presented,and an interpretation of the importance of the pore pressure calculation in an undrained DEM simulation is proposed.

2. Discrete element method

The DEM was first proposed by Cundall and Strack (1979) to simulate micromechanical behavior of an assembly of granular particles based on Newton’s second law and a force-displacement relationship.In this paper,calculation of the motion of each particle and simulation of the interaction between contacting particles adopt a two-dimensional (2D) DEM that is similar to the DISC proposed by Bathurst (1985). Fig. 1 shows a flow chart of this approach. In each computational cycle, a force-displacement relationship is first applied to the contact overlaps to calculate the contact force between each pair of particles in contact with each other. Afterwards, the Newton’s second law is employed to conduct a motion computation for each particle, which results in new displacements and contact overlaps for the particles in the assembly and thereby forms a computational loop for the following cycles.

Fig. 2. Linear contact interaction between two particles in contact.

In this model,surface interpenetration is introduced to present a contact between a pair of contacting particles. The inter-granular contact force between them is determined by considering a linear contact interaction, as schematically shown in Fig. 2. Springs with normal and tangential stiffnesses knand ktare employed to compute the normal and shear contact forces fnand ftbetween particles A and B,given the displacements Δnand Δtin the normal and tangential directions, respectively.

To maintain the static equilibrium and dissipate the kinetic energy of the granular system, both friction damping and viscous damping are incorporated in this model as introduced by Cundall and Strack (1979). The former is applied to limit the inter-particle shear contact force, which is switched on only when sliding between particles occurs. The latter can be further extended into contact damping and global damping, both of which are modeled by dashpots. Normal contact damping (Dn) and tangential contact damping(Dt)operate on each contact point and can be represented as

where vnand vtare the normal and tangential relative velocities of a pair of contacting particles, respectively; cnand ctare the corresponding contact damping coefficients which are proportional to the normal and tangential contact stiffnesses knand ktin the following manner:

where β is a coefficient of proportionality.

Different from contact damping,global damping is achieved by connecting each particle to a fixed reference through envisaging a dashpot. It acts on each particle’s velocity via the coefficients cmand cI,which are related to the mass m and moment of inertia I of a particle, respectively:

where α is a coefficient of proportionality, and ω is applied to amplify the impact of the rotational damping component.

In this study, the time step Δt proposed by Strack and Cundall(1978) is employed, which is a fraction of the critical time step Δtcthat can be estimated by considering a single degree of freedom mass-spring system:

where Cfracis the fraction,which can be assigned as 0.1 as indicated by Strack and Cundall(1978);mminis the lowest value of all particle masses in a granular system; and kmaxis the highest value of normal or tangential contact stiffness (knor kt) in the granular system.

3. Constant volume method

The constant volume method is in essence a DEM applied to dry particles. However, by maintaining a constant total volume of a granular assembly and assuming that both the solid particles and the fluid are incompressible, an undrained condition can be simulated and the undrained response of fully saturated granular media can be obtained. Traditionally, in a 2D case, a constant volume can be preserved by applying equal values but opposite signs of the constant strain rate at the boundaries of a granular assembly:

Fig. 3. Flow network for a set of polygons: (a) Fluid flow network, and (b) Conduit pipe.

Fig. 4. Pore structure evolution example: (a) Initial pore structure, and (b) Pore structure in the next cycle.

4. Individual pore-based fluid-coupled DEM

4.1. Fluid-particle coupling

The individual pore-based fluid-coupled DEM proposed by Zhang (2018) is used in this study. This method includes three components: the particle motion calculations, the pore fluid movement computations, and the particle-fluid coupling computations. Calculation of the motion of particles in this method uses the DEM that is similar to the DISK proposed by Bathurst(1985)as introduced above. The pore fluid movement computations in this method inherited the idea proposed by Hakuno and Tarumi(1988)and Olivera(2004).In addition,a“pore groups”scheme is added to overcome the difficulty of preserving fluid mass balance during dramatic pore structure changes in the strain-softening phase.Specifically,in this method,pore pressures are generated due to the deformation of the pores in such a way that the voids in the granular system are treated as being fully saturated with fluid which is treated as visco-elastic media. Thus, the pore pressure change Δuiin a void i in the granular system during any computational step can be written as

where Bfis the bulk modulus of the fluid media,Viis the volume of void i in the previous step, and ΔViis the volume change of void i from the previous step. ΔViin this equation comes from both volume change due to particles movements ΔVpi and the input or output of the fluids from the interconnected pores.The former can be directly computed by tracking the geometry of a pore through two successive cycles. Each pore is detected by connecting the centers of contacting disc particles, which forms polygons enclosing the pore, as indicated in Fig. 3a. The latter can be solved using the Hagen-Poiseuille theory considering that an interconnected flow network can be constructed by connecting the centroid of gravity of each pore(see Fig.3a).Hence,the volumetric flow rate q between any pair of interconnected pores is proportional to their pore pressure difference:

where u1and u2are the pore pressures of the two interconnected pores,μ is the dynamic viscosity of the fluid media,and d and L are respectively the diameter and length of the conduit pipe that connects any two interconnected pores (see Fig. 3b). Hence, an expression for volume change in void i,ΔVi,within time step Δt can be represented as

where Δqjrepresents the volumetric flow rate of the fluids in conduit j,and N is the total number of pores in the assembly.Next,combining Eqs.(9)-(11),a system of ordinary differential equations can be obtained as

Fig.5. Flow network constructed by conduit pipes used to determine global coefficient of permeability k.

The above system of first-order differential equations can be solved using a suite of codes “RKSUITE” that was built based on the Runge-Kutta method (Brankin et al., 1991). Afterwards, uisolved from Eq. (12) can be transformed into the corresponding pore pressure forces for the immediately nearby particles via multiplying it by the associated affecting area. For each particle,adding all the pore pressure forces from the surrounding pores to its resultant force derived from motion calculations will result in a new resultant force, which will participate in the forcedisplacement computation of DEM in the following cycle.

4.2. “Pore groups” scheme

To successfully conduct individual pore-based computations using the above equations,it is essential to maintain the fluid mass balance, particularly during dramatic pore structure changes.During the strain-softening phase,pore coalescence due to contact loss and pore subdivision due to contact creation occur simultaneously, which affects the composition of groups of pores. The“pore groups” scheme proposed by Zhang (2018) continuously tracks pore space modifications in any two successive computational steps.It identifies the sub-volumes containing pore groups in the current time step and one-to-one maps them into the subvolumes uniquely identified in the previous cycle. Therefore, a fluid mass balance is accurately guaranteed during any dramatic microstructure changes. Fig. 4 demonstrates a simplified scenario of this scheme.Two pores I and II are initially formed by 15 particles labeled sequentially from 1 to 15 during a time step, as shown in Fig.4a.In the next cycle(see Fig.4b),due to particle movements,a pore structure with three pores III,IV and V results from a contact disintegration between particles 6 and 7, and two newly created contacts between particles 3 and 4 and between particles 6 and 12.Since pores III,IV and V are newly derived from pores I and II,their ΔVivalues cannot be calculated individually.Thus,they have to be grouped together and their volume summation is compared with the volume summation of pores I and II to compute ΔViso that a fluid mass balance is guaranteed. In this case, pores III, IV and V share the same volumetric strain. The “pore groups” scheme exhibits its advantage when modifications of pore space involve several contact creations and disintegrations taking place simultaneously that affect groups of pores, which is common during a strain-softening phase.

4.3. Selection of conduit diameter

To effectively use the above fluid flow coupling technique in Zhang’s model, a reasonable choice of conduit diameter d is necessary.As indicated by Olivera(2004)and Zhang(2018),this can be achieved by constructing a relationship between the conduit diameter d and the global coefficient of permeability k on the basis of Darcy’s law. Fig. 5 shows schematically a simulation conducted on a rectangular assembly aimed at testing the global coefficient of permeability k.The assembly has a thickness of D,which is equal to the average diameter of the particles in the assembly.Conduit pipes connecting interconnected pores introduced in Fig. 3 are represented as double thick lines in the rectangular assembly.Initially,a pressure gradient is assigned uniformly in the vertical direction throughout the rectangular assembly, which induces fluid flow through conduit pipes from the top boundary to the bottom boundary. Subsequently, the fluid flow system is brought to a steady state in which the amounts of fluid flow into the top boundary and out of the bottom boundary are equal. Then the volumetric flow rate into the rectangular assembly Qinat the top boundary can be represented by the summation of the volumetric flow rates into the individual conduit pipes on the top boundary qn:

According to Darcy’s law, Qincan also be represented using the coefficient of permeability k,the hydraulic gradient i,and the crosssectional area A:

Thus,based on Eqs.(13)and(14),a relationship between k and qncan be established:

Since qnis related to the conduit diameter d given by the Hagen-Poiseuille theory (Eq. (10)), a relationship between k and d can be constructed by repeating the same simulation but using different values of d. Fig. 6 shows the resultant relationship between k and d from these simulations and the regime of d for different types of soils based on their range of k values given by Das (2010). It can be seen in Fig. 6 that the global coefficient of permeability k is very sensitive to the changes in conduit diameter d when d is low, and a small increase in d can induce a significant increase in the global coefficient of permeability k.It is worthy of note that the relationship between k and d presented in Fig. 6 is exclusive to the network constructed in this study.Based on Fig.6,two values of d=80 μm and 320 μm will be used in the following simulations to compare the responses of granular media subject to undrained shearing given by the two methods and to study the effect of permeability on undrained responses.

Fig.6. Relationship between global coefficient of permeability k and conduit diameter d.

Fig. 7. Granular assembly samples with different void ratios before undrained shearing: (a) Sample D, and (b) Sample A.

5. Micromechanical descriptors

A comparison study of micromechanical responses of a granular material subjected to shearing given by the two methods was performed in this study. Three descriptors that are exclusive to micromechanics and closely related to the structure and fabric of a granular material are adopted for analysis. They are the average coordination number γ, the contact normal anisotropy parameter an,and the normal contact force anisotropy parameter af. In the following, a brief introduction of them is presented.

5.1. Average coordination number γ

The average coordination number γ is the average number of immediate neighboring particles that each particle is in contact with in an assembly.It can be mathematically formulated using two times the number of physical contacts M divided by the number of particles N in an assembly:

In this study,only the particles that carry load are considered in computing the average coordination number (mechanical coordination number), and particles with zero or one contact with their immediate neighboring particles(named floaters in this study)are excluded from the calculation. Based on the studies performed by Athanasiou-Grivas and Harr (1982) and Smith et al. (1929), the average coordination number is strongly correlated with the void ratio, and generally the higher the average coordination number,the lower the void ratio.

5.2. Contact normal anisotropy parameter an

The contact normal anisotropy parameter andescribes the intensity of the contact normals in the principal directions of contactnormal anisotropy in a granular system. It is a parameter taken from the truncated Fourier series approximation of the contact normal distribution function S（θ） proposed by Konishi (1978) and Rothenburg (1980):

Table 1 Physico-mechanical properties of the particles and fluid.

Table 2 Parameters for the undrained testing.

Table 3 Summary of initial physical properties of the tests.

where anand bnare the second- and fourth-order coefficients of contact normal anisotropy, respectively; and θaand θbare the second- and fourth-order principal directions of the contact normal anisotropy in a granular system, respectively. Typically,the larger the value of an, the higher the level of contact normal anisotropy: a larger difference of intensity of the contact normals falls into the major and minor principal stress directions.

5.3. Normal contact force anisotropy parameter af

Similar to the anintroduced above, the normal contact force anisotropy parameter afrepresents the extent to which the normal contact forces fall into the principal directions of the normal contact force anisotropy.The parameter afis from the truncated Fourier series approximation of the normal contact force distribution functionproposed by Rothenburg (1980):

Fig. 8. Comparisons of stress-strain responses using the two methods: (a) Samples A and B, and (b) Samples C and D.

6. Summary of the test program

Fig.9. Comparisons of pore pressure variations using the two methods:(a)Samples A and B, and (b) Samples C and D.

A comparison of the undrained behavior of a granular material given by the two aforementioned methods was conducted in this study.Four samples A,B,C and D were subjected to both methods of analysis to make the comparison.The void ratio of the four samples varies significantly sothat bothstrain-softening and strain-hardening types ofresponse can beachieved.The four samples were prepared by initially generating a sample(sample D)with about 10,000 particles(see Fig. 7a of sample D and its particle size distribution). In this original sample, a great many floaters were formed by alternately increasing and decreasing the particle diameters during the initial sample preparation process.Once the initial sample D was prepared,samples A, B and C were created by removing different amounts of floaters from the pore space of sample D while maintaining the main structure.Specifically,one floater per pore was removed to produce sample C,two floaters per pore were removed to form sample B,and all the floaters inthe pore space were removed to create sample A(see Fig.7b).Inthisfloater removal process,the particlesize distributionof samples A, B and C may vary slightly from the original sample D;however, it will not change significantly due to the randomness of floater formation and deletion. Afterwards, an equilibrium state is introduced in all four samples,which induces small variations of the number of floaters in the samples. Different from the traditional method of preparing a loose sample using a large inter-particle friction value during the process of sample generation, this floater generation-removal method allows the preparation of a granular assembly with a much higher void ratio to guarantee the exhibition of strain-softening and a static liquefaction type of behavior when it is subjected to shearing.

Particles in the four samples are assumed to be quartzite(Franklin and Dusseault, 1991), and their physico-mechanical properties are listed in Table 1. The fluid in the simulations using the individual pore-based fluid-coupled DEM is taken to be water at 20°C, and its physico-mechanical properties are also listed in Table 1.The conduit diameter d is selected as 80 μm in this group of simulations to represent the hydraulic conductivity in a granular assembly with fine components based on Fig. 6. All four samples are prepared with an initial confining stress σ11equal to 100 kPa.Simulations using the constant volume method are conducted under a constant strain rate mode,which applies the same values of constant strain rate at the boundary but with different signs as indicated by Eq. (8). Simulations using the fluid-coupled DEM are performed under an undrained biaxial compression mode, which maintains the average stress σ11constant, and increases the average value of stress σ22by applying a constant strain rateat the boundary. The relevant parameters for conducting both simulations are presented in Table 2. Table 3 shows both macroscopic and microscopic initial conditions of the four samples.

Fig.10. Comparisons of the stress paths using the two methods: (a)Samples A and B,and (b) Samples C and D.

Fig.11. Comparisons of the stress ratio variations using the two methods:(a)Samples A and B, and (b) Samples C and D.

7. Simulation results

7.1. Comparison of the macroscopic undrained behaviors

Comparisons of the macroscopic undrained behaviors with respect to the stress-strain relationship, the pore pressure variations, the stress path, and the stress ratio changes of the four samples subjected to shearing given by the two methods are presented in Figs.8-11,respectively.It can be seen from these figures that for both methods, samples A and B both exhibit a strainsoftening type of behavior, while samples C and D both present a strain-hardening type of behavior. For ease of comparison in the following, the results of samples A and B are presented together and are shown in Figs. 8a, 9a, 10a and 11a, while the results of samples C and D are grouped together and are presented in Figs.8b,9b,10b and 11b.

7.1.1. Stress-strain response

Fig. 8 presents the stress-strain response of the four samples given by the two methods. It can be seen in the figure that the overall stress-strain responses derived by the two methods are similar for samples A and B. However, for both samples, the peak strength values obtained using the individual pore-based fluidcoupled DEM are higher than those given by the constant volume method.In addition,sample A,which has a higher void ratio,shows a larger discrepancy between the peak strength values. At the end of the strain-softening stage,the shear stresses obtained using the individual pore-based fluid-coupled DEM fall below those obtained using the constant volume method for both samples A and B.Thus,the residual strengths obtained using the individual pore-based fluid-coupled DEM are lower than those obtained using the constant volume method in the steady-state.

Unlike the responses for samples A and B,there is a high level of consistency between the stress-strain relationships obtained for samples C and D using the two methods. Especially at the early to intermediate shearing stage, shear stresses given by the two methods are almost the same for both samples C and D.Afterwards,along with a continuous progress of strain-hardening, the shear stresses given by the individual pore-based fluid-coupled DEM increase at lower rates than those derived using the constant volume method.Above about 2%shear strain in sample C and 5%shear strain in sample D,the shear stresses given by the constant volume method become noticeably higher than those given by the fluidcoupled DEM, and this trend continues to the end of simulations.

Fig.13. Comparisons of the number of floaters nf using the two methods: (a)Samples A and B, and (b) Samples C and D.

From this comparison,it is recognized that there are noticeable discrepancies between the shear stress values given by the two methods for the loose samples at the early and intermediate shearing stages, which are not exhibited by the dense samples.However, beyond the intermediate shearing phase, all samples exhibit the same response as the shear stress obtained using the individual pore-based fluid-coupled DEM becomes lower than that obtained using the constant volume method. Therefore, the constant volume method results in a higher value of shear strength at a large strain.

7.1.2. Pore pressure variations

Fig. 14. Comparisons of the contact normal anisotropy parameter an using the two methods: (a) Samples A and B, and (b) Samples C and D.

A comparison of the variations of the pore pressure of the four samples obtained using the two methods are presented in Fig.9.In this comparison, pore pressure changes given by the individual pore-based fluid-coupled DEM are calculated directly from the simulation, while pore pressure changes obtained using the constant volume method are indirectly evaluated by subtracting the current effective confining stress value from the initial effective confining stress at the beginning of shearing as indicated by Shafipour and Soroush (2008), Yimsiri and Soga (2010), and Keishing and Hanley (2018). It can be seen in Fig. 9 that the pore pressure variations for all four samples obtained using both methods match their corresponding stress-strain response presented in Fig.8.Both samples A and B exhibit an initial agreement of pore pressure generations during the early shearing stage, followed by a gradual deviation between the two methods with continued shearing.Beyond about 1%shear strain for sample A and 2%shear strain for sample B,the pore pressures obtained using the individual pore-based fluid-coupled DEM become higher than those given by the constant volume method. The differences between the pore pressure values given by the two methods further enlarge with continued shearing for both samples A and B, until large strains are reached and the processes of strain-softening have been completed. The difference between the pore pressure values given by the two methods for samples A and B then stays at an almost constant value, indicating the achievement of a steadystate.

Fig. 15. Comparisons of the normal contact force anisotropy parameter af using the two methods: (a) Samples A and B, and (b) Samples C and D.

Table 4 Execution time for the four samples using the two methods.

Different from samples A and B, discrepancies of pore pressure obtained using the two methods are much smaller for samples C and D, particularly at the early to intermediate stages, they are nearly superimposed. Thereafter, from about 2% shear strain for sample C and 5%shear strain for sample D,the pore pressures obtained using the individual pore-based fluid-coupled DEM exhibit markedly higher values than those given by the constant volume method until the end of simulation.Sample C reached a steady-statewithanalmost constant pore pressure value at a large strain,while the pore pressure of sample D was still increasing by the end of simulation,indicating that a strain-hardening behavior continued.

Based on above comparison of pore pressure variations given by the two methods for the four samples, it is shown that individual pore-based fluid-coupled DEM produces higher values of pore pressuresatlarge strains forbothlooseanddensesamples,althoughatthe early stage of shearing, pore pressure values given by the two methods are very close to each other. From samples A-D,the strain where a deviation of pore pressure values given by the two methods starts to be exhibited becomes higher and higher,which agrees with the corresponding stress-strain responses shown in Fig.8.

7.1.3. Stress paths

The stress paths obtained using the two methods generally exhibit similar patterns (see Fig. 10); however, obvious discrepancies can be seen when the sample is loose. For samples A and B, the stress paths given by the constant volume method are located below the stress paths given by the individual pore-based fluid-coupled DEM throughout the course of the simulations. But the residual strengths and the mean effective stresses at the steady-state given by the constant volume method are higher for both samples A and B. On the other hand, stress paths for samples C and D obtained using the two methods superimpose except that results obtained from the constant volume method exhibit higher shear strengths and larger mean effective stresses at large strains. Therefore, based on Fig. 10, in sequence from samples A-D, the stress paths computed by the individual porebased fluid-coupled DEM gradually change from exhibiting higher stress levels than those given by the constant volume method to exhibiting responses close to those given by the constant volume method. However, for all four samples, the shear stress and mean effective stress values at the end of simulations are higher for the constant volume method.

Fig.16. Comparisons of stress-strain responses using the two methods(d=320 μm):(a) Samples A and B, and (b) Samples C and D.

Fig.17. Comparisons of pore pressure variations using the two methods(d=320 μm):(a) Samples A and B, and (b) Samples C and D.

7.1.4. Stress ratio variations

The stress ratio q/p′variations for the four samples given by the two methods are shown in Fig.11. It can be seen in the figure that similar to the stress path responses shown in Fig. 10, distinctive differences between the stress ratios are presented when the sample is loose. At the initial shearing stage, samples A and B exhibit similar stress ratio responses.Then for both samples A and B, the stress ratio obtained using the individual pore-based fluidcoupled DEM increases at a much higher rate than that obtained using the constant volume method. The difference between the stress ratios given by the two methods continues to enlarge until a steady-state is achieved at a large strain, where the difference remains constant. Contrastingly, the stress ratio responses for samples C and D given by the two methods almost merge with each other throughout the whole course of simulation,which agree with their associated stress path changes presented in Fig.10. For both samples C and D, the stress ratios computed using both methods are nearly the same in all states.

Fig. 18. Comparisons of the stress paths using the two methods (d = 320 μm): (a)Samples A and B, and (b) Samples C and D.

From this comparison, we see that the stress ratios, q/p′, at the steady state given by the two methods are very different when the void ratio of the sample is high. However, there is a tendency of reduction in the differences between the effects of the two simulation methods on the steady state stress ratio q/p′as the void ratio of the sample decreases.

7.2. Comparison of the microscopic responses

Comparisons of the responses of the three micromechanical descriptors in addition to the number of floaters as given by the two methods for the four samples are displayed in Figs. 12-15,respectively. Fig.12 shows the responses of the average coordination number γ that resulted from the two methods.Fig.13 presents the variations of the number of floaters nfgiven by the two methods. Fig.14 shows the changes in the contact normal anisotropy parameter anobtained from the two methods, and Fig. 15 shows a comparison of variations of the normal contact force anisotropy parameter afobtained using the two methods. Again,for ease of comparison, results for samples A and B are presented together, while results for samples C and D are shown together.

Fig.19. Comparisons of the stress ratio q/p′ using the two methods (d = 320 μm): (a)Samples A and B, and (b) Samples C and D.

Fig. 20. Comparisons of the average coordination number γ using the two methods(d = 320 μm): (a) Samples A and B, and (b) Samples C and D.

7.2.1. Comparison of the microscopic behaviors of samples A and B

The overall micromechanical responses of samples A and B given by the two methods are similar for the four descriptors.At the initial shearing stage,the values of all four descriptors given by the two methods are very close to each other for both samples. Subsequently, values of all descriptors computed using the two methods start to diverge. The average coordination number γobtained using the individual pore-based fluid-coupled DEM falls below that given by the constant volume method, while the number of floaters nf,the contact normal anisotropy parameter an,and normal contact force anisotropy parameter afcomputed using the individual pore-based fluid-coupled DEM all fall above those given by the constant volume method.In addition,sample A shows a larger difference than sample B between all four descriptors given by the two methods. With continued shearing, the difference in values given by the two methods continues to increase for all four descriptors, until a steady-state is achieved when the differences given by the two methods for all four descriptors stay constant.

Fig.21. Comparisons of the number of floaters nf using the two methods(d=320 μm):(a) Samples A and B, and (b) Samples C and D.

The differences after the initial shearing stage given by the two methods for all four descriptors for samples of A and B are due to differences in the excess pore pressure dissipation rates given by the two methods. The constant volume method simulates a scenario in which the excess pore pressures of the individual pores equalize instantly, whereas in the individual pore-based fluidcoupled DEM employed in this study, the excess pore pressure dissipation rate depends on the permeability of the granular assembly. In this group of simulations, due to the presence of fine particles in the granular assembly, a relatively low conduit diameter is assigned to the model to simulate a flow network with relatively low connectivity.Therefore,the equalization of the excess pore pressure of individual pores takes some time. It is shown in Fig.9a that the pore pressure in both samples A and B increases to a high value after the initial shearing stage. Thus, the local high excess pore pressure generated cannot dissipate instantly when the connectivity of the flow network is relatively low. This induces a large number of contact disintegrations in both directions during the strain-softening phase compared to those given by the constant volume method.

Fig. 22. Comparisons of the contact normal anisotropy parameter an using the two methods (d = 320 μm): (a) Samples A and B, and (b) Samples C and D.

The above suggestion that a larger number of contact losses occur after the initial shearing stage given by the individual porebased fluid-coupled DEM model is supported by the corresponding responses of the micromechanical descriptors at this stage.It is shown in Figs.12a,13a and 14a and 15a that after an initial shearing stage, both samples A and B exhibited a more rapid decrease of γand more rapid increases of anand afwhen computed using the individual pore-based fluid-coupled DEM than that when the constant volume method was used. These responses indicate a larger number of contact disintegrations occurred at this stage when the individual pore-based fluid-coupled DEM was used.Specifically,there is a larger difference between the contact normal intensities and normal contact force intensities falling into the major and minor principal stress directions(in this case,the vertical and horizontal directions) for the individual pore-based fluidcoupled DEM as indicated by the higher values of anand afthan that of the constant volume method.Therefore,a larger number of contact disintegration occurred in the horizontal direction for the individual pore-based fluid-coupled DEM than that of the constant volume method.This indicates that its relatively low pore pressure dissipation rate induces a local high pore pressure, which pushes the surrounding particles away mainly in the horizontal direction with the continuous application of a vertical load. Additionally, in sample A,whose pore pressure is higher than sample B,the impact of a local high pore pressure pushing its surrounding particles away is greater; therefore, it increases the level of discrepancy between the results given by the two methods.

7.2.2. Comparison of the microscopic behaviors of samples C and D

Fig. 23. Comparisons of the normal contact force anisotropy parameter af using the two methods (d = 320 μm): (a) Samples A and B, and (b) Samples C and D.

The micromechanical descriptors of samples C and D exhibit very different responses from those obtained in samples A and B for both methods.For sample C,all four descriptors computed using the two methods are much closer than those shown by samples A and B.For sample D,all of the four descriptors computed using the two methods merge. The reason for the consistency of the micromechanical responses given by the two methods in samples C and D is that they are both dense samples. Hence, when using the individual pore-based fluid-coupled DEM, negative pore pressures are generated under shearing (see Fig. 9b), which causes the microstructure to exhibit a tendency to dilate. The negative pore pressure generated after the initial shearing stage in samples C and D results in contact creation in both directions,which is shown by continuously increasing γ and a much slower increasing rate of anand af.In comparison,when dense samples C and D undergo shearing using the constant volume method,there is also a tendency of the microstructures to dilate such as what usually occurs when a dense sample is sheared under a drained condition.Since the total volume is preserved in the constant volume method,the dilation of the microstructure is restricted,thus more contacts are created in both directions, which is shown by a continuously increasing γ and a much lower rate of increase of anand af. Therefore, the negative pore pressure generated using the individual pore-based fluid-coupled DEM and the preservation of the volume given by the constant volume method both induced the same results of contact creation in both directions.

Fig.24. Contour of pore pressure distribution of sample A(a)at peak strength(d=80 μm),(b)at peak strength(d=320 μm),(c)at large strain(d=80 μm),and(d)at large strain(d = 320 μm).

7.3. Comparison of computational efficiencies of the two methods

For both methods, the execution time required for each of the four samples to achieve about 10% shear strain is different. Both methods are capable of generating detailed micromechanical information on each particle, such as variations of inter-granular contact orientations, normal contact force values and orientations,and tangential contact force values and orientations.Besides,the individual pore-based fluid-coupled DEM model can also generate micro-pore information such as the shape and micro-pore pressure values; therefore, the sample size and particle packing intensity affect the length of computation time. For instance,sample D, whose particle packing intensity is the highest among the four samples, contains the largest number of inter-granular contacts and interconnected pores, and hence, requires the longest computation time. Table 4 shows a comparison of the execution times for the four samples given by the two methods. A computer with an Intel i7-4770 K processor running at a frequency of 3.6 GHz was used for the simulations. It can be seen from the table that the execution time when using the individual pore-based fluid-coupled DEM is about 15 times that of using the constant volume method. Nevertheless, it is worthy of note that the very long execution time in the current version of the individual porebased fluid-coupled DEM code is due to the recreation of a polygon list at every time step without taking into account the current polygon list.Mapping is then provided between the current list and the list at the previous step. However, the execution time will be dramatically reduced if only an update of the list is implemented.There is being attempted in on-going work by the authors.

Fig. 25. Micro-pore pressure equalization process in sample A with d = 80 μm.

8. Discussion

8.1. Effect of permeability on undrained responses

From the above comparisons of macroscopic and microscopic responses,it is recognized that very different responses resulted for the two loose samples A and B for the two methods, while the responses of the two dense samples C and D derived by the two methods are much closer to each other. For samples A and B, the very different responses given by the two methods are explained by the different excess pore pressure dissipation rates generated in the two methods.The constant volume method simulates a scenario in which the excess pore pressures of the individual pores equalize instantly,while the excess pore pressure dissipation rate generated in the individual pore-based fluid-coupled DEM depends on the permeability of the granular assembly. On the other hand, for samples C and D,the much closer behavior achieved using the two methods is due to the resulting contact creation in both directions,although for the individual pore-based fluid-coupled DEM, it is because of the negative pore pressure generation, while with regard to the constant volume method,it is due to the preservation of the volume during shearing.

In an attempt to support the above interpretation,particularly to present the impact of permeability (the conduit diameter) on undrained response of the two loose samples A and B, a series of simulations using the individual pore-based fluid-coupled DEM but with a much higher value of the coefficient of permeability k were conducted and the results were compared with those from the constant volume method presented above. The much higher coefficient of permeability k is achieved by selecting a much higher value of conduit diameter d based on the relationship between k and d presented in Fig. 6. As an example, a conduit diameter of 320 μm is selected to represent the permeability k for the clean coarse sand. Macroscopic and microscopic responses for the four samples obtained using the individual pore-based fluid-coupled DEM with a conduit diameter equal to 320 μm and their comparison with those given by the constant volume method are presented in Figs.16-23, respectively. It can be seen in these figures that when the conduit diameter increases to 320 μm, both the macroscopic and microscopic responses given by the two methods for the four samples become much close to each other especially for samples A and B, although small discrepancies are still exhibited.This is because when the conduit diameter increases to 320 μm,the process of excess pore pressure equalization becomes much faster,which is close to the instant equalization given by the constant volume method, hence, behaviors obtained from using the two methods are much closer to each other.

8.2. Contour of local pore pressure distribution

The above interpretation of the impact of permeability on the process of pore pressure equalization, and therefore on the undrained behavior,can also be supported by considering the contour of local pore pressure distribution. As an example, the contour of the local pore pressure distribution of sample A given by the individual pore-based fluid-coupled DEM with conduit diameters of both 80 μm and 320 μm are presented in Fig. 24. Fig. 24a and b demonstrate the contours of local pore pressure distribution of sample A at peak strength given by the two conduit diameters.It is shown in this figure that both of them display higher micro-pore pressure values in the major principal stress direction near the top and bottom boundaries, while lower micro-pore pressure values in the minor principal stress direction occur near the left and right boundaries. Nevertheless, the ranges of micro-pore pressure values presented in Fig.24a and b are very different,a much wider range of micro-pore pressure distribution with 65 kPa and-25 kPa as maximum and minimum micro-pore pressure values, respectively,are displayed in Fig.24a compared to a discrepancy between the maximum and minimum micro-pore pressures of only 1 kPa between each other as shown in Fig.24b. The very different spans of micro-pore pressure values at the peak strength shown in Fig. 24a and b support the characterization of impact of permeability on the rates of excess pore pressure equalization presented above. When the conduit diameter is lower, micro-pore pressure variances are higher among pores in different locations in the sample. A random selection of five pores at different locations in Fig.24a and a test of their pore pressure equalization process were performed by fixing the boundary position and letting the fluids flow freely within the assembly.Fig.25 shows the results from this test, which indicate that the changing rate of all five micro-pore pressures was initially rapid, and then gradually reduced and approached zero at about 50,000 cycles where the five micro-pore pressure values merged together.

At the steady state (see Fig. 24c and d), the same pattern that higher micro-pore pressure values are shown in the major principal stress direction while lower micro-pore pressure values occur in the minor principal stress direction is exhibited as that at the peak strength state. Furthermore, micro-pore pressure values change evenly from the boundary to interior in both principal stress directions given by both conduit diameters. The only difference between them is that a much smaller range of pore pressure values occurs within the whole sample when the conduit diameter is increased to 320 μm.Nevertheless,a conduit diameter of 320 μm is still not large enough for the pressures of pores in any location of the assembly to be equalized instantaneously, which explains the remaining small difference in responses of the sample given by the two methods when the conduit diameter is increased to 320 μm.Therefore, the constant volume method may be capable of simulating the undrained behavior of a loose granular material with very coarse grains; however, it is not applicable when fine particle components are present due to its instant equalization of pore pressures.

9. Conclusions

A comparison of the undrained behaviors of granular media simulated using the individual pore-based fluid-coupled DEM and the constant volume method is presented in this paper. Four discshaped granular assemblies with fine particle components were prepared with a wide range of initial void ratios for the simulations.Both macroscopic and microscopic responses of the four samples obtained using the two methods were studied.It was found that the two loose samples exhibited higher values of peak strength and lower values of residual strength when simulated using the individual pore-based fluid-coupled DEM than that when the constant volume method was used.For the two dense samples adopted,both the macroscopic and microscopic responses given by the two methods were very close to each other at the early and intermediate shearing stages. This tendency increased as the void ratio of the sample decreased.

The reasons that different behaviors are predicted by the two methods for the two loose samples while similar behaviors are predicted by the two methods for the two dense samples are clarified by a study of the micromechanical responses of the four granular assemblies. It is demonstrated that the differences in behavior given by the two methods for the loose samples are due to the difference in their excess pore pressure dissipation rates. The constant volume method describes a scenario of instant equalization of the excess pore pressure, while the process of excess pore pressure equalization in the individual pore-based fluid-coupled DEM depends on the coefficient of permeability of the granular assembly. Therefore,the constant volume method may be capable of simulating the undrained behaviors of loose granular material with very coarse grains, which has a high permeability and thus rapid pore pressure equalization. However, it is not applicable when fine particles are present and pore pressures equalize more slowly.For the dense samples,the similar behavior at the early and intermediate shearing stages given by the two methods occurs because both methods induce the same response of contact creation in both directions. In the individual pore-based fluid-coupled DEM, this is due to the generation of a negative pore pressure;while in the constant volume method,it is due to the preservation of a constant volume. Hence, using both methods, a very close agreement in the undrained responses is induced.The deficiency of the constant volume method for a very dense sample as indicated by Hanley et al. (2013) was not noticed because the two dense samples adopted in this study are only medium dense indicated by their void ratio and average coordination number values.

Declaration of competing interest

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

Acknowledgments

The authors would like to express their acknowledgments to Dr.Timothy Topper for his insightful comments regarding our English during the preparation of this paper.