Experiments and 3D DEM of Triaxial Compression Tests under Special Consideration of Particle Stiffness

Discrete element modelling is commonly used for particle-scale modelling of granular or particulate materials. Developing a DEM model requires the determination of a number of micro-structural parameters, including the particle contact stiffness and the particle-particle friction. These parameters cannot easily be measured in the laboratory or directly related to measurable, physical material parameters. Therefore, a calibration process is typically used to determine the values for use in simulations of physical systems. This paper focuses on how to define the particle stiffness for the discrete element modelling in order to perform realistic simulations of granular materials in the case of linear contact model. For that, laboratory tests and numerical discrete element modelling of triaxial compression tests have been carried out on two different non-cohesive soils i.e. poorly graded fine sand and gap graded coarse sand. The results of experimental tests are used to calibrate the numerical model. It is found that the numerical results are qualitatively and quantitatively in good agreement with the laboratory tests results. Moreover, the results show that the stress dependent of soil behaviour can be reproduced well by assigning the particle stiffness as a function of the particle size particularly for gap graded soil.

1. Introduction

Numerical modelling methods, such as finite element or finite differential methods are generally used to investigate geotechnical engineering problems, in order to assess for instance the response of soil subjected to imposed loads and/or changes in boundary conditions. For the continuum numerical models such as the finite element and the finite differential methods, the variables such as displacement and stress are assumed to vary continuously, despite the evident discontinuous nature of the soil. The continuum numerical models are primarily based on the mathematical modelling of the observed phenomena at macroscopic scale, but do not reproduce the local discontinuous nature of the soil material. However, these local discontinuities play a major role in the behaviour of granular materials, Belheine et al. (2009) [1] . These local discontinuities induce special features such as anisotropy, micro-fractures or local instabilities, which are not easy to capture or understand by means of continuum modelling, Hu et al. (2010), Arthur and Menzies (1972) [2] [3] .

Therefore, the discrete element method (DEM) is an alternative approach, which considers the discrete nature of granular materials, and provides new insight as far as the constitutive model for soil material concerns. The Particle Flow Code (PFC) is a DEM Code, which is based on the work of Cundall (1971), Cundall and Hart (1979), Cundall and Hart (1992), Cundall (2001) [4] [5] [6] [7] [8] . In PFC the mechanical behaviour of soil material was described in terms of the movement of each particle and the inter-particle forces acting on each contact particle. The discrete element method has been proven to be a very useful tool in obtaining complete qualitative information level with simple assumptions and few parameters at the microscopic scale, Belheine et al. (2009) and Bagi (2005) [1] [9] . DEM is used over a wide range of applications which include construction and building materials, Coetzee and Els (2009), Harkness et al. (2011), Stahl and Konietzky (2011), Indraratna et al. (2009), Lu and McDowell (2010), Suhr and Six (2017) [10] - [15] , geotechnical applications, Belheine et al. (2009), Chen et al. (2012), Van Lysebetten et al. (2014) [1] [16] [17] , mining, Grima et al. (2010, 2011) [18] [19] , post-harvest, Boac et al. (2014) [20] , soil-structure interaction, Coetzee (2009, 2010) [21] [22] [23] , mixing and milling, Alian et al. (2015), Cleary (2015) [24] [25] [26] and amongst others. The discrete element method (DEM) has become the method of choice for researching and engineering to validate and optimise the design related to granular material or to assess the phenomena on grain scale.

Before DEM modelling can be carried out with confidence, an accurate set of input particle parameters is needed Coetzee (2016) [27] . The material macro properties (real material parameters e.g. unit weight, cohesion, angle of internal friction) and the DEM micro parameters (e.g. particle parameters, particle normal stiffness kn, particle shear stiffness ks, particle coefficient of friction μ) have to be distinguished. In the literature, micro parameters are often not measured and the values are assumed without specific justification. How the parameter values were obtained is often not mentioned and whether they were measured is not clear. Calibration procedures are required to obtain reliable micro parameters for the granular materials. Two approaches are generally applied for the determination of the DEM parameters, Marigo and Sitt (2015) [28] . The first approach is the most used and consists of a calibration procedure where either in-situ or laboratory experiments are performed to measure a specific macro property. The experiment is then numerically reproduced by following the laboratory or field setup and procedures as closely as possible. The micro parameter values are then changed iteratively until the predicted macro response matches the measured numerical result. A possible problem with this approach is that the macro response of the numerical experiment can be influenced by more than one parameter. This means that there would be no unique solution since more than one combination of the parameter values may result in the same macro behaviour, Marigo and Sitt (2015), Coetzee (2016, 2017) [27] [28] [29] .

The following authors used the first calibration approach. Belheine et al. (2009) [1] carried out a 3D discrete element modelling of granular material. The DEM parameters were calibrated using the triaxial test on Labenne sand. They concluded that the peak stress of the numerical sand sample does not only depend on local friction parameters but also on the rolling resistance; and the deformation response depends strongly on local friction. Plassiard et al. (2009) [30] proposed a calibration approach for spherical particles for discrete element modeling considering the transfer of the moment between the spherical particles. They concluded that the incremental response is well described by elastoplasticity law with a single mechanism, and a non-associative flow rule. Harkness et al. (2016) [11] investigated the effects of confining pressure ranging from 15 to 200 kPa on the behaviour of scaled railway ballast in triaxial tests using discrete element modelling (DEM). They introduced a contact law modelling damage at the contacts between particles. They concluded that simulation results show excellent agreement with the laboratory data for monotonic triaxial tests. The DEM simulation of high pressure triaxial tests on crushable sand using the octahedral shear stress as fracture criterion showed also that particle crushing is essential to reproduce the realistic behaviour of sand in particular the volumetric contraction in high-pressure shear tests, de Bono and McDowell (2014) [31] . However, McDowell and Li (2016) [32] made use of tetrahedral clumps with breakable asperities to model scaled railway ballast using the discrete element method, and concluded that the strength affects the macroscopic shear strength at both high and low confining pressures, while the effects of the number of asperities diminishes with increasing confining pressure due to asperity breakage. Zhao et al. (2011) [33] used the discrete element method to simulate granular assembly behaviour with different initial conditions under three different loading conditions e.g. plane strain, conventional triaxial compression, and direct shear. They concluded that there is no a unique critical state line for specimens with different initial void ratios under different loading conditions. They showed that the anisotropy is shown to have significant effects on the soil behaviour and is related to the non-uniqueness of the critical state line. Stahl and Konietzky (2011) [12] developed calibration procedure under consideration of grain shape, grain size and relative density in order to reproduce the stress strain behaviour of stiff granular material under different loading conditions. Coetzee et al. (2009) [21] and Lommen et al. (2014) [34] showed also that the shear box results were influenced by the contact stiffness and friction coefficient. Li et al. [35] made use of triaxial compression tests and a calibration procedure based on a response surface method to determine the normal stiffness, the shear stiffness and the friction coefficient of rock fill material by means of Particle Flow Code (PFC). The results of their investigation showed that the normal stiffness, the tangent stiffness and the friction coefficient of rockfill materials will slightly increase with increase of confining in triaxial compression tests. Thus, these results confirm the stress dependent behaviour of the rockfill material. Based on DEM results, the particle friction coefficient and stiffness were determined from energy principles and direct shear tests. Derakhshani et al. (2015) [36] used spherical particles to model quartz sand. They determined both the particle sliding friction and rolling friction by modelling the flow of sand through a sandglass. They showed that two independent parameters were needed to be measured in order to determine a unique set of values for the two unknown particle sliding friction and rolling friction. For this purpose, they used the angle of repose that formed as the sand flowed through the sandglass as well as the discharging time. The results indicated that the coefficient of rolling and sliding of sand particles are 0.30 and 0.52, respectively.

The second approach for the determination of the input parameter values of the particle is to directly measure the values on the particle level. Some of these parameters are easy to measure while others are very challenging, depending on the particle size. Several approaches were made in literature, but they were all applied to particles in the millimetre and above size range, Marigo and Sitt (2015), Coetzee (2016) [27] [28] . Even if the micro parameter values can be directly measured, it does not necessarily mean that the DEM model would show the same level of accuracy on a macro level. This approach would only be accurate if the shape of the particles, the size of the particles and the particles size distribution are modelled accurately and if the contact model is an accurate representation of the contact behaviour, Barrios et al. (2013) [37] . In practice, it is difficult to accurately model the particle size and shape when large industrial scale systems are modelled. The particle size often has to be increased and the particle shape cannot be accurately modelled due to computational limitations. The advantage of this direct measurement approach is that the resulting micro properties are not dependent on the contact model or the specific DEM code used, Coetzee (2016), González-Montellano (2011) [27] [38] . Some researchers have tried to experimentally measure the micro parameters. Vu-Quoc et al. (2000) [39] measured the coefficient of restitution in soybeans using drop tests and Gonzalez-Montellano et al. (2012) [40] measured the micro parameters such as particle density, modulus of elasticity, particle-wall coefficient of restitution, particle-particle coefficient of restitution, and the particle-wall coefficient of friction of glass beads, maize grains and olives, and validated the procedure by modelling silo discharge, Gonzalez-Montellano et al. (2012) [41] . Paulick et al. (2015) [42] developed a technique to examine the contact behaviour between two particles and to measure the particle contact stiffness while eliminating the possible container deformation. This was however only applied to glass beads with a particle diameter from 0.8 to 3.0 mm and a maximum compression force of 80 N.

The literature review above reveals that the behaviour of granular material can be predicted numerically with the discrete element method. This method potentially uses material or particle properties to describe the behaviour of a particle and its interactions with other particles or walls. One characteristic particle parameter is the particle contact stiffness, whose determination was not still performed by means of a definitive comprehensive procedure. Therefore, this paper focuses on the determination of the particle stiffness for granular soil. For that, triaxial compression tests have been carried out on two granular soils in the laboratory. Then the results of these triaxial tests in terms of the mobilized angle of internal friction have been calibrated using the first calibration approach described above. Hereby, the commercial DEM code, Particle Flow Code in Three (3) Dimensions PFC3D has been used. The comparison between the experimental and numerical results shows that the macroscopic behaviour of real granular material can be quantitatively and qualitatively reproduced by means of a discrete element method, i.e. Particle Flow Code (PFC3D). Furthermore, analysis accounting for particle stiffness depending on particle size is carried out. The results of this analysis show best comparison between the experimental and numerical results.

2. Experiments

2.1. Soil Materials Tested

Two (2) granular soils have been used for the present study. The curve of the grain size distribution obtained from sieve analysis in accordance with German code DIN 18123 [43] , and the geotechnical properties for the soils tested are shown in Figure 1 and Table 1, respectively. The minimum porosity nmin and the maximum porosity nmax have been determined in accordance with German code DIN 18126 [44] . The soil A1 is poorly graded fine to medium sand, whereas the soil E2 is gap graded coarse sand.

2.2. Procedure and Results of Triaxial Tests in Laboratory

The experimental consolidated drained (CD) triaxial compression tests on soils A1 and E2 have been carried out in accordance with the German code DIN 18137 [45] . The soil specimens are reconstituted in a cylinder of 5 cm diameter

Table 1. Properties of soils used for triaxial compression tests.

Figure 1. Grain size distribution curves of investigated soils A1 and E2.

and 12 cm height. A relative density of D = 0.75 is reached. The applied relative density is defined as $D=\left({n}_{\mathrm{max}}-n\right)/\left({n}_{\mathrm{max}}-{n}_{\mathrm{min}}\right)$ . Here, nmax is the maximum porosity, nmin is the minimum porosity and n is the actual porosity of the soil. Water content of 10% and 7% is achieved for the soil material A1 and E2, respectively.

Thereafter, the specimens are frozen at a temperature of −20˚C and therefore have sufficient strength to be installed in the triaxial cells. The applied water content is the optimum one for which the grains of the soil get just wet so that at “frozen state” a solid bond grains-ice occurs. Higher water content would lead to an unfavourable volumetric increase during the freezing process of the sample.

After unfreezing, the saturation of each soil sample is performed with a back pressure of 6 bar until a degree of saturation Sr > 90% is reached. The saturation of the specimen limits the air trapped in the pore structure of the soil sample. As result, the air trapped moved into the compression liquid i.e. water and does not falsify the test result through compressed air in the soil sample. After the saturation phase follows the consolidation phase. Hereby, the initial stress state before shearing the soil sample is achieved. The applied confining pressures are 50, 100 and 200 kPa. Each sample is subjected to the above mentioned constant confining pressure and additionally to a constant shear rate. The resulting axial stress is measured by means of a load cell. Since the soil specimen consists of non-cohesive soils, the shearing of the soil sample is carried out with a constant shear rate of 0.08 mm/min in accordance with the German code DIN 18137 [45] in order to avoid undrained shear strain due to possible excess pore water pressure in the soil specimen.

The results of the experimental investigations of triaxial compression test on soil A1 and E2 are presented in Figure 2, which shows non-linear stress-strain behaviour of soils investigated, Ahlinhan et al. (2016) [46] . The higher the confining stress (σ2) is, the higher is the deviatoric stress (σ1 − σ2) at the peak state and post peak state. After the peak state, the critical state is reached and the deviatoric stress remains approximately constant (see Figure 2(b)). However, the samples of the soil A1 show some hardening effect in the plastic zone (see Figure 2(a)). Soil E2 shows a maximum dilatancy for an axial strain of 1%, which corresponds to the peak stress. The dilatancy is less pronounced for the soil A1. This could be explained by the gap-gradation of the soil E2, which shows more instability and therefore more dilatancy. The mechanical instability of the soil material starts at the point of the maximum dilatancy, which corresponds to the deformation at the peak stress, Belheine et al. (2009) [1] .

The normal vertical stress σ1 and the confining stress σ2 are related to the mobilized angle of internal friction φmob as follows:

$\mathrm{sin}{\phi }_{mob}=\left({\sigma }_{1}-{\sigma }_{2}\right)/\left({\sigma }_{1}+{\sigma }_{2}\right)$ . (1)

Figure 2. Results of experimental triaxial test for a relative density D = 0.75. (a) Deviatoric stress vs axial strain for soil Material A1; (b) Deviatoric stress vs axial strain for soil Material E2; (c) Volumetric strain vs axial strain for soil Material A1; (d) Volumetric strain vs axial strain for soil Material E2.

Equation (1) is valid for non-cohesive soils such as soil material A1 and E2.

Figure 3 shows the mobilized angle of internal friction for the soils A1 and E2. The mobilized angle of internal friction increases almost linear with the strain up to the peak. This reflects the linear elastic behaviour of the soil materials. Beyond the peak, the curve of the mobilized internal friction tends towards an asymptote, and the mobilized angle of internal friction becomes approximately constant.

3. Numerical Modelling

3.1. Discrete Element Method and PFC3D

The numerical modelling is performed by means of a discrete element method (DEM) with Particle Flow Code (PFC3D) Itasca (2005) [47] using the most widely applied linear elastic model. Hereby, the soil grains are idealized as spheres, Belheine et al. (2009), Derakhshani et al. (2015), Tom Woerden et al. (2004) [1] [36] [48] . The material behaviour is simulated by a linear contact stiffness constitutive law (normal stiffness kn, shear stiffness ks) and the sliding law in accordance with Coulomb law (friction coefficient μ) as shown in Figure 4. The maximum shear force ${F}_{\mathrm{max}}^{s}$ is related to the normal force Fn as follows ${F}_{\mathrm{max}}^{s}=\mu \cdot {F}^{n}$ .

For a run time, the contact stiffness is calculated at each contact, using two linear springs in series. These two springs are either particle-particle springs or particle-wall springs depending on the specific contact. Particles are allowed to overlap at the contact points. When particles have contact, the contact force is calculated as function of relative displacement and stiffness governed by a constitutive contact model.

The force-displacement law is described in PFC for both particle-particle and ball-wall contacts. For particle-particle contact, the relevant equations are presented for the case of two spherical particles, labelled A and B in Figure 4. For particle-wall contact, the relevant equations are also presented for the case of a spherical particle and a wall, labelled p and w, respectively, in Figure 5. In both cases, ${U}^{n}$ denotes overlap.

(a) (b)

Figure 3. Mobilized friction angle, (a) for soil material A1; (b) for soil material E2.

Figure 4. Shear and normal contact model in particle flow code, Itasca (2008) [47] .

Figure 5. Notation used to describe (a) particle-particle contact (b) particle-wall contact, Itasca (2008) [47] .

For particle-particle contact, the unit normal, ${n}_{i}$ , that defines the contact plane is given by

${n}_{i}=\left({x}_{i}^{\left[B\right]}-{x}_{i}^{\left[A\right]}\right)/d$ . (2)

Here, ${x}_{i}^{\left[A\right]}$ and ${x}_{i}^{\left[B\right]}$ are the position vectors of the centres of balls A and B, and d is the distance between the particles centres:

$d=|{x}_{i}^{\left[B\right]}-{x}_{i}^{\left[A\right]}|=\sqrt{\left({x}_{i}^{\left[B\right]}-{x}_{i}^{\left[A\right]}\right)\left({x}_{i}^{\left[B\right]}-{x}_{i}^{\left[A\right]}\right)}$ . (3)

The overlap ${U}^{n}$ , defined to be the relative contact displacement in the normal direction, is given by

${U}^{n}={r}^{\left[A\right]}+{r}^{\left[B\right]}-d$ (Particle − Particle) (4)

${U}^{n}={r}^{\left[b\right]}-d$ (Particle − Wall) (5)

The contact force vector ${F}_{i}$ (which represents the action of particle A on particle B for particle-particle contact, and represents the action of the particle on the wall for particle-wall contact) can be decomposed into normal and shear components with respect to the contact plane as:

${F}_{i}={F}_{i}^{n}+{F}_{i}^{s}$ (6)

Here, ${F}_{i}^{n}$ and ${F}_{i}^{s}$ denote the normal and shear component vectors, respectively.

The normal contact force vector for the applied linear elastic contact is calculated by

${F}^{n}={K}^{n}{U}^{n}$ (7)

where ${K}^{n}$ is the resultant normal stiffness at the contact. The value of ${K}^{n}$ is determined by the linear contact stiffness model by:

${K}^{n}=\left({k}_{n}^{A}\cdot {k}_{n}^{B}\right)/\left({k}_{n}^{A}+{k}_{n}^{B}\right)$ (8)

${K}^{n}=\left({k}_{s}^{A}\cdot {k}_{s}^{B}\right)/\left({k}_{s}^{A}+{k}_{s}^{B}\right)$ (9)

Here, ${k}_{n}^{A}$ and ${k}_{s}^{A}$ are the normal and tangential stiffness for the particle A, respectively. ${k}_{n}^{B}$ and ${k}_{s}^{B}$ are the normal and tangential stiffness for the particle B, respectively.

The frictional sliding starts at the contact point if the contact forces ${F}^{s}$ and ${F}^{n}$ satisfy the frictional Mohr-Coulomb equation:

$‖{F}^{s}‖-\mu \cdot ‖{F}^{n}‖\le 0$ (10)

with μ the particle-particle friction coefficient.

The values of kn, ks and μ are determined via the calibration on triaxial tests in laboratory.

For the consideration of the dissipation of energy of particles during motion and interaction, PFC offers two types of damping: global damping applied to particles and viscous damping applied to contacts. Global damping applies a damping force with magnitude proportional to unbalanced force to each particle. This is usually used in modelling quasi-static problems. Viscous damping adds a normal and shear dashpot to each contact. This damping acts proportional to relative velocities of two particles. Hence, damping forces are added to normal and shear contact force.

When the problem has reached a quasi-static state, damping is not any more considered in the numerical model. The default value for the global damping coefficient of 0.7 is applied. For dynamic simulations, local damping is deactivated and only viscous damping is used, Itasca (2008) [47] .

To generate numerical particle packing, numerous methods exist, Bagi (2005) [9] . The more common methods used are:

1) “Defined Position Method”, hereby particles are generated at given position. For that, there are numerous appropriate algorithms to create the packing depending on desired density. The disadvantage of this method is the irregular lower density of the packing.

2) “Fill and Expand Method” is a typical approach to generate the required number of particles into a domain of interest. Then the particles diameters are gradually increased until the desired relative density is reached. This technique may lead to high stress for particles, however this technique recommended for triaxial test due to the high stress level for such test.

3) “Gravitation Deposition Method”, hereby the particles with defined size are generated without overlapping in a random distribution inside a desired volume, which is larger than the target volume. Then the particles fall down under the effect of the acceleration of gravity into the domain or into a container, and find their equilibrium position. This method is time consuming.

The “Fill and expand” method is applied to generate the numerical particles assembly for the present numerical analysis, De Bono and McDowell (2014), Itasca (2008) [31] [47] .

For the numerical simulation of the triaxial test, the particles have been generated in accordance with the grain size distribution curve of soils presented in Figure 1. For example, the grain size distribution curve of soil material E2 is split in four parts as shown in Figure 6. The parts 0 - 1 and 1 - 2 are considered as fine fraction, whereas the parts 2 - 3 and 3 - 4 are coarse one. The minimum and the maximum grain size of each part are entered as input in PFC3D for the generation of the numerical sample with respect to the grain size distribution.

In reality, the soil material consists of different grain sizes, which are depicted by the curve of the grain size distribution. Since the grain size distribution curve is considered for the present analysis, the question of how to assign the particle stiffness in function of the particle size needs to be addressed in Section 3.2.

3.2. Approach of Particle Stiffness Depending on Particle Size

As shown in Section 1, the particle stiffness of fine granular material is different from the particle stiffness of coarse material. The assignment of the particle stiffness can be addressed from two points of view. Considering a real homogeneous soil material with similar grain sizes, it can be assumed that each grain of the same material could exhibit the same stiffness. Therefore, for a first assumption it can be considered that the particle stiffness is the same for both fine particles and coarse fraction, i.e. kn = constant and ks = constant for a DEM model. However, from micromechanical point of view, the interaction i.e. force-displacement between particles should rather be considered. It can be believed that the particle stiffness depends on the particle size, therefore the particle stiffness has to be assigned in function of the particle diameter. This is not generally considered in the existing DEM modelling using linear elastic contact

Figure 6. Splitting of the grain size distribution of soil E2 for numerical particle generation.

model. However, the particle stiffness is a function of particle diameter for the Hertz model with non-linear elastic contact law. Therefore, the approach of particle size dependent stiffness does not need to be considered, when the Hertz contact model is applied.

In Hakuno and Tarumi (1988) [49] [50] , the resulting normal stiffness Kn and tangential stiffness Ks between two particles of different size were determined. Hereby, a contact theory for an elastic cylinder was applied. When two particles of radii r1, r2 with the same Young’s modulus E and Poison’s ratio ν, are subjected to compression force q per unit depth as shown in Figure 7, the displacement δ, and the overlap Un between the two particles can expressed as follows:

$\delta =\left[\frac{2q\left(1-{\nu }^{2}\right)}{\text{π}E}\right]\left[\frac{3}{2}+\mathrm{ln}\left(\frac{4{r}_{1}}{{U}^{n}}\right)+\mathrm{ln}\left(\frac{4{r}_{1}}{{U}^{n}}\right)\right]$ (11)

${\left({U}^{n}\right)}^{2}=\left[\frac{8{r}_{1}{r}_{2}}{\text{π}\left({r}_{1}+{r}_{2}\right)}\right]\left[\frac{q\left(1-{\nu }^{2}\right)}{\Pi E}\right]$ (12)

Then, Kn can be calculated for a particle of radius r as follows:

${K}^{n}=\frac{q}{\delta }=\frac{\text{π}E}{2\left(1-{\nu }^{2}\right)\left(1.5+2\mathrm{ln}\left(\frac{4r}{{U}^{n}}\right)\right)}$ (13)

Considering two small particles with radius r1 with the resulting stiffness ${K}_{1}^{n}$ and two larger particles of radius r2 with the resulting stiffness ${K}_{2}^{n}$ , such as r2 = λ·r1, λ > 1 (Figure 7), the ratio ${K}_{2}^{n}/{K}_{1}^{n}$ can be expressed as follows:

$\frac{{K}_{2}^{n}}{{K}_{1}^{n}}=\frac{1.5+2\mathrm{ln}\left(\frac{4{r}_{1}}{\sqrt{\frac{4{r}_{2}\left(1-{\nu }^{2}\right)q}{\text{π}E}}}\right)}{1.5+2\mathrm{ln}\left(\frac{4{r}_{2}}{\sqrt{\frac{4{r}_{1}\left(1-{\nu }^{2}\right)q}{\text{π}E}}}\right)}$ (14)

Figure 7. Two particles models with different radii.

Or

$\frac{{K}_{2}^{n}}{{K}_{1}^{n}}=\frac{1.5+2\mathrm{ln}\xi }{1.5+2\mathrm{ln}\left({\lambda }^{1/2}\xi \right)}$ (15)

with

$\xi =\frac{4{r}_{1}}{\sqrt{\frac{4{r}_{2}\left(1-{\nu }^{2}\right)q}{\text{π}E}}}$ (16)

Table 2 shows the parameters for sand grain as given in Hakuno and Tarumi (1988) [49] [50] .

The evaluation results of the particle stiffness based on the equations 14 to 16 for particle diameters between 0.06 mm and 3 mm are presented in Figure 8 and Figure 9. The stiffness ratio can be determined from Figure 8 and Figure 9. For instance, the first grain part of the fine fraction for soil E2 is characterized by particle size between d1 = 0.06 mm and d2 = 0.1 mm. Then, the average radius $d=\left({d}_{1}+{d}_{2}\right)/2=0.08$ mm and the average radius r = 0.00004 m. Based on the Figure 8, ξ = 88.7 for the average radius of 0.00004 m. Since the ratio of radius can be calculated for each part of grain size distribution curve, then the stiffness ratio ${K}_{2}^{n}/{K}_{1}^{n}$ can be derived from Figure 9. Table 3 presents for example the dependence of the particle stiffness on the particle size for the soil material E2.

3.3. Numerical Simulation of Triaxial Test

The determination of relevant particle parameters i.e. the normal contact stiffness kn, the shear contact stiffness ks, and the friction contact coefficient µ for the used non-cohesive soils has been carried out by means of the discrete element model with PFC3D. The model has been calibrated with the experimental behaviour of real material shown in Figure 2 and Figure 3.

For the numerical simulation of the triaxial test the particles are generated using the fill and expand method inside a cylinder with diameter b at least ten times the maximum particle diameter and a height h about two times the cylinder diameter. For a realistic simulation of reliable triaxial tests by means of the discrete element method, the container size shall be large enough to limit the boundary effect such as silo effect. Pohl (2005) [51] stated that 15 particles at least shall be located on the container length in order to reproduce the soil material parameters (angle of internal friction, void ratio). In practice, the numerical sample size shall be larger than 10 times the maximum grain size, Biarez and

Table 2. Parameters for sand grain.

Figure 8. ξ-value in function of particle radius.

Figure 9. Stiffness ratio ( ${K}_{2}^{n}/{K}_{1}^{n}$ ) as a function of ratio of particle radius (r2/r1).

Table 3. Dependence of particle stiffness on particle size for soil E2 and A1.

Hicher (1994) [52] . Considering these requirements, the numerical sample of sand A1 contained 24922 particles (Figure 10), which can be run with the

Figure 10. Numerical triaxial sample for sand A1.

hardware with 4 × 3.00 GHz and 2.00 GB RAM. This quantity of particles is larger than the applied particle number for many of the similar simulation for cemented particles De Bono et al. (2014), Utili and Nova (2008), Camusso and Barla (2009) [53] [54] [55] or crushable soils, De Bono and McDowell (2014), Bolton et al. (2008), Harireche and McDowell (2002) [31] [56] [57] .

The top and bottom of cylinder are closed with two walls, which simulate the loading plates of the triaxial cell. The stiffness of the cylinder is set to 10% of the stiffness of the particle in order to simulate the latex membrane of triaxial cell, Tom Woerden et al. (2004) [48] .

First the isotropic stress state is achieved by means of the servo-control-method. The top and bottom plates are moved vertically to each other, until the required isotropic stress for each numerical sample in each triaxial cell was reached (50, 100 or 200 kPa). This corresponds to the consolidation phase of the triaxial test in the laboratory. Thereafter the numerical sample is also loaded in a strain-controlled mode by specifying the opposite velocities for the frictionless top and bottom walls. The calculation is run until the critical state of the numerical sample is reached.

A calibration is performed using the first calibration approach, by fitting the numerical particle behaviour to the real soil behaviour in the laboratory triaxial tests presented Figure 2 and Figure 3. During the calibration process the model parameters such as the normal stiffness kn and the shear stiffness ks are varied in order to match the stress-strain curve of real material for strain lower than 1%, i.e. the elasticity modulus, while the particle-particle friction coefficient is kept constant. Thereafter, the particle-particle friction coefficient µ was varied to adjust the peak stress.

4. Results and Discussions

The numerical and experimental mobilized friction angle and the volumetric strain versus the axial strain for soil A1 and E2 for different contact stiffness are shown in Figure 11 and Figure 12, respectively. Table 4 summarizes the calibration parameter for the soil A1 and E2.

For the axial strain lower than 2.5%, the higher the contact stiffness is, the higher is the modulus of elasticity. The peak state is reached for an axial strain of 5%, thereafter the numerical soil shows a plastic behaviour. Hence, the non-linear stress-strain behaviour of the soil can be indeed reproduced by means of the discrete element modelling. A good comparison between experimental and numerical mobilized friction angle can be observed for kn = ks = 105 N/m for soil A1 (Figure 11). Hereby, there is only a slight difference between the two numerical approaches with and without the particle diameter dependent of stiffness. This can be explained by the uniformity of the grading curve of the soil material A1 and the relative small range between maximum and minimum particle diameter. However, for gap graded soil E2 the analysis results show clear difference between the approach of particle size independent stiffness and particle size dependent stiffness (see Figure 12). The approach of grain size dependent stiffness reproduces better the stress dependent of soil behaviour. Therefore, it can be assumed that the stress dependent of soil behaviour can be better taken into account for the approach of grain size dependent grain stiffness.

However, the volumetric strain shows that the discrete element assembly exhibits more contraction and more dilation volume than the real soil. This result is in agreement with the results reported by Tom Woerden et al. (2004) [48] , McDowell and Li (2016) [32] . Similar results are also reported by de Bono and McDowell (2014) [31] based on the DEM simulations, and by Bolton et al. (2008) [55] based on the observations in realistic laboratory triaxial conditions.

The results of the discrete element modelling of the large shear test in Coetzee (2016) [27] show that the numerical contraction and dilation are not very comparable with the experimental contraction and dilation.

In order to best fit the experimental and the numerical volumetric strain, Belheine et al. (2008) [1] introduced in their DEM-model additional parameters such as rolling stiffness coefficient and non-dimensional plastic coefficient of the contacts. However, general validated procedures to measure these additional

Table 4. Calibration parameter of soil A1 and E2.

Figure 11. Experimental versus numerical results for soil material A1. (a) Mobilized friction angle with axial strain for σ2 =50 kPa; (b) Mobilized friction angle with axial strain for σ2 =100 kPa; (c) Mobilized friction angle with axial strain for σ2 =200 kPa; (d) Volumetric strain with axial strain for σ2 =50 kPa; (e) Volumetric strain with axial strain for σ2 =100 kPa; (f) Volumetric strain with axial strain for σ2 =200 kPa.

coefficients are not yet well established. As the soil grains have been modelled by means of spheres, a high coefficient of friction needs to be applied in order to cover the rotation resistance of the real soil grains (Table 4). Similar high coefficient of friction was applied for the discrete modelling of the soil grain by means of sphere without clump model, Tom Woerden et al. (2004) [48] .

The simulation and laboratory results show that the higher the confining pressure is, the smaller is the mobilized friction. These results are comparable with the results of Kozicki et al. (2012) [58] . In reality, soil particles are not spherical,

Figure 12. Experimental versus numerical results for soil material E1. (a) Mobilized friction angle with axial strain forσ2 =50 kPa; (b) Mobilized friction angle with axial strain for σ2 =100 kPa; (c) Mobilized friction angle with axial strain for σ2 =200 kPa; (d) Volumetric strain with axial strain for σ2 =50 kPa; (e) Volumetric strain with axial strain for σ2 =100 kPa; (f) Volumetric strain with axial strain for σ2 =200 kPa.

they are rather characterized by a polygonal edge form, which resists rolling. This may lead to crushing grains under triaxial loading condition.

Figure 13 and Figure 14 show the experimental and numerical Mohr’s-circles as well the failure envelope. The failure envelope for the peak stresses represents the straight line passing through the origin, as the investigated soil materials are non-cohesive.

The angle of internal friction is the angle between the straight line of the

(a) (b)

Figure 13. Mohr’s circle representation, stress failure envelope and angle of internal friction for soil A1. (a) Experimental results; (b) Numerical results.

(a) (b)

Figure 14. Mohr’s circle representation, stress failure envelope and angle of internal friction for soil E2. (a) Experimental results; (b) Numerical results.

failure envelope and the abscissa axis. The experimental and numerical angle of internal friction for soil A1 is equal to 31.74˚ and 31.22˚, respectively (Figure 13). For the soil E2 the experimental and numerical angle of internal friction is equal to 32.95˚ and 33.08˚, respectively (Figure 14). Hence, a good agreement between the experimental and numerical results can be stated. Similar angle of internal friction (31˚ for initial void ratio of 0.74) is obtained by means of DEM triaxial test using clusters of spheres, Kozicki et al. (2012) [58] .

5. Conclusions

For the design, analysis and optimization of geotechnical structures, e.g. foundations, dams, dykes, and filters, discrete element method (DEM) can be applied. Shear tests such as triaxial compression tests have been carried out to assess the shear behaviour of the soil, to determine the shear parameters and to calibrate the input particle parameters for the DEM.

Experimental and numerical triaxial compression tests have been carried out on two non-cohesive soils materials A1 (poorly graded soil) and E2 (gap graded soil).

The results of the experimental tests have been used to calibrate the discrete element model. It is found that the shear behaviour of the two soil materials investigated can be qualitatively and quantitatively reproduced by means of a Discrete Element Method (DEM) such as Particle Flow Code (PFC3D). However, the volumetric strain is overestimated and can be explained by the applied spherical particle. Furthermore, the particle size dependent particle stiffness has been analysed. The results show that the particle size dependent particle stiffness influences negligibly the shear behaviour of the poorly graded soil (A1). In contrast, the influence of particle size dependent particle stiffness on the shear behaviour of the gap graded soil (E2) is definitive. Hence, the stress dependent soil behaviour can be better reproduced, when the particle stiffness is assigned as a function of particle size as presented in this paper.

Moreover, the calibrated particles parameters can be used for modelling and analysis of practical issues e.g. suffusion piping, erosion, ballast for road constructions, foundations, stability of slope, etc.

Conflicts of Interest

The authors declare no conflicts of interest.

Cite this paper

Ahlinhan, M. , Houehanou, E. , Koube, M. , Doko, V. , Alaye, Q. , Sungura, N. and Adjovi, E. (2018) Experiments and 3D DEM of Triaxial Compression Tests under Special Consideration of Particle Stiffness. Geomaterials, 8, 39-62. doi: 10.4236/gm.2018.84004.