Abstract
The reliability of discrete element method (DEM) numerical simulations is significantly dependent on the particle-scale parameters and boundary conditions. To verify the DEM models, two series of biaxial compression tests on ellipse-shaped steel rods are used. The comparisons on the stress-strain relationship, strength, and deformation pattern of experiments and simulations indicate that the DEM models are able to capture the key macro- and micromechanical behavior of inherently anisotropic granular materials with high fidelity. By using the validated DEM models, the boundary effects on the macrodeformation, strain localization, and nonuniformity of stress distribution inside the specimens are investigated using two rigid boundaries and one flexible boundary. The results demonstrate that the boundary condition plays a significant role on the stress-strain relationship and strength of granular materials with inherent fabric anisotropy if the stresses are calculated by the force applied on the wall. However, the responses of the particle assembly measured inside the specimens are almost the same with little influence from the boundary conditions. The peak friction angle obtained from the compression tests with flexible boundary represents the real friction angle of particle assembly. Due to the weak lateral constraints, the degree of stress nonuniformity under flexible boundary is higher than that under rigid boundary.
1. Introduction
The natural granular materials such as sands and gravels universally have the characteristics of anisotropy due to deposition under gravity or compaction. A number of studies in the bearing capacity of shallow foundations [1–3] and slope stability [4, 5] demonstrated that the deformation and strength anisotropy of the granular materials played a significant role on the geotechnical engineering.
The mechanical behavior of granular materials with inherent fabric anisotropy has been investigated using almost all the available laboratory testing methods such as triaxial compression tests [6, 7], direct shear tests [8, 9], plane strain compression tests [6, 10, 11], and hollow cylinder torsion shear tests [10, 12]. All of these experimental results indicate that the deformation and strength of inherently anisotropic granular materials are significantly dependent on the direction of applied stresses with respect to the bedding plane. In order to correlate these macrodeformation behaviors to the evolution of fabric characteristics, various new testing technologies including microstructural observation of thin sections fixed by resin [13], X-ray CT [14, 15], and stereophotogrammetry [16] have been used. However, these methods are too expensive or even impossible to capture the particle-scale quantities during the whole process of deformation.
Instead of making efforts on the particle-scale fabric measurement of real 3D laboratory experiments, the biaxial compression tests were conducted using two-dimensional rod assemblages [19, 20]. In the tests conducted by Konishi et al. [19], the photoelastic rods with oval cross-section were used to investigate the inherent anisotropy and shear strength. Their test results indicated that the deformation behavior of these 2D rods resembled that of real granular materials to a great extent. However, compared to the 3D laboratory experiments, it is much easier to catch the evolution of the fabric characteristics during the deformation process of the specimen.
The discrete element method (DEM) is capable of providing the detailed information about particle movement, rotation, and interaction between particles. A large number of numerical simulations for the biaxial/triaxial compression tests [18, 21–25] and direct/simple shear tests [18, 22, 26–28] have demonstrated that DEM is a powerful tool to study the microdeformation mechanism of granular materials. However, these DEM models differ greatly in the simulation of particle shape and boundary conditions, which have great effects on the macro- and particle-scale responses of granular materials.
The present paper aims at simulating the biaxial compression tests of ellipse-shaped steel rod assembly with high fidelity. The DEM model is validated by comparing the macro- and particle-scale responses of laboratory experiments and numerical simulations for two series of biaxial compression tests. The effects of boundary conditions on the stress-strain relationship, strength, strain localization, and stress nonuniformity are investigated.
2. Validation of Discrete Element Models
2.1. Biaxial Compression Experiments
Two series of biaxial compression tests on ellipse-shaped steel rod assembly are used to validate the DEM models in this paper. The biaxial compression test equipment was developed by the second author [17]. Its structure diagram was shown in Figure 1. A rectangular sample container , 240 mm in height and 120 mm in width, was constituted by the top plate , bottom plate , and two side plates . The base was supported by the vertical loading platform of a conventional triaxial compression apparatus and the component labeled as was the reaction frame. During the shearing, the vertical deformation of the sample was controlled by the vertical movement of the loading platform , while the top plate was kept immovable. It should be pointed out that the base together with both side plates and bottom plate moved upward at the same speed as the movement of loading platform in this equipment, which was not common for the compression tests. The vertical pressure applied on the sample was measured by the force gauge . The force gauge and was used to measure the confined pressure applied on the left and right side platens, respectively. Each of the left and right side plates together with the base of the frame was installed two displacement sensors, and totally six displacement sensors, denoted by DT1 to DT6, were used.
The materials tested were the ellipse-shaped steel rods with a uniform aspect ratio (the ratio of the minor axis length [] to the major axis length []) of 1 : 2 and a length of 40 mm. The aggregate of the specimen was made by mixing three kinds of rods with their major axis length of 4 mm, 2 mm, and 1 mm. And their mass ratio was controlled to be 8 : 2 : 1.
To investigate the loading direction-dependent responses of the rod assembly, the specimens with various tilting angles, denoted by , were fabricated as Figure 2 shown. The tilting angle is defined as the angle between the bedding plane and the plane of the major principal stress. One black rectangular frame was used to contain the rod assembly, whose inside dimensions were 240 mm high, 120 mm wide, and 50 mm long. To fix the black rectangular frame at a prescribed tilting angle, one transparent organic glass with marked lines and holes was used. The specimen with the tilting angle of was fabricated as follows. Firstly, according to the required tilting angle , the horizontal black rectangular frame was rotated clockwise by the angle of and fixed on the organic glass using bolts. Then the mixed iron rods were placed into the frame layer by layer by hand while keeping the major axis of rods horizontal. When the frame was filled with iron rods, small shaking was applied for 1 minute to uniform the rod assembly. After that, the frame was removed from the organic glass and returned back to the horizontal direction by rotating counterclockwise by . Finally the rod assembly was pushed horizontally to the rectangular sample container of the biaxial compression equipment using an organic glass plate, which has the same inside width and height as the frame and the rectangular sample container . Till now the specimen with the tilting angle of was prepared and ready for biaxial compression tests. Two series of biaxial compression tests by changing tilting angles and confining pressures were conducted.
(a)
(b)
2.2. Discrete Element Model
The DEM simulation package PPDEM developed by Fu and Dafalias [18, 22] was used in this study. As described in the papers of Fu and Dafalias [18, 22], the PPDEM is capable of characterizing any noncircular particle shapes by using “polyarc” element. The initial fabric anisotropy of the specimen can be well represented by modeling the deposition process under gravity. In addition, local quantities such as local stress, strain, particle orientation, rotation, and void ratio can be measured conveniently by defining a polygon-shaped “mask”, whose vertex is attached to a particle.
To simulate the biaxial compression experiments of the rod assembly described above, three particle sizes with the major axis length of 4 mm, 2 mm, and 1 mm, respectively, were produced with their number ratio of 1 : 1 : 2, which was the same as the tested rod assembly. The biaxial compression specimens with various tilting angles were produced using the same method as described by Fu and Dafalias [18]. As Figure 3(a) shows, a “master pack” of 30000 particles was fabricated firstly by particle pluviation, whose bedding plane is horizontal. Then the “master pack” was rotated counterclockwise by the tilting angle of , and the biaxial compression specimens were “trimmed” horizontally out of the master pack. Around 10000 particles were included in the “trimmed” specimen with the initial size of 240 mm in height () and 120 mm in width (). The initial void ratios of all the specimens with different tilting angles were .
(a)
(b)
When the specimen was fabricated, four rigid walls were applied as the boundary of the specimen. The loading in the numerical simulations was controlled to be the same as that in laboratory experiments, which was shown in Figure 3(b). After the specimen was consolidated isotropically at the required confining pressure of , the shearing began. In the vertical direction, the bottom wall moved upward at a specified rate while the top wall was kept immovable. The two end walls were free to move in the horizontal direction. For the left and right lateral walls, the horizontal confining pressure of was maintained constant by the wall servocontrol. However, both lateral walls moved vertically at the same speed as the bottom wall.
As described by Fu et al. [29], the overlap-area contact law was adopted for the interparticle behavior in PPDEM. The research conducted by Mirghasemi et al. [30] has demonstrated that the format of contact laws has minor effects on the macroscopic behavior of particle assemblage as long as the model parameters are appropriately selected. Thus the contact model for the tested steel rods has not been measured, and the overlap-area contact law is used for the numerical simulations. By conducting parameter sensitivity analysis, it is found that two parameters, the interparticle friction angle and the friction angle between particle and wall, have significant effects on the macromechanical behavior of particle assembly. These two parameters are chosen to be 30° and 10°, respectively, by comparing the stress-strain relationship between experiments and simulations for two series of tests varied in tilting angles and confining pressures.
2.3. Verification of Discrete Element Model
Two series of biaxial compression tests are used to validate the discrete element models. One is the tests with the tilting angle of varying from 0° to 90° with interval of 15° while keeping the same confining pressure kPa; the other is conducted by changing the confining pressures at the same tilting angle . Figures 4(a)–4(d) compare the evolution of the stress ratio , volumetric strain , and deformation pattern of laboratory experiments and numerical simulations for the first series of tests with , 30°, 60°, and 90°, respectively. The principal stresses of and are calculated using the same method as the laboratory experiments. The is obtained by dividing the average vertical force of two end walls by the specimen width, and the is calculated by dividing the horizontal force of two lateral walls by the specimen height.
| (a) ; |
| (b) ; |
| (c) ; |
| (d) ; |
As Figures 4(a)–4(d) show, the key direction-related mechanical behavior of granular materials with inherent fabric anisotropy can be captured by numerical simulations with high fidelity, although the initial shear modulus of all simulations is a little bit higher than that of experiments. The response of the granular materials is significantly dependent on the loading direction for both simulations and experiments. For and , the principal stress ratio reaches a peak followed by strain softening. As the tilting angle increases, the strain softening is weakened. For and , the development of the principal stress ratio tends to be strain hardening, which progressively approaches plateaus and then remains constant. With the continuation of the deformation, the specimen contracts firstly and then dilates. The dilation is reduced with the increase of the tilting angle . It should be pointed out that these results are qualitatively similar to the plane strain test results obtained by Oda et al. [6] and Tatsuoka et al. [10]. In addition, the deformation pattern of specimens at typical states of A to G visually looks similar, which shows that both experiments and simulates should possess the same particle-scale deformation mechanism.
Figure 5 gives the comparison of the peak friction angle with respect to the tilting angle . The peak friction angle corresponds to the maximum principal stress ratio in Figures 4(a)–4(d), which is calculated through the curve of the principal stress ratio versus axial strain by assuming zero cohesion. It can be seen that the evolution of the peak friction angle with respect to the tilting angle follows the same tendency for both experiments and simulations. As the tilting angle increases, the peak friction angle decreases. The biggest difference of between simulations and experiments is about 2°, which happens at .
Figure 6 shows the comparison of experiments and simulations for the second series of tests, in which three different confining pressures are applied for the tilting angle . For , the stress-strain relationship shows the characteristics of strain softening under three different confining pressures. The effects of confining pressure can be modeled. With the increase of confining pressure, the peak strength reduces. The specimen contracts followed by dilation. The maximum volumetric contraction increases as the confining pressure increases.
3. Investigation of Boundary Effects
It should be pointed out that the two lateral platens move vertically at the same speed as the bottom platen in the above biaxial compression tests and simulations, which is not common for compression tests. In the following, this mode of boundary control is denoted by Rigid boundary A. To investigate the effects of boundary condition, two other boundaries are used as Figure 7 shows. One boundary, denoted by Rigid boundary B, is the same as Rigid boundary A except that the two lateral walls are free to move in the vertical direction. The other boundary, denoted by Flexible boundary, resembles the conventional triaxial compression tests. The top and bottom boundaries are simulated by the rigid walls. The two lateral boundaries are flexible like membrane. The confining pressure is directly applied on particles as described by Fu and Dafalias [18].
(a)
(b)
3.1. Boundary Effects on the Stress-Strain Relationship and Strength
Figures 8(a)–8(d) compare the development of the principal stress ratio and volumetric strain with axial strain under the above three different boundary conditions for the tilting angle , 30°, 60°, and 90°, respectively. All the principal stresses of and are calculated by dividing the force applied on the wall by the relevant specimen size, except that the is directly applied on particles under Flexible boundary. The averaged specimen width, dividing the specimen volume by the height, is used to calculate under Flexible boundary.
(a)
(b)
(c)
(d)
It can be seen that the stress-strain relationship and strength are dependent on the boundary condition. Much stronger strain softening happens under Flexible boundary compared to the other two rigid boundaries. For and 90°, the development of stress ratio shows strain hardening characteristics under Rigid boundary A and Rigid boundary B, while the marked drop of still occurs after the peak under Flexible boundary. The reason for the strong strain softening under Flexible boundary may be due to the lateral bulging of the specimen at large deformation. The maximum peak friction angles and the initial shear modulus are achieved under Rigid boundary A, and the corresponding minimum values are obtained under Flexible boundary for any tilting angle . The effects of boundary condition on the development of volumetric strain are not as significant as on the stress ratio . At the initial stage of volume contraction, the volumetric strain is almost the same for three different boundary conditions. As the deformation continues, some minor differences are observed.
The DEM package PPDEM used is capable of measuring the average local stresses in any domain inside the specimens by defining a mask as described by Fu and Dafalias [18]. To obtain “real” average stresses inside the specimens, the mask is defined firstly as Figure 9 shows. The mask defined under two rigid boundaries is a little bit smaller than the whole specimen and the same mask is used during the shearing process. However, under Flexible boundary, due to the distortion of the specimen, the particles bulged outside are not included in the mask, and the mask is changed at each 4% axial strain. The average stresses inside the mask are calculated.
The development of the “real” stress ratio measured by the above masks is shown in Figures 10(a)–10(d) for tilting angle , 30°, 60°, and 90°, respectively. It can be seen that the “real” stress-strain relationship measured inside the specimen is almost the same except that some minor difference happens for . The minor difference of the among different boundary conditions may be due to the different masks used. However, reviewing the stress-strain relationship presented in Figure 8, in which the stresses are calculated by the force applied on the wall, the stress ratio is affected significantly by the boundary conditions.
(a)
(b)
(c)
(d)
To investigate the boundary effects on the strength, Figure 11 gives the peak friction angle calculated on the basis of Figures 8 and 10, which are indicated with “by wall” and “by mask”, respectively, for three different boundary conditions. It can be found that the peak friction angle decreases as the tilting angle increases for all cases. The peak friction angle calculated “by wall” under Rigid boundary A is the maximum, which is almost 1.5° higher than that under Rigid boundary B and 3.5° higher than that under Flexible boundary. However, when the “real” stresses inside the specimen are calculated by “mask”, the difference of the peak friction angle among three different boundary conditions is less than 1° as far as the same tilting angle is considered. In addition, one interesting phenomenon found in Figure 11 is that, under Flexible boundary, the peak friction angle obtained by “wall” is very close to the real value obtained by mask. The difference between them is less than 1°, which indicates that the peak friction angle obtained usually by triaxial compression tests can represent the true strength of the granular materials.
3.2. Boundary Effects on the Strain Localization
The strain localization happens in almost all kinds of granular materials during the shearing process. The shear localization pattern is very complex and may be affected by many parameters, among which the effect of the boundary condition is significant [10, 16]. To study the effects of boundary conditions on the strain localization, Figures 12 and 13 show the grid deformation and the particle rotation contours at the axial strain of 15% under Rigid boundary B and Flexible boundary for the tilting angle , 30°, 60°, and 90°, respectively. The grid was originally “painted” on the consolidated specimen. Considering that the specimens have almost the same deformation pattern under two rigid boundaries, only the results under Rigid boundary B are presented in Figure 12.
(a)
(b)
(a)
(b)
It can be seen that the grid deformation and the particle rotation contour give the same shear localization region. The primary shear plane is dependent on the boundary condition. Under Rigid boundary B, the primary shear plane extends from corner to corner due to the strong constraint of the four rigid walls. However, under Flexible boundary, the constraint of the boundary is much weaker, and large bulging of the particles is found. At any tilting angle , the primary shear plane extends upward from left to right, which was denoted as Type-b failure plane in Tatsuoka et al. [10]. And the shear localization mainly focuses in the middle part of the specimen. However, under Rigid boundary B, the primary shear plane produced is significantly dependent on the direction of loading. Different types of shear planes are found when the tilting angle changes from 0° to 90° as Figure 12 shows. For and 60°, the shear plane extending from the left-up corner to the right-down corner is dominant, which was denoted as Type-a failure plane in Tatsuoka et al. [10]. For , the primary shear plane is Type-b mode. The X-type shear plane happens for .
3.3. Boundary Effects on the Stress Nonuniformity Inside the Specimen
To investigate the boundary effects on the stress nonuniformity inside the specimen, three different masks, denoted by mask-up, mask-mid, and mask-bot, respectively, are defined in the upper middle, central middle, and bottom middle part of the specimen as shown in Figure 14. Each mask contains over 1000 particles. The averaged stresses in these masks are calculated. Given typical examples of and 90°, Figures 15 and 16 show the evolution of stress ratio measured in different masks with axial strain under the conditions of Rigid boundary B and Flexible boundary. The curves of Figure 10 under the same conditions are plotted for comparison, which are indicated with “mask-whole.”
(a)
(b)
(a)
(b)
It can be found that the boundary conditions affect the stress distribution inside the specimen. The stresses in the central middle part of the specimen are higher than those in the upper and bottom parts. And the stresses denoted by “mask-whole” can represent the average stresses of the upper middle, central middle, and bottom middle masks. In addition, the degree of stress nonuniformity under Flexible boundary is higher than that under Rigid boundary B. The reason for the high stress nonuniformity under Flexible boundary can be explained as follows. Under Flexible boundary, the lateral constraints are weak. More particles in the middle part of the specimen extrude and they cannot transfer the vertical stresses efficiently. Thus high forces concentrate in the central middle particles.
4. Conclusion
The DEM numerical simulations are a very promising tool to investigate the macro- and micromechanical behavior of granular materials. However, its reliability is significantly dependent on the particle-scale parameters and boundary conditions. In this paper, two series of biaxial compression tests varied in bedding plane inclination angles, and confining pressures are conducted on the ellipse-shaped steel rod assembly. The DEM models are validated by comparing the stress-strain relationship, strength, and deformation pattern of experiments and simulations. On this basis, three different boundary conditions are applied to investigate the boundary effects on the macrodeformation, strain localization, and the nonuniformity of stress distribution inside the specimen. The main conclusions can be summarized as follows.(1)The stress-strain relationship and strength measured by the force on the wall are significantly dependent on the boundary conditions. The peak friction angle obtained under rigid boundary is higher than that under flexible boundary.(2)The boundary condition has minor effects on the mechanical behavior of particle assembly inside the specimen. The peak friction angle obtained under flexible boundary is the closest to the real friction angle of granular materials. (3)The strain localization pattern of granular materials with inherent fabric anisotropy is dictated by the boundary condition and the bedding plane inclination angle. Under the rigid boundary conditions, various types of shear planes are observed for different loading directions.(4)The boundary condition affects the stress distribution inside the specimen. The degree of stress nonuniformity under flexible boundary is higher than that under rigid boundary due to the weak lateral constrains under flexible boundary.
Acknowledgments
The present work is financially supported by the Beijing Natural Science Foundation (Contract no. 8133053) and the National Natural Science Foundation of China (Contract nos. 51079075 and 51009002).