Numerical simulations of granular soil using elliptical particles

Computers and Geotechnies 16 (1994) 153-169 © 1994 Elsevier Science Limited Printed in Great Britain. All rights reserved 0266-352X/94/$7.00 ELSEVIER

NUMERICAL SIMULATIONS OF GRANULAR SOIL USING ELLIPTICAL PARTICLES

Tang-Tat Ng Department of Civil Engineering University of New Mexico, Albuquerque, NM 87131.

ABSTRACT Two-dimensional numerical experiments on a random array of circular and elliptical particles were performed using the discrete element method (DEM) to study the behavior of idealized granular soil. A new program called ELLIPSE2 was developed and used to show the stress-strain behavior of granular soft under monotonic drained compressive loading and cyclic constant volume loading at various shear strains (10~ to 103). In these numerical simulations, the granular soil was represented by a 2-D random array of elastic, rough quartz particle. The random array consists of a mixture of circular and elliptical particles. The numerical results were compared with available experimental data on actual sands with good agreement. These results are extremely encouraging and clearly show the potential of a program such as ELLIPSE2 for studying granular soil behavior.

INTRODUCTION

The DEM simulations using assemblies of round particles (discs or spheres) have shown capabilities for studying the behavior of granular media. Many researchers use the DEM codes to simulate the behavior of real sands. While the shear strength of the assemblage agrees well with experimental observations on glass beads, it is low compared with real sands. To increase the shear strength of the assemblage, dense arrays or artificially prohibiting particle rotation were used. One of the many possible reasons for the low shear strength behavior is the particle shape (roundness of the particles). One expects that by changing the shape of the particles from circular to elliptical, there will be better agreement with the data for sands. Therefore, elliptical particles have been introduced recently [Ting, 1991~; Ting and Corkum, 19922; Rothenburg and Bathurst, 19913, 19924; Ng, 19923]. Positive results have been shown for monotonic loading. The goal of this work is to show that better cyclic simulations can also be obtained using 153

154

elliptical particles. This paper briefly describes the newly developed numerical code ELLIPSE2 (ellipse in 2-D) and the corresponding contact detection scheme. A specimen consisting of circular and elliptical particles was isotropically compressed and then sheared under monotonic drained and cyclic constant volume (undrained) loading.

THE DEM AND COMPUTER PROGRAM ELLIPSE2 The DEM employs an explicit finite difference scheme to simulate the interactions between rigid bodies in relatively dense assemblages. It was originally developed by Cundall [19716] to analyze problems in rock mechanics.

In addition, the method has been used to

investigate flow in granular media and process engineering, as well as in geotechnical problems. The proceedings of the 1st U.S. and 2nd International DEM Conferences [Mustoe et al., 19897, WiUiams and Mustoe, 1993s] are excellent references. The DEM has also been linked with the boundary element method and the finite element method to analyze stresses around tunnel openings and to simulate cone penetration tests [Lorig and Brady, 19849; CundaU and Hart, 1989t°; Huang et al., 1992H]. The methodology and detailed descriptions of the DEM are available in the literature [CundaU and Strack, 1979~2; Rothenburg and Bathurst, 1989~3; Ng, 1989"; Ting et al., 1989~s]. Two-dimensional and three-dimensional simulations using circular particles have demonstrated the ability to simulate soil behavior qualitatively, and quantitatively in some cases [e.g. Cundall and Strack, 1979~2; CundaU, 1988t6; Rothenburg and Bathurst, 1989~3; Ng, 198914; Ng and Dobry, 1991a~; 1992~s]. Recently, 2-D monotonic drained simulations using elliptical particles have been reported [Rothenburg and Bathurst, 19914, 1992s; Ng, 19926: Ting and Corkum, 19923]. When the non-spherical shape of sand particles has been approximated by an ellipse, the simulations show a better agreement with actual sands.

Program ELLIPSE2 Program ELLIPSE2 was developed based on the program CON-BAL2-1 [Ng and Dobry, 1991bt9]. It uses a non-linear Hertzian normal contact law and a linear normal-force dependent

155

tangential contact relationship.

The normal and tangential stiffness of the contact are,

respectively, I(~ = 2G,r/(l - ~,,) and Kt = 4G,r/(2 - J,,), where r = [3(1 - p,)NR/SG,] ',3 is the radius of the contact area, G,, ~,,and R are the shear modulus, Poisson's ratio, and radius of the spheres, respectively, and N is the normal force. The contact slides when the tangential force, T, reaches/~,N (where p, is the friction coefficient), hence IT[ < /~,N. Particle rotation is allowed while a periodic space is used to eliminate boundary effects. The major effort in implementing elliptical particles is the contact detection scheme which is discussed in the next section.

Contact Detection Between Two Ellinses

Contact detection schemes for two ellipses have been reported elsewhere [Rothenburg and Bathurst, 19913; Ting and Corkum, 19922]. However, numerical difficulties such as identical roots, complex roots, and numerical inaccuracies may occur for some orientations. Furthermore, schemes in 3-D are not as straightforward as in 2-D since the intersection of two ellipsoids is a curve. The author develops a detection scheme which does not encounter the difficulties mentioned above and the same technique can be applied to 3-D models. A similar method has been used by Ting (19932°). Figure 1 shows three contact scenarios for two ellipses: (a) no contact, Co) "simple" contact with one or two intersection points, (e) "complex" contact with more than two intersection points. Case (c) is not considered because it is physically incorrect and denotes numerical simulation error only (excessive overlap). The following methodology addresses cases (a) and Co). As shown in Fig. 2, case (b) is used to illustrate the method. The contact point, P, is defined as the mid point of the line segment jointing point 1 to point 2. Point 1 (or 2) is the point on one ellipse which creates the smallest potential with respect to the other ellipse. The potential of a point at (x,y) with respect to an ellipse is defined as:

f ( x , y ) = a z x 2 + a z x y + a 3 Y 2 + a 4 x + a s y + ae

with a~ = real number, i = 1, 2, 3 . . . . 6.

156

(a) No Contact

(b) Simple Contact

(c) Complex Contact

Figure 1. Three Scenarios Between Two Ellipses

Yl

alx2+a2xy+a3Y2+a4x+asy+a6=O

Icontact

~-~~o:o) /

'

r

X

Figure 2. Detection of Contact Between Two Ellipsoids

157

If fix,y) = 0, point P is on the ellipse (case b with one intersection point); if fix,y) > 0, point P is outside the ellipse (case a); if fix,y) < 0, point P is inside the ellipse (case b with two intersection points). The points 1 and 2, as shown in Fig. 2, are obtained by minimizing fix,y).

The

corresponding radii of curvature are obtained from basic analytical geometry. The problem then reduces to the contact problem of two circles with two radii of curvature. A detailed description of the mathematics can be found in Ng (1992 s) and Ng and Lin (199Yl). For two 3-D ellipsoids, fix,y) will change to f(x,y,z) but the same procedure can be used.

NUMERICAL SIMULATIONS OF GRANULAR SOIL USING ELUPSE2

Monotonic drained and cyclic constant volume experiments were performed on a random array of quartz elliptical and circular particles. The properties of these quartz particles are: shear modulus, G, = 28,957 MPa, Poisson's ratio, v, = 0.15, and friction coefficient,/z, = 0.5. Based on a preliminary study of an array using only 158 identical elliptical particles, a large fluctuation in the stress-strain curve is observed and more than 68 ~ of the computation effort was spent on contact detection. Therefore, a mixture of elliptical and circular particles is used to reduce the computational time. The array has 227 particles, of which 54 are elliptical particles with aspect ratios of 1.5. The remaining 173 particles are circular with a radius equal to the minor length of the ellipse. Although the majority of the particles are circles, the simulations show better agreement with sands than two-sized specimens with only circular particles. This array was consolidated isotropically to a porosity of 0.158 under a consolidation pressure, ~0, of 294 kPa. The configuration of this specimen after consolidation is shown in Fig. 3. Figure 4 shows some micromechanical statistics of the specimen, including: (a) A histogram of the distribution of the contact force angle in polar coordinates, Co) The distribution of varticle coordination number, (c) A histogram showing the distribution of the normalized mobilized angle (T//~0, (d) A histogram of the distribution of the normalized contact force (the contact force is normalized with the maximum contact force). In addition to these plots, other useful fabric statistics (Rothenburg and Bathurst, 19924) can also be generated at any loading step. This ability to extract all information from DEM

158

/

\ ‘\

(al

(cl

PARTICLE

CONTACT

ORIENTAION

NORMAL

DIRFICTION

lb)

PARTICLE

Cd)

--_ - i -- /-

/’

ORIENTAION

CONTACT

Figure 3. Geometry Configuration and Micromechanical of An Elliptical Specimen under An Isotropic Pressure = 294 kPa)

ANGLE

ANGLE

Statistics Consolidation

159

SO

. , . 2 ( } - - T . . . . . . . ,,

I

/'" /

(b)

-\(a)

"

, ,I,D"-'~""" x

"""~

>,-

zu zs

.

.

.

.

.

.

.

.

.

.

.

.

W

\

"

......-

mm.v,~.N~

o wn~ i,

/ /

0

I

CONTACT FORCE ANGLE

CONTACTS PER PARTICLE 30

30

: >.-

U Z W

o w n~

15

]0

S

0

l

(c)

'

(d) >-

Q W

b_

i,

0 0.0

0.5 I .0 NORMALIZEO MOBILIZE INTERPARTICLE FRICTION ANGLE

0 0.0

0 .S

NORMALIZEO

CONTACT FORCES

Figure 4. Micromechanical Statistics of An Assemblage of Elliptical Particles ( (To = 294 kPa)

1,0

160

particle simulations is extremely important and provides direct insight into what is taking place to the microstructure evolution during loading. This information can be used to interpret the macroscopic response of the particulate medium based on micromechanical features, and may lead to physically correct models to predict macroscopic behavior. Monotonic Drained Simulations

Two "CID" compressive tests, where o2 is kept constant while ol increases during shear, have been performed on two simulated specimens with identical porosity (0.158) and different initial confining pressur-~ (Oo = 294 and 588 kPa). The corresponding curves of ck = sin4[(o~ - oz)/(ol + o2)] and volumetric strain versus axial strain (el) are plotted in Fig. 5 for these two runs. Some interesting observations can be made: (1)

The curves for these two tests are similar to sands, with more dilation for the one with the lower initial confining pressure. The expected results show that the initial consolidation pressure affects the dilation behavior and the initial modulus but not the maximum shear strength;

(2)

The strain (ep) where the maximum q/p value occurs in these simulations is about 1.25% which is closer to the typical range for sands than simulations which use only circular particles. In the simulations using circular particles, with the same material properties and the contact law used in this paper, ep always occurs below 0.5 %

[Cundall, 198816; Ng, 198914; Ng and Dobry, 1991b19]; (3)

The maximum friction angle for these two tests is about 30°, which is comparable with the low end of the typical values found for sands [Lambe and Whitman, 196922]. A denser specimen should yield a higher friction angle.

Cyclic Constant Volume Simulations A very small cyclic shear strain % ,, 10-6 was applied to the specimen (Oo = 294 kPa) used in the monotonic drained simulations and the maximum shear modulus (Gmax) was calculated. Applying different 7c, Fig. 6 shows the variation of the normalized shear modulus, G/Gmax, and the damping ratio, h, calculated for the first cycle of constant volume loading. Also shown in Fig. 6 is the experimental bands for sands reported by Seed et al. [1984~]. The agreement is good. In contrast, simulations with circular particles gave a higher damping ratio

161

35 (a)

A

Go = 294 kPa

+

~

25

~ 2o ""

15

om

10 I

-05 I

0

1

Ep

I

I

2

3

4

1.5 (b)

o_. 0.s

J

0

-0.5

0

'

'

'

1

2 Axial Strain (%)

3

4

Figure 5. Results of Drained Monotonic Compressive Tests (CID) For Two Different Confining Pressures

162

1 . 0 ---~

L~ 0.5

band

(Seed

for s a n d s / e t al., 1984)

"

~

~i

0 Numedcal Simulation 0.0 10-4

10-3

Cyclic s t r a i n

10-2

amplitude

i0 -I

7c (%)

z0 ......................... i ........................ !........ - / b ) o

-

N = I cycle

i /

0 °P--I

bl)

E xp e r i ."m..e.n.t.a.1..... ~ for s a n d s . ~ _

I0

...........

........ band

/

o~ c~

o

T 10-4

( ( i ii~IIT I0-)

i i l illll)

i I s llsl

10-2

I0 -I

Cyclic strain amplitude 7c (%)

Figure 6. Shear Modulii and Damping Ratios versus Cyclic Shear Strain for An Assembly of Elliptical Particles with ~o ---294 kPa

163

than the experimental values for sands when % > 0.05 % [Ng, 198914]. The same specimen was also subjected to cyclic strain-controlled loading beyond one cycle for four different values of % (0.01, 0.03, 0.05 and 0.1%), without volume change. The numerical simulations are compared with the results of cyclic constant volume experiments on one type of sand (Finn and Bhatia, 1980~). This type of numerical simulation is also similar to an undrained cyclic torsional test. In these simulations, the difference (¢o- ¢2) can be interpreted as a pore water pressure change under the reasonable assumptions of incompressible fluid and a unique effective stress path. The "pore pressure ratio" is then defined as r, ffi I - ¢2/¢7o. A similar interpretation has been used by Finn and Bhatia (198024). When r, is smaller than 0.8, the average inertia force is g e n ~ l l y less than 0.8% of the magnitude of average contact force, which shows that the simulations are still quasi-static. However, the inertia force increases abruptly to about 40% as r, approaches 0.93 as a result of the significantly reduced contact forces. The numerical results, when r, is larger than 0.8, are also shown here. The results of a cyclic test with % = 0.05 % are shown in Fig. 7.

The calculated

hysteresis loops are shown in Fig. 7a, which shows a decrease of shear stress as the number of cycles increases. The shear stress decreases much faster in the first cycle. This shear stress decrease, or cyclic modulus degradation, is similar to experimental results obtained with cyclic triaxial strain-controlled testing [Li et al, 198825] and cyclic torsional strain-controlled testing on saturated sands [Ishibashi, 198826]. The graph of shear stress versus mean stress for this simulation, shown in Fig. To, moves progressively toward the origin. We can identify an effective stress failure envelope both in the compression and extension regions in Fig. 7b. These results are similar to cyclic triaxial strain controlled experiments on saturated sands [Li et al., 198825; Vucetic and Dobry, 198827]. As shown in Fig. 7c, the pore pressure ratio versus the number of cycles is similar to that observed in cyclic strain-controlled laboratory tests on sands [Finn et al., 197128; Cho et al., 197629]. In laboratory observations of cyclic tests, the rate of pore pressure buildup per cycle decreases with the number of cycles. The pore pressure fluctuates during each cycle with a net pore pressure increase after each cycle. The pore pressure buildup eventually approaches the consolidation pressure, where Au = ~o, r= --- 1 and thus "initial liquefaction" (also called the "o - 0 condition') occurs.

The same trend is also observed in Fig. 7c.

Figure 7d shows the trend of decreasing average coordination number, Cn, as the number of cycles increases. A very interesting aspect of this calculated curve is that the coordination

164

L =

~

Z

'"

/

(ed~l) $$OJ|S IllgqS

JeqwnN uo.qeu!pjooo

(gd~) sseJ|S mgLIS

" -

165

o

0.8 ~

2

..,~'~"

" o.o~

0.°f;

0

5

10

15

20

25

30

Number of Cycles N

Figure 8.

Pore PressureBuildup as a Function of Number of Cycles for A

Simulated Specimen(227 particles,2 Sizes, O'o = 294 kPa)

1

~

/

.75

~

0

.

2

%

UJ

._~ =

.5

I¢1 n

.25

/ I

0 Q.

0

f

I

I

I

I

I

I

I

f

2

4

6

8

10

12

14

16

18

20

Number of Cycles N ~gure9.

Pore Pressure Generation In A Sand Subjected to Cyclic Strain-Controlled Constant Volume Test (from Finn and Bhatia, 1980)

166

number fluctuates with a small net decrease before the specimen reaches the state of initial liquefaction at which point a large decrease and violent fluctuation occur. This indicates that most of the particles within the specimen are almost floating. Generally speaking, the particles form new contacts when 7° increases during a load cycle (specimen dilates), while contacts tend to be destroyed as % decreases. Figure 8 illustrates that the pore pressure ratio, ru, increases significantly as both % and the number of cycles increase. The rate of pore pressure buildup is faster when 7o is greater, but the rate of pore pressure buildup for constant % diminishes with the number of cycles. Comparing the simulations with cyclic constant volume experiments on a sand (as shown in Fig. 9) reveals similar qualitative trends.

CONCLUSIONS The DEM simulations using an assembly of circular and elliptical particles show a better representation of the behavior of granular soils under monotonic and cyclic loading than using only circular particles. The results of monotonic simulations, shown in Fig. 5, yield a better agreement with sands for the maximum friction angle and the strain where the peak strength occurs, than the results of using assemblies of only circular particles. Moreover, the damping behavior (see Fig. 7) of cyclic constant volume simulations has better agreement with sands than simulations using only circular particles. These results are very promising and show that one can use elliptical particles to study the behavior of cohesionless soils. Currently, the program ELLIPSE2 are being used to study the effect of particle shape, and inherent (geometry induced) anisotropy during specimen generation.

Finally, the 3-D

program ELLIPSE3 was recently completed at UNM and numerical simulations using the program are currently underway. It is expected that the 3-D program ELLIPSE3 will reveal and explain the behavior of cohesionless soils in a more convincing fashion.

167

ACKNOWLEDGMENTS The initial portion of this study was supported by the National Science Foundation under Grant MSS-9105988 and the cognizant program director is Dr. Mehmet T. Tumay. This support is gratefully acknowledged. The supercomputer time provided by CRAY Inc. is also greatly appreciated. The author also thanks the reviewer for his valuable comments.

REFERENCES .

Ting, J. "An Ellipse-based Micromechanical Model for Angular Granular Materials," Proc. ASCE Enfineerino Mechanics Special~ Conference, Columbus, OH, (1991) 1214-1218.

. Ting, J.M. and Corkum, B.T. "Computational Laboratory for Discrete Element Geomechanics," J. Conmutine in Civil Ene incerine, ASCE, Vol. 6., No. 2., (1992) 129-146. 3.

Rothenburg, L. and Bathurst, R.J. "Numerical Simulation of Idealized Granular Assemblies with Plane Elliptical Particles," Comvuters and Geotechnics, Vol. 11 (1991) 315-329.

4.

Rothenburg, L. and Bathurst, R.J. "Micromechanical Features of Granular Assemblies," Geotechnique, Vol. 42, (1992) 79-85.

.

Ng, T.-T. "Numerical Simulations of Granular Soil Using Elliptical Particles," Microstructural Characterization in Constitutive Modeling of Metals and Granular Media,

The ASME Summer Mech. and Materials Conferences, Tempe, AR, (1992) 95-118. 6.

Cundall, P.A. "A Computer Model for Simulating Progressive, Large-scale Movements in Block Rock Systems," Proc. Svmv. Int. Soc. Rock Mech., Nancy H, Art. 8. (1971).

. Mustoe, G.G.W., Henriksen, M. and Huffelmaier, H.P. Proc. 1st U.S. Conference on Discrete Element Methods, Golden, CO, (1989); also Engineering Computations, Vol. 9 No. 2, April (1992). 8. Williams, J. and Mustoe, G.G.W. Prec. 2nd International Conference on Discrete Element Methods~ Boston, MA, March 18-19, (1993). .

Lorig, L.I. and Brady, B.H.G. "A Hybrid Computational Scheme for Excavation and Support Design in Jointed Rock Media," Design and Performance of Under, round ~ , IRSM/BGS, Cambridge, England, (1984) 105-112.

10. CundaU, P.A. and Hart, R.D. "Numerical Modeling of Discontinua," Proc. 1st US Conf, on Discrete Element Methods, Golden, Co,(1989); also Engineering Computations, Vol. 9 No. 2, April (1992).

168

11. Huang, A.B., Ma, M.Y. and Lee, J.S. "Effects of Stress History on Cone Penetration Tests in Sand: A Micromechanical Study," Microstructural Characterization in Constitutive Modelin~ of Metals and Granular Media, The ASME Summer Mech. and Materials Conferences, Tempe, AZ, (1992) 69-80. 12. Cundall, P.A. and Strack, O.D.L. "A Discrete Numerical Model for Granular Assemblies," C,eotechnique, Vol. 29, No. 1, (1979) 47-65. 13. Rothenburg, L. and Bathurst, J.R. "Analytical Study of Induced Anisotropy in IdealiTed Granular Materials, Geotechniaue, Vol. 39, No. 4, (1989) 601-614. 14. Ng, T.-T. "Numerical Simulation of Granular Soils under Monotonic and Cyclic Loading: A Particulate Mechanics Approach," PJa,]~DJ~gllali~, Department of Civil Engineering, Rensselaer Polytechnic Institute, Troy, NY (1989). 15. Ting, J.M., Corkum, B.T., Kauffman, C.R., and Greco, C. "Discrete Numerical Model for Soil Mechanics," J. Geotechnical Engineering, ASCE, Vol. 115, No. 3., (1989) 379-398. 16. Cundall, P.A. "Computer Simulations of Dense Sphere Assemblies," in Microntechani¢$ o_fGranular Materials, Satake and Jenkins, eds. (1988) 113-123. 17. big, T.-T. and Dobry, R. "Numerical Experiments Using the Discrete Element Method," Proc. ASCE Entineering Mechanics Soecialtv Conference, Columbus, OH, (1991a) 1234-1238. 18. Ng, T.-T. and Dobry, R. "A Nonlinear Model for Soil Mechanics," Int. J. for Numerical and Analytical Methods in Geomechanics, Vol. 16, No. 4, (1992) 247-264. 19. Ng, T.-T. and Dobry, R. "CONBAL - Simulated Granular Material using Quartz Spheres with the Discrete Element Method," Report to NSF, Rensselaer Polytechnic Institute, Troy, NY, (1991b). 20. Ting, J.M. "A Robust Algorithm For Ellipse-Based Discrete Element Modelling of Granular Materials," to appear in Computers and Geotechnies, (1993). 21. Ng, T.-T. and Lin, X. "Simulations of Fabric Anisotropy on Natural Deposited Granular Soil Using the Discrete Element Method," submitted for possible publication (1993). 22. Lambe, T.W. and Whitman, R.V. Soil MechaMcs, John Wiley & Sons (1969). 23. Seed, H.B., Wong, R.T., Idriss, I.M., and Tokimatsu, K. "Moduli and Damping Factors for Dynamic Analyses of Cohesionless Soils," Ret~ort UCB/EERC-84/14, Univ. of Calif., Berkeley, CA, (1984). 24. Finn, W.D.L. and Bhatia, S.K. "Endochronic Theory of Sand Liquefaction," Proc. 7th World Conference on Earthauake Engineering, Istanbul, Turkey, (1980), 149-156. 25. Li, X.S., Chan, C.K. and Shen, C.K. "An Automated Triaxial Testing System," Advanced Triaxial Testin~ of Soil and Rock. ASTM, STP 977, (1988) 95-106.

169

26. Ishibashi, I. (1988) Personal Communication. 27. Vucvtic, M. and Dobry, R. "Cyclic Tfiaxial Strain-Controlled Testing of Liquefiable Sands," Advanced Triaxial Testin¢ of Soil and Rock, ASTM, STP 977, (1988) 475-485. 28. Finn, W.D.L., Picketing, D.J. and Bransby, P.L. "Stress Strain Relation for Sand in Simple Shear," Seismic Problems in GeotechnicalEngineering, ASCE National Convention, Denver, Co, (1971). 29. Cho, Y., Rizzo, P.C. and Humphries, W.K. "Saturated Sand and Cyclic Dynamic Tests," Liquefaction Problems in Geotechnical Engineering, ASCE Annual Convention and Exposition, Philadelphia, PA, September, (1976) 285-312. Received 5 August 1992; revised version received 20 December 1992; accepted 2 July 1993