Abstract
In this study, a 3D simple shear model using DEM is built based on the boundary condition of an NGI-type bidirectional simple shear apparatus. Stack of rings used as lateral constraints in a bidirectional simple shear test is modelled by layers of clumps which is possible to be moved by particles; different contact types and parameters are used to model the sand-loading caps, sand-latex membrane, and sand-sand contacts. A simple shear test using the bidirectional simple shear apparatus is performed for the calibration of the 3D DEM simple shear model. By analyzing the simulation results, the following can be concluded. (1) Rings generated by clumps can provide an accurate boundary condition, effective in computation since no contact force is needed for a clump. (2) In the simulation, the orientation of average contact force changed dramatically during shear. It is in the vertical direction (90°) before shear and changes to 45° at 40% shear strain. No shear band is observed which is consistent with the test, and particles move uniformly. (3) In the simulation, the degree of noncoaxiality is the greatest at the beginning of shear, and it is decreased during shear. However, the degree of noncoaxiality is still large at 20% shear strain where there is a 10° difference between the rotation angle of principal stress and principal strain increment.
1. Introduction
The NGI-type simple shear test is mostly used in geotechnical testing of soils, and its one characteristic is using a stack of rings as a lateral constraint for samples. The stack of rings is usually made of frictionless Teflon coated rings, and the rings have different lateral displacements during shear. When two neighboring rings have a large displacement, the sample will be shear to failure at this plane. This characteristic gives the possibility to explore the shear failure plane and gives a more reliable shear behavior, which is different from another widely used test-direct shear test which has a fixed shear plane.
In past decades, the discrete element method (DEM) has been used by researchers to study the mechanical behaviors of granular materials from mesoscope. Different from the finite element method (FEM), DEM can be applied when large displacement occurs in the model. Compared with laboratory tests, numerical simulation could provide more testing options and useful results in different ways. For example, some complex stress paths, such as shape figure-8 oval stress path, which cannot be simulated by unidirectional simple shear tests, can be simulated by numerical simulations. In addition, lateral stress that cannot be measured in simple shear tests can be easily determined in simulated tests. Knowing its advantages in enriching laboratory testing possibilities, many researchers used simulated laboratory tests to study the physical and mechanical behaviors of materials [1, 2].
However, due to the complexity and difficulty of building a reliable simple shear model using DEM, especially defining the contact model parameters and boundary conditions in a DEM model, many previous studies cannot accurately reflect laboratory tests.
In previous DEM studies, the agreement of experimental results and simulated result is widely used to validate DEM models [3–5]. Asadzadeh and Soroush [6] studied the shear behavior of glass bead with drained simple shear test. Parameters, such as the friction between boundaries and particles and the friction between particles, are selected based on relevant studies. Young’s modulus is selected from the model that gives the best fit of shear stress-strain behavior of laboratory tests. Results showed that Coetzee and Els [7] studied the calibration process of a blade moving through granular material. Two laboratory tests, shear tests, and compression tests are used to determine the internal friction angle and stiffness of granular material, respectively. Results show that shear behavior is related to particle friction coefficient and the particle stiffness, and compression behavior is related to particle stiffness. Experimental tests are also used to get certain particle contact parameters. Bernhardt et al. [5] studied the shear behavior of steel spheres under drained simple shear test. Dondi et al. [8] combined triaxial test and DEM simulation to calibrate particle contact parameters with consideration of particle shape effect. The inner-particle friction value is determined by the method described by Cavarretta et al. [9]. Ball-wall friction values are measured from tilt tests. By comparing the shear stress-strain and shear strain-volumetric strain behavior of DEM with laboratory tests, it is concluded that increasing particle numbers in DEM will give a better approximation of shear stress-strain and shear strain-volumetric strain behavior of laboratory tests.
Based on the above studies, the model parameters for a simple shear test simulation were well calibrated. However, the boundary condition of an NGI-type simple apparatus is complex and difficult to model in commercial DEM software, such as PFC 3D in which a wall (boundary condition) cannot be moved by particles.
In previous studies, the stack of rings in a simple shear model is usually modelled as two rotating walls in 2D models and a stack of moving short hollow cylindrical walls in 3D models. Gutierrez and Muftah [3] used left and right frictionless wall to model the rings in a 2D model, and the walls are confined by semirigid walls and minor principal stress. During shear, the two walls rotate at a constant angular velocity. Qian et al. [4] used left and right walls to model the rings in a 2D model, and the walls rotate during shear. Asadzadeh and Soroush [6] modelled the rings as small rigid walls in a 3D model, and the walls have linearly varied rate of shearing during shear in which the top ring has a zero velocity and the bottom ring has a velocity of 0.1 mm/s. Bernhardt et al. [5] modelled the rings as 10 layers of walls in a PFC 3D model, and the bottom wall has a fixed velocity. Zhang et al. [10] used cylindrical sidewalls to model the rings in a PFC 3D model. It should be noted that ring movement methods are not described in detail by Bernhardt et al. [5] and Zhang et al. [10]. Since the walls in PFC 3D cannot be passively moved by particles, it is reasonable to assume that ring movement in the two studies is similar to the method described by Asadzadeh and Soroush [6].
All the above studies indicated that DEM simulation is useful to study the simple shear test with acceptable results consistency. However, the boundary conditions in these DEM models were different from those in laboratory tests since the velocity of different walls was set to constant values. In such boundary condition, the failure plane may also be changed because the specimen was subjected to not only shear force caused by the top and bottom plate but also radial forces caused by rings.
Due to the above problems in existing studies, this study aims at developing a 3D distinct element model for simple shear tests that can better reflect its boundary condition and mechanical behavior of tested material. The model is built using PFC 3D based on the VDDCSS—the first commercially available bidirectional simple shear apparatus as shown in Figure 1. More details of this testing apparatus are described by Li et al. [11–13]. Different from previous models, the steel rings are simulated by “clump” elements in PFC 3D, which can be passively moved by particles within the specimen.
(a)
(b)
(c)
2. DEM Model
In previous studies, the linear contact model [7, 12], linear rolling resistance model [3, 9], and Hertz–Mindlin nonlinear model [5, 6, 14] are widely used to simulate granular materials in DEM simulation of simple shear tests. Among them, the linear contact model is the simplest contact model in DEM. The linear rolling resistance model is used in this study which is built based on the linear model since the shape effect can be partly considered in the linear rolling resistance model [3].
In the linear model, the contact interface cannot resist rotation, and it is divided into two linear components and two dashpot components for a ball-ball contact. Two linear components act perpendicular to each other, and two dashpot components are parallel with two linear components, as shown in Figure 2. Linear elastic and frictional behaviors are considered by the linear component, and viscous behavior is considered by the dashpot component, as shown in Figure 2. The linear elastic is simulated by two perpendicular linear springs. The one in the direction connecting the center of two particles has a constant normal stiffness, , and the one perpendicular to it has a constant shear stiffness, . The dashpot force is simulated by two perpendicular dashpots indicated by the normal critical-damping ratio and shear critical-damping ratio . The one in the direction connecting the center of two particles has a constant normal critical-damping ratio, and the one perpendicular to it has a shear critical-damping ratio. The two red lines denote notional surfaces, and the distance between the two notional surfaces is called surface gap, . Contact between two particles will be activated when the surface gap is equal to or smaller than zero.
Many micromechanical behaviors can be reflected by considering rolling resistances, such as the steric effect caused by particle roughness or nonsphericity at the contact point. The linear rolling resistance model has the behavior similar to that in linear model, expect the linearly increased internal moment. The internal moment increases till a maximum value equal to current normal force multiplied by a rolling resistance coefficient and effective contact radius is reached [15].
In the VDDCSS, the specimen is formed within a latex membrane, and a stack of low-friction Teflon coated rings, 1.16 mm high each, is placed outside the testing specimen. In the 3D DEM model, the specimen is formed within layers of clump rings, each being 1.16 mm high. The clump ring is considered as a rigid wall, and it is fixed in Z axis while being free in XY plane. There is no contact between two clump rings, replicating the low friction of the Teflon coated rings, as shown in Figure 3. It should be noted that the wall-type boundaries with fixed velocity were used in previous studies because additional forces act on neighboring particle that push the particles with the fixed velocity. As a result, the particle movement is predetermined, and the measured stress is interfered by additional forces. This is why the clump with free movement in XY plane is used in this study.
(a)
(b)
In the VDDCSS, the top cap and bottom cap that are in contact with specimen during a test have several ring-shape upheavals, as shown in Figure 4. Those are used to increase the friction between specimen and loading facilities. In the DEM model, top and bottom planes are used to model the top cap and bottom cap, and the effect of upheaval is considered in the contact parameters.
In the VDDCSS, a cylindrical specimen with 70 mm diameter and 17.5 mm height (before consolidation) is tested. In the DEM model, a cylindrical specimen with 70 mm diameter is formed with particles with the same grading curve as experiment, as shown in Figure 5 [9]. In trial 3D DEM tests, a desired porosity cannot fully fill the cavity (simulated mould), so a smaller porosity is used to generate particles. Particles are generated by a given porosity of 0.37 in a box with 70 mm width and 22 mm height using radius expansion method. 14061 generated particles are then mobilized and balanced by gravity, as shown in Figure 6.
(a)
(b)
(c)
In the VDDCSS, the specimen is consolidated by moving downwards the top cap till target vertical stress is reached. In the DEM model, the specimen is consolidated by moving downwards the top plane till target vertical stress is reached. Servo wall function is used to move the wall with continuous calculation of current vertical stress. When the calculated vertical stress is close to the target vertical stress within a given tolerance value, the servo system is stopped.
In a VDDCSS drained simple shear test, a specimen is sheared by moving bottom plate in the horizontal direction while the top plate is used to maintain a fixed vertical stress by moving up and down. In a 3D DEM model, a specimen is sheared moving bottom plane in the horizontal direction while the top plate is used to maintain a fixed vertical stress, as shown in Figure 7. To maintain a fixed vertical stress in a drained test, the top plane moves up and down by wall servo function and continuously calculate the current vertical stress on the top plate.
(a)
(b)
The stress is used to control the drained test by maintaining a constant vertical load while adjusting the position of top plane. Five measurement balls are created to measure stresses within defined volume. A measurement ball with a radius equal to 0.35 times of sample height is created in the center of the sample, and four measurement balls with the same size are created next to the center one along X and Y directions, as shown in Figure 8. The measurement balls are in the center of the sample. This is due to the fact that previous DEM studies showed that stresses in the center of the sample are uniform [6], while other parts have large boundary effect. Shear stresses and normal stresses are measured by this method. Shear strain is calculated by dividing shear displacement by sample height. Timestep is fixed at . Trial tests indicated that the model only requires reasonable computation time and is suitable for further simulations.
3. Model Parameters
The calibration used the experimental results from a drained simple shear test. Leighton Buzzard fraction B was tested via the dry funnel method [11]. The specimen has a void ratio of 0.30 and vertical stress of 200 kPa after consolidation. The measured shear stress-strain behavior and shear strain-volumetric strain relation are used to calibrate the parameters in DEM.
In the 3D DEM model, three types of contacts are used, which are ball-ball contact, ball-pebble contact, and ball-facet contact. Ball-ball contact is used to define the contact between sand particles, and it is the main part in the calibration. Three parameters are separately studied, which are friction coefficient, Young’s modulus, and rolling resistance coefficient. Values that give the best presentation of tested shear stress-strain behavior and shear strain-volumetric strain relation are used in following studies, and the details of the calibration are omitted.
The ball-pebble contact is used to define the contact between sand and clump rings which has a small friction coefficient. Asadzadeh and Soroush [6] used the friction coefficient of 0 for similar contact. However, since the friction between membrane and sand is not zero in the laboratory test, a small friction coefficient of 0.1 is used in this study. The ball-facet model is used to define the contact between sand particles and loading walls (top plate and bottom plate) which has a large friction coefficient. A large friction coefficient equal to 0.9 is used in this study, as shown in Figure 9. In the model, local damp ratio is 0.7, normal and shear critical damp ratios are 0.2, and timestep is fixed to . Details of parameters are listed in Table 1.
Figure 10(a) shows the shear stress and shear strain relations in the simple shear tests and DEM simulation, and Figure 10(b) shows the corresponding shear strain and volumetric strain relations. In Figure 10(a), it can be seen that the modelled shear stress and shear strain behavior is in a good agreement with that in the simple shear test, and the maximum shear stress ratios are 0.53 and 0.51, respectively. In Figure 10(b), it is shown that the modelled shear strain and volumetric relation is close to that in the simple shear test. In addition, the phase transformation is well captured in the simulation at a corresponding shear strain in the simple shear test. In simple shear test, the specimen transformed from contraction to dilation at 15% shear strain, while the transformation is observed at 20% shear strain in DEM simulation. In calibration, increasing Young’s modulus can effectively improve the agreement of that in DEM simulation and the simple shear test. However, the increased Young’s modulus led to a higher simulated shear stress compared with the simple shear test. Considering the neglected shape effect in DEM simulation, the DEM model gives a reasonable agreement in terms of phase transformation and strain changes, since it is difficult to reproduce experimental volumetric strain in a DEM simulation [9].
(a)
(b)
4. Results
Figure 11 shows the contact forces at a block with y = −0.002 m to y = 0.002 m (along shear direction) during different stages of a simulation; the thickness in the figure indicates the magnitude of contact force. Figure 11(a) shows the contact force network after 200 kPa consolidation stress is reached. It can be seen that large contact forces are uniformly distributed along vertical direction which is also the direction of major principal stress. Figure 11(b) shows the contact force network at the shear strain of 40%; it shows that large contact forces rotate towards the shearing direction. Figure 12 shows the qualitative distribution of average contact force after consolidation, at 5% shear strain, and at 40% shear strain. A bin shows that the magnitude of average contact force increases by 10° in the rose diagram. It is clearly indicated that, after consolidation, large contact forces are mainly near 90°, while they change to 45° at the shear strain of 40%.
(a)
(b)
(a)
(b)
Figure 13 shows the movement of particles during simulation. In Figure 13(a), particles after consolidation are colored in red, green, and blue which indicate different x positions. Figure 13(b) shows the updated location of particles at 40% shear strain. On the cutting plane, particles in the same color show linear increased displacement from top to bottom, and particles in different colors are generally parallel to each other. Since there is no obvious buckling of columns, it can be concluded that no shear band is developed during shear at 40% shear strain which is in consistence with test observations. Buckled column was reported by Gutierrez and Muftah [3] in a 2D simple shear DEM simulation, in which a higher rolling resistance coefficient of 0.5 was used. A higher rolling resistance coefficient increases the bending stiffness of a column and, at the same time, increases the instability related to buckling. A lower rolling resistance coefficient allows particles to move independently and is unable to form a column due to lack of rigidity; as a result, colored column will not buckle locally. In this study, the simulated Leighton Buzzard sand is a subrounded particle, so a small rolling resistance coefficient of 0.2 is used which leads to the uniform movement of particles.
(a)
(b)
In the DEM simulation, the complete stress state can be measured via measurement function, and it is transformed to principal stresses for further analysis, as shown in Figure 14. In addition, rotation of principal stress and principal strain increment is presented in Figure 15, where is the rotation angle of principal stress and is the rotation angle of principal axes of strain increment. The major principal stress increases dramatically in the first 15% shear strain, and then the value is almost stable during shear. The intermediate principal stress and minor principal stress have a gradual increasement during shear. The rotation angle of principal stress has a dramatic increasement (from 0° to 35°) in the first 15% shear strain, followed by gradual increasement, while the rotation angle of principal strain increment has a gradual increasement from 43° to 45° in the first 15% shear strain and kept constant at 45°. It can be concluded that the difference of the two rotation angles is the greatest at the beginning of shear (43°), and the difference is then decreased during shear. The difference indicates the noncoaxiality, meaning the degree of noncoaxiality is the greatest at the beginning of shear. It should be noted that, different from previous studies on noncoaxiality in which the degree of noncoaxiality is close to 0 at 10–20% shear strain [9, 16, 17], the degree of noncoaxiality in this study is still large at that range. This is due to the face that this study used samples (in test and simulation) with larger contractive behavior and larger shear strain at peak shear stress. It can be concluded that the degree of noncoaxiality should be carefully treated for different cases.
5. Conclusion
In this study, a 3D simple shear model using DEM is built based on the boundary condition of an NGI-type bidirectional simple shear apparatus. Stack of rings used as lateral constraints in a bidirectional simple shear test is modelled by layers of clumps which is possible to be moved by particles, and different contact types and parameters are used to model the sand-loading caps, sand-latex membrane, and sand-sand contacts. A simple shear test using the bidirectional simple shear apparatus is performed for the calibration of the 3D DEM simple shear model. By analyzing the simulation results, the following can be concluded:(1)Rings generated by clumps can provide an accurate boundary condition and are effective in computation since no contact force is needed for a clump.(2)In the simulation, the orientation of average contact force changed dramatically during shear. It is in the vertical direction (90°) before shear and changes to 45° at 40% shear strain. No shear band is observed which is in consistence with the test, and particles move uniformly.(3)In the simulation, the degree of noncoaxiality is the greatest at the beginning of shear, and it is decreased during shear. However, the degree of noncoaxiality is still large at 20% shear strain where there is a 10° difference between the rotation angle of principal stress and principal strain increment.
Abbreviations
| : | Vertical stress |
| : | Vertical stress after consolidation |
| : | Normalized vertical stress |
| : | Shear stress |
| : | Normalized shear stress |
| (%): | Shear strain |
| (%): | Volumetric strain |
| : | Constant normal stiffness |
| : | Constant shear stiffness |
| : | Friction coefficient |
| : | Surface gap |
| : | Normal critical-damping ratio |
| : | Shear critical-damping ratio. |
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (NSFC Contract no. 51708040), Natural Science Foundation of Shaanxi Province (Project Code: 2019JQ-683), Fundamental Research Funds for the Central Universities of Ministry of Education of China (Grant no. 300102210216), and Zhejiang Natural Science Foundation (Project Code: LQ19E090003).