- Email: [email protected]
Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model
Computers and Geotechnics 67 (2015) 83–93
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model Antonio Salazar, Esteban Sáez ⇑, Gislaine Pardo Department of Structural and Geotechnical Engineering, Pontificia Universidad Católica de Chile, Vicuña Mackenna 4860, Macul, Santiago, Chile
a r t i c l e
i n f o
Article history: Received 2 October 2014 Received in revised form 20 January 2015 Accepted 23 February 2015
Keywords: DEM Direct shear test Micro-parameters Rolling friction model Real grain size distribution
a b s t r a c t In this research, a direct shear test was modeled using the 3D Discrete Element Method (DEM). The results are compared against experimental data. The real sand was modeled as about 70,000 spheres, with sizes consistent with the real grain size distribution, using a rolling friction model to include the sand’s grain shape. It is shown that, the stress path can be appropriately reproduced, but the dilatancy were very difficult to replicate. Finally, specimens with up-scaled material were tested. Although the unscaled material was the most representative, the response of the particles that were up-scaled twice was slightly less accurate, but with a considerable decrease of computation time. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction In this paper, a 3D Discrete Element Method (DEM) model of the direct shear test was developed to reproduce the response of coarse sand in order to explore the ability of this modeling approach to reproduce this widely used test considering a realistic grading distribution. Several features of this method were explored, such as the ability to monitor and analyze displacements and contact forces of all particles composing the material. This laboratory test is widely used in geotechnical engineering, because of its simplicity and easy analysis of the results. The test uses a box with two halves, where the lower half is moved horizontally while the upper half remains fixed. A confinement force is applied to the top wall, which is free to move vertically. The registered reaction force depends on the properties of the shear band resulting from the induced relative displacement between both rigid half-boxes (Fig. 1). Several studies have shown that it is possible to qualitatively reproduce laboratory results with the discrete element technique using 2D models, despite the planar formulation limitations [1,2]. Cui and O’sullivan [3] were probably the first to perform a comparison of physical test data and direct shear test DEM models using perfect metallic spheres. Later, Yan and Ji [4] compare their DEM results with real soil direct shear tests performed on irregular ⇑ Corresponding author. E-mail addresses: [email protected] (A. Salazar), [email protected] (E. Sáez), [email protected] (G. Pardo). http://dx.doi.org/10.1016/j.compgeo.2015.02.017 0266-352X/Ó 2015 Elsevier Ltd. All rights reserved.
limestone rubbles. In their approach, they used clumps to simulate the shape of the real rubbles. In both studies, a good match between numerical results and real test data was achieved. Most recently, Härtl and Ooi [5] used the Jenike direct shear test to investigate the particle shape and microscopic friction on the bulk friction. They used single sphere particles and a clump of two overlapping spheres to investigate the shape effect. Their results show that the particle interlocking has a greater effect than the packing density on the bulk friction. In the present research a 3D approach was chosen, because the 2D models fail to accurately predict the peak and residual friction angles under shear flows [6], due to out-of-plane migration of the particles, among other effects. In recent years, the microprocessor industry has made little improvement in the core’s clock speed, giving advantage to multicore architectures [7,8] and the number of processors per machine will continue to increase. This change encourages the use of highperformance parallel software to fully take advantage of this technology. The DEM is especially suitable for this type of code because its elements have independent degrees of freedom, which makes it simple to parallelize. Harnessing this benefit, the coarse sand was initially modeled without size scaling. Nevertheless, to reduce the computational cost and to allow the use of a large number of particles, simple individual spheres were used with a rolling friction model representing the particle shape effect [9]. The effect of the model’s micro-mechanical parameters were studied with a sensitivity analysis. Two specimens with different initial void ratios considered, in order to highlight the effect of the system’s compactness on the response. We present the
84
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Sand grains photos Fig. 1. Direct shear test.
response of the calibrated model as compared to the laboratory results, which were obtained from previous stages of this research [10,11]. Also, a brief analysis of the chain force and evidence of outof-plane particles migration is shown. Despite advances in computing technology and numerical algorithms, computational cost is still the major limitation of DEM. Due to this issue, a scaled material was investigated and also compared with the experimental data. Due to laboratory restrictions, several down-scaling methodologies exist to test representative soil samples using small grain sizes. Among them, the parallel graduation method [12] was found suitable and able to capture the original soil response [13–15] with a finer gradation curve. In this paper, we use the parallel graduation method, but in the other direction (up-scaling), to develop specimens with particles twice and four times bigger than the original material. The scaled specimens were tested with the same calibrated parameters as the unscaled specimen.
(b) Grain size distribution of the real and simulated material Fig. 2. Grain samples and grain size distribution.
2. Real sand and DEM model description 2.1. Tested material
2.2. Brief description of DEM
The selected material was an angular non-cohesive coarse sand. Its main properties are listed in Table 1. A series of direct shear tests were performed on this sand under three different confinement pressures: 160 kPa, 80 kPa and 40 kPa. The laboratory test was instrumented with force and displacement transducers. Fig. 2a shows a series of photos, taken with a microscope, of the sand’s grains and processed with Balu toolbox for Matlab [16] to extract their apparent boundaries. Under a visual analysis, particles of all sizes have a sphericity of about 0.8 and a roundness of about 0.5, hence a subrounded shape [17]. Fig. 2b shows the sand grain size distribution. According to the Unified Soil Classification System (USCS), it corresponds to a poorly graded sand (SP) and all the material is retained in four sieves. The sand has one percent of fines (diameter less than 0:075 mm) and most of the mass is retained in sieve N 16, which are the biggest grains (between 2.16 mm and 1.18 mm of diameter).
The DEM is a numerical technique used to simulate the behavior of granular aggregates, jointed rock masses and chemical powders, among others. The main feature of this methodology is that the material is simulated as an assemblage of individual particles that interact with each other. This approach allows to capture the relative movements and rotations of the particles without the need of a sophisticated constitutive model. Also, it can achieve large displacements because each particle has its own independent degrees of freedom. For granular applications, the dynamic equilibrium condition uses the Eq. (1) (Newton’s second law) and Eq. (2), where M i and Ii are the mass and the moment of inertia of each particle i; Fji and Mji are the contact forces applied by particle j on i and the moment induced by this force; Fe and Me are the resultant of the body-force field over the i’th particle (only gravity in this case) and the torque because of this force.
€i ¼ Mi x
Table 1 Sand properties.
X ji F þ Fe
ð1Þ
j
Property
Value 3
Density: c [kg/m ] Void ratio: e0 Porosity: n0 Angle of internal friction: /0 Cohesion: c0 [kPa]
1578 0.679 0.405 35 8
€i ¼ Ii H
X ji M þ Me
ð2Þ
j
The time is divided into intervals, and the displacements are calculated in each step by the double integration of Eqs. (1) and (2). This time step must be small enough to avoid numerical instabilities. In
85
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
this research, the maximum time step was selected as the minimum of the Hertz or the Rayleigh’s time step [18]. The DEM software used to solve the problem was LIGGGHTS [19] which is a 3D Open Source Code with several granular features already implemented and ready to run in parallel. Further information can be found in the developer’s website.
Fij
ei Centroid j
2.3. Inter- granular contact modeling Eq. (3) shows contact model adopted in this investigation, which is based in Hertz’s work [20]. In this equation, E is the equivalent Young modulus (Eq. (4)), R is the equivalent particle radius (Eq. (5)) and dn is the overlap of the contact. mi and Ri in Eqs. (4) and (5) are the Poisson ratio and the radius of the particle i which is in contact with particle j.
pffiffiffi 3 4 Fn ¼ E R d2n 3
ð3Þ
1 m2j 1 m2i 1 ¼ þ Ej Ei E 1 R
¼
ð4Þ
1 1 þ Ri Rj
ð5Þ
The magnitude of the tangential force is given by the Mindlin’s theory [21] and limited by the Coulomb friction criteria shown in Eq. (6), where l is the friction coefficient. If this limit is exceeded, relative sliding takes place and the tangential force remains constant.
F tmax 6 ljF n j
ð6Þ
Since the particles being considered are perfect spheres, the Ai’s et al. [22] Elastic–plastic spring-dashpot rolling friction model was considered to take into account the shape effect between the sand’s grains. This model adds a torque (Mr ) in the dynamical equilibrium of the particles, which includes a spring torque (Mkr ) controlled by the rolling viscous damping coefficient (g) and a viscous torque (M dr ) controlled by the user by the rolling friction coefficient (lr ). The model implemented is shown in Eqs. (7)–(10).
M r ¼ M kr þ M dr
ð7Þ
The spring torque (Mkr ) is the product of the rotation of the particle (Dhr ) and the rolling stiffness (kr ). This torque is limited by Eq. (9), which is controlled by lr , similar to the Coulomb’s friction model. The rolling spring stiffness is defined in Eq. (10), where hm r is the full movilization rolling angle.
ej
Centroid i
−Fij Fig. 3. Real contact.
2.4. DEM material The granular material is simulated as an assembly of spheres with a radius according to the grain size distribution of the tested sand. Fig. 2b shows the comparison between the grain size of the real and simulated material. The model does not include the onepercent of fine sand because it would increase the computational cost considerably and it does not have any significant influence on the material’s macroscopic behavior. The percentage of sand mass retained in each of the sieves of the ASTM gradation test was simulated with particles whose diameters were linearly distributed between the opening size of the sieve where they mass was retained and the opening size of the preceding sieve, as can be seen in Fig. 4a. The two exception is that the finest part was not modeled and the diameter of the coarsest particles was limited to 1.6 mm. Fig. 4b shows that the amount of the modeled mass accumulated per sieve is slightly greater than the real material, because the mass of the finest portion that is not included is distributed among the other sieves. 2.5. Boundary conditions and fabric generation All the walls, except the top one, were modeled as rigid boundaries. This rigid boundary condition does not have inertia and its velocity is imposed to control its position during the simulation. Walls interact with the particles, following the Hertz’s contact model, to calculate the repulsive forces. The top boundary is modeled as a servo-controlled wall, in which its speed is related to the current stress by the algorithm shown in Eq. (12). The process is controlled by the constant kp ; ki and kd , which multiply the error between the current force in the wall and the target value (Det ), P the cumulative error ( t Det ) and the error derivative, respectively.
X Det DetDt Det þ kd Dt t
M kr;tþDt ¼ M kr;t kr Dhr
ð8Þ
V wall ¼ kp Det þ ki
jMkr;tþDt j 6 lr RF n
ð9Þ
The specimen is composed of approximately 70,000 particles. All the particles are randomly generated inside the box and are divided in three layers of 0.8 cm each. After the first layer is created, the particles fall freely until the static equilibrium is reached due to the gravity; then the process is repeated with the other two layers. After the settlement is completed, the servo-controlled wall is activated to confine the soil. Finally, a finite velocity is imposed on the lower half in order to reach a horizontal target displacement of 8 mm. The model, shown at the beginning of the shear process, is shown in Fig. 5. In order to get a loose or a dense fabric (eoi ), different coefficients of friction (li ) were used in the settlement process and during the shearing modeling. If lower values of li are used during generation, a dense fabric is obtained. The chosen Poisson’s ratio coefficient is a typical value used to simulate sand, while E was chosen to have a realistic contact stiffness. During shearing, the
kr ¼
lr RF n
ð10Þ
hm r
The viscous torque (M dr ) is controlled by g (Eq. (11)) and is only active while the rolling torque is not fully mobilized, assisting particle stabilization more than energy dissipation [22].
( Mdr;tþDt
¼
pffiffiffiffiffiffi 2g Ikr h_r
if jMkr;tþDt j < lr RF n
0
if jMkr;tþDt j ¼ lr RF n
ð11Þ
The reason for using this kind of model is because the contact force between two particles does not always pass through the centroid of the particles (an implicit assumption when spheres are used), as depicted in Fig. 3.
ð12Þ
86
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93 Table 2 Tested sand properties. Property
Specimen A (dense) 3
Specimen B (loose)
cs [kg/m ]
2650
2650
m li gi lri
5:0 108 0:256 0 0 0 0:489
5:0 108 0:256 0:3 0 0 0:594
E [kg/m2]
e0i obtained
in Table 2. The difference between these two specimens is the initial coefficient of friction li and consequently, the obtained void ratio e0i . The unit solid weight cs , used for both specimens, was obtained from laboratory measurements. The Young Modulus E and the Poisson’s ratio m are required to calculate the contacts’ stiffness and were selected among values proposed in the literature [1–4,6], but were small enough to avoid a time step that was too small. Fig. 6a–c shows the results using different combinations of the coefficient of friction l, the coefficient of rolling viscous damping g and the coefficient of rolling friction lr during the shearing phase. The shear stress, used in these comparisons, was calculated according to Eq. (13), where T is the horizontal reaction force of the upper half of the box (Fig. 1) and A is the effective shear area which depends directly on the horizontal displacement of the lower half (Dh ).
(a) Particle number histogram
rs ¼
(b) Mass histogram Fig. 4. Mass and particle number histograms.
z y
x
Fig. 5. Developed DEM model.
value of l; lr and g were modified to highlight their effect on the macroscopic response (Table 2). No rolling friction is considered during specimen generation.
3. Sensitivity analysis Two specimens were modeled and tested in order to explore the effects of the micro-mechanical parameters on the macroscopic response. The properties used for the model construction are listed
T AðDh Þ
ð13Þ
According to Fig. 6a, the shear stress comparison shows that the friction coefficient l has an influence for small displacements, and it tends to increase the peak value by approximately 30% when varied from 0.3 to 0.6. In the studied cases, the large-strain stress value is approximately the same for horizontal displacements larger than 4–5 mm, indicating that by varying only the Coulomb’s friction coefficient it is not possible to control the response of the model at large strains. Regarding volume evolution, the friction coefficient l increases the specimen dilatation without considerably affecting the path’s shape (Fig. 6a). Indeed, due to increase of particle friction, a looser configuration could be reached at large strains. Fig. 6b shows the sensibility on the rolling viscous damping g. In this case, a direct effect of this parameter on the stress path or on the displacement responses is not as clear, but a slight increment in the peak value is detected. The little effect of this parameter indicates that during most of the shear process the viscous torque is fully mobilized (Eq. (11)) and that the particles are rolling-sliding. Fig. 6c shows that the most important effect of the coefficient of rolling friction lr is the increase of the shear stress reached at large displacements. This parameter also increases the peak value, but in a less pronounced manner than l. The displacement comparison in Fig. 6c shows that lr increases the large-strain volume of the specimens. This influence cannot be noticed at a low-strain, where all the curves seem to follow the same path, highly influenced by the initial material fabric. As lr takes into account the particle’s shape, large values of this parameter can be associated with grains that are less round. In this case, it is reasonable to obtain a supplementary geometric interlocking at large strain, associated with a large apparent macroscopic angle of friction. Due to the shape effect, a looser critical state is also expected. Nevertheless, in general terms, the laboratory displacement path was very difficult to replicate. All three figures show that the main difference between the specimen A and B is that A gives a larger peak at low strain (a
87
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) μ influence (μr =0.2 & η=0.2).
(c) μr influence (μ=0.4 & η=0.5).
(b) η influence (μ=0.4 & μr =0.5).
Fig. 6. Sensitivity analysis of the DEM parameters.
typical response of a relatively dense sand), while the response for specimen B correspond to the usual stress–strain path of loose sand. Also, the shear stress comparison shows that specimen A has a higher initial slope (low-strain stiffness), in agreement with the experimental response. This slope can be controlled by the stiffness of the contacts and the initial void ratio. Specimen B has a lower large-strain vertical displacement in comparison to specimen A, which is also a typical experimental result for loose sands. The behavior difference was expected because of the different initial void ratio obtained for both simulated specimens. Compared to the experimental maximum (emax = 0.74) and minimum (emin = 0.44) void ratios, specimen A has a relative density of 84% classifying it as a very-dense material while specimen B is a medium-dense material with a relative density of 49%.
Table 3 Calibrated material parameters. Property
Specimen C
E [kg/m2]
lri l g lr
5:0 109 0:256 0:3 0 0 0:3 0:1 0:4
Obtained: c [kg/m3] Obtained: e0i
1612 0:643
m li etai
88
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
4. Calibrated model against experimental data A third specimen was created with the properties shown in Table 3, a higher contact stiffness was used in order to obtain a higher low-strain stiffness, while using a looser fabric. A greater void ratio was obtained, even when the same li was used, because higher contact stiffness avoids larger overlaps. The void ratio
obtained is slightly lower than the experimental value, which is expected because of the different shape of the real grains and the simulated particles [23]. With a confinement of 160 kPa, the model can appropriately reproduce the shear stress path during the whole test (stress calibration in Fig. 7). For the 80 kPa and 40 kPa confinements, the model does not reproduce the peak value, but the residual stress
Fig. 7. Calibrated material response.
(a) Initial force chain
(b) Final force chain Fig. 8. Initial and final chain force.
89
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Polar contact force distribution at the beginning of the shear test
(a) Initial velocity field
(b) Polar contact force distribution at the end of the shear test
(b) Final velocity field
Fig. 9. Polar distribution of the contact forces at the beginning and the end of the shear test on the shear plane.
is accurately reproduced for all the confinements. This might be due to a higher density of the sand in the laboratory tests compared to the DEM models. The dilatancy of the simulated specimen increases as the confinement stress decreases, in agreement with laboratory results (displacement calibration in Fig. 7). Additionally, the shapes of the curves are very similar. Nevertheless, the model fails to accurately predict the final sample’s volume. Indeed, this large strain volume is highly controlled by the particles’ shape, one of the major limitations of the method due to the use of spheres. Thus, the rolling friction model used to approximately include the particles’ shape effect, is not able to fully reproduce the large strain state of the reference soil. In the case of 160 kPa test, the final laboratory void ratio is 0.72; while the simulation gives a value of 0.703. The initial difference was reduced by the increase of the model’s dilatancy, but difference still remains, due to the shape of the sand grains. Fig. 8 shows the initial and final network of contact forces. The forces are evenly distributed at the beginning of the test (Fig. 8a), which is expected due to the initial vertical compression. During the test, a diagonal soil strut was formed, as can be seen in Fig. 8b. This might be due to the high geometric interlocking of
Fig. 10. In-plane velocity field (imposed displacement plane).
the material used, which prevents the particles of the central part of the sample to relatively slide. The aforementioned is evidenced when orientations and amplitude of the contact forces are analyzed (Fig. 9). Fig. 9a shows the distribution of the contact force at the beginning of the shear test, after the settlement process, while Fig. 9b shows the contact force distribution at the end of the test. At the beginning, the vertical contact forces are almost the double the horizontal contact forces, because of the vertical confinement imposed to the specimen. At the end of the test, the predominant contact forces are no longer vertical, but almost have the same orientation as the strut (approximately 19:5 ). The magnitude of the force in this orientation (between 0 and 36 ) is about four times larger than the initial vertical forces, while the magnitude of forces in the normal direction to the strut (between 90 and 126 ) is about one third of the initial vertical magnitude, highlighting the stress-state rotation during the shear test. In Figs. 10 and 11 the test was divided into small regions of 3 1.2 mm and the particles’ velocities were averaged in each of these regions in order to reduce the number of arrows and to represent the velocity as a continuous field. When the mean radius of the grain size distribution used is 0.5 mm, each of these zones represent the average value of about 360 particles. Fig. 10 shows a lateral view of the initial and final velocity field of the particles over
90
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Initial velocity field Fig. 12. Scaled material real grain size distribution.
(b) Final velocity field Fig. 11. Out-of-plane velocity field.
the imposed displacement plane. At the beginning of the test, all particles have similar velocities, which means that the specimen behaves as a continuous material (Fig. 10a). When large strains are reached, the bottom particles have a higher speed than top ones, hence a sliding surface dividing the specimen in two halves appears (Fig. 10b). A shear band of about 2.5 times the mean diameter of the particles is depicted in this figure. The particles do not follow a clear path in the out-of-plane direction (normal to the imposed displacement plane), also the magnitude of the average speed in this plane is about 30 times less than the in-plane speed. Nevertheless, the out-of-plane component of the speed indicates that there was some particle migration across this direction, justifying the choice of a 3D model (Fig. 11b). At the beginning of the test, the particles have moved almost completely in the vertical direction, with a very small out-of-plane component (Fig. 11a). It is interesting to note that out-of-plane migration particle movement is not symmetric across the observed section, probably due to the random generation of particle size and location among the sample. This asymmetry increases for large strains.
the simulations [13]. The purpose of using the up-scaled material is to investigate the decrease in the computation time and the accuracy degradation with the change in the ratio between the specimen size and the maximum particle diameter (Dspecimen =Dmax ). Two specimens with twice and four times bigger particles were studied. The grain size distribution is shown in Fig. 12. The Dspecimen =Dmax ratio of the twice scaled material is 7.5, while the ratio of the four times scaled material is 3.75. Although the Dspecimen =Dmax ratio of the four scaled material is out of the recommended range (Dspecimen =Dmax [24]), this scaled factor was tested for completeness. The time-computing comparison shown in Table 4 was done using only one core, to avoid the performance decrease when using a large number of processors with low data. Also, due to the large amount of time it would take to run the whole test in a single core, and because the purpose of this comparison is to illustrate computing time reduction, each of these simulations consists in a tenth of the required steps to complete the shear and the settlement process. All the simulations were carried out using the same timestep, hence the three specimens include the same number of steps and the same total simulated time. When scaling two-times, the computation time is reduce by 94% (20 times), while the four times scaled is reduce by 99.5% (200 times). This is because the reduction in the number of particles (86.8% and 98.7%) leads to an important decrease in number of contacts (96.4% and 99.8%), which defines the size of the DEM model. It is clear that even with a small coarse-graining scaling, the computation time decreases dramatically, allowing the development of large models with a reasonable running time. Because large fluctuations appear when larger particles are considered, both coarse-grained specimens were tested five times with the same set of parameters shown in Table 3 and an average was calculated (Figs. 13 and 14). The fluctuation in the results of
Table 4 Computational cost comparison. Unscaled
5. Parallel gradation In order to cut down the running time, the parallel gradation technique has been explored to reduce the number of particles in
Time (min) Number of particles Number of neighbors
Two-times scaled
Four-times scaled
Value
[%]
Value
[%]
Value
[%]
1903 64,944 2,754,758
100 100 100
107 8610 100,596
5.6 13.2 3.6
9.3 845 4033
0.5 1.3 0.15
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
91
(a) Twice scaled average - 160 [kPa]
(b) Twice scaled average - 80 [kPa]
(c) Twice scaled average - 40 [kPa]
Fig. 13. Shear test results with two-times scaled material.
the scaled materials are explained because when using the scaled sample, the number of contacts is significantly reduced, hence any change affects proportionally more in the macroscopic response than the unscaled case where a large number of contacts exist (Fig. 14). The paths of the scaled material, which was scaled two times, have a low deviation from the average (Fig. 13), while it increases greatly for the specimen that has been scaled four times (Fig. 14). This high fluctuation in the stress response is not as clear in the displacement results, because the servo-controlled boundary implementation responds more to a trend than to a single value. As the algorithm relates the stress with the velocities, and not directly to the displacements, the change in the vertical stress must be held for some time steps in order to see changes in the displacements.
When averages are compared with laboratory and unscaled material results (Figs. 13 and 14), very similar stress paths are obtained for both scales; however, the little increase in the stress path of the material which was scaled four times evidences that the twice scaled is lightly more accurate. The displacement comparison shows that the material which was scaled four times has a higher displacement than the twice scaled material, which is explained by the increase of the particle size. The material that was scaled four times is an inconvenient choice to reproduce the shear test, because of the high fluctuation of the results, despite the huge decrease of the computational cost. However, the satisfactory enough simulation of the stress paths and the huge decrease in running time make this strategy a suitable choice for problems where the ratio Dspecimen =Dmax is big enough to avoid boundaries effects.
92
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Four times scaled average - 160 [kPa]
(b) Four times scaled average - 80[kPa]
(c) Four times scaled average - 40[kPa]
Fig. 14. Shear test results with four times bigger material.
6. Conclusions When a discontinuous methodology such as DEM is used, the initial fabric density is very difficult to control compared to other continuum mechanics based methodologies. Specifically, in the implemented direct shear model, it was very difficult to find the correct combination of friction parameters and sphere generation to satisfactorily reproduce the initial laboratory relative density. The rolling friction model was used to represent the grain’s shape effect. It is concluded that the grain’s shape is essential to properly represent the material response under large strains. In general terms, the shape of the curves obtained with the DEM model approximately matches the experimental data. Indeed, the stress response of the model is very similar to the laboratory results, but the displacement paths are similar only qualitatively, especially at high horizontal displacements.
The particle force network clearly illustrates a shear zone and the formation of a strut. Additionally, the 3D model that was developed allows to point out an out-of-plane particle’s migration. Hence, the use of direct bi-dimensional discontinuous models under macroscopic plane-strain condition models must be done carefully. Because of the accuracy of the results compared to experimental data, the unscaled material seems to be the best model to reproduce the shear test. However, when the computational cost is relevant, a parallel grading curve might be considered. The results show a 94% (20 times) reduction in the running time with the model scaled twice, and just a little increment in the stress path and the displacement response. With the material that was scaled four times, the time reduction attained was 99.5% (200 times), but with a large fluctuation in the results, which makes this model unreliable. This leads us to conclude that the maximum up-scale
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
value must be studied for each problem probably in terms of the ratio Dspecimen =Dmax . The small loss in accuracy is totally compensated by the huge time reduction when using coarse-grained materials, based on the parallel grading technique. Hence, this methodology, theoretically, allows us to make bigger and more complex models with reasonable computing time. Nevertheless, the size increase technique raises the fluctuation of results and must be repeated several times to identify the tendency. Acknowledgement This work was possible due to a grant from the Fondo Nacional de Desarrollo Científico y Tecnólgico FONDECYT, Grant No. 11100157. References [1] Thornton C, Zhang L. Numerical simulations of the direct shear test. Chem Eng Technol 2003;26(2):153–6. [2] Masson S, Martinez J. Micromechanical analysis of the shear behavior of a granular material. J Eng Mech 2001;127(10):1007–16. [3] Cui L, O’sullivan C. Exploring the macro-and micro-scale response of an idealised granular material in the direct shear apparatus. Geotechnique 2006;56(7):455–68. [4] Yan Y, Ji S. Discrete element modeling of direct shear tests for a granular material. Int J Numer Anal Methods Geomech 2010;34(9):978–90. [5] Härtl J, Ooi JY. Numerical investigation of particle shape and particle friction on limiting bulk friction in direct shear tests and comparison with experiments. Powder Technol 2011;212(1):231–9. [6] Fleischmann PM, Drugan JW. Quantitative comparison of two-dimensional and three-dimensional discrete-element simulations of nominally twodimensional shear flow. Int J Geomech 2013;13(3):205–12. [7] Geer D. Chip makers turn to multicore processors. Computer 2005;38(5):11–3. [8] Borkar S, Chien AA. The future of microprocessors. Commun ACM 2011;54(5):67–77.
93
[9] Wensrich C, Katterfeld A. Rolling friction as a technique for modelling particle shape in DEM. Powder Technol 2012;217:409–17. [10] Pardo G, Sáez E. Experimental and numerical study of the arching effect in coarse sand. Comput Geotech 2014(57):75–84. [11] Pardo G. Estudio experimental y numérico del efecto arco. Master’s thesis, Pontificia Universidad Católica de Chile, Santiago, Chile; 2013. [12] Lowe J. Shear strength of coarse embankment dam materials. In: Proc, 8th int congress on large dams, vol. 3; 1964. p. 745–61. [13] Verdugo R, de la Hoz K. Strength and stiffness of coarse granular soils. In: Soil stress–strain behavior: measurement, modeling and analysis. Springer; 2007. p. 243–52. [14] Verdugo R, de la Hoz K. Caracterización geomecánica de suelos granulares gruesos. Revista Internacional de Desastres Naturales, Accidentes e Infraestructura Civil 2006;6(2). [15] Ramamurthy T, Gupta K. Response paper to how ought one to determine soil parameters to be used in the design of earth and rockfill dams. In: Proceedings of Indian geotechnical conference. New Delhi, India: Indian Geotechnical; 1986. p. 15–9. [16] Mery D. BALU: a toolbox Matlab for computer vision. pattern recognition and image processing; 2011.. [17] Mitchell JK, Soga K, et al. Fundamentals of soil behavior. New York: Wiley; 1976. [18] Li Y, Xu Y, Thornton C. A comparison of discrete element simulations and experiments for sandpiles composed of spherical particles. Powder Technol 2005;160(3):219–28. [19] Kloss C, Goniva C, Hager A, Amberger S, Pirker S. Models, algorithms and validation for opensource DEM and CFD–DEM. Int J Prog Comput Fluid Dyn 2012;12(2):140–52. [20] Hertz H. Ueber die berührung fester elastischer körper. J reine angew Math 1882;92:156–71. [21] Mindlin R. Compliance of elastic bodies in contact. J Appl Mech 1949;16:259–68. [22] Ai J, Chen J-F, Rotter JM, Ooi JY. Assessment of rolling resistance models in discrete element simulations. Powder Technol 2011;206(3):269–82. [23] Biarez J, Hicher P-Y. Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires. Rev francaise genie civ 1997;1(4):607–31. [24] Holtz W, Gibbs H. Triaxial shear tests on pervious gravelly soils. J Soil Mech Found Eng Div ASCE 1956;82(1):1–22.
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model Antonio Salazar, Esteban Sáez ⇑, Gislaine Pardo Department of Structural and Geotechnical Engineering, Pontificia Universidad Católica de Chile, Vicuña Mackenna 4860, Macul, Santiago, Chile
a r t i c l e
i n f o
Article history: Received 2 October 2014 Received in revised form 20 January 2015 Accepted 23 February 2015
Keywords: DEM Direct shear test Micro-parameters Rolling friction model Real grain size distribution
a b s t r a c t In this research, a direct shear test was modeled using the 3D Discrete Element Method (DEM). The results are compared against experimental data. The real sand was modeled as about 70,000 spheres, with sizes consistent with the real grain size distribution, using a rolling friction model to include the sand’s grain shape. It is shown that, the stress path can be appropriately reproduced, but the dilatancy were very difficult to replicate. Finally, specimens with up-scaled material were tested. Although the unscaled material was the most representative, the response of the particles that were up-scaled twice was slightly less accurate, but with a considerable decrease of computation time. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction In this paper, a 3D Discrete Element Method (DEM) model of the direct shear test was developed to reproduce the response of coarse sand in order to explore the ability of this modeling approach to reproduce this widely used test considering a realistic grading distribution. Several features of this method were explored, such as the ability to monitor and analyze displacements and contact forces of all particles composing the material. This laboratory test is widely used in geotechnical engineering, because of its simplicity and easy analysis of the results. The test uses a box with two halves, where the lower half is moved horizontally while the upper half remains fixed. A confinement force is applied to the top wall, which is free to move vertically. The registered reaction force depends on the properties of the shear band resulting from the induced relative displacement between both rigid half-boxes (Fig. 1). Several studies have shown that it is possible to qualitatively reproduce laboratory results with the discrete element technique using 2D models, despite the planar formulation limitations [1,2]. Cui and O’sullivan [3] were probably the first to perform a comparison of physical test data and direct shear test DEM models using perfect metallic spheres. Later, Yan and Ji [4] compare their DEM results with real soil direct shear tests performed on irregular ⇑ Corresponding author. E-mail addresses: [email protected] (A. Salazar), [email protected] (E. Sáez), [email protected] (G. Pardo). http://dx.doi.org/10.1016/j.compgeo.2015.02.017 0266-352X/Ó 2015 Elsevier Ltd. All rights reserved.
limestone rubbles. In their approach, they used clumps to simulate the shape of the real rubbles. In both studies, a good match between numerical results and real test data was achieved. Most recently, Härtl and Ooi [5] used the Jenike direct shear test to investigate the particle shape and microscopic friction on the bulk friction. They used single sphere particles and a clump of two overlapping spheres to investigate the shape effect. Their results show that the particle interlocking has a greater effect than the packing density on the bulk friction. In the present research a 3D approach was chosen, because the 2D models fail to accurately predict the peak and residual friction angles under shear flows [6], due to out-of-plane migration of the particles, among other effects. In recent years, the microprocessor industry has made little improvement in the core’s clock speed, giving advantage to multicore architectures [7,8] and the number of processors per machine will continue to increase. This change encourages the use of highperformance parallel software to fully take advantage of this technology. The DEM is especially suitable for this type of code because its elements have independent degrees of freedom, which makes it simple to parallelize. Harnessing this benefit, the coarse sand was initially modeled without size scaling. Nevertheless, to reduce the computational cost and to allow the use of a large number of particles, simple individual spheres were used with a rolling friction model representing the particle shape effect [9]. The effect of the model’s micro-mechanical parameters were studied with a sensitivity analysis. Two specimens with different initial void ratios considered, in order to highlight the effect of the system’s compactness on the response. We present the
84
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Sand grains photos Fig. 1. Direct shear test.
response of the calibrated model as compared to the laboratory results, which were obtained from previous stages of this research [10,11]. Also, a brief analysis of the chain force and evidence of outof-plane particles migration is shown. Despite advances in computing technology and numerical algorithms, computational cost is still the major limitation of DEM. Due to this issue, a scaled material was investigated and also compared with the experimental data. Due to laboratory restrictions, several down-scaling methodologies exist to test representative soil samples using small grain sizes. Among them, the parallel graduation method [12] was found suitable and able to capture the original soil response [13–15] with a finer gradation curve. In this paper, we use the parallel graduation method, but in the other direction (up-scaling), to develop specimens with particles twice and four times bigger than the original material. The scaled specimens were tested with the same calibrated parameters as the unscaled specimen.
(b) Grain size distribution of the real and simulated material Fig. 2. Grain samples and grain size distribution.
2. Real sand and DEM model description 2.1. Tested material
2.2. Brief description of DEM
The selected material was an angular non-cohesive coarse sand. Its main properties are listed in Table 1. A series of direct shear tests were performed on this sand under three different confinement pressures: 160 kPa, 80 kPa and 40 kPa. The laboratory test was instrumented with force and displacement transducers. Fig. 2a shows a series of photos, taken with a microscope, of the sand’s grains and processed with Balu toolbox for Matlab [16] to extract their apparent boundaries. Under a visual analysis, particles of all sizes have a sphericity of about 0.8 and a roundness of about 0.5, hence a subrounded shape [17]. Fig. 2b shows the sand grain size distribution. According to the Unified Soil Classification System (USCS), it corresponds to a poorly graded sand (SP) and all the material is retained in four sieves. The sand has one percent of fines (diameter less than 0:075 mm) and most of the mass is retained in sieve N 16, which are the biggest grains (between 2.16 mm and 1.18 mm of diameter).
The DEM is a numerical technique used to simulate the behavior of granular aggregates, jointed rock masses and chemical powders, among others. The main feature of this methodology is that the material is simulated as an assemblage of individual particles that interact with each other. This approach allows to capture the relative movements and rotations of the particles without the need of a sophisticated constitutive model. Also, it can achieve large displacements because each particle has its own independent degrees of freedom. For granular applications, the dynamic equilibrium condition uses the Eq. (1) (Newton’s second law) and Eq. (2), where M i and Ii are the mass and the moment of inertia of each particle i; Fji and Mji are the contact forces applied by particle j on i and the moment induced by this force; Fe and Me are the resultant of the body-force field over the i’th particle (only gravity in this case) and the torque because of this force.
€i ¼ Mi x
Table 1 Sand properties.
X ji F þ Fe
ð1Þ
j
Property
Value 3
Density: c [kg/m ] Void ratio: e0 Porosity: n0 Angle of internal friction: /0 Cohesion: c0 [kPa]
1578 0.679 0.405 35 8
€i ¼ Ii H
X ji M þ Me
ð2Þ
j
The time is divided into intervals, and the displacements are calculated in each step by the double integration of Eqs. (1) and (2). This time step must be small enough to avoid numerical instabilities. In
85
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
this research, the maximum time step was selected as the minimum of the Hertz or the Rayleigh’s time step [18]. The DEM software used to solve the problem was LIGGGHTS [19] which is a 3D Open Source Code with several granular features already implemented and ready to run in parallel. Further information can be found in the developer’s website.
Fij
ei Centroid j
2.3. Inter- granular contact modeling Eq. (3) shows contact model adopted in this investigation, which is based in Hertz’s work [20]. In this equation, E is the equivalent Young modulus (Eq. (4)), R is the equivalent particle radius (Eq. (5)) and dn is the overlap of the contact. mi and Ri in Eqs. (4) and (5) are the Poisson ratio and the radius of the particle i which is in contact with particle j.
pffiffiffi 3 4 Fn ¼ E R d2n 3
ð3Þ
1 m2j 1 m2i 1 ¼ þ Ej Ei E 1 R
¼
ð4Þ
1 1 þ Ri Rj
ð5Þ
The magnitude of the tangential force is given by the Mindlin’s theory [21] and limited by the Coulomb friction criteria shown in Eq. (6), where l is the friction coefficient. If this limit is exceeded, relative sliding takes place and the tangential force remains constant.
F tmax 6 ljF n j
ð6Þ
Since the particles being considered are perfect spheres, the Ai’s et al. [22] Elastic–plastic spring-dashpot rolling friction model was considered to take into account the shape effect between the sand’s grains. This model adds a torque (Mr ) in the dynamical equilibrium of the particles, which includes a spring torque (Mkr ) controlled by the rolling viscous damping coefficient (g) and a viscous torque (M dr ) controlled by the user by the rolling friction coefficient (lr ). The model implemented is shown in Eqs. (7)–(10).
M r ¼ M kr þ M dr
ð7Þ
The spring torque (Mkr ) is the product of the rotation of the particle (Dhr ) and the rolling stiffness (kr ). This torque is limited by Eq. (9), which is controlled by lr , similar to the Coulomb’s friction model. The rolling spring stiffness is defined in Eq. (10), where hm r is the full movilization rolling angle.
ej
Centroid i
−Fij Fig. 3. Real contact.
2.4. DEM material The granular material is simulated as an assembly of spheres with a radius according to the grain size distribution of the tested sand. Fig. 2b shows the comparison between the grain size of the real and simulated material. The model does not include the onepercent of fine sand because it would increase the computational cost considerably and it does not have any significant influence on the material’s macroscopic behavior. The percentage of sand mass retained in each of the sieves of the ASTM gradation test was simulated with particles whose diameters were linearly distributed between the opening size of the sieve where they mass was retained and the opening size of the preceding sieve, as can be seen in Fig. 4a. The two exception is that the finest part was not modeled and the diameter of the coarsest particles was limited to 1.6 mm. Fig. 4b shows that the amount of the modeled mass accumulated per sieve is slightly greater than the real material, because the mass of the finest portion that is not included is distributed among the other sieves. 2.5. Boundary conditions and fabric generation All the walls, except the top one, were modeled as rigid boundaries. This rigid boundary condition does not have inertia and its velocity is imposed to control its position during the simulation. Walls interact with the particles, following the Hertz’s contact model, to calculate the repulsive forces. The top boundary is modeled as a servo-controlled wall, in which its speed is related to the current stress by the algorithm shown in Eq. (12). The process is controlled by the constant kp ; ki and kd , which multiply the error between the current force in the wall and the target value (Det ), P the cumulative error ( t Det ) and the error derivative, respectively.
X Det DetDt Det þ kd Dt t
M kr;tþDt ¼ M kr;t kr Dhr
ð8Þ
V wall ¼ kp Det þ ki
jMkr;tþDt j 6 lr RF n
ð9Þ
The specimen is composed of approximately 70,000 particles. All the particles are randomly generated inside the box and are divided in three layers of 0.8 cm each. After the first layer is created, the particles fall freely until the static equilibrium is reached due to the gravity; then the process is repeated with the other two layers. After the settlement is completed, the servo-controlled wall is activated to confine the soil. Finally, a finite velocity is imposed on the lower half in order to reach a horizontal target displacement of 8 mm. The model, shown at the beginning of the shear process, is shown in Fig. 5. In order to get a loose or a dense fabric (eoi ), different coefficients of friction (li ) were used in the settlement process and during the shearing modeling. If lower values of li are used during generation, a dense fabric is obtained. The chosen Poisson’s ratio coefficient is a typical value used to simulate sand, while E was chosen to have a realistic contact stiffness. During shearing, the
kr ¼
lr RF n
ð10Þ
hm r
The viscous torque (M dr ) is controlled by g (Eq. (11)) and is only active while the rolling torque is not fully mobilized, assisting particle stabilization more than energy dissipation [22].
( Mdr;tþDt
¼
pffiffiffiffiffiffi 2g Ikr h_r
if jMkr;tþDt j < lr RF n
0
if jMkr;tþDt j ¼ lr RF n
ð11Þ
The reason for using this kind of model is because the contact force between two particles does not always pass through the centroid of the particles (an implicit assumption when spheres are used), as depicted in Fig. 3.
ð12Þ
86
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93 Table 2 Tested sand properties. Property
Specimen A (dense) 3
Specimen B (loose)
cs [kg/m ]
2650
2650
m li gi lri
5:0 108 0:256 0 0 0 0:489
5:0 108 0:256 0:3 0 0 0:594
E [kg/m2]
e0i obtained
in Table 2. The difference between these two specimens is the initial coefficient of friction li and consequently, the obtained void ratio e0i . The unit solid weight cs , used for both specimens, was obtained from laboratory measurements. The Young Modulus E and the Poisson’s ratio m are required to calculate the contacts’ stiffness and were selected among values proposed in the literature [1–4,6], but were small enough to avoid a time step that was too small. Fig. 6a–c shows the results using different combinations of the coefficient of friction l, the coefficient of rolling viscous damping g and the coefficient of rolling friction lr during the shearing phase. The shear stress, used in these comparisons, was calculated according to Eq. (13), where T is the horizontal reaction force of the upper half of the box (Fig. 1) and A is the effective shear area which depends directly on the horizontal displacement of the lower half (Dh ).
(a) Particle number histogram
rs ¼
(b) Mass histogram Fig. 4. Mass and particle number histograms.
z y
x
Fig. 5. Developed DEM model.
value of l; lr and g were modified to highlight their effect on the macroscopic response (Table 2). No rolling friction is considered during specimen generation.
3. Sensitivity analysis Two specimens were modeled and tested in order to explore the effects of the micro-mechanical parameters on the macroscopic response. The properties used for the model construction are listed
T AðDh Þ
ð13Þ
According to Fig. 6a, the shear stress comparison shows that the friction coefficient l has an influence for small displacements, and it tends to increase the peak value by approximately 30% when varied from 0.3 to 0.6. In the studied cases, the large-strain stress value is approximately the same for horizontal displacements larger than 4–5 mm, indicating that by varying only the Coulomb’s friction coefficient it is not possible to control the response of the model at large strains. Regarding volume evolution, the friction coefficient l increases the specimen dilatation without considerably affecting the path’s shape (Fig. 6a). Indeed, due to increase of particle friction, a looser configuration could be reached at large strains. Fig. 6b shows the sensibility on the rolling viscous damping g. In this case, a direct effect of this parameter on the stress path or on the displacement responses is not as clear, but a slight increment in the peak value is detected. The little effect of this parameter indicates that during most of the shear process the viscous torque is fully mobilized (Eq. (11)) and that the particles are rolling-sliding. Fig. 6c shows that the most important effect of the coefficient of rolling friction lr is the increase of the shear stress reached at large displacements. This parameter also increases the peak value, but in a less pronounced manner than l. The displacement comparison in Fig. 6c shows that lr increases the large-strain volume of the specimens. This influence cannot be noticed at a low-strain, where all the curves seem to follow the same path, highly influenced by the initial material fabric. As lr takes into account the particle’s shape, large values of this parameter can be associated with grains that are less round. In this case, it is reasonable to obtain a supplementary geometric interlocking at large strain, associated with a large apparent macroscopic angle of friction. Due to the shape effect, a looser critical state is also expected. Nevertheless, in general terms, the laboratory displacement path was very difficult to replicate. All three figures show that the main difference between the specimen A and B is that A gives a larger peak at low strain (a
87
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) μ influence (μr =0.2 & η=0.2).
(c) μr influence (μ=0.4 & η=0.5).
(b) η influence (μ=0.4 & μr =0.5).
Fig. 6. Sensitivity analysis of the DEM parameters.
typical response of a relatively dense sand), while the response for specimen B correspond to the usual stress–strain path of loose sand. Also, the shear stress comparison shows that specimen A has a higher initial slope (low-strain stiffness), in agreement with the experimental response. This slope can be controlled by the stiffness of the contacts and the initial void ratio. Specimen B has a lower large-strain vertical displacement in comparison to specimen A, which is also a typical experimental result for loose sands. The behavior difference was expected because of the different initial void ratio obtained for both simulated specimens. Compared to the experimental maximum (emax = 0.74) and minimum (emin = 0.44) void ratios, specimen A has a relative density of 84% classifying it as a very-dense material while specimen B is a medium-dense material with a relative density of 49%.
Table 3 Calibrated material parameters. Property
Specimen C
E [kg/m2]
lri l g lr
5:0 109 0:256 0:3 0 0 0:3 0:1 0:4
Obtained: c [kg/m3] Obtained: e0i
1612 0:643
m li etai
88
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
4. Calibrated model against experimental data A third specimen was created with the properties shown in Table 3, a higher contact stiffness was used in order to obtain a higher low-strain stiffness, while using a looser fabric. A greater void ratio was obtained, even when the same li was used, because higher contact stiffness avoids larger overlaps. The void ratio
obtained is slightly lower than the experimental value, which is expected because of the different shape of the real grains and the simulated particles [23]. With a confinement of 160 kPa, the model can appropriately reproduce the shear stress path during the whole test (stress calibration in Fig. 7). For the 80 kPa and 40 kPa confinements, the model does not reproduce the peak value, but the residual stress
Fig. 7. Calibrated material response.
(a) Initial force chain
(b) Final force chain Fig. 8. Initial and final chain force.
89
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Polar contact force distribution at the beginning of the shear test
(a) Initial velocity field
(b) Polar contact force distribution at the end of the shear test
(b) Final velocity field
Fig. 9. Polar distribution of the contact forces at the beginning and the end of the shear test on the shear plane.
is accurately reproduced for all the confinements. This might be due to a higher density of the sand in the laboratory tests compared to the DEM models. The dilatancy of the simulated specimen increases as the confinement stress decreases, in agreement with laboratory results (displacement calibration in Fig. 7). Additionally, the shapes of the curves are very similar. Nevertheless, the model fails to accurately predict the final sample’s volume. Indeed, this large strain volume is highly controlled by the particles’ shape, one of the major limitations of the method due to the use of spheres. Thus, the rolling friction model used to approximately include the particles’ shape effect, is not able to fully reproduce the large strain state of the reference soil. In the case of 160 kPa test, the final laboratory void ratio is 0.72; while the simulation gives a value of 0.703. The initial difference was reduced by the increase of the model’s dilatancy, but difference still remains, due to the shape of the sand grains. Fig. 8 shows the initial and final network of contact forces. The forces are evenly distributed at the beginning of the test (Fig. 8a), which is expected due to the initial vertical compression. During the test, a diagonal soil strut was formed, as can be seen in Fig. 8b. This might be due to the high geometric interlocking of
Fig. 10. In-plane velocity field (imposed displacement plane).
the material used, which prevents the particles of the central part of the sample to relatively slide. The aforementioned is evidenced when orientations and amplitude of the contact forces are analyzed (Fig. 9). Fig. 9a shows the distribution of the contact force at the beginning of the shear test, after the settlement process, while Fig. 9b shows the contact force distribution at the end of the test. At the beginning, the vertical contact forces are almost the double the horizontal contact forces, because of the vertical confinement imposed to the specimen. At the end of the test, the predominant contact forces are no longer vertical, but almost have the same orientation as the strut (approximately 19:5 ). The magnitude of the force in this orientation (between 0 and 36 ) is about four times larger than the initial vertical forces, while the magnitude of forces in the normal direction to the strut (between 90 and 126 ) is about one third of the initial vertical magnitude, highlighting the stress-state rotation during the shear test. In Figs. 10 and 11 the test was divided into small regions of 3 1.2 mm and the particles’ velocities were averaged in each of these regions in order to reduce the number of arrows and to represent the velocity as a continuous field. When the mean radius of the grain size distribution used is 0.5 mm, each of these zones represent the average value of about 360 particles. Fig. 10 shows a lateral view of the initial and final velocity field of the particles over
90
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Initial velocity field Fig. 12. Scaled material real grain size distribution.
(b) Final velocity field Fig. 11. Out-of-plane velocity field.
the imposed displacement plane. At the beginning of the test, all particles have similar velocities, which means that the specimen behaves as a continuous material (Fig. 10a). When large strains are reached, the bottom particles have a higher speed than top ones, hence a sliding surface dividing the specimen in two halves appears (Fig. 10b). A shear band of about 2.5 times the mean diameter of the particles is depicted in this figure. The particles do not follow a clear path in the out-of-plane direction (normal to the imposed displacement plane), also the magnitude of the average speed in this plane is about 30 times less than the in-plane speed. Nevertheless, the out-of-plane component of the speed indicates that there was some particle migration across this direction, justifying the choice of a 3D model (Fig. 11b). At the beginning of the test, the particles have moved almost completely in the vertical direction, with a very small out-of-plane component (Fig. 11a). It is interesting to note that out-of-plane migration particle movement is not symmetric across the observed section, probably due to the random generation of particle size and location among the sample. This asymmetry increases for large strains.
the simulations [13]. The purpose of using the up-scaled material is to investigate the decrease in the computation time and the accuracy degradation with the change in the ratio between the specimen size and the maximum particle diameter (Dspecimen =Dmax ). Two specimens with twice and four times bigger particles were studied. The grain size distribution is shown in Fig. 12. The Dspecimen =Dmax ratio of the twice scaled material is 7.5, while the ratio of the four times scaled material is 3.75. Although the Dspecimen =Dmax ratio of the four scaled material is out of the recommended range (Dspecimen =Dmax [24]), this scaled factor was tested for completeness. The time-computing comparison shown in Table 4 was done using only one core, to avoid the performance decrease when using a large number of processors with low data. Also, due to the large amount of time it would take to run the whole test in a single core, and because the purpose of this comparison is to illustrate computing time reduction, each of these simulations consists in a tenth of the required steps to complete the shear and the settlement process. All the simulations were carried out using the same timestep, hence the three specimens include the same number of steps and the same total simulated time. When scaling two-times, the computation time is reduce by 94% (20 times), while the four times scaled is reduce by 99.5% (200 times). This is because the reduction in the number of particles (86.8% and 98.7%) leads to an important decrease in number of contacts (96.4% and 99.8%), which defines the size of the DEM model. It is clear that even with a small coarse-graining scaling, the computation time decreases dramatically, allowing the development of large models with a reasonable running time. Because large fluctuations appear when larger particles are considered, both coarse-grained specimens were tested five times with the same set of parameters shown in Table 3 and an average was calculated (Figs. 13 and 14). The fluctuation in the results of
Table 4 Computational cost comparison. Unscaled
5. Parallel gradation In order to cut down the running time, the parallel gradation technique has been explored to reduce the number of particles in
Time (min) Number of particles Number of neighbors
Two-times scaled
Four-times scaled
Value
[%]
Value
[%]
Value
[%]
1903 64,944 2,754,758
100 100 100
107 8610 100,596
5.6 13.2 3.6
9.3 845 4033
0.5 1.3 0.15
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
91
(a) Twice scaled average - 160 [kPa]
(b) Twice scaled average - 80 [kPa]
(c) Twice scaled average - 40 [kPa]
Fig. 13. Shear test results with two-times scaled material.
the scaled materials are explained because when using the scaled sample, the number of contacts is significantly reduced, hence any change affects proportionally more in the macroscopic response than the unscaled case where a large number of contacts exist (Fig. 14). The paths of the scaled material, which was scaled two times, have a low deviation from the average (Fig. 13), while it increases greatly for the specimen that has been scaled four times (Fig. 14). This high fluctuation in the stress response is not as clear in the displacement results, because the servo-controlled boundary implementation responds more to a trend than to a single value. As the algorithm relates the stress with the velocities, and not directly to the displacements, the change in the vertical stress must be held for some time steps in order to see changes in the displacements.
When averages are compared with laboratory and unscaled material results (Figs. 13 and 14), very similar stress paths are obtained for both scales; however, the little increase in the stress path of the material which was scaled four times evidences that the twice scaled is lightly more accurate. The displacement comparison shows that the material which was scaled four times has a higher displacement than the twice scaled material, which is explained by the increase of the particle size. The material that was scaled four times is an inconvenient choice to reproduce the shear test, because of the high fluctuation of the results, despite the huge decrease of the computational cost. However, the satisfactory enough simulation of the stress paths and the huge decrease in running time make this strategy a suitable choice for problems where the ratio Dspecimen =Dmax is big enough to avoid boundaries effects.
92
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
(a) Four times scaled average - 160 [kPa]
(b) Four times scaled average - 80[kPa]
(c) Four times scaled average - 40[kPa]
Fig. 14. Shear test results with four times bigger material.
6. Conclusions When a discontinuous methodology such as DEM is used, the initial fabric density is very difficult to control compared to other continuum mechanics based methodologies. Specifically, in the implemented direct shear model, it was very difficult to find the correct combination of friction parameters and sphere generation to satisfactorily reproduce the initial laboratory relative density. The rolling friction model was used to represent the grain’s shape effect. It is concluded that the grain’s shape is essential to properly represent the material response under large strains. In general terms, the shape of the curves obtained with the DEM model approximately matches the experimental data. Indeed, the stress response of the model is very similar to the laboratory results, but the displacement paths are similar only qualitatively, especially at high horizontal displacements.
The particle force network clearly illustrates a shear zone and the formation of a strut. Additionally, the 3D model that was developed allows to point out an out-of-plane particle’s migration. Hence, the use of direct bi-dimensional discontinuous models under macroscopic plane-strain condition models must be done carefully. Because of the accuracy of the results compared to experimental data, the unscaled material seems to be the best model to reproduce the shear test. However, when the computational cost is relevant, a parallel grading curve might be considered. The results show a 94% (20 times) reduction in the running time with the model scaled twice, and just a little increment in the stress path and the displacement response. With the material that was scaled four times, the time reduction attained was 99.5% (200 times), but with a large fluctuation in the results, which makes this model unreliable. This leads us to conclude that the maximum up-scale
A. Salazar et al. / Computers and Geotechnics 67 (2015) 83–93
value must be studied for each problem probably in terms of the ratio Dspecimen =Dmax . The small loss in accuracy is totally compensated by the huge time reduction when using coarse-grained materials, based on the parallel grading technique. Hence, this methodology, theoretically, allows us to make bigger and more complex models with reasonable computing time. Nevertheless, the size increase technique raises the fluctuation of results and must be repeated several times to identify the tendency. Acknowledgement This work was possible due to a grant from the Fondo Nacional de Desarrollo Científico y Tecnólgico FONDECYT, Grant No. 11100157. References [1] Thornton C, Zhang L. Numerical simulations of the direct shear test. Chem Eng Technol 2003;26(2):153–6. [2] Masson S, Martinez J. Micromechanical analysis of the shear behavior of a granular material. J Eng Mech 2001;127(10):1007–16. [3] Cui L, O’sullivan C. Exploring the macro-and micro-scale response of an idealised granular material in the direct shear apparatus. Geotechnique 2006;56(7):455–68. [4] Yan Y, Ji S. Discrete element modeling of direct shear tests for a granular material. Int J Numer Anal Methods Geomech 2010;34(9):978–90. [5] Härtl J, Ooi JY. Numerical investigation of particle shape and particle friction on limiting bulk friction in direct shear tests and comparison with experiments. Powder Technol 2011;212(1):231–9. [6] Fleischmann PM, Drugan JW. Quantitative comparison of two-dimensional and three-dimensional discrete-element simulations of nominally twodimensional shear flow. Int J Geomech 2013;13(3):205–12. [7] Geer D. Chip makers turn to multicore processors. Computer 2005;38(5):11–3. [8] Borkar S, Chien AA. The future of microprocessors. Commun ACM 2011;54(5):67–77.
93
[9] Wensrich C, Katterfeld A. Rolling friction as a technique for modelling particle shape in DEM. Powder Technol 2012;217:409–17. [10] Pardo G, Sáez E. Experimental and numerical study of the arching effect in coarse sand. Comput Geotech 2014(57):75–84. [11] Pardo G. Estudio experimental y numérico del efecto arco. Master’s thesis, Pontificia Universidad Católica de Chile, Santiago, Chile; 2013. [12] Lowe J. Shear strength of coarse embankment dam materials. In: Proc, 8th int congress on large dams, vol. 3; 1964. p. 745–61. [13] Verdugo R, de la Hoz K. Strength and stiffness of coarse granular soils. In: Soil stress–strain behavior: measurement, modeling and analysis. Springer; 2007. p. 243–52. [14] Verdugo R, de la Hoz K. Caracterización geomecánica de suelos granulares gruesos. Revista Internacional de Desastres Naturales, Accidentes e Infraestructura Civil 2006;6(2). [15] Ramamurthy T, Gupta K. Response paper to how ought one to determine soil parameters to be used in the design of earth and rockfill dams. In: Proceedings of Indian geotechnical conference. New Delhi, India: Indian Geotechnical; 1986. p. 15–9. [16] Mery D. BALU: a toolbox Matlab for computer vision. pattern recognition and image processing; 2011.