^{1}

^{*}

^{1}

^{1}

^{2}

^{3}

^{4}

^{1}

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.

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) [

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) [

Before DEM modelling can be carried out with confidence, an accurate set of input particle parameters is needed Coetzee (2016) [_{n}, particle shear stiffness k_{s}, 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) [

The following authors used the first calibration approach. Belheine et al. (2009) [

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) [

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 PFC^{3D} 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 (PFC^{3D}). 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.

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 [_{min} and the maximum porosity n_{max} have been determined in accordance with German code DIN 18126 [

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 [

Property | Soil material | |
---|---|---|

A1 | E2 | |

Unit weight of grains γ_{s} [t/m^{3}] | 2.65 | 2.65 |

Minimum porosity n_{min} | 0.40 | 0.27 |

Maximum porosity n_{max} | 0.52 | 0.40 |

Coefficient of uniformity C_{u} = d_{60}/d_{10} | 2.10 | 13.9 |

Coefficient of curvature C c = d 30 2 / ( d 10 × d 60 ) | 1.00 | 6.70 |

and 12 cm height. A relative density of D = 0.75 is reached. The applied relative density is defined as D = ( n max − n ) / ( n max − n min ) . Here, n_{max} is the maximum porosity, n_{min} 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 S_{r} > 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 [

The results of the experimental investigations of triaxial compression test on soil A1 and E2 are presented in _{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

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

sin φ m o b = ( σ 1 − σ 2 ) / ( σ 1 + σ 2 ) . (1)

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

The numerical modelling is performed by means of a discrete element method (DEM) with Particle Flow Code (PFC^{3D}) Itasca (2005) [_{n}, shear stiffness k_{s}) and the sliding law in accordance with Coulomb law (friction coefficient μ) as shown in ^{n} as follows F max s = μ ⋅ 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

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

n i = ( x i [ B ] − x i [ A ] ) / d . (2)

Here, x i [ A ] and x i [ B ] are the position vectors of the centres of balls A and B, and d is the distance between the particles centres:

d = | x i [ B ] − x i [ A ] | = ( x i [ B ] − x i [ A ] ) ( x i [ B ] − x i [ A ] ) . (3)

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

U n = r [ A ] + r [ B ] − d (Particle − Particle) (4)

U n = r [ b ] − 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 = ( k n A ⋅ k n B ) / ( k n A + k n B ) (8)

K n = ( k s A ⋅ k s B ) / ( k s A + k s B ) (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 ‖ − μ ⋅ ‖ F n ‖ ≤ 0 (10)

with μ the particle-particle friction coefficient.

The values of k_{n}, k_{s} 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) [

To generate numerical particle packing, numerous methods exist, Bagi (2005) [

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) [

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 ^{3D} 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.

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. k_{n} = constant and k_{s} = 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

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) [^{n} and tangential stiffness K^{s} between two particles of different size were determined. Hereby, a contact theory for an elastic cylinder was applied. When two particles of radii r_{1}, r_{2} with the same Young’s modulus E and Poison’s ratio ν, are subjected to compression force q per unit depth as shown in ^{n} between the two particles can expressed as follows:

δ = [ 2 q ( 1 − ν 2 ) π E ] [ 3 2 + ln ( 4 r 1 U n ) + ln ( 4 r 1 U n ) ] (11)

( U n ) 2 = [ 8 r 1 r 2 π ( r 1 + r 2 ) ] [ q ( 1 − ν 2 ) Π E ] (12)

Then, K^{n} can be calculated for a particle of radius r as follows:

K n = q δ = π E 2 ( 1 − ν 2 ) ( 1.5 + 2 ln ( 4 r U n ) ) (13)

Considering two small particles with radius r_{1} with the resulting stiffness K 1 n and two larger particles of radius r_{2} with the resulting stiffness K 2 n , such as r_{2} = λ·r_{1}, λ > 1 (

K 2 n K 1 n = 1.5 + 2 ln ( 4 r 1 4 r 2 ( 1 − ν 2 ) q π E ) 1.5 + 2 ln ( 4 r 2 4 r 1 ( 1 − ν 2 ) q π E ) (14)

Or

K 2 n K 1 n = 1.5 + 2 ln ξ 1.5 + 2 ln ( λ 1 / 2 ξ ) (15)

with

ξ = 4 r 1 4 r 2 ( 1 − ν 2 ) q π E (16)

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 _{1} = 0.06 mm and d_{2} = 0.1 mm. Then, the average radius d = ( d 1 + d 2 ) / 2 = 0.08 mm and the average radius r = 0.00004 m. Based on the

The determination of relevant particle parameters i.e. the normal contact stiffness k_{n}, the shear contact stiffness k_{s}, and the friction contact coefficient µ for the used non-cohesive soils has been carried out by means of the discrete element model with PFC^{3D}. The model has been calibrated with the experimental behaviour of real material shown in

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) [

Parameter | Values |
---|---|

Unit weight γ | 26.50 kN/m^{3} |

Young’s modulus E | 13.93∙10^{+6} kN/m² |

Poison’s ratio ν | 0.30 |

Soil E2 | Fine fraction | Coarse fraction | ||
---|---|---|---|---|

Part 1 | Part 2 | Part 3 | Part 4 | |

Grain fraction (mm) | 0.06 to 0.1 | 0.1 to 0.8 | 0.8 to 1 | 1 to 3 |

Average radius (mm) | 0.04 | 0.225 | 0.45 | 1 |

Radius ratio r_{i}/r_{1} | 1 | 5.625 | 11.25 | 25 |

Normal stiffness | k_{n} | 0.85k_{n} | 0.81k_{n} | 0.76k_{n} |

Soil A1 | Fine fraction | Coarse fraction | ||

Grain fraction (mm) | 0.06 to 0.1 | 0.1 to 0.2 | 0.2 to 0.6 | |

Average radius (mm) | 0.04 | 0.075 | 0.2 | |

Radius ratio r_{i}/r_{1} | 1 | 1.875 | 5 | |

Normal stiffness | k_{n} | 0.93k_{n} | 0.85k_{n} |

Hicher (1994) [

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) [

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) [

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 _{n} and the shear stiffness k_{s} 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.

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

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 k_{n} = k_{s} = 10^{5} N/m for soil A1 (

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) [

The results of the discrete element modelling of the large shear test in Coetzee (2016) [

In order to best fit the experimental and the numerical volumetric strain, Belheine et al. (2008) [

Soil | k_{n}[N/m] | k_{s} [N/m] | µ_{e} [-] |
---|---|---|---|

A1 | 10^{5} | 10^{5} | 2 |

k_{n}_{1} = 10^{5} k_{n}_{2} = 0.93∙10^{5} k_{n}_{3} = 0.85∙10^{5} | k_{s}_{1} = 10^{5} k_{s}_{2} = 0.93∙10^{5} k_{s}_{3} = 0.85∙10^{5} | 2 | |

E2 | 10^{5} | 10^{5} | 5 |

k_{n}_{1} = 10^{5} k_{n}_{2} = 0.85∙10^{5} k_{n}_{3} = 0.81∙10^{5}, k_{n}_{3} = 0.76∙10^{5} | k_{s}_{1} = 10^{5} k_{s}_{2} = 0.85∙10^{5} k_{s}_{3} = 0.81∙10^{5}, k_{s}_{3} = 0.76∙10^{5} | 5 |

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 (

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) [

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

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

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 (

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 (PFC^{3D}). 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.

The authors declare no conflicts of interest regarding the publication of this paper.

Ahlinhan, M.F., Houehanou, E., Koube, M.B., 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. https://doi.org/10.4236/gm.2018.84004