Large strain finite element analysis of sand: model, algorithm and application to numerical simulation of tire–sand interaction

Computers and Structures 74 (2000) 253±265

www.elsevier.com/locate/compstruc

Large strain ®nite element analysis of sand: model, algorithm and application to numerical simulation of tire± sand interaction C.H. Liu a,*, 1, J.Y. Wong b, H.A. Mang c a MARC Analysis Research Corporation, 260 Sheridan Avenue, Suite 309, Palo Alto, CA 94306, USA Transport Technology Research Laboratory, Department of Mechanical and Aerospace Engineering, Carleton University, Ottawa K1S 5B6, Canada c Institute for Strength of Materials, Technical University of Vienna, A-1040 Vienna, Austria

b

Received 29 January 1998; accepted 23 January 1999

Abstract Based on a review of several existing models for the elastic responses of geomaterials, a new nonlinear elastic law is suggested. A fully implicit return mapping algorithm for non-associative elastoplastic models of the geomaterials is presented within the regime of large strains. This algorithm is formulated in the principal axes. Algorithmic tangent moduli are derived. A modi®ed critical state model in conjunction with the new nonlinear elastic law is implemented into the general purpose ®nite element software MARC, using the proposed algorithm. Furthermore, the tire±sand interaction is numerically simulated, employing the modi®ed model and the algorithm. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Finite element; Large strain; Plasticity; Sand

1. Introduction The motivation of this work is to develop a computational method for predicting traction performance of tires moving on sand. Compared to existing analytical and empirical methods (see [1] and the references therein), computational methods, such as the ®nite element method, permit modeling machine±terrain interaction in a very detailed manner without introducing many simplifying assumptions. With the widespread

* Corresponding author. Tel.: +1-650-329-6800; fax: +1650-323-5892. 1 Formerly, Transport Technology Research Laboratory, Department of Mechanical and Aerospace Engineering, Carleton University, Ottawa K1S 5B6, Canada.

availability of high speed computing facilities and of powerful engineering analysis software, coupled with recent advances in soil mechanics and computational methods, computational terramechanics has attracted a great deal of interest as a means for improving the design of construction equipment and o€-road vehicles [2±9] and for reaching a better understanding of the soil compaction process [10±14]. The early work on numerical simulations of machine±terrain interaction can be traced to Perumpral et al. [2] in 1971. Idealizations such as elasticity and isotropy of the soil, and a prescribed contact area were used to simplify the analysis. In their work, the stress distributions on the tire±soil interface, obtained experimentally by Onafeko and Reece [15], were used as input data to compute soil response. Yong et al. [3,4]

0045-7949/00/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 0 4 9 - 8

254

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

adopted a nonlinear elastic constitutive relation for soil and considered both rigid and deformable wheels. To study the deformation in the soil under the action of a tire, either load or displacement boundary conditions must be speci®ed on the wheel±soil interface at the outset. Consequently, the analysis by Yong et al. [3,4] appears to be suitable only for predicting the soil deformation (compaction) caused by the wheel load, if the boundary conditions on the wheel±soil interface are known. In recent years, research on numerical simulations of soil±machine interaction has advanced considerably by describing the soil as an elastoplastic medium and by using improved contact algorithm [5± 14]. The work by Aubel et al. [6±9] on the rolling process of wheels on soft soil is characterized by the pressure-dependent Drucker±Prager plasticity model for soil taking into account its frictional behavior, and the contact algorithm available in the general purpose ®nite element program ABAQUS. A weakness of this method is that the Drucker±Prager model is not able to represent soil compaction behavior when the soil is subjected to hydrostatic pressure. Moreover, there is no experimental veri®cation for the numerically obtained predictions. Applications of critical state soil mechanics to the study of soil deformation under the action of agricultural implements and vehicles have also been reported recently in the literature [10,11]. The mechanical behavior of geomaterials strongly depends on their physical conditions, such as on their saturation and density, and on the process the material being formed. Because of the high complexity and variability, a large number of constitutive models for geomaterials have been reported (see, e.g., [16,17] and references therein), depending on di€erent mathematical idealizations. Pressure-dependent elastoplastic models are considered to be the most appropriate representations of the behavior of geomaterials because of the existence of inelastic deformations of the materials subjected to shear stress or/and hydrostatic pressure. Among these models, the modi®ed Cam-clay model [18] and the cap model [19] are two of the most commonly used models in geotechnical engineering. The deformation within an elastoplastic model can be divided into an elastic and an inelastic part. Both are nonlinear in nature. Even though the elastic deformation is generally small compared to the inelastic components, it is normally not negligible, particularly for unloading/reloading process. The choice of di€erent elastic models will a€ect not only the accuracy of the results, but also the stability of algorithms. Various elastic models for geomaterials have been presented (e.g., [20±25]). Several of these models will be discussed in this paper. After a model is established, its implementation becomes the critical issue in order to achieve a reasonably good stability and accuracy in ®nite element

analysis. Because of the mathematical complexities of models of geomaterials, such as the nonlinear elastic law and the pressure-dependent hardening/softening law, resulting from the complex nature of geomaterials, implementation of such models is not trivial. Research along this direction has been reported in recent years (e.g., [26±28]). It appears that most of the implementations were performed for the case of small deformation and/or by means of explicit integration of the constitutive equations of elastoplasticity. Only very recently, applications of algorithmic framework for ®nite deformation plasticity [29,30] to the Cam-clay model were reported by Simo and Meschke [31] and Meschke et al. [32]. As far as the speci®c problem of tire±sand interaction is concerned, the mechanical behavior of sand and its deformation characteristics, such as plasticity and localized failure resulting from normal and shear stresses on the vehicle±sand interface, govern the traction mechanism of a moving vehicle. Therefore, it is important to have a reasonable mechanical representation of sand. Moreover, since large sand deformations occur, it is necessary to formulate and implement the sand model within the regime of large strains. Unfortunately, all studies on machine±terrain interaction mentioned above ([2±14]) did not take the e€ect of large deformations into account. It is the objective of this paper to report our research on the modeling and the associated algorithm for sand and on applications to simulation of tire± sand interaction. Nonlinear elastic constitutive equations for elastoplastic models of geomaterials will be addressed. A general approach for implementing the elastoplastic models of geomaterials within the regime of ®nite deformation, including a fully implicit integration approach and the algorithmic tangent moduli, will be presented. A modi®ed critical state model in conjunction with a new nonlinear elastic law has been implemented into the general purpose ®nite element program MARC, within the proposed algorithmic framework. Material parameters for the sand will be calibrated by using experimental data from hydrostatic and triaxial tests. The applicability of the model to represent sand behavior will be demonstrated by means of numerical analyses of the tests at di€erent load levels. Numerical results obtained from simulations of the interaction between the wheel and the sand, including the traction performance of the wheeldrawbar pull at di€erent slips, stress distributions in the sand and on the wheel±sand interface, and sand compaction, will be presented and compared to available experimental data. The physical problem simulated in this paper is the one of tires moving on dry sand. Consequently, pore pressure and time-dependent phenomena of geomaterials are not considered here.

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

2. Remarks on the modeling elastic behavior The critical state soil model proposed by Roscoe and Burland [18] employs a nonlinear relationship between the hydrostatic pressure p and the elastic volumetric strain Eev , resulting from the assumption of a linear relationship between the speci®c volume v (ˆ 1 ‡ e,e is the void ratio) and the logarithmic hydrostatic pressure, i.e., v ˆ N ÿ k ln… ÿ p †,

…1†

where N is a material parameter de®ned as the speci®c volume at unit compressive pressure, k is the swelling/ recompression index. The elastic bulk modulus K is obtained as Kˆ

dp dp v ˆ ÿ p: ˆ e dEv dv=v k

…2†

The elastic shear modulus G is either chosen as a constant or calculated as 3…1 ÿ 2n † Gˆ K 2…1 ‡ n †

…3†

by assuming a constant value of Poisson's ratio n. The well-known linear relationship (1) has been widely incorporated into elastoplastic constitutive equations of geomaterials. However, a number of shortcomings of this de®nition was reported by Butter®eld [20] and Hashiguchi [21]. Undesirable characteristics include physically meaningless results at very large stress levels and unsatisfactory predictions for many types of soils, among others. From computational point of view, the physically meaningless characteristics may result in numerical problems [23,32]. Moreover, for the case of large deformations, a consistent algorithmic treatment of this nonlinear elastic law is dicult [21,23]. To overcome the diculties associated with Eq. (1), Butter®eld [20] and Hashiguchi [21] suggested a linear relationship between ln v and ln…ÿp†, i.e., ln v ˆ Q ÿ r ln… ÿ p †,

…4†

where Q and r are material parameters. Eq. (4) leads to p Kˆÿ : r

…5†

To complete the de®nition of the isotropic elastic law, either a constant value of G or of n can be assumed. There is no guarantee that the two elastic models (Eqs. (1) and (4)) preserve the energy for any closedloop stress or strain path. In fact, if a constant value of n is used, the value of G is a function of K, and in

255

turn, a function of p. The principle of conservation of energy is generally not satis®ed [22]. As an alternative to the linear relationships between v and ln…ÿp† or between ln v and ln…ÿp†, a linear relationship between p and ln v can be used. It is interesting to observe that di€erentiation of this linear relationship will lead to a constant value of the bulk modulus and, therefore, to a constant value of the shear modulus when using the same procedure for determining G as described previously. Constant elastic parameters is simple and also commonly used in computational geomechanics. However, they are lacking not only physical background but also mathematical consistency once large deformations occur [33]. The principle of conservation of energy requires that, mechanical energy is neither generated or dissipated in a closed-loop stress or strain path. Based on this fact, Lade and Nelson [22] derived the following expression for Young's modulus: " E ˆ Mpa

I1 pa

2

1 ‡ n J2 ‡6 1 ÿ 2n p2a

#l ,

…6†

where I1 is the ®rst stress invariant, J2 is the second deviatoric stress invariant, M and l are material parameters, and pa is the atmospheric pressure. A comparison of results obtained by this expression with experimental results from triaxial tests on sand for various stress paths demonstrated the ability of Eq. (6) to capture the nonlinear elastic behavior of geomaterials [22,34]. In a previous paper by the ®rst two authors [24], a nonlinear elastic law  1=1ÿl1 p ˆ m1 Eev ,

 1=1ÿl2 k s kˆ 2m2 k ee k

…7†

was suggested. In Eq. (7), k s k is the norm of deviatoric stress tensor s; k ee k is the norm of elastic deviatoric strain tensor ee ; m1 , l1 , m2 and l2 are material parameters. It may be interpreted as a hyperelastic constitutive relation, which automatically satis®es the principle of conservation of energy [22]. The explicit form of the strain energy function can be obtained by integrating Eq. (7). Di€erentiation of Eq. (7) yields   m1 l 1 dp ˆ K dEev ˆ dEev , p 1 ÿ l1 …8†   m2 e l2 e d k s kˆ 2G d k e kˆ 2 dke k: ksk 1 ÿ l2 The two sets of material parameters …m1 ,l1 † and …m2 ,l2 † can be determined directly by ®tting numerical results to experimentally obtained unloading/reloading volumetric strain±stress and deviatoric strain±stress curves, respectively.

256

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

Fig. 1. Typical experimental results of unloading tests: (a) hydrostatic behavior and (b) shear behavior.

The form of K (or G ) in Eq. (8) is similar to the form of E proposed by Lade and Nelson in Eq. (6), obtained by using the principle of conservation of energy. However, the former provides more ¯exibility. It allows consideration of changes of Poisson's ratio. This is important since the void ratio may change considerably because of large deformations. In this paper, new nonlinear elastic laws will be used. They have the following form:   Ee p ˆ ÿK0 l1 ln 1 ÿ v , l1   k ee k k s kˆ ÿ2G0 l2 ln 1 ÿ : l2

…9†

In Eq. (9), K0 , l1 , G0 and l2 are material parameters. Di€erentiation of Eq. (9) yields the explicit forms of K and G:

dp ˆ K dEev ˆ

K0 dEe 1 ÿ Eev =l1 v

2G0 d k s kˆ 2G d k e kˆ d k ee k : 1ÿ k ee k =l2

…10†

e

Hence, K0 and G0 can be interpreted as the bulk modulus and the shear modulus, de®ned on undeformed con®guration …Eev ˆk ee kˆ 0†, respectively. Eq. (9) is found to be superior to Eq. (7) as regards the ®tting of experimental results. Moreover, as shown in Fig. 1, all material parameters can be interpreted physically. This facilitates the calibration of these material parameters.

3. Algorithm for large strain plasticity 3.1. Multiplicative elastoplasticity and exponential return mapping algorithm This subsection contains a concise outline of an algorithmic framework for ®nite deformation plasticity which serves as the basis for implementing constitutive models for geomaterials at large strains. For a more elaborate description of the underlying theory, see Ref. [29]. The elastic domain is assumed to be de®ned by f …tt,q† < 0. f …tt ,q† ˆ 0 determines a convex yield surface in the stress space. This function has the same form as the one used in the case of small deformation theory. t is the Kirchho€ stress tensor in the current con®guration. q is a stress-like internal variable characterizing plastic hardening of materials. For simplicity, it is assumed to be a scalar function. In contrast to the additive decomposition of the strain tensor in the theory of elastoplasticity for small strains, in the ®nite deformation regime the deformation gradient F is decomposed into an elastic part Fe and a plastic part Fp in a multiplicative manner: F ˆ Fe Fp . The stress response is obtained by virtue of the hyper elastic constitutive relation, characterized by a free energy function c ˆ c…be ,x†, where be ˆ Fe FeT is the elastic part of the left Cauchy-Green deformation tensor and x is a strain-like internal variable, conjugate to q, i.e., q ˆ ÿ@ c=@ x. Exploiting the principle of maximum plastic dissipation, the local evolution equations, i.e., the associative plastic ¯ow rule and the hardening law are obtained as [29] @f e bÇ ˆ lbe ‡ be lT ÿ 2g be , @t

…11†

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

@f x_ ˆ g , @q

…12†

where g is the plastic consistency parameter and ÿ1 l ˆ FÇ F is the spatial velocity gradient. Eqs. (11) and (12) are the basis for extension of the theory of associative plasticity to the one of non-associative plasticity within the large strain regime. For the latter theory, the Eqs. (11) and (12) take the form @g e bÇ ˆ lbe ‡ be lT ÿ 2g be , @t

…13†

@g x_ ˆ g , @q

…14†

x ˆ xn ‡ Dg

@g : @q

257

…20†

b is the principal Kirchho€ stress vector containing bA …A ˆ 1,2,3†. Note that the materials are assumed to be isotropic. Eqs. (19) and (20) and the discrete KuhnTucker conditions: Dgr0,

…21†

f…b b,q†R0,

…22†

Dgf…b b,q† ˆ 0,

…23†

where g ˆ g…tt ,q† is the potential of plastic ¯ow which has an identical form as the respective potential in the small deformation theory. Within the time interval tE‰tn ,tn‡1 Š, an algorithmic approximation to the local evolution Eqs. (13) and (14), characterized by a hyperelastic trial step, followed by return mapping constructed on the basis of an exponential approximation to the ¯ow rule in the updated con®guration, leads to

represent a form identical to the one for small deformation theory. However, depending on the form of c, the evolution equations of the elastoplastic model in stress space are generally not the same as the ones for the small deformation theory. The algorithmic tangent moduli, resulting from linearization consistent with the return mapping algorithm, are obtained as

  @ g…t ,q † etr e b ˆ exp ÿ 2Dg b , @t

cˆ

x ˆ xn ‡ Dg

@ g…t ,q † , @q

…15†

…16†

where Dt ˆ t ÿ tn ,Dg ˆ gDt and betr ˆ fben f T , with f as the relative deformation gradient de®ned in the deformed con®guration at tn . For simplicity, the subscript t has been omitted in this paper. Making use of the spectral decomposition of be and t , i.e., be ˆ

3 X ÿ e
2 …A † lA m ,

…17†

Aˆ1

tˆ

3 X bA m…A † ,

…18†

Aˆ1

where m…A† ˆ n…A† n…A† , with n…A† denoting the principal directions of be in the current con®guration and leA and bA …A ˆ 1,2,3† are the elastic part of the principal stretches and the principal Kirchho€ stresses, respectively, Eqs. (15) and (16) may be recast in terms of the logarithmic principal strains EeA ˆ ln …leA † as [29] E e ˆ E etr ÿ Dg

@g , @b

…19†

3 X 3 X Aˆ1 Bˆ1

‡2

tr…A † aep

mtr…B † AB m

3 X bA ctr…A † ,

…24†

Aˆ1

etr tr…A† where aep , which AB ˆ @ bA =@ EB . The explicit form of c is independent of the speci®c plasticity model, was given by Simo and Taylor [35] for ®nite deformation elasticity.

3.2. Implicit integration An outline of the implicit integration algorithm for solving the elastoplastic evolution Eqs. (19)±(23) is given below. Step 1. Initialize the state variables at k ˆ 0, where k is the counter of the local iteration steps in the return mapping process, by setting the elastic strains E e…0† ˆ E etr , the hardening variable x…0† ˆ xn and the plastic consistency parameter Dg…0† ˆ 0. Step 2. Evaluate the stress b …k† , the hardening parameter q…k† by means of the hyperelastic constitutive relation, i.e., … †

bk ˆ

… †

@c k , @ Ee

…25†

258

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

q…k † ˆ ÿ

… †

@c k : @x

…26†

For the speci®c elastic model de®ned by the Eqs. (9) and (10), the stresses are calculated using Eq. (9). Step 3. Evaluate 8 …k † > > > b …k † ÿ @ c > < @ Ee r…1k † ˆ > > …k † @ c…k † > > :q ‡ @x r…2k † ˆ

f

…k †



Step 6. Update the variables by means of (



9 > > > > = , > > > > ;

…27†

 …k † E e …k † ÿ E etr ‡ Dg…k † h , … † x k ÿ xn

  … † ˆ f b k ,q…k † ,

…28†

…29†

… †

db bk dq…k †

dEE e…k † … † dx k

)

h
…k † i ÿ ˆ ÿB…3k † r2…k † ‡ B…1k † r1…k † ‡ d Dg…k † h , "(



B…1k †

ˆ

… †

db bk dq…k †

)

…35†

# ‡

r…1k †

,

…36†

i.e.,
ÿ Dg…k‡1† ˆ Dg…k † ‡ d Dg…k † ,

…37†

  … † … † b …k‡1† ˆ b k ‡ d b k ,

…38†

  q…k‡1† ˆ q…k † ‡ d q…k † ,

…39†

E e…k‡1† ˆ E e…k † ‡ dEEe…k † ,

…40†

… †

…41†

…k†

with h ˆ ‰@ g…k† =@ b , ÿ @ g…k† =@ qŠT . Step 4. Check the convergence according to the following criteria: k r…k† k r2…k† k RTOL2 and 1 k RTOL1, …k† f RTOL3. The iterations will be terminated if convergence is achieved. Step 5. Perform return mapping iteration. For a ®nite incremental step of d…Dg…k† †, linearization of the Eqs. (27)±(29), and setting r1…k‡1† ˆ 0,r2…k‡1† ˆ 0 and f …k‡1† ˆ 0, leads to [23]   ÿ … †
f …k † ÿ h…k † B3…k † r2…k † ‡ B…1k † r1…k † d Dg k ˆ , …30† …k † h…k † B…3k † h where h h…k † ˆ @ f

B3…k †

=@ b ,@ f

…k †

i =@ q ,

…31†

iÿ1 h ˆ B…1k † ‡ Dg…k † B2…k † ,

B1…k † ˆ

B2…k †

…k †

8 … † … † > @ 2c k @ 2c k > > > < @ Ee@ Ee @ Ee@ x

…32† 9ÿ1 > > > > =

… † … † > > > @ 2c k @ 2c k > > > > > ÿ :ÿ ; e @ x@ E @ x@ x

8 > @ 2 g…k † @ 2 g…k † > > > < @ b@ b @ b@ q ˆ > @ 2 g…k † @ 2 g…k † > > > ÿ :ÿ @ q@ b @ q@ q

9 > > > > = : > > > > ;

Set k ˆ k ‡ 1 and proceed to Step 3.

3.3. Algorithmic tangent moduli The algorithmic tangent moduli are obtained by means of consistent linearization of the stress±strain relations Eqs. (25) and (26), the ¯ow rule Eq. (19) and of the hardening law Eq. (20), respectively. Taking the total di€erential of Ee and x according to Eqs. (19) and (20), respectively, yields 

dEE e dx



 ˆ

dEEetr 0



 ÿ d…Dg†h ÿ DgB2

 db b : dq

…42†

Substituting the constitutive relations (resulting from Eqs. (25) and (26)) 

,

… †

x…k‡1† ˆ x k ‡ dx k :

…33†

dEE e dx



 ˆ B1

db b dq

 …43†

into Eq. (42) leads to  …34†

db b dq



 ˆ B3

dEE etr 0



 ÿ d…Dg†h :

…44†

Making use of Eq. (44) and of the consistency condition df ˆ 0 yields

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

259

Fig. 2. Critical state model: meridional section and critical state line.

d…Dg† ˆ

  hB3 dEEetr : hB3 h 0

Inserting Eq. (45) into Eq. (44) results in    etr  db b dEE , ˆ Aep dq 0

…45†

…46† Fig. 3. Hydrostatic test of the Ottawa sand: model predictions and experimental results.

with Aep ˆ B3 ÿ

B3 h h B3 : hB3 h

…47†

The resulting algorithmic elastoplastic tangent operator aep ‰33Š , de®ned in principal axes, relating the di€erentials of the principal Kirchho€ stresses to the di€erenep tials of the principal strains is obtained as aep ij ˆ Aij , where i ˆ 1,2,3 and j ˆ 1,2,3. In consequence of the non-associative plastic ¯ow rule, aep is generally nonsymmetric. Finally, the algorithmic moduli c are obtained from Eq. (24). 4. Numerical simulations of tire±sand interaction 4.1. Critical state model A modi®ed critical state model in conjunction with the proposed nonlinear elastic law (i.e., Eqs. (9) and (10)) will be used to represent the mechanical behavior of sand in numerical simulations of the interaction of the tire with the sand. In contrast to the original formulation proposed by Roscoe and Burland [18], the ellipsoidal yield surface (see Fig. 2) for the critical state model is described as  s 2  k s k2 qÿT f…p, k s k ,q† ˆ ‡ p‡ 2 M2   q‡T ÿ R0: …48† 2 Eq. (48) is a convex function which guarantees the numerical stability of the algorithm. In Eq. (48), M

denotes the slope of the critical state line in the k ÿs k ÿp plane; the hardening parameter q represents the consolidation pressure; and T denotes the hydrostatic strength of the sand. The evolution of the plastic strains is governed by the associative ¯ow rule. Because of the nature of sand materials, the hardening, however, is characterized by a non-associative rule, i.e., q ˆ H…x†

with x ˆ Epv ,

…49†

where Epv is the plastic part of volumetric strain. The explicit form of H…x† depends on the speci®c material considered. It can be determined by means of experimental volumetric stress±strain curves. In this work, the following form is used:   x q ˆ H…x† ˆ q0 ÿ D ln 1 ÿ , …50† W where q0 is the initial consolidation pressure; D and W are hardening parameters. 4.2. Calibration of material parameters The material parameters for Ottawa Sand were calibrated by experimental data, obtained from hydrostatic and triaxial tests at 100, 150 and 200 kPa con®ning pressures. They are K0 ˆ 11000 kPa, l1 ˆ 0:005, G0 ˆ 9750 kPa, l2 ˆ 0:050, M ˆ 1:130, D ˆ 1060 kPa, W ˆ 0:046, and q0 ˆ qres ˆ T ˆ 22:246 kPa. qres denotes a lower limit for q for the case of pronounced softening.

260

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

and the experimental results further con®rms the usefulness of the proposed modi®ed critical state model for representation of the mechanical behavior of Ottawa sand within the considered range of con®ning pressure. 4.3. Numerical simulations of tire±sand interaction

Fig. 4. Triaxial test of the Ottawa sand: model predictions and experimental results.

Fig. 3 shows the results from two hydrostatic tests and the numerical predictions. The proposed critical state soil model for representation of the mechanical behavior of Ottawa sand under isotropic compression is quite satisfactory. In particular, the very good agreement between experimentally obtained unloading/ reloading behavior and model predictions demonstrates the usefulness of the proposed nonlinear elastic constitutive equations. Please note that the hyperelastic model does not capture the hysteresis observed in the experiments during the unloading/reloading process. A comparison of the experimental data of triaxial tests of the Ottawa sand with model predictions is shown in Fig. 4. For each con®ning pressures, the test was repeated twice to ensure the reliability of the obtained results. The agreement between the numerical

The traction performance of two wheels, with di€erent dimensions, on two di€erent types of sand was simulated using the modi®ed critical state soil model described previously and the ®nite element program MARC. To evaluate the capability of the proposed numerical method, the predicted performance of the two wheels was compared with two sets of available experimental data. One set of experimental data was obtained by Jastremski [36] with a pneumatic tire of 12.5/75R20 at a high in¯ation pressure on the Ottawa sand. The other set of experimental data was obtained by Onafeko and Reece [15] with a rigid wheel of diameter of 1245 mm and a width of 305 mm on loose sand. At an in¯ation pressure of 345 kPa, the pneumatic tire of 12.5/75R20, with a normal load of 9.23 kN, exhibited a relatively small de¯ection. Consequently, in the simulation it was represented by a rigid wheel. As the ratio of the tire width to sinkage is relatively large on the Ottawa sand, the problem was treated as two dimensional. The diameter of the tire is 977 mm. In the simulation of the tire±sand interaction, the sand is represented by a block of 10000 mm in length and 2500 mm in depth. It is discretized by 1067 bilinear elements. Plane strain conditions are assumed. Horizontal and vertical boundary conditions are applied at both sides and at the bottom of the block, respectively. The ®nite element mesh is shown in Fig. 5. At ®rst, a vertical load of 9.23 kN is applied to the rigid wheel representing the tire 12.5/75R20. Then, the

Fig. 5. Finite element mesh for the simulations of the wheel±sand interaction.

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

Fig. 6. Numerically predicted drawbar pull as a function of horizontal moving distance of a tire 12.5/75R20 from standstill, at an in¯ation pressure of 345 kPa and a normal load of 9.23 kN, and at four di€erent slips on the Ottawa sand.

wheel is moved horizontally at given slips. Assume that the horizontal moving velocity, the angular velocity and the radius of the wheel are V, o and r, respectively. The wheel slips is de®ned as …or ÿ V †=or. The friction law is de®ned as t ˆ mN…1 ÿ eÿj=k †,

…51†

where N is the normal pressure, t is the shear stress, j is the relative slipping distance between the wheel and the sand on their interface, and m and k are parameters chosen, based on experimental data, as 0.32 and 38.0 mm, respectively. Fig. 6 shows numerically predicted drawbar pull as

Fig. 7. Comparison of predicted and measured traction performance of a tire 12.5/75R20, at an in¯ation pressure of 345 kPa and a normal load of 9.23 kN, on the Ottawa sand (measured data from [36]).

261

Fig. 8. Numerically predicted drawbar pull as a function of horizontal moving distance from standstill, for a rigid wheel with a diameter of 1245 mm and a normal load of 9.28 kN, at ®ve di€erent slips on a loose sand.

a function of the moving distance of the wheel from standstill, at four di€erent slips, namely, at 5, 20, 35 and 50%. It is seen that there is a considerable transient period before steady-state operating conditions are reached. A comparison of the predicted and the measured traction performance under steady-state operating conditions is illustrated in Fig. 7. Compared to experimental results [36], the predicted drawbar pull appears to be lower at lower ranges of slip, and higher at higher ranges of slip. To further demonstrate the capability of the proposed method based on modi®ed critical state soil mechanics and on the ®nite element method, the traction performance of a rigid wheel, with a diameter of 1245 mm and a width of 305 mm, on loose sand was simulated and compared with experimental data reported in [15]. In addition to the drawbar performance, experimentally obtained stress distributions along the wheel-sand interface were compared with corresponding numerical results. In the numerical simulation, a vertical load of 9.28 kN was applied to the rigid wheel. Then, the wheel is moved horizontally at ®ve slips, namely, 3.1, 7.1, 12.1, 17.1 and 22.1%. The sand was represented by the ®nite element mesh shown in Fig. 5. The boundary conditions were the same as mentioned previously. Material parameters for the loose sand used in the experiments were not given in [15]. In the simulation, the material parameters for the sand were selected to achieve approximately the same amount of sinkage as in the experiments. The parameters used were K0 ˆ 11000 kPa, l1 ˆ 0:005, G0 ˆ 9750 kPa, l2 ˆ 0:050, M ˆ 1:130,D ˆ 1,060 kPa, W ˆ 0:260, and q0 ˆ qres ˆ

262

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

Fig. 9. Comparison of predicted and measured traction performance of a rigid wheel, with a diameter of 1245 mm and a normal load of 9.28 kN on a loose sand (measured data from [15]).

T ˆ 22:246 kPa. Compared to the respective parameters of the Ottawa Sand mentioned previously, this set of material parameters results in a higher degree of compressibility, which is typical for the mechanical behavior of loose sand. The friction law along the wheel-sand interface is given by Eq. (51), with m ˆ 0:57 and k ˆ 32 mm. Fig. 8 shows the numerically predicted drawbar pull as a function of the horizontal moving distance of the rigid wheel from standstill, at ®ve di€erent slips. Again, prior to reaching a steady-state operating condition, there is a considerable transient period. A comparison of the numerically obtained drawbar pull±slip relation with experimental data reported in [15] is shown in Fig. 9. Close agreement is observed. Both the predicted and the measured distributions of the contact pressure and the friction force along the wheel±sand interface, at 3.1 and 22.1% slips, are shown in Fig. 10. In contrast to experimental observations, the predicted sinkage for di€erent slips does not di€er very much. Hence, the contact pressure for both 3.1 and 22.1% slips are close. The friction force at 22.1% slip is, however, considerably higher than the one at 3.1% slip. This follows from the longer slipping distance between the wheel and the sand on the interface. The di€erence between the predicted and the measured friction force shows that the commonly used friction law de®ned by Eq. (51), resulting from shear tests, does not adequately describe the rolling±slipping process. Fig. 11 shows the deformed pattern of the sand in the vicinity of the moving rigid wheel and the distribution of the vertical stress at 3.1% slip. It can be seen that residual

Fig. 10. Comparison of the predicted and measured stress distributions on the rim of a rigid wheel on loose sand: (a) at 3.1% slip; (b) at 22.1% slip (measured data from [15]).

stresses in the sand remain after the wheel passes by. The predicted maximum vertical pressure on the wheel-sand interface is around 130 kPa. The soil compaction under the moving wheel at 3.1% slip is depicted in Fig. 12. It is shown that the volumetric plastic strain, which indicates the compaction resulting from the loading of the moving wheel, decreases with the increase of depth. The major compaction occurs in the layer immediately under the sand surface.

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

263

Fig. 11. Deformed pattern of the sand in the vicinity of a rigid wheel, with a diameter of 1245 mm and a normal load of 9.28 kN; distribution of vertical stress (kPa) at 3.1% slip.

5. Summary and conclusions A general approach for implementing non-associative constitutive models for geomaterials within an ecient algorithmic framework for ®nite deformation plasticity has been proposed. It is characterized by a fully implicit integration and algorithmic tangent moduli, resulting in a quadratic rate of convergence in global Newton's iterations. Based on a review of the existing models describing the elastic behavior of geomaterials, a new nonlinear elastic law was suggested. The essential features of this model are satisfaction of the principle of conservation of energy, ¯exibility with

respect to consideration of the change of Poisson's ratio and the physically meaningful interpretations of all model parameters. A modi®ed critical state model in conjunction with the new nonlinear elastic law was implemented into the ®nite element software MARC using the described algorithmic approach. Moreover, numerical simulations of the tire±sand interaction by the proposed model were reported. It is expected that consideration of large strains yields more accurate numerical results. The usefulness of the proposed constitutive model to represent the material behavior of the Ottawa sand was demonstrated experimentally. The numerically obtained results for the

Fig. 12. Deformed pattern of the sand in the vicinity of a rigid wheel, with a diameter of 1245 mm and a normal load of 9.28 kN; distribution of plastic volumetric strain at 3.1% slip.

264

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

tire±sand interaction were compared to experimental data. Reasonably close agreement between the predicted and the measured drawbar performance was observed. Qualitatively agreement between the predicted and the measured stress distributions on the wheel±soil interface was obtained. Because of the complex nature of the problem considered, further improvements of the models are needed. The contact pressure and friction force along the tire±soil interface and, consequently, the traction performance of the tire, mainly depend on the friction law used. A more accurate representation of the shearing behavior for this kind of rolling±slipping process is needed to obtain more realistic results. Moreover, local failure as well as signi®cant amount of ¯ow of sand under moving tires cannot be modeled adequately by the conventional ®nite element method. This is the reason why the numerical results do not show signi®cant di€erence in wheel sinkage at di€erent slips. Further studies along these directions should be carried out to improve the numerical simulations. Finally, consideration of deformable tires and three dimensional cases are also topics of further research.

[7]

[8]

[9]

[10]

[11]

[12]

[13]

Acknowledgements This research is supported by a research grant from the Natural Sciences and Engineering Research Council of Canada (OGP0005590), awarded to Professor J.Y. Wong, Carleton University. The assistance given by Mr. Zeljko Knezevic, Graduate Research Assistant, Transport Technology Research Laboratory, Carleton University, in experimental work is appreciated.

References [1] Wong JY. Terramechanics and o€-road vehicles. New York: Elsevier Science, 1989. [2] Perumpral JV, Liljedal JB, Perlo€ WH. A numerical method for predicting the stress distribution and soil deformation under a tractor wheel. Journal of Terramechanics 1971;8(1):9±22. [3] Yong RN, Fattah EA. Prediction of wheel-soil interaction and performance using ®nite element method. Journal of Terramechanics 1976;13(4):227±40. [4] Yong RN, Fattah EA, Boosinsuk P. Analysis and prediction of tyre-soil interaction and performance using ®nite elements. Journal of Terramechanics 1978;15(1):43±63. [5] Jarzebowski A, Uchiyama K, Hashiguchi K, Ueno M, Nohse Y. Numerical analysis of soil±wheel interaction. In: Proceedings of 11th International Conference of ISTVS, Lake Tahoe, NV, USA. 1993. p. 340±9. [6] Aubel T. FEM simulation of the interaction between

[14]

[15]

[16]

[17]

[18]

[19]

[20] [21]

[22]

[23]

elastic tire and soft soil. In: Proceedings of 11th International Conference of ISTVS, Lake Tahoe, NV, USA. 1993. p. 791±802. Aubel T. The interaction between the rolling tyre and the soft soilÐFEM simulation by VENUS and validation. In: Proceedings of 6th European ISTVS Conference. Austria: Vienna, 1994. p. 169±88. Fervers CW. FE-simulation of tyre-pro®le-e€ects on terrain mobility of vehicles. In: Proceedings of 6th European ISTVS Conference. Austria: Vienna, 1994. p. 618±33. Schmid IC. Interaction of Vehicle and terrain results from 10 years research at IKK. Journal of Terramechanics 1995;32(1):3±26. Chi L, Tessier S. Finite element analysis of soil compaction reduction with high ¯otation tires. In: Proceedings of the 5th North American ISTVS Conference/ Workshop, Saskatoon, SK, Canada. 1995. p. 167±76. Foster WA, Johnson CE, Raper RL, Shoop SA. Soil deformation and stress analysis under a rolling wheel. In: Proceedings of the 5th North American ISTVS Conference/Workshop, Saskatoon, SK, Canada. 1995. p. 194±203. Chi L, Kushwaha RL. Finite element analysis of forces on a plane soil cutting blade. Canadian Agricultural Engineering 1989;31:135±40. Chi L, Tessier S, LagueÈ C. Finite element prediction of soil compaction induced by various running gears. Transaction of the ASAE 1993;36(3):629±36. Chi L, Tessier S, LagueÈ C. Finite element prediction of soil compaction by liquid manure spreaders. Transaction of the ASAE 1993;36(3):637±44. Onafeko O, Reece AR. Soil stress and deformation beneath rigid wheel. Journal of Terramechanics 1967;4(1):59±80. Chen WF, Mizuno E. Nonlinear analysis in soil mechanicsÐtheory and implementation. Amsterdam: Elsevier, 1990. Desai CS, Siriwardane HJ. Constitutive laws for engineering materials with emphasis on geologic materials. Englewood Cli€s, New Jersey: Prentice-Hall, 1984. Roscoe KH, Burland JB. On the generalized stress strain behavior of ``wet'' clay. In: Heyman J, Leckie FA, editors. Engineering Plasticity. Cambridge: Cambridge University Press, 1968. p. 535±609. Di Maggio FL, Sandler IS. Material model for granular soils. Journal of Applied Mechanics Division [special issue]. ASCE 1971;97:935±50. Butter®eld K. A natural compression law for soils (an advance on e-log, p ' ). Geotechnique 1979;29:469±80. Hashiguchi K. On the linear relations of V-ln p and ln, v- ln p for isotropic consolidation of soils. International Journal for Numerical and Analytical Methods in Geomechanics 1995;19:367±76. Lade PV, Nelson RB. Modeling the elastic behavior of granular materials. International Journal for Numerical and Analytical Methods in Geomechanics 1987;11:521± 42. Liu CH. Traction mechanism of automobile tires on snowÐan investigation by means of the ®nite element

C.H. Liu et al. / Computers and Structures 74 (2000) 253±265

[24]

[25] [26]

[27]

[28]

[29]

method. Ph.D thesis, Institute for Strength of Materials, Technical University of Vienna, 1994. Liu CH, Wong JY. Finite element implementation of a critical state soil model at large strains. In: Meguid SA, editor. Mechanics in Design. Engineering Mechanics and Design Laboratory. Toronto, Canada: University of Toronto, 1996. p. 127±36. Liu CH, Wong JY. Numerical simulation of tire±soil interaction based on critical state soil mechanics. Journal of Terramechanics 1997;33:209±21. Sandler IS, Rubin D. An algorithm and a modular subroutine for the cap model. International Journal for Numerical and Analytical Methods in Geomechanics 1979;3:177±86. Borja RI, Lee SR. Cam-clay plasticity. Part I: Implicit integration of elasto-plastic constitutive relations. Computer Methods in Applied Mechanics and Engineering 1990;78:49±72. Borja RI. Cam-clay plasticity. Part II: Implicit integration of constitutive equation based on a nonlinear elastic elastic stress predictor. Computer Methods in Applied Mechanics and Engineering 1991;88:225±40. Simo JC. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping

[30] [31] [32] [33]

[34] [35]

[36]

265

schemes of the in®nitesimal theory. Computer Methods in Applied Mechanics and Engineering 1992;99:61±112. Bathe KJ. Finite element procedures. Englewood Cli€s, New Jersey: Prentice-Hall, 1996. Simo JC, Meschke G. A new class of algorithms for classical plasticity extended to ®nite strains. Application to geomaterials. Computational Mechanics 1993;11:253±78. Meschke G, Liu CH, Mang HA. Large strain ®nite element analysis of snow. Journal of Engineering Mechanics, ASCE 1996;122:591±602. Simo JC, Pister KS. Remarks on constitutive equations for ®nite deformation problems: Computational implications. Computer Methods in Applied Mechanics and Engineering 1984;46:201±15. Panos P, Sun Y. Fine Ottawa Sand: Experimental behavior and theoretical predictions. Journal of Geotechnical Engineering 1992;118:1906±23. Simo JC, Taylor RL. Quasi-incompressible ®nite elasticity in principal stretches. Continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering 1991;85:273±310. Jastremski JA. Further developments in the modeling of the performance of pneumatic tires over deformable terrain. M.Eng. thesis, Carleton University, 1994.