DEM simulation of drained triaxial tests for glass-beads

Powder Technology 364 (2020) 123–134

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

DEM simulation of drained triaxial tests for glass-beads Lassaad Hazzar a,⁎, Mathieu Nuth b, Mohamed Chekired c a b c

Conception des aménagements de production, hydraulique et géotechnique, Hydro-Québec, Montreal, Quebec H2L4P5, Canada Department of Civil Engineering, Université de Sherbrooke, Sherbrooke, Quebec J1K2R1, Canada Institut de recherche d'Hydro-Québec (IREQ), Hydro-Québec, Varennes, Quebec J3X1S1, Canada

a r t i c l e

i n f o

Article history: Received 23 April 2018 Received in revised form 1 August 2019 Accepted 30 September 2019 Available online 22 January 2020 Keywords: Triaxial test DEM Glass-beads Stress-strain Volumetric strain Micro parameter Macro parameter

a b s t r a c t The triaxial compression test is the most common laboratory test used to measure the shearing mechanical properties of granular soils. In this paper, a series of numerical drained triaxial compression tests were performed in laboratory with different confining pressures on calibrated glass-beads ranging from 1.25 mm to 2.0 mm in diameter. The confining pressures used for the triaxial tests were 50, 100, and 200 kPa. A model was developed using the Discrete Element Method (DEM) and computed by Particle Flow Code in three dimensions (PFC3D) for the simulation of the triaxial compression tests. The model specimen was an assembly of spherical particles, which were defined by a set of micro parameters, idealized using a stiffness contact model. The influence of several micro parameters as the friction coefficient, normal stiffness and shear stiffness on the stress-strain behavior of the beads was analyzed. Numerical test results show that the micro parameters of the model can greatly influence the macroscopic behavior of the glass-beads. The calibration process was validated by comparing to the experiment results. The relationships between micro parameters and macro parameters were built for further glass-beads analysis. It was showed also that the DEM can accurately capture the very small-strain behavior of the glass-beads and to predict the elastic modulus, Poisson's ratio, angle of friction and dilatation angle. © 2020 Elsevier B.V. All rights reserved.

1. Introduction The glass-beads are a typical granular material characterized by an assembly of particles and a mechanical behavior depending on the size of these particles, their arrangement and particle-to-particle friction. When deformations occur in these granular materials, the external forces may cause internal fabric changes, caused by particles sliding, rolling, and interlocking. Those changes will produce a different response of the material behavior. Understanding of such material response is very important in the design of structures such as retaining walls, foundations systems, and dams, because the analyses of these systems are based on the strength and deformation behavior of the material beneath or adjacent to them. The shearing mechanical properties are the important parameters to characterize the behavior of granular material (glass-beads) and the triaxial test is commonly used to determine these parameters. Some laboratories as well as numerical triaxial tests on glass-beads have been performed in the last few decades. The limited conclusions of conducted research revealed an ambiguity which leads to further investigations. Moreover, relationships between micro parameters and macro parameters for glass-beads are still elusive. In order to investigate the effect of particle size, Kolbuszewski and Frederick [1] used ballotini ⁎ Corresponding author. E-mail address: [email protected] (L. Hazzar).

https://doi.org/10.1016/j.powtec.2019.09.095 0032-5910/© 2020 Elsevier B.V. All rights reserved.

(glass-beads) in their triaxial tests and concluded that the friction angle and the dilatancy increase with the increment of particle size. Three years later, shear tests on soda-lime silica glass-beads were carried out in both solid and hollow cylindrical triaxial devices using either stress-controlled or strain-controlled lateral boundaries [2]. The specimens were sheared in triaxial compression and extension, maintaining the mean stress constant, and volume strains were monitored. The results indicated that the mold used in preparing the specimen created a strong but localized ordering of particles at the boundary and that this type of initial fabric had an important effect on the volumetric behavior of the hollow cylinder specimens. Martinez [3] and Roussel [4] have conducted further, drained triaxial compression tests on glass-beads. Tests results showed that, as the confining pressure increases, the dilatancy tendency of the beads decreases producing a decrease in the peak stress ratio. Additional, it has been noticed that the surface roughness and the particle size significantly affect the stress-strain and volumetric strain behavior of the beads. Recently, Amirpour et al. [5] performed undrained triaxial compression tests on polydisperse spherical glass-beads samples with different particle-size distributions. The experimental results also confirm that shear strength of tested samples were affected by their particle size distribution. Numerical simulation using the discrete element method (DEM, [6]) is a very realistic tool for discrete granular material such as glass-beads, which can displace independently from one another and interact only at contact points. The DEM was introduced by Cundall [7] for the analysis

124

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

of the rock-mechanics problems and later was applied to soil mechanics by Cundall and Strack [6]. Since then, the DEM has been used over a range of applications to model the behaviour of different materials as it allows the micro parameters of the model to be adjusted based on the macro behaviour of the material. Lenaerts et al. [8] used the DEM to simulate the grain straw separation with bendable straw particles and calibrated the micro parameters to match the realistic properties of the straw. Boac et al. [9] provided a summary for the applications of the DEM in simulating the post-harvest operations of wheat, corn, rice, and rapeseed. Nan et al. [10] established the relationship between the micro and macro properties of rod-like particles using the DEM. The DEM is a powerful tool to model granular materials behaviours, such as stress-strain patterns, monitor changes in volume, and strength of soil [11]. In addition to this, the DEM can capture better particle behavior at the micro scale than continuum models especially for glass-beads [12] than that obtained from FEM method [13–15]. It is of great vital to determine which micro parameters should be chosen to achieve the macro behavior from real triaxial tests [16,17]. Much of the previous research works were conducted using 2D modelling, which cannot capture the real stress state of triaxial compression test [18,19]. However, 3D DEM modelling provides a realistic stress state [12,20–22]. In view of the above mentioned issues, this paper presents and discusses the results of a series of 3D DEM analyses in order to show how well a triaxial compression model based on a local discontinuous description, can reproduce the macro behavior of glass-beads in a quantitative manner. These simulations are compared with a series of triaxial laboratory tests performed on glass-beads in the Geotechnical laboratory of the Sherbrooke University (Quebec, Canada). To calibrate the micro properties of the numerical loose and dense beads according to the glass-beads laboratory tests, a numerical testing procedure for the drained triaxial tests was developed. This calibration procedure is a crucial issue. The three main input micro parameters in the numerical model are: the normal contact stiffness, shear contact stiffness and local friction between beads. Several parametric studies were done to evaluate the influence of each of the above parameters on the output macro parameters such as the modulus of elasticity, the Poisson's ratio, the friction angle, the dilatancy angle, the peak and the post-peak regime. It was found that the elastic response of glass-beads is mainly governed by the normal contact stiffness or the ratio of normal to shear contact stiffness's. The friction angle and the dilatancy angle are dependent on inter-particle friction coefficient. In addition, the observed deformation depends strongly on local friction between beads.

When the saturation phase of the test is completed, the maximum back pressures was hold constant and the cell pressure was increased until the difference between the cell pressure and the back pressure equals the desired consolidation pressure. Three samples are tested at different consolidation pressures of 50, 100 and 200 kPa. After final consolidation phase, each sample was sheared under drained condition (CD test). During shearing, the cell pressure σc was kept constant (the specimen is confined by compressing the water in the cell) while advancing the axial load piston downward against the sample cap with strain rate of 0.1 mm/min. It is assumed that no shear stresses occur at the end plates, σc and the axial stress: σa ≈ σc + (F/A) can be taken as the minor σ3 and major σ1 principal stresses, respectively. Fig. 2 shows the state of stress on a cylindrical specimen. The following equations are used to analyze the results of a DTC test:

2. Triaxial tests on glass-beads

Deviator stress q ¼ σ 1 −σ 3

ð1Þ

2.1. Experimental setup

Axial strain ε a ¼ Δh=h0

ð2Þ

A series of drained triaxial compression (DTC) tests were performed on glass-beads to study their shear strength and volumetric behavior. The proposed cylindrical specimen has a diameter of 38 mm and a height of 76 mm. The apparatus used for performing these tests is ELDYN (GDS Enterprise level Dynamic triaxial testing system) as shown in Fig. 1. ELDYN system can utilize to carry out traditional triaxial tests such as UU, CU, CD and monitor the changes in the volumetric strain of the specimen during the test. In these tests, the glass-beads were formed inside a split mold with the help of funnel to control of the height and gain dense samples. Air pluviation preparation method was used to prepare the specimens. Each layer of beads was compacted to achieve the required bulk density (corresponding to a porosity n of 0.34). The formed sample was then placed within the triaxial cell. After assembling the traixial cell and fill it with water, samples were saturated using back pressure. Samples saturation was verified using the pore-pressure parameter B (B = Δu/Δσ). They were saturated under back pressures between 600 and 750 kPa until the B reaches 0.91 +/− 5 [5].

Volumetric strain ε v ¼ ΔV=V 0

ð3Þ

Fig. 1. Triaxial test apparatus used to measure stress-strain of glass-beads.

where: Δh and h0 are the change in height and initial height of the specimen, respectively; and, ΔV and V0 are the change in volume and initial volume of the specimen, respectively. All tests were terminated when strain to failure reached the determined value of 20% axial strain. Data acquisition was used to measure the load, shear stress, volume change and axial deformation. 2.2. Characteristics of the tested glass-beads The glass-beads used in the experimental program were Whitehouse polydisperse reference standards, which are typically spherical shape, clear color, and originally commissioned by the European Community Bureau of Reference (BCR). The principal technical specifications of these tested glass-beads are summarized in Table 1. The DTC tests were carried out for uniform specimens of glass-beads ranging from 1.25–1.40 mm in particle diameter.

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

125

2.3. Experimental results The stress-strain and volumetric strain versus nominal axial strain curves of the glass-beads experienced four stages during the whole process are presented in Fig. 3a and b, respectively. In the stress-strain curves of the four different tests, after the peak is reached, a pronounced peak softening is observed; once the peak stress is reached, the deviator stress q level out smoothly until the critical state condition is reached (Fig. 3a). However, the volumetric strain versus axial strain curves of the four different tests show a continue increase in the volume change throughout the test (Fig. 3b). It is observed that the volumetric dilation is very small. While with the higher confining pressure, the compressed volumetric strain becomes smaller. 3. DEM model simulating the compression triaxial 3.1. Model generation

Fig. 2. Loading conditions in a DTC test.

Table 1 Technical specification of considered glass-beads (SODA-LIME GLASS). Mechanical properties Young modules Rigidity modulus Poisson's ratio

Physical properties 6.89 × 104 N/mm 2.96 × 104 N/mm 0.21

Specific gravity Coefficient of fraction Static Coefficient of fraction dynamic

2.43–2.49 g/cm 0.9–1.0

A 3D DEM numerical model was developed to simulate the triaxial tests mentioned above. Glass-beads were represented by an assembly of spherical particles that move independently of one another and interact only at the contact points. To prepare a model specimen, a cylinder (38 mm × 76 mm) was created using the three-dimensional particle flow code PFC3D [23] cylindrical wall to represent the rubber membrane used in the tests. The cylinder was filled with the model particles generated initially with random orientations (Fig. 4). The DEM model is created in an effort to closely resemble the physical cylinder of the laboratory test. To eliminate the anisotropy of using the uniform particles, the radius of the particles elements were uniformly distributed between Rmin (0.625 mm) and Rmax (0.700 mm). The number of balls generated, N, was 53,097. In PFC3D, the void ratio, e, can be calculated using the total volume, V of the sample minus the total volume of the particles, as follows: 1−

0.7–0.8

e¼

  PN 4 R þ Rmax 3 π min 1 3 2 V

Fig. 3. Experimental (a) deviator stress, and (b) volumetric change vs axial strain of glass-beads under different confining pressures σc (50 kPa, 100 kPa, 200 kPa).

ð4Þ

126

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

ΔF s ¼ K s ∙ΔU s

ð6Þ

where: − − − − −

Fn is the total normal force transferred by the whole system; ΔFs is the incremental shear force; Un is the normal displacement; ΔUs is the incremental shear displacement; Kn is the equivalent normal stiffness of the system of the two particles or springs; − Ks is the equivalent shear stiffness of the system of the two particles or springs;

The composite of equivalent normal stiffness Kn and shear stiffness Ks of the two particles in contact can be expressed as follows [23]: A

Kn ¼

A kn

B

ð7Þ

B

ð8Þ

þ kn

A

Ks ¼

B

kn :kn

B

ks :ks A ks

þ ks

where kAn and kBn are the normal stiffness's and kAs and kBs are the shear stiffness's of the two particles or springs, A and B. These stiffnesses can be expressed as follows: Fig. 4. Numerical model of triaxial test. A;B

kn ¼ 4ERA;B The particles were cycled and allowed to settle and reached equilibrium where the unbalanced force was less than the 1e−3 Newton. 3.2. Constitutive model In order to simulate the interaction between glass-beads particles, a contact stiffness model was adopted in this paper. This contact model is the simplest form of contact for calibration, available in PFC3D [23], and acts like glue between two spheres, connecting them together. Three of the simple contact models, which are shear and normal stiffness, static and sliding friction, and inter-particle cohesion/adhesion, can be utilized. The elastic relationship between the contact force and relative displacement between particles can be provided by the stiffness model [23]. In this model, the cohesive bond connecting two particles A and B of radii RA and RB can be treated as two springs connected in series (Fig. 5). The force–displacement relationships for the whole system and for each particle can be written in the following form: Fn ¼ Kn ∙ Un

ð5Þ

A;B

ks

¼ 4GRA;B

ð9Þ ð10Þ

where E and G are the Young's modulus and shear modulus, respectively, of the material of particles. Introducing the relationships (9) and (10) into the formula (7) and (8), respectively, we obtain the expression for the equivalent stiffness's Kn and Ks in the following forms: K n ¼ 4E

RA :RB ~ ¼ 2ER RA þ RB

ð11Þ

K s ¼ 4G

RA :RB ~ ¼ 2GR RA þ RB

ð12Þ

~ is the harmonic mean of the radii RA and RB. where R In many-particle system, using the average of the radii of all the par~ the equivalent ticles in the specimen 〈R〉 in Eq. (11) or (12) instead of R, tangential and normal contact stiffness components are obtained in the

Fig. 5. Schematics of contact stiffness model: (a) Details of bonded particles, (b) Details of contact model, and (c) Details of normal and shear bond strength.

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

following form: K n ¼ 2EhRi

ð13Þ

K s ¼ 2GhRi

ð14Þ

The average radius is calculated considering all the particles in the assembly by:

127

observing the effect on the macro parameters to find some general trends that could be used for creating a material with certain characteristics. These ideas were then used to find suitable micro parameters to model glass-beads for which triaxial test results were available. 4. Procedure for calibrating micromechanical parameters 4.1. Calibration parameters

p 1 X R hRi ¼ Np i¼1 i

N

ð15Þ

where Np is the total number of particles. When the tensile and/or shear stress at a contact exceeds the tensile and shear strengths of the bond, the bond breaks and separation and/or frictional sliding can occur. The slip condition at each contact point is expressed as: Fr ¼ μ ∙ Fs

ð16Þ

where, − μ is the particle coefficient of friction between the contacting bodies; − Fr is the friction force. 3.3. Modelling procedure The triaxial test was performed in two principal stages: isotropic compression and deviatoric loading. After the DEM specimen was generated and the material properties were assigned referring to the adopted contact-stiffness model, particles were subjected to isotropic compression with a low strain rate over a large number of small timesteps. During isotropic compression, a servo-control mechanism [24] was introduced to maintain the desired confining pressure (σc). Three different σc were considered, i.e., σc = 50, 100, and 200 kPa. The detailed implementation of the code for a stress-controlled wall in PFC3D is given in Dai [25]. The assemblies were considered to be equilibrated when the tolerance of stress obtained from the walls was less than 0.5%, and the ratio of the mean static unbalanced force to the mean contact force was less than 10−5. Fig. 4 shows the final state of the assembly for σc of 200 kPa. After isotropic compression, the assemblies were subjected to vertical compression by a downward displacement of the top wall at a constant velocity, whereas a constant σc acted on lateral walls enforced by the servocontrol mechanism. To ensure quasi-static shearing, the shear rate must be sufficiently slow such that kinetic energy generated during shearing is negligible. This can be formulated in terms of an experiential parameter I, which is defined as the ratio of the mean unbalanced force to the mean contact force in the particle system. For quasi-static shear, I should be less than 10−5 [26]. The constant velocity acting on the top wall was set to 2 mm/s in the simulation, making I less than 10−5 for the entire test. During the loading process: stresses were obtained using the forces on the DEM walls and outer boundaries of the specimen; the strains were calculated directly from wall displacements, the stressstrain and volumetric strain versus axial strain curves were plotted; and the macro parameters (modulus of elasticity, Poisson's ratio, friction angle and dilatancy angle) were deduced using fish functions. Fig. 6 shows the diagram steps for the modelling procedure of the triaxial test.

Given that many microscopic parameters cannot be obtained directly from laboratory experiments, numerical calibration is required [27]. The set of DEM input parameters required to capture the macro-scale response of the glass-beads was selected using the existing knowledge in the literature and methodological trial and error to match the DEM response with the laboratory measurements. The micro parameters are adjusted to simulate the stress-strain volumetric strain versus axial strain curves of glass-beads. According to the previous research results [6,16,28–30], the DEM can be used to study the behavior of real geomaterials, under the appropriate controlled conditions including micro parameters, boundary conditions and loading conditions. During the calibration procedure, the micro parameters introduced in a proposed contact model are selected by comparing the experimental results with the numerical ones. This procedure allows firstly to estimate the values of the micro parameters (kn, kn/ks, and μ) and subsequently to reproduce, respectively, the macroscopic behavior of glass-beads characterized by the modulus of elasticity E, the Poisson's ratio ν, the friction angle ϕ, and the dilatancy angle ψ. These characteristics are the features of the stress–strain and volumetric strain versus axial strain curves deduced from the triaxial tests (see Fig. 7). The Fig. 7 shows four (04) curves: 1. Blue Curve: Deviator stress q versus axial strain εa for dense soil; 2. Green Curve: Deviator stress q versus axial strain εa for loose soil;

3.4. Determining the macro parameters The simulations of the triaxial tests were made with certain values for micro parameters. These micro parameters are normal contact stiffness kn, shear contact stiffness ks and inter-particle friction angle μ. The macroscopic quantities that result from the tests or were applied as boundary conditions, were the three stresses in directions perpendicular to the cell walls σ1 and σ2 = σ3 and the three strains ε1 and ε2 = ε3. Simulations were made by adjusting certain micro parameters and

Fig. 6. Diagram describing the steps of the triaxial test modelling.

128

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

Fig. 7. Typical behavior obtained with triaxial tests.

3. Brown Curve: volumetric strain εv versus axial strain εa for dense soil; 4. Pink Curve: volumetric strain εv versus axial strain εa for loose soil; 4.1.1. Prediction of elasticity modulus E. The method used to determine modulus of elasticity, E, was described by Yang et al. [31]. Shear strength, τs was determined at 15% strain. The axial strain at 50% of the deviatoric stress was found and this point is highlighted on the curve. A line was drawn through the origin and this point. The slope of this line is found by regression and was used to determine E. The criterion for the prediction of E is illustrated in Fig. 8a. 4.1.2. Prediction of Poisson's ratio ν Poisson's ratio ν is a measure of the ratio of transverse strain to longitudinal strain resulting from a change in normal stress under compression [32]. When ν approaches 0.5 in isotropic materials it indicates the

incompressibility of the material whereas when it approaches 0, it indicates the ease with which a material can be compressed. The change in the water level in the burette on the triaxial test apparatus was equivalent to the volume change of the specimen as the test was confined by water. The ν was calculated by the following equation: ν¼

εa −εv 2ε a

ð17Þ

4.1.3. Prediction of friction angle ϕ The shear strength of soils is usually represented using MohrCoulomb theory represented by Mohr circles. Those circles represent the state of stresses of a soil specimen in the plane that contains major stress, σ1, and minor stress, σ3. If Mohr circles of different specimens of the same material subjected to different confining stresses are

Fig. 8. Criteria for calculating of: (a) Young's modulus, E, and (b) friction angle, ϕ.

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

129

Fig. 9. Effect of the particle friction coefficient (μ) on the numerical macro performance of glass-beads: (a) Stress-strain, and (b) Volumetric change versus axial strain under confining pressure 200 kPa.

drawn together, then the friction angle, ϕ, of the material can be estimated from the slope of the line tangent to the circles, known as failure envelope (Fig. 8b). However, depending on the engineering problem under consideration, either the peak friction angle or the constant volume friction angle is needed. Those angles can be calculated using the Mohr-Coulomb failure criterion (for granular material):  σ1 ¼ σ3

1 þ sinϕ 1− sinϕ

Hence, the peak friction angle can be calculated as: ϕp ¼ sin−1



 ð19Þ p

where the values of σ1 and σ3 are taken at the peak. Similarly, the constant volume friction angle can be calculated as:

 ð18Þ

σ 1 −σ 3 σ1 þ σ3

ϕcv ¼ sin−1



σ 1 −σ 3 σ1 þ σ3

 ð20Þ cv

Fig. 10. Effect of the ratio kn/ks on the numerical macro performance of glass-beads: (a) Stress-strain, and (b) Volumetric change versus axial strain for loose test under confining pressure 200 kPa.

130

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

where the values of σ1 and σ3 are taken at the constant volume stage. 4.1.4. Prediction of dilatancy angle ψ Dilatancy can be defined as the volume change associated with the application of shear stresses. An increase in volume, or expansion, is known as positive dilation, while a decrease in volume, or contraction, is known as negative dilation. The amount of dilatancy that a granular material can experience is dependent on the particle interlocking, which depends on the fabric of the material. The dilatancy angle ψ can be estimated from the volumetric strain versus axial strain curve of a material subjected to CTC with the following expression [33]: ψ ¼ sin−1

    .   εv εv εv ¼ sin−1 − 2þ γ εa εa

ð21Þ

particle friction coefficient (see Fig. 9a). The same results can be observed from simulation research conducted by Cui and O'Sullivan [34] and experimental results obtained by Abriak and Caron [35]. Fig. 9b shows that the volumetric strain increases as the friction coefficient rises. In fact, particles are more difficult to move with higher friction coefficient, so the specimen volume changes much more easily than system with lower friction coefficient. It is important to mention that the friction coefficient between boundary walls and particles is another important parameter, which influences the global shear behavior of materials. However, this parameter was set to zero in many references in order to avoid the influence of boundary walls. Wu et al. [36] highlighted the influence of friction coefficient between boundary walls and particles in the deformation of samples. Their research showed that the friction coefficient has a slight influence on the volumetric strain.

4.2. Calibration methodology During the calibration procedure, the microscopic parameters that were introduced in the model (k n, k n/ks, and μ) were selected by comparing the experimental results with the numerical ones. However, some of the parameters can have a cross effect on the macroscopic ones (for example μ on the maximum σpeak). Thus, a calibration methodology is proposed to overcome this cross dependency. This methodology is: 1. The normal contact stiffness kn and or the ratio kn/ks were varied to match the deformation modulus and Poisson's ratio of the real material while all other parameters of the test were kept constant. It is worth noting that contact stiffness has a significant effect on the deformation characteristics. Young's modulus and the Poisson's ratio were determined from an initial small loading with a maximum axial strain of 1%. 2. Using the same assembly, a series of simulation tests were carried out with the appropriate elastic local parameters. The inter-particle friction coefficient was varied to adjust the dilatancy curve to that of the real material, while all other parameters were kept constant. 5. Effect of micro parameters on the macro performance of glassbeads With the set up as described above, there are the two linear springs in a contact between particles, one for the normal stiffness given by the contact Kn and one for the shear stiffness given by Ks. The normal force is proportional to the normal stiffness and the overlap of two particles. The shear force is proportional to the shear stiffness and the relative rotation between two particles. The effects of the micro parameters, represented by the inter-particle friction coefficient μ and the ratio of normal to shear stiffness kn/ks were discussed in the following subsections. 5.1. Effect of inter-particle friction coefficient μ The Friction coefficient between particles μ is a physical parameter dependent on the mineral, water content and chemical composition etc. and very hard to measure. This parameter plays an important role in DEM analysis. The impact of inter-particle friction coefficient μ on the global shear behavior of glass-beads was investigated. The ratio of normal to shear stiffness k n /k s was set to be 1.0. For a confining pressure σ c of 200 kPa, Fig. 9 shows the influence of friction coefficient corresponding to 0.1, 0.2, 0.3, and 0.4 on the global deviator stress q and volumetric strain ε v versus axial strain εa. It can be observed that the peak stress and the residual deviatoric stress increase with increasing inter-

Fig. 11. Effect of local properties at particle contact on macro-properties of glass-beads for a confining pressure of 200 kPa: (a) Effect of the ratio kn/ks on Young's modulus E and Poisson's ratio ν, and (b) Effect of friction coefficient μ on friction angles ϕp and ϕcv and dilatancy angle ψ.

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

131

Fig. 12. Comparisons of (a) deviatoric stress and (b) volumetric strain between experimental results and numerical results under different confining pressures.

5.2. Effect of normal to shear stiffness ratio kn/ks A series of simulations were performed for varying normal to shear stiffness ratio kn/ks. The friction coefficient was kept constant at 0.2 and normal to shear stiffness ratio were set to 1, 2, 5, and 10.0. For a confining pressure of 200 kPa, Fig. 10 shows the influence of normal to shear stiffness ratio on the global deviator stress q and volumetric strain εv versus axial strain εa. Fig. 10a shows that the elastic stiffness increases significantly with the increasing normal to shear stiffness ratio kn/ks (i.e. normal stiffness kn) of glass beads. In addition, the peak strength increased slightly, thus the angle of internal friction increased slightly. On the other hand, Fig. 10b shows that the volumetric strain increases very slightly as the ratio kn/ks rises.

5.3. Relationships between micro and macro parameters The initial Young's modulus E is calculated by the elastic stage of stress strain curve (see Fig. 4 and Fig. 5a) and the Poisson's ratio ν, which plays a significant role in the influence on the triaxial tests, is calculated according to Eq. (4) and Fig. 4. Fig. 11a shows the relationship between the elastic macroscopic parameters (E and ν) and the normal to shear stiffness ratio kn/ks for a confining pressure of 200 kPa. It can be seen that the Young's modulus and the Poisson's ratio decrease with the ratio k n /ks in a non-linear manner. A little ratio of kn/ks represents high ability of tangential deformation resistance for particles in contact, and the resemble effect is to increase the lateral deformation for the same axial deformation, so Poisson's ratio is upper. The resistance in lateral

Table 2 Microscopic parameters of the discrete element model.

generates more stable force network in the axial direction, which makes the sample stiffer in the axial direction. Considering that, the cohesion of glass beads was zero, internal friction angle at the peak (ϕ p ) and at the residual shear stress (ϕcv) can be calculated according to Eqs. (6) and (7), respectively. On the other side, the angle of dilatancy (ψ) is deduced based on the Eq. (8). The variations of these angles (ϕp, ϕcv, and ψ) with different friction coefficient μ under the confining stress of 200 kPa are shown in Fig. 11b. It is demonstrated in Fig. 11b that the three angles (peak friction angle, residual friction angle, and dilatancy angle) increase with inter-particle friction coefficient in a non-linear manner. Particle contact with a larger inter-particle friction coefficient exhibits stronger resistance to sliding, so only a larger force is able to force particles to slide. Meanwhile, a larger rotation would happen for particle contact with larger interparticle friction before sliding happens, which causes a higher dilatancy. A suitable particle friction coefficient can be chosen for further simulations when the desired peak and residual strengths are matched with the laboratory results under different confining stresses. 6. Comparison between the numerical model and the experimental results During the numerical triaxial test, the model parameters (interparticle friction coefficient μ, normal to shear stiffness ratio kn/ks, and normal contact stiffness kn), were chosen by comparing the experimental results with the numerical ones. Laboratory tests are available Table 3 Macroscopic parameters values of the glass-beads for each experimental test. Parameters

Parameters

Mean value

Normal contact stiffness kn (Pa) Ratio kn/ks Inter-particle friction coefficient μ Damping constant

18.2 × 108 4.70 0.13 0.7

Young's modulus E (MPa) Poisson's ratio ν Peak friction angle ϕp (°) Constant volume friction angle ϕcv (°) Dilatancy angle ψ (°)

Tests 50

100

200

28.49 0.25 27.75 23.13 14.30

53.65 0.25 26.16 21.44 13.95

95.84 0.28 25.63 21.54 12.55

132

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

Table 4 Macroscopic parameters values of the glass-beads for each numerical test. Parameters

Tests

Young's modulus E (MPa) Poisson's ratio ν Peak friction angle ϕp (°) Constant volume friction angle ϕcv (°) Dilatancy angle ψ (°)

50

100

200

27.42 0.245 26.40 19.83 14.7

56.98 0.26 24.90 19.30 13.83

97.71 0.29 25.79 22.31 11.92

under different confining pressures (50, 100, and 200 kPa). In the last section, several attempts were made to evaluate the influence of the input microscopic parameters on the deformation and the characteristic strength.

For each of the three micro parameters (kn, ks, and μ), their effects on the macro scale have been described. The results suggested that a glass-beads material with a certain elastic stiffness and angles of internal friction and dilatancy could be produced by adjusting the micro parameters. This meant that a certain value of the internal friction could be reached by adjusting the friction coefficient μ and then the desired elastic stiffness could be reached by adjusting the micro normal to shear stiffness ratio kn/ks. The experiment curves were fitted with DEM simulation results shown in Fig. 12 with the parameters reported in Table 2 for confining pressures of 50, 100 and 200 kPa. Fig. 12a shows the deviatoric stress via axial strain curves obtained from numerical tests are basically consistent with the results obtained from the laboratory for each confining pressure. The glass-beads is in the elastic stage when axial strain is less than 3%, which the

Fig. 13. Variation of: (a) deformation modulus, (b) Poisson's ratio, (c) friction angle, and (d) dilatancy angle, versus confining pressure.

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

relationship between stress and strain is linear. The stress-strain curve is parabolic relationship which the model is in the elasticplastic transition stage with the gradual increase in axial strain. It demonstrates that the cracks have emerged which lead to the shear strength decreases, relatively. There is a critical state when the value of shear strength achieves the maximum that the cracks penetrate through the shear zone. It is observed that the results of the laboratory and numerical simulations almost have the same peak shear strengths. However, the strain-softening phenomenon occurs, however, in post peak regime when confining pressure is 50 kPa or 100 kPa, that the residual strength is less than that in laboratory. There are several reasons that may explain this difference. First, the shear strength of the specimen is mainly based on friction after the shear surface formed. The confining wall of the model is an unchangeable cylinder which is different from the rubber material that in laboratory. It leads to different normal stress on the shear surface under the same confining stress. Second, as compressing progressed, particle friction coefficient may change due to particle rotation and re-arrangement in the shear surface in the DEM analysis, which will directly affect the residual shear strength. Fig. 12b shows the volumetric strains for experimental and numerical tests under different confining pressures. It can be seen that the numerical volumetric strains were very close from experiments up to 5%. Above this value, the numerical volumetric strain was systematically larger than experimental tests. Overall, the results indicate that the non-linear stress–strain behavior of glass-beads including dilatancy is covered by the numerical model. This means that the model is predictive. Tables 3 and 4 show the macroscopic parameters for experimental and numerical, respectively, glass-beads. It is observed that the deformation modulus increases with increasing confining pressure, and the experimental and numerical parameters are much closed. In Fig. 13a and b, the variation of the deformation modulus, Poisson's ratio versus the confining pressure are plotted. The deformation modulus depends strongly on the confining pressure, unlike the Poisson's coefficient, which does not vary much with confinement. Fig. 13c and d gives the evolution of the macroscopic friction angles (peak friction angle and constant volume friction) and the dilatancy angle versus confining pressure. The friction angles decreases with increasing confining pressure (50–100 kPa) by approximately 1.5 deg., then a slight increase for 200 kPa, which is a phenomenon also observed in physical samples. The same observation is noted for the dilatancy angle. 7. Discussions of results The objective of the above calibration process was to create a numerical glass-beads sample, which would reproduce the mechanical response of SODA-LIME glass-beads (Table 1) and could directly be used to solve the boundary problem under consideration. The calibration was carried out using triaxial laboratory test results. Simulations are done under the same test conditions as experiments (porosity, boundary condition and loading). Using the calibration procedure with DEM proposed in Section 3, we could calibrate the numerical model of SODA-LIME glassbeads. Several parametric studies were performed. The results presented in Fig. 11 are the results of this methodology. It is worth noting that the macroscopic response in the elastic region is mainly influenced by the normal to shear stiffness ratio parameter. The initial modulus E is independent of the microscopic friction coefficient μ. The Poisson's ratio decreases when the ratio between normal stiffness and shear stiffness decreases. The same observation was made in Collop et al. [37]. A coupled effect between the inter-particle friction coefficient, and the normal to shear stiffness ratio on the peak and post-peak regime was also

133

found. The inter-particle friction coefficient governs the friction and dilatancy curves. It is observed that the failure stress increases with increasing confining pressure. The initial stiffness is also found to increase with an increase in the confining pressure. Poisson's ratio varies slightly with the confining pressure: this dependence could be due to the variation of the number of contacts. The peak friction angle does not vary much in the stress range studied, as in Fig. 13c. The variation lies between 26.40 and 25.80 for confining pressure from 50 to 200 kPa. The variation of dilatancy angle with confining pressure is displayed in Fig. 13d. Once again, this quantity does not vary much in the studied stress range and it falls in the range measured in physical material. The numerical and experimental results are in good qualitative and quantitative agreement; as seen in Fig. 12. 8. Conclusions DEM simulations were carried out to simulate glass-beads triaxial testing. To reach this goal, a calibration procedure is necessary. Some conclusions can be summarized as follows: − A calibration methodology for numerical triaxial tests of a glassbeads material is proposed to be followed. Firstly, choose an appropriate stiffness ratio so that we can have required Poisson's ratio. Secondly, with stiffness ratio fixed, we start to search for contact modulus by comparing Young's modulus with experiments. Finally, we try to find a particle friction coefficient so that the peak friction angle and dilatancy angle are in agreement with experiments. − In spite of its simplicity, the DEM model is capable to reproduce the most important macroscopic properties of glassbeads materials without necessitating a description of the granular structure perfectly. In terms of the deformation properties, a good agreement between the experimental and the numerical material is obtained (the chosen microscopic values of kn = 18.2 × 108 Pa and kn/ks = 4.7 give macroscopic values E = 97.7 MPa and ν = 0.29 for a numerical material which corresponds well with the macroscopic values E = 95.8 MPa and ν = 0.28 for the physical material). A non-linear relationship between the elastic parameters at the micro and macro scales is found. − The corresponding shear strength parameters agree well. The friction angle (25.79°) of the numerical sand sample, which was determined by a linear failure envelope, is slightly bigger than that of the physical material (25.63°). When calibrating the numerical material, a considerable influence of the interparticle friction coefficient on the volumetric strain was found. This article focuses on the study of shear behavior of homogeneous glass-beads. However, these granular materials are often very heterogeneous as they are composed of particles of different sizes. In the next step, the shear behavior of heterogeneous glass-beads will be explored with the numerical tool developed in this paper. References [1] J. Kolbuszewski, M.R. Frederick, The significance of particle shape and size on the mechanical behavior of granular materials, Proc. Euro. Conf. Soil Mech. Found. Eng. Paper 9 (1963) 253–263. [2] I. Ishibashi, J.T. Jenkins, J.W. Choi, L. Clifton, I.V. Parker, The influence of boundaries on the volumetric behavior of solid and hollow cylindrical specimens of glassbeads, Soils Found. 36 (2) (1996) 45–55. [3] B.N. Martinez, Strength properties of granular materials. MSc. Thesis, Dep. Civil Environ. Eng., Louisiana State University and Agricultural and Mechanical College, USA, 2003. [4] L.E. Roussel, Experimental investigation of stick-slip behavior in granular materials. MSc. Thesis, Dep. Civil Environ. Eng., Louisiana State University and Agricultural and Mechanical College, USA, 2005. [5] A. Samaneh, M.N. Hussien, M. Karray, M. Chekirad, V. Roubtsova, Laboratory investigation on shear strength of granular material considering the effect of particle size distribution, Proc. 66th Can. Geo. Conf., 2013.

134

L. Hazzar et al. / Powder Technology 364 (2020) 123–134

[6] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique 29 (1) (1979) 47–65. [7] P.A. Cundall, A computer model for simulating progressive large-scale movements in block rock mechanics, Proc. Symp. Int. Soc. Rock Mech (1971) 2–8. [8] Lenaerts, B., Aertsen, T., Tijskens, E., De Ketelaere, B., Ramon, H., De Baerdemaeker, J.: Simulation of grain-straw separation by discrete element modeling with bendable straw particles. Comput. Electron. Agric. 101, 24–33 (2014). [9] J.M. Boac, R.P.K. Ambrose, E. Mark, M.E. Casada, R.G. Maghirang, Applications of discrete element method in modeling of grain postharvest operations, Food Eng. R 6 (4) (2014) 128–149. [10] W. Nan, Y. Wang, Y. Liu, H. Tang, DEM simulation of the packing of rod-like particles, Advan. Pow. Tech 26 (2) (2015) 527–536. [11] T.G. Sitharam, J.S. Vinod, L. Rothenburg, Shear behavior of glass beads, Proc. Inter. Conf. Micro. Gran. Med., Powder and Grains (2005) 257–260. [12] S.G. Bardenhagen, J.U. Brackbill, D.L. Sulsky, A numerical study of stress distribution in sheared granular material in two dimensions, Phys. Rev. E 62 (2000) 3882–3890. [13] J.E. Andrade, C.F. Avila, S.A. Hall, N. Lenoir, G. Viggiani, Multiscale modeling and characterization of granular matter: from grain kinematics to continuum mechanics, J. Mech. Phys. Solids 59 (2011) 237–250. [14] C. Fang, W. Wu, Y.J. Lin, Y.H. Chen, T. Doanh, On sintering of tiny glass beads in oscillating diametric compressions, Granul. Matter 18 (4) (2016) 1–18. [15] Y. Zheng, S. Deng, L. Vantuan, S. Yujie, Micro-mechanism of the intermediate principal stress effect on the strength of granular materials, Soil Behav. Geomech 236 (2014) 465–475. [16] L.T. Shao, G. Liu, F.T. Zeng, Recognition of the stress-strain curve based on the local deformation measurement of soil specimens in the triaxial test, Geotech. Test. J. 39 (4) (2016) 658–672. [17] T. Hiroshi, I. Kazuyoshi, Numerical simulation of triaxial test using two and three dimensional DEM, Proc. Powders and Grains (2001) 177–180. [18] D. Darius Markauskas, R. Kačianauskas, Compacting of particles for biaxial compression test by the discrete element method, J. Civ. Eng. Manag. 12 (2) (2006) 153–161. [19] N. Belheine, J.P. Plassiard, Numerical simulation of drained triaxial test using 3D discrete element modelling, Comput. Geotech. 36 (2009) 320–331. [20] L. Hazzar, N. Mathieu, 3D distinct element modeling of mechanical lab tests for unsaturated granular materials, Proc. CSCE Annual Conf., Regina, Saskatchewan, Canada, Paper GEN-234, 2015. [21] K. Cheng, Y. Wang, Y.B. Mo, Q. Yang, Y. Guo, Determination of micro-parameters in DEM through the whole surface deformation measurement in the triaxial test, Proc. 7th Inter. Conf. Dis. Elem. Meth (2017) 1429–1439.

[22] J.M. Ting, M. Khwaja, L.R. Meachum, J.D. Rowell, An ellipse-based discrete element model for granular materials, Num. Analy. Meth. Geomech 17 (9) (1993) 603–623. [23] Itasca Consulting Group, Particle flow code (PFC3D), Version 5, vol. 0, Minneapolis, MN USA, 2013. [24] A. Di Renzo, F.P. Di Maio, Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes, Chem. Eng. Sc 59 (3) (2004) 525–541. [25] B.B. Dai, Micromechanical investigation of the behavior of granular materials. Phd. Thesis, University of Hong Kong, China, 2010. [26] S. Masson, J. Martinez, Micromechanical analysis of the shear behavior of a granular material, J. Eng. Mech. 127 (10) (2001) 1007–1016. [27] T. Koyama, L. Jing, Effects of model scale and particle size on micro-mechanical properties and failure processes of rocks–a particle mechanics approach, Eng. Anal. Bound. Elem. 31 (5) (2007) 458–472. [28] C. Thornton, Numerical simulations of deviatoric shear deformation of granular media, Géotechnique 50 (1) (2000) 43–53. [29] T.G. Sitharam, S.V. Dinesh, N. Shimizu, Micromechanical modeling of monotonic shear behaviour of granular media using three dimensional DEM, Int. J. Numer. Anal. Methods Geomech. 26 (2002) 1167–1189. [30] T.G. Sitharam, J.S. Vinod, Critical state behaviour of granular materials from isotropic and rebounded paths: DEM simulations, Granul. Matter 11 (2009) 33–42. [31] G. Yang, T. Yu, H. Liu, Numerical simulation of undrained triaxial test using 3D discrete element modelling, Geotech. Special Publ. ASCE. (2011) 222. [32] A.M. Mouazen, H. Ramon, J. De Baerdemaeker, SW-soil and water: effects of bulk density and moisture content on selected mechanical properties of sandy loam soil, Biosyst. Eng. 83 (2) (2002) 217–224. [33] T. Schanz, P.A. Venneer, Angles of friction and dilatancy of sand, Géotechnique. 46 (1) (1996) 145–151. [34] L. Cui, C. O'Sullivan, Development of a mixed boundary environment for axisymmetric DEM analyses: proc, 5th Int. Conf. Micromech. Granul. Media (2005) 1–5. [35] N.E. Abriak, J.-F. Caron, Experimental study of shear in granular media, Adv. Powder Technol. 17 (2006) 297–318. [36] K. Wu, P. Pizette, F. Becquart, S. Rémond, N.E. Abriak, W. Xu, S. Liu, Experimental and numerical study of cylindrical triaxial test on mono-sized glass beads under quasistatic loading condition, Adv. Powder Technol. 28 (1) (2017) 155–166. [37] A.C. Collop, G.R. McDowell, Y.W. Lee, Modelling dilation in an idealised asphalt mixture using discrete element modelling, Granul. Matter 8 (2006) 175–184.