- Email: [email protected]
Finite element simulation of localization in granular materials by micropolar continuum approach
Available online at www.sciencedirect.com
Computers and Geotechnics 35 (2008) 548–562 www.elsevier.com/locate/compgeo
Finite element simulation of localization in granular materials by micropolar continuum approach Haydar Arslan a
a,*
, Stein Sture
b
Exxon Mobil Upstream Research Company, P.O. Box 2189, Greenway Plaza-3, Room 704-A, Houston, TX, United States b University of Colorado, 26 UCB, Boulder, Colorado 80309-0026, United States Received 23 March 2007; received in revised form 13 October 2007; accepted 24 October 2007 Available online 8 January 2008
Abstract There are two different approaches in continuum mechanics: a continuum considered as solid state – classical continuum mechanics – or a continuum with the microstructure. In order to explain the influences of microstructures on the behavior of materials, higher order effects such as moments that reflect the grain geometry and kinematics (e.g. rotations) must be incorporated, which in turn lead to a basic asymmetry of shear stress and micropolar (or Cosserat) effects. Cosserat continuum incorporates the characteristic length that relates the couple stresses and micro-curvatures at the constitutive model. Finite Element solutions of localization analysis have been carried out in Cosserat continuum to investigate length scale-size effect on the overall behavior of granular materials. The formulations of Cosserat elastoplasticity have been implemented into the nonlinear commercial finite element code using the user defined element utility. The mesh dependency effect in FE solution is reduced. The numerical model was used to predict the load–displacement behavior of a shallow footing and failure behavior of a slope. A proper treatment of the localization phenomenon in granular materials has been observed and simulated in Finite Element analysis. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Localization; Length scale; Size effect; Micropolar; Granular materials
1. Introduction Strain localization in metals, composite materials and geomaterials has been studied since the early 1900s. The strain localization phenomenon has received considerable research attention lately, and much of recent activities have focused on developing mathematical and numerical models for capturing or simulating the experimentally observed strain localization phenomena. The theory of strain localization is an analytical framework to analyze shear bands in materials. It was first proposed by Hadamard [1] and later developed by Thomas [2], Hill [3], Mandel [4], and Rice [5]. It was applied to predict the existence and orientation of shear bands within various *
Corresponding author. Tel.: +1 720 231 9239. E-mail addresses: [email protected] (H. Arslan), stein.sture@ colorado.edu (S. Sture). 0266-352X/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2007.10.006
types of materials, including elastoplastic soils and rocks [6–14]. Rudnicki and Rice [6] established that the elastoplastic Drucker–Prager model predicts localization during triaxial compression of pressure-sensitive materials. Vardoulakis [10] extended the bifurcation analysis of Hill [15], Hill and Hutchinson [16] to sands. By using a constitutive model of the deformation theory of plasticity, he analytically derived the empirical relation that Arthur et al. [17] proposed for defining the inclination of shear bands. Desrues [18] extensively and accurately measured the direction of shear bands during plane-strain compression of sand specimens. However, he could not describe or explain his experimental results by using the relationships developed by Arthur et al. [17] and Vardoulakis [10]. Finite Element analysis studies related to the strain localization process using conventional continuum mechanics principles have shown that shear band formation is mesh
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
dependent. The mesh dependent solutions might lead to a vanishing shear band thickness as the mesh size is refined in the finite element solution. Tejchman [19], Tejchman et al. [20], Maier [21], Manzari [22], Voyiadis et al. [23], Alsaleh et al. [24,25], Alshibli et al. [26], Arslan [27] confirmed that the shear band formation is mesh independent in Cosserat continuum analysis due to the regularization of the deformations by additional (rotational) degrees of freedom. 2. Background on localization analysis solution techniques Coulomb [28] first introduced the concept of slip planes in frictional materials. The original concept of Coulomb has since evolved and reached substantial maturity after the development and implementation of the theory of plasticity as it applies to soils. Roscoe [29] developed a slip plane theory for pressure sensitive materials such as soils, where the subsequently termed Roscoe criterion is based on kinematics rather than static or failure conditions. The Mohr–Coulomb criterion, which is only based on statics would explain the orientation of the failure plane but does not provide any information related to geometry or kinematics. The classical Mohr–Coulomb criterion as developed for two dimensions, calculates the orientation of shear band as: hc ¼ ð45 þ /=2Þ
ð1Þ
where / is the internal friction angle and hc is the angle between the major principal direction and shear band (Fig. 1). It should be noted that the equation does not include any variables related to failure, except failure plane orientation based on static concepts. Roscoe [29] introduced a kinematical assumption for the orientation of the shear band as: hR ¼ ð45 þ w=2Þ
Fig. 1. Mohr–Coulomb solution to shear band orientation.
ð2Þ
549
where hR is the angle between the shear band and the direction of the major principle stress and w is the angle of dilatancy, which is calculated by: w¼
e_ v e_ q
ð3Þ
where e_ v is volumetric strain rate and e_ q is the shear strain rate. Arthur et al. [17] proposed an empirical equation for shear band orientation based on plane strain experimental results: 1 1 hA ¼ 45 þ ð/ þ wÞ ¼ ðhC þ hR Þ ð4Þ 4 2 Alshibli [30] have conducted plane strain experiments on different sands and the experimental results verified the Eq. (4) proposed by Arthur et al. [17]. The equation was theoretically confirmed by Vermer [31]. By using plasticity concepts, Mandel [4] derived the critical hardening parameter for a granular material subjected to loading in plane strain. The critical hardening parameter is associated with the possible existence of discontinuous bifurcation modes and possible loss of uniqueness. The orientation of possible band of localized deformation is determined by the values of a critical hardening parameter Hcr. Mandel did not derive expressions for the bifurcation directions. Vermeer [31] rederived Mandel’s equation that includes the orientation of the shear band as follows: tan2 hV ¼
2 þ ðsin / þ sin wÞ 2 ðsin / þ sin wÞ
ð5Þ
As can bee seen, Arthur’s and Vermeer’s expressions are in good agreement. Vermeer derived the critical hardening parameter as: H cr 1 2 ðsin / sin wÞ ¼ 2G 16ð1 tÞ
ð6Þ
Vardoulakis [8] performed theoretical and experimental investigations for the shear band orientation by using the Mohr–Coulomb failure criterion. He found an expression for the shear band orientation similar to those developed by Arthur and Vermeer. Drucker–Prager and Mohr–Coulomb plasticity are fundamentally different failure theories. Formulations and the terms that used for the yield functions are different. However, the terms are based on friction and cohesion and they both predict similar results. Ottosen and Runesson [32] showed that the localization angle coincides with the failure angle of the yield condition for the Mohr–Coulomb yield criterion based on associatedplastic flow concepts. The main disadvantage of the many different classical Mohr–Coulomb criteria is that the failure plane is independent of the intermediate principal stress. Furthermore, it does not explain the material behavior after the material failure initiation. The Drucker–Prager yield condition, which takes into account all three principal stresses generalizes the
550
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Mohr–Coulomb failure condition. Thus, the stress dependent localization analysis may be investigated by this model. However, the Drucker–Prager model does not provide failure mode information. Hill [15], introduced material bifurcation concepts and processes to plasticity. Later on, Rudnicki and Rice [6], Ranecki and Bruhns [33], Runesson and Mroz [35], Rynesson et al. [34], Ottosen and Runesson [36] considered the localization as a discontinuous bifurcation problem. The bifurcation analysis considers a necessary condition for discontinuous bifurcation, which is expressed as a singularity of the localization tensor which is similar to the elastic acoustic tensor [55]. The localization condition signals the potential formation of a strain discontinuity. From a mathematical point of view, these complicating features are related to the socalled loss of ellipticity of the governing differential equation. The boundary value problem becomes ill-posed, i.e., it does not have a unique solution with continuous dependence on the given data. Thus, a discontinuity approach is necessary to develop a continuum formulation that follows
the failure behavior correctly. The readers are referred to chapter three for detail explanations. The initial thickness of the localization band depends on the material microstructure. Mathematically, narrow zones of highly concentrated evolving strain can be represented in many different ways. There are mainly three types of discontinuity in continuum mechanics as shown in Figs. 2 and 3. The first type of discontinuity incorporates strong discontinuities, i.e., jumps in displacements across a discontinuity curve (in two dimensions) or discontinuity surface (in three dimensions). In the simplest case, the strong discontinuity can be considered as stress-free, which necessarily leads to a stress singularity at the discontinuity line. In physical terms, the strong discontinuity corresponds to a sharp crack. The strong discontinuity approach has been used as a solution technique for quasi-static solid mechanics problems. It exhibits jumps in the displacement field across a line in two dimensional problems. The line can be named as slip line or discontinuity line. Chakrabarty [38] introduced the strong discontinuity approach to a von
Fig. 2. Types of discontinuity [37].
Fig. 3. Classification of discontinuities for localization.
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Mises-based plasticity model to explain the observed shear bands in metals. Later on, Simo et al. [39], Oliver [40,41], Larson et al. [42], and Steinmann [43] continued to use the concept of a strong discontinuity as a tool for capturing localization in solid mechanics. Another possible kinematic description represents the region of localized deformation by a band of a small but of finite thickness, separated from the remaining part of the body by two weak discontinuities. Since the displacement is continuous, the strain components in the plane tangential to the discontinuity surface must remain continuous as well, and only the out-of-plane components can have a jump. In physical terms, the band between the weak discontinuities corresponds to a damage process zone with an almost constant density of microdefects. Instead of lumping all the inelastic effects into a surface, it is possible to distribute them uniformly across the width of a band of a finite thickness. As Fig. 3 illustrates, a weak discontinuity is characterized by discontinuous but bounded strain fields. Since the displacement is continuous, the strain components in the plane tangential to the discontinuity surface must remain continuous as well, and only the out-of-plane components can exhibit a jump. Runesson et al. [34], Peric [14], Willam et al. [44], Iordache and Willam [45], and Arslan et al. [55] have used the weak discontinuity concept for bifurcation solutions in elasto-plastic materials. Finally, the most regular description uses a continuously differentiable displacement field, and the strain field remains continuous. Strain localization is manifested by high strains in a narrow band, with a continuous transition to much lower strain levels in the surrounding parts of the body. In physical terms, this corresponds to a damage process zone with a continuously varying concentration of defects. As the experimental and numerical studies have shown, the shear band has a finite thickness inside granular materials. The formation of weak discontinuities characterized by continuous displacements but discontinuous strains, which concentrate into a band of finite width, thereby this approach is more suitable to capture and calculate the thickness of the shear band. The weak discontinuity has been used in this study to model and capture localization. It must be recognized that the predictions of localization conditions is not only dependent on the chosen discontinuity approach but also to the constitutive description. The underlying physical definition of the behavior of granular materials must be considered for a well representative formulation. When the applied stress reaches the reduced strength, softening starts and the stress decreases. Consequently, the material outside the weaker region must unload elastically, because its strength limit has not been reached. This leads to the conclusion that the size of the softening region cannot exceed the size of the region having minimum strength. In the present one-dimensional setting, loss of ellipticity occurs, when the tangent modulus ceases to be
551
positive. From the numerical point of view, ill-posedness is manifested by pathological sensitivity of the results to the size of finite elements. The limit solution represents a displacement jump at the center, with zero strains everywhere else. Plasticity models have been used to capture discontinuities. As Jirasek [37] explained in detail, the strong discontinuity can be formulated as a traction-separation law, which complements the stress–strain law valid for the continuous part of the body. Instead, models with localization bands bounded by weak discontinuities can be considered as simple regularizations of models with strong discontinuities. Full regularization of the localization problem can be achieved by a proper generalization of the continuum theory. Various techniques have been implemented to provide a physically acceptable solution. Some of these techniques impose restrictions on the constitutive moduli in the postlocalization regime, while others artificially restrict the size of the finite element relative to the localization zone. The former is based on enriching the continuum with non-conventional constitutive relations in such a way that an internal or characteristic length scale is introduced. Non-local theories of this kind include the Cosserat continuum theory [46], higher gradient theories, and the integral or gradient theory. Localization analysis in Cosserat continua has typically focused on analytical and finite element simulation [20– 22,26,47–49,55]. Main research interests have been to illustrate the applicability of the higher order micro-continuum theory to elastic–plastic solid mechanics analysis. The progress in implementation of Cosserat continuum in material modeling has led to more accurate and realistic constitutive relations for granular materials. Instead of going into detailed formulation, the nature and progress in Cosserat continuum approach will be summarized. Additional material parameters in Micropolar continuum and different solution approaches by different researchers will be reviewed to clarify the main underlying concepts of the Cosserat continuum. 3. Method of analysis The Cosserat continuum belongs to the larger class of generalized continua which introduce intrinsic length scales in continuum mechanics via higher order gradients, additional degrees of freedom in fully non-local constitutive equations [50–53,58]. Eringen [52] defined a Cosserat medium as a continuous collection of particles that behave like rigid bodies. Accordingly, translational and rotational degrees of freedom describe displacement and rotation of an underlying microstructure. Cosserat theory incorporates a local rotation of points (Fig. 4) as well as the translation assumed in classical elasticity. N and n are the socalled directors that show the orientations of particles in the reference and the actual configuration, respectively. The rotation of the joint from
552
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
moments created by couple stresses. The infinitesimal strain tensor encompasses the displacement gradient and micro-rotation. Thus, the strain tensor can be written as: e ¼ ru e x
ð9Þ
The micro-curvature strain tensor that is conjugate to the couple stress can be derived from the gradient of the independent micro-rotation field as: j ¼ rx Fig. 4. Definition of the quantities used in the Cosserat theory.
particle i to j depends on the movement of the system and Additionally therefore equals the continuum rotation x. the particles can rotate freely and thus add an independent rotation x* to the total rotation x. The stress in Cosserat continua incorporates a couple stress l (a torque per unit area) as well as the traditional force stress (force per unit area) (Fig. 5). As can be seen in Fig. 5, stress components in Cosserat continuum is similar to free body diagram in civil engineering. There are normal and shear forces and couple stresses can be thought as flexure. The kinematic of Cosserat point can be explained by FBD solution. It must be noted that body couples and body moments are ignored in this research. While classical continuum is appropriate for most materials, there are a number of elements that are better represented with rotational degree of freedom in addition to translational. If a material is composed of relatively rigid micro-scale structures that can rotate independently, a micropolar or Cosserat continuum might be more appropriate. The motion of each point in a micropolar continuum is described both by displacement and a rotation. The stress and strain components of the Cosserat continuum can be written as: e ¼ ½exx ; eyy ; ezz ; exy ; eyx ; jxz lc ; jyz lc T r ¼ ½rxx ; ryy ; rzz ; rxy ; ryx ; mxz lc ; myz lc
ð7Þ T
ð8Þ
Note that the stress tensor is no longer symmetric. The nonsymmetric component is needed to balance additional
ð10Þ
As explained in detail in the previous sections, Cosserat continuum is a generalization of classical continuum approach. Thus, similar to classical continua the differential operator matrix L for two dimensions must be determined. De Borst [49] explained in detail that the L operator as: 3 2 o 0 0 ox 7 6 o 6 0 oy 0 7 7 6 60 0 0 7 7 6 7 6 o 7 6 0 1 ð11Þ L¼6 ox 7 6 o 0 þ1 7 7 6 oy 7 6 6 0 0 lc o 7 4 ox 5 o 0 0 lc oy As can be seen in Eq. 11, the differential operator includes the length scale. Stress–strain relations for Cosserat continuum is developed by using a thermomechanical approach [54] r_ ¼ D_e
ð12Þ
where D is stiffness matrix and for the Cosserat continuum can be defined as: 3 2 2Ga1 2Ga2 2Ga2 0 0 0 0 7 6 0 0 0 0 7 6 2Ga2 2Ga1 2Ga2 7 6 6 2Ga2 2Ga2 2Ga1 0 0 0 0 7 7 6 7 6 D¼6 0 0 0 G þ Gc G Gc 0 0 7 7 6 6 0 0 0 G Gc G þ Gc 0 0 7 7 6 7 6 0 0 0 0 2G 0 5 4 0 0 0 0 0 0 0 2G
Fig. 5. Stress components in Cosserat continuum.
ð13Þ
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
where G is shear modulus, a1 = (1 v)/(1 2v), and a2v/ (1 2v). The elastic stiffness matrix includes additional material parameter Gc that is defined as Cosserat shear modulus. The Cosserat shear modulus has a great effect on the non-symmetric stress state. The parameter lc has the dimension of length and therefore a length scale will be introduced in continuum description. The details of FE formulation can be found in Arslan [27] and Arslan et al. [59,60]. The stress rate, plastic flow, and consistency condition are defined as: uu D0 Duw e_ e_ p r_ 0 _ s_ ¼ ¼ D0 : ð_c c_ p Þ ¼ : ð14Þ j_ j_ p Dwu Dww l_ 0 0 u e_ p _ ¼ k_ m ¼ k_ oQ=or ¼ km ð15Þ c_ p ¼ j_ p oQ=ol mw t oF oF F_ ¼ ð16Þ : s_ þ k_ ¼ nu : r_ þ nw : l_ H p k_ os ok where nu is the gradient of the yield function with respect to stress, nw is the gradient of the yield function with respect to couple stress; mu and mw are the gradient of the plastic potential with regard to s and l, k_ is the plastic multiplier. n : D0 : e_ k_ ¼ hp
ð17Þ
where hp = Hp + nt:D0:m The component of the elastoplastic tangent operator Dep can be written as: Duu ep
¼
Duu 0
1 Duu : mu nu : Duu 0 hp 0
1 uu D : mu nu : Duu 0 hp 0 1 uu Dwu D : mu nu : Duu ep ¼ 0 hp 0 1 uu ww Dww D : mu nu : Duu ep ¼ D0 0 hp 0 Duw ep ¼
ð18Þ
where hp reduces to: ww u w w hp ¼ H p þ nu : Duu 0 : m þ n : D0 : m
ð19Þ
where the hardening parameter Hp is defined as: oF ð20Þ ok There are many yield criteria for different types of materials. As the response of granular materials is pressure sensitive, hydrostatic state stress must be taken into account for a realistic yield criterion. The Drucker–Prager yield criterion [15] is a modification of the von Mises yield criterion that accounts for the influence of the hydrostatic stress component. Since the early developments in the 1950s and 1960s by Drucker and Prager, this yield criterion has become an established framework for modeling the inelastic behavior of soils. Nowadays, it has earned a broad interest and acceptance in geomechanics. The yield surface
Hp ¼
553
is a circular cone that overcomes numerical difficulties of the Mohr–Coulomb yield criteria. As it is powerful and capable in simulating the constitutive behavior of granular materials, a Drucker–Prager yield criteria is used for elastoplastic implementation of Cosserat continuum. A Drucker–Prager type yield function is used for the elastoplasticity of Cosserat continua. The yield function can be defined as: pffiffiffiffiffi F ¼ J 2 þ aJ 1 b J 1 ¼ trðrÞ; ! 1 s s 1 rs þ rt s : s þ 2l:l J2 ¼ ð21Þ with ss ¼ 2 2 lc The calibration parameters a and b may be expressed as: 6c cos / b ¼ pffiffiffi 3ð3 sin /Þ
2 sin / a ¼ pffiffiffi 3ð3 sin /Þ
ð22Þ
The plastic potential Q may be expressed the same as the yield function F with different calibration parameters replaces a and b for A and B. Thus, the gradients of the yield and plastic potential function become: u ( p1ffiffiffiffi ss þ aI ) 2 oF=or n oF 2 J2 ¼ ¼ n¼ ¼ ð23Þ 1ffiffiffiffi l w p os oF=ol n 2 J 2 l2c u ( p1ffiffiffiffi ss þ AI ) 2 oQ=or m oQ 2 J2 ¼ ¼ ð24Þ ¼ m¼ l 1 w pffiffiffiffi os oQ=ol m 2 J 2 2
lc
The equivalent plastic strain for isotropic hardening can be defined as: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c_ ¼ e_ p : e_ p þ l2c j_ p : j_ p ð25Þ Thus, the hardening–softening parameter may be calculated as: oF oa ob ocp ¼ I 1 ðrÞ ð26Þ Hp ¼ ok ocp ocp ok here cp is the equivalent plastic strain. The parameters are calibrated by Alshibli’s [30] plane strain experiments on Ottowa sand. The finite element implementation for the discretized equations was done using the nonlinear commercial finite element program ABAQUS. This program does not provide an element with material rotation; therefore a User Element Subroutine (UEL) is needed to solve the system of the finite element equations within the Micropolar framework. A 4-noded isoparametric element with four integration points was used. However, a selective reduced integration technique was used to avoid any possible volumetric locking during the softening regime. In this sense, full integration was used for all the state variables and only reduced integration technique was used for the volumetric strains. The finite element program ABAQUS [57] uses Newton Raphson iteration technique to fulfill the static equilibrium
554
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
equations and the load–displacement increments are updated using an implicit integration scheme within the standard version of ABAQUS. The problem in hand is a mixed control problem (load–displacement control) and all the internal state variables (such as stresses, plastic work, void ratio, etc.) are updated within the UEL using the explicit forward Euler integration scheme. Thereafter the ABAQUS postprocessor is used to show the analysis results. 4. Geotechnical engineering applications The purpose of this section is to evaluate the performance of the proposed Cosserat finite element implementation. In that study the effect of particle size and confinement on the orientation and thickness of the shear band was investigated. The plane strain numerical model is shown in Fig. 6. Instead of applying the displacement from the bottom of specimen as in plane strain tests, the vertical displacement is applied to the top of the specimen to avoid any numerical problems. The elastic material properties that are used in the Finite element analysis are summarized in Table 1. The length scale is in the range of 5–15 times the medium soil grain size which is related with the shear band thickness obtained from experimental studies. The length scales are ranging from 40% to 80% of the maximum size of a typical finite element in the finite element mesh. As Fig. 6 illustrates, the finite element results are in good agreement with the experiments. The model captures both thickness and orientation of localization quite well. The prediction of the stress–strain curves showed that the model fits the experimental results of granular materials under plane strain. The summary of FE and experimental results are given in Table 2. The finite element model is used for two common geotechnical engineering problems which are both plane strain.
Table 1 Material Parameters for the granular material Young’s modulus, E (kPa)
72,000
Poisson’s ratio, m Length scale (mm) Cosserat shear modulus, Gc Friction angle Cohesion (kPa) Dilatancy angle
0.26 Changes 0.5 * G 40° 5.0 15°
The first example will be the simulation of behavior of the granular material beneath a shallow foundation. A strip footing was assumed to be on the ground surface. The footing was assumed to be very stiff (steel material was assumed) and the soil stratum modeled as the foundation soil was assumed to be 10.0 m width and 5.0 m thick to avoid boundary effects. The loading was performed through pushing the footing towards the foundation soil a total of 25.0 mm. The second example is the stability of a slope with a footing resting on its crest. The footing is elastic and 100 times stiffer than the soil. The footing is subjected to a prescribed displacement at its center. The example assumes plane strain and Poisson’s ratio m = 0.3. The same material properties have been used in this study as given in Table 1. 5. Slope stability There are many situations where footings are constructed on sloping surfaces or adjacent to a slope crest, such as footings for bridge abutments on sloping embankments. It is well known that many parameters affect the stability of soil. The impact of length scale (or particle size as a physical meaning) on the failure of a slope will be tested with the finite element implementation of Cosserat continuum. Two different length scales have been used to illustrate the effect of the length scale. Mesh sizes and other
Fig. 6. Finite Element Experiment comparison of deviatoric stress–strain behavior of granular materials under plane strain condition.
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
555
Table 2 Summary of the FE-experiment comparison Soil type
Confining pressure (kPa)
d50 (mm)
F-75 Ottawa F-75 Ottawa Coarse Slica Sand Coarse Slica Sand
15 100 15 100
0.22 0.22 1.60 1.60
Shear band thickness (mm)
Shear band inclination
Experimental
FE
Experimental
FE
2.97 2.91 17.33 17.00
3.00 2.90 17.50 17.50
51.6 53.7 51.4 53.2
52.0 54.0 52.0 52.0
Fig. 7. Finite element mesh.
parameters are fixed and the length scale effect on the failure mechanism of slope will be compared. Fig. 7 illustrates the finite element mesh used in the analysis of slope failure model. As Figs. 8 and 9 illustrate, the thickness of localized zone is increasing with the increase of length scale. Fig. 10 shows the displacement vectors obtained by the Cosserat finite element solution. It can be seen that the largest vertical and horizontal displacements take place on the top surface and failure surface. The micro-rotation distribution is shown in Fig. 11a. As the figure shows, the Cosserat rotations are noticeable only in the shear zone. Fig. 11b illustrates the load–displacement plots for three
different length scales. As the figure shows, increasing length scale gives higher load responses. It is the key objective to relate the discrete level parameters to continuum parameters (stresses and strains). The finite element simulations in elastoplastic Cosserat elements enriched the classical continuum solutions by accommodating the microstructure to macro structure. 6. Footing punching problem The bearing capacity of a surface footing on granular soils has been a challenging problem approached experimentally, analytically and numerically by numerous
Fig. 8. Equivalent plastic strain distribution of slope failure (lc = 1.0 mm).
Fig. 9. Equivalent plastic strain distribution of slope failure (lc = 2.0 mm).
556
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Fig. 10. Velocity field.
Fig. 11a. Micro-rotation distribution of slope failure (lc = 2.0 mm).
Fig. 11b. Reaction force–displacement relation.
researchers. In the present example the behavior of granular soil under a rectangular footing a well-centered vertical displacement is investigated. Dimensions of the model are chosen large enough to eliminate the boundary condition effects. Bearing resistance is mobilized by applying equal displacement increments to the nodes beneath one half of the footing. As discussed by Arslan [27], displacement control provides numerical convenience, better stability, fewer iterations and represents physical reality (Fig. 12). The footing problem is solved for two different length scales. The first example is the finite element solution of the footing punching problem for the length scale, lc = 1.0 mm. Contours of resultant vertical stresses, shear
stresses and plastic strain viewed in plan are shown in Figs. 13–15, 16a, 16b, 17–19. The vertical stress distribution is uniform over most of the footing (Fig. 13). However, due to singularity at the footing edge, a high value of shear stress was predicted near the edge. The shear stress distribution is given in Fig. 14. As the figure shows, shear stress is localized at the corner of the footing. The granular foundation failure surface is not mobilized simultaneously, but in a progressive manner, so called ‘progressive failure’ [56]. The wedge Hill-type mechanism accompanying failure of the smooth footing is seen in Fig. 15. The deepest part of the failure mechanism observed for the smooth square footing takes place at 1.25 times the width of the footing.
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
557
Fig. 12. Footing punching problem.
Fig. 13. Vertical stress distributions inside the granular soil under footing.
Fig. 14. Shear stress distributions inside the granular soil under footing.
The kinematic mechanism at failure observed in the finite element analysis for rectangular footings is illustrated in Fig. 15, indicating a double wedge Hill-type mechanism (Fig. 20). Considering the geometry of the kinematic mechanisms accompanying failure in the finite element analyses
of rectangular footings various observations were made: plane strain conditions prevailed at approximately 1.25B from the ends of the footings (B is the width of the footing). Micro-rotation fields for this analysis are shown in Fig. 16a. The micro-rotational profile along the shear band
558
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Fig. 15. Equivalent plastic strain distributions inside the granular soil under footing.
Fig. 16a. Micro-rotation distributions inside the granular soil under footing.
Fig. 16b. Micro-rotational profile along the shear band.
is shown in Fig. 16b. As the figure shows, maximum microrotation occurs at the center of the shear band. Displacement vectors from the finite element analyses for a smooth square footing is shown in Fig. 17 and clearly indicates the different direction of soil movement directly
below the footing. As shown in the figures, localized plastic zones initiated at footing edges, and then further penetrated vertically in a narrow zone under the edge. The results of the numerical model showed the foundation soil failed with a localized shear failure mode. Unlike classical continuum, the incremental boundary value problem remains well posed in the micropolar solution and analysis can be continued until a collapse load analysis is attained or the vertical displacement of the footing exceeds a prescribed limit. Figs. 17 and 18 illustrate the finite element results of the footing problem for lc = 1.5 mm. The deepest part of the failure mechanism observed for the smooth square footing increased from 1.25 to 1.50 times the width of the footing. The different length scales provide a comparison to assess the particle size effect on the bearing capacity of the soil. The finite element results illustrates that the larger grain size leads to a deeper failure mechanism. 7. Length scale-size effect study on footing punching problem The effect of length scale on the bearing capacity of shallow foundations will be studied next. To demonstrate
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
559
Fig. 17. Velocity field.
Fig. 18. Vertical stress distributions inside the granular soil under footing.
Fig. 19. Equivalent plastic strain distributions inside the granular soil under footing.
significance of the length scale, three different cases are considered. The length scale effect on the analysis of different footing will be examined for three different scales and
three different footing dimensions (Fig. 21). In order to construct a length scale-bearing capacity relation, other parameters are fixed where only the length scales and the
560
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Fig. 20. Double wedge Hill-type failure mechanism [56].
Fig. 21. Footing and domain dimensions for three different cases.
size of the domain are varied. Footing dimensions and the size of domain changed at the same ratio as shown in Fig. 21. The example can be explained as: A. Analyzing the size effect for a footing by changing the dimensions of footing having the same material properties and loading condition. The footing sizes studied are 1.0, 10.0, and 100.0 m. The ratio of soil media dimension to footing dimension remains constant at 20:1 for all three examples. B. Length scale effect for a footing with the same dimensions but different length scales. The results are presented as plots of reaction force-vertical displacement curves. Fig. 22 shows load displacement curves of the foundations due to vertical displacement controlled loading of the footing for three different length scales and footing width (L) equals to 1. A 12 mm punching displacement has been applied. As Fig. 22 illustrates, an increase of the length scale from 0.1 mm to 10 mm leads to a clear increase of stiffness. While peak reaction forces of the footing for the length scale of 0.1 is 520 kN/m, it is 680 and 800 kN/m for the length scales 1.0 and 10.0 mm respectively. Virtually the same effect of length scale is shown by the Figs. 23 and 24. As can be seen, length scale has a great effect on the stiffness and peak reaction forces of the foundation. Similar to the plane strain example, the increase of the size of the footing leads to a decrease of peak reaction
Fig. 22. Footing reaction force–displacement curve (L = 1.0 m.).
forces. This can be due to higher stress rotations and higher micro-rotations effect for the larger dimensions. The resultant reaction force of the footing increases first, shows a pronounced peak, drops later and reaches slowly a residual state. The maximum reaction force of the footing increases with increasing length scales. The material becomes more ductile after the peak with increasing length scale. A footing on soil produces a strain field, which varies from relatively high strains immediately beneath the footing to infinitesimal values at large distances beneath and laterally away from the footing. Thus, there is a higher strain gradients due to higher stress concentrations at the corner of footing. Micropolar continua enables to take into account the
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
561
The absence of intrinsic length scale leaves the width of localized zone unspecified. Thus, the finite element solution of the classical continuum is mesh sensitive. To solve the mesh sensitivity and to investigate the effect of microstructure, Cosserat continuum theory is implemented into Drucker Prager criterion and subsequently used in finite element analysis. The finite element analyses were carried out using the commercial nonlinear finite element code ABAQUS. The Cosserat element, with the additional degrees of freedom was implemented in ABAQUS using the UEL interface. This formulation was used in order to analyze selected geotechnical problems or configuration. Some of the major findings of this research are: Fig. 23. Footing reaction force–displacement curve (L = 10.0 m.).
The numerical model was used to predict the load–displacement behavior of a shallow footing and failure behavior of a slope. A Hill-type failure mechanism accompanying failure of the smooth footing was observed. The change of length scale affects the failure mechanism. The bearing capacity (maximum reaction force) of the footing increases with increasing length scales. The material becomes more ductile after the peak with increasing length scale. The thickness of the localized zone in the slope failure increases with increasing length scale. FE calculations show that localization in granular materials can realistically be determined with Cosserat elastoplastic constitutive law. Acknowledgements
Fig. 24. Footing reaction force–displacement curve (L = 100 m.).
higher strain gradient and effect of length scale would be observed clearly under such a higher strain gradient. The bearing capacity of soil to support a footing depends on the properties of soil such as density shear strength also mechanical properties such as particle sizes. Length scale has significantly affected the calculated bearing capacity of a foundation. This analysis clearly demonstrates the effect of the micropolar approach in the numerical simulation of bearing capacity and the significance of length scale in the analysis of bearing capacity. It must be noted that the mesh independent results can be obtained within enhanced continuum if the mesh size length scale ratio is small. For large problems, however, it is not possible without remeshing to achieve this due to the limited calculation time. On the other side, if too large length scale is artificially assumed in calculations, the FE results would not be realistic since both strength and ductility change with length scale. 8. Conclusions The basic deficiency of a classical continuum with regard to microstructure phenomena is the lack of a length scale.
The authors would like to thank Prof. Kaspar Willam for his guidance and help during the study. The comments and suggestions of anonymous reviewers that helped to improve the quality of the paper are also greatly appreciated. The authors gratefully acknowledge the financial support provided by NASA under Contract No. NCC8-242. References [1] Hadamard J. Propagation des Ondes et les Equations d’Hydrodynamique. New York: Chelsea; 1949 [reprinted 1903]. [2] Thomas TY. Plastic flow and fracture in solids. New York (NY): Academic Press; 1961. [3] Hill R. Acceleration waves in solids. J Mech Phys Solids 1962;10: 1–16. [4] Mandel J. Conditions de Stabilite et Postulat de Drucker. In: Kravtchenko, Siireys, editors. Proceedings of the IUTAM symposium on rheology and soil mechanics. Berlin: Springer; 1964. p. 58–68. [5] Rice J. The localization of deformation. In: Koiter W, editor. Theoretical and applied mechanics. North Publishing Company; 1976. [6] Rudnicki JW, Rice JR. Conditions for the localization of deformation in pressure sensitive dilatant materials. J Mech Phys Solids 1975;23:371–94. [7] Rice JR, Rudnicki JW. A note on some features of the theory of localization of deformation. Int J Solids Struct 1980;16:597–605. [8] Vardoulakis I. Equilibrium bifurcation of granular earth bodies. Adv Anal Geotech Instabilities 1978;13:65–120.
562
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
[9] Vardoulakis I. Bifurcation analysis of the triaxial test on sand samples. Acta Mech 1979;32:35–54. [10] Vardoulakis I. Shear band inclination and shear modulus of sand in biaxial tests. Int J Numer Anal Meth Geomech 1980;4:103–19. [11] Vardoulakis I. Constitutive properties of dry sand observable in the triaxial test. Acta Mech 1981;38:219–39. [12] Vardoulakis I. Shear-band and liquefaction in granular materials on the basis of Cosserat continuum theory. Ingenieur-Archiv 1989;59: 106–13. [13] Molenkamp F. Comparison of frictional material models with respect to shear band initiation. Geotechnique 1985;35:127–43. [14] Peric D. Localized deformation and failure analysis of pressure sensitive granular materials. Ph.D. dissertation, University of Colorado-Boulder; 1990. [15] Hill R. A general theory of uniquness and stability in elastic–plastic solids. J Mech Phys Solids 1958;6:36–249. [16] Hill R, Hutcinson KW. Bifurcation phenomena in the plane strain test. J Mech Phys Solids 1975;23:239–64. [17] Arthur J, Dunstan T, Al-Ani Q, Asadi A. Plastic deformation and failure in granular media. Geotechnique 1977;27:53–74. [18] Desrues J. Localization de la deformation plastique dans les materiaux granulaires. These d’etat, Universite’ de Grenoble; 1984. p. 180–85. [19] Tejchman J. Modelling of shear localisation and autogeneous dynamic effects in granular bodies. Publication series of the Institute of Soil and Rock Mechanics, vol. 140. University Karlsruhe; 1997. [20] Tejchman J, Herle I, Wehr J. FE-Studies on the influence of initial void ratio, pressure level, and mean grain diameter on shear localization. Int J Numer Anal Meth Geomech 1999;23:2045–74. [21] Maier T. Nonlocal modeling of softening in hypoplasticity. Comput Geotech 2003;30:599–610. [22] Manzari MT. Application of micropolar plasticity to post failure analysis in geomechanics. Int J Numer Anal Meth Geomech 2004;28:1011–32. [23] Voyiadjis JZ, Alsaleh MI, Alshibli KA. Evolving internal length scales in plastic strain localization for granular materials. Int J Plasticity 2005;21:2000–24. [24] Alsaleh MI, Voyiadjis GZ, Alshibli KA. Modeling strain localization in granular materials using micropolar theory: mathematical formulations. Int J Numer Anal Meth Geomech 2006;30(15):1501–24. [25] Alsaleh MI, Alshibli KA, Voyiadjis GZ. Influence of micro-material heterogeneity on strain localization in granular materials. ASCE: Int J Geomech 2006;6(4):248–59. [26] Alshibli KA, Alsaleh MI, Voyiadjis GZ. Modeling strain localization in granular materials using micropolar theory: numerical implementation and verification. Int J Numer Anal Meth Geomech 2006;30(15):1525–44. [27] Arslan H. Localization analysis of granular materials in Cosserat Elastoplasticity –Formulation and finite element implementation. PhD dissertation, University of Colorado-Boulder; 2006. [28] Coulomb C.A. Essai sur une application des re`gles des maximis et minimis a` quelques proble`mes de statique relatifs a` l’architecture. Me`m. Acad. Roy. Pre`s. Divers Savants 7, Paris; 1776. [29] Roscoe KH. Tenth rankine lecture: the influence of strains in soil mechanics. Geotechnique 1970;20:129–70. [30] Alshibli K.A. Localized deformations in granular materials. Ph.D. dissertation, University of Colorado at Boulder, CO; 1995. [31] Vermeer P.A. A simple shear-band analysis using compliances. In: IUTAM conference on deformation and failure of granular materials, Delft; 1982. p. 493–9. [32] Ottosen N, Runesson K. Acceleration waves in elastoplasticity. Int J Solids Struct 1991;28:135–59. [33] Raniecki B, Bruhns OT. Bounds to bifurcation stresses in solids with non-associated flow law at finite strain. J Mech Phys Solids 1981;29:153–72.
[34] Ottosen NS, Runesson K, Peric D. Discontinuous bifurcations of elastic–plastic solutions at plane stress and plane strain. Int J Plasticity 1991;7:99–121. [35] Runesson K, Mroz Z. A note on non associated plastic flow rules. Int J Plasticity 1989;5:639–58. [36] Ottosen N, Runesson K. Properties of discontinuous bifurcation solutions in elasto-plasticity. Int J Solids Struct 1991;27:401–21. [37] Jira´sek M. Modeling of localized inelastic deformation, Short Course LID, Prague; 19–23 September 2005. [38] Chakrabarty J. Theory of plasticity. New York: McGraw-Hill; 1987. [39] Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by softening relations in rate-independent solids. Comp Mech 1993;12:277–96. [40] Oliver J. Continuum modelling of strong discontinuities in solid mechanics. In: Owen DRJ, Onate E, editors. Computational plasticity. Fundamentals and applications, vol. 1. Pineridge Press; 1995. p. 455–79. [41] Oliver J. Continuum modelling of strong discontinuities in solid mechanics using damage models. Comput Mech 1995;17(1/2): 49–61. [42] Larsson R, Runesson K, Sture S. Embedded localization band in undrained soil based on regularized strong discontinuity theory and finite element analysis. Int J Solids Struct 1996;33(20–22):3081–101. [43] Steinmann P. A finite element formulation for strong discontinuities in fluid saturated porous media. Mech Cohes-Frict Mater 1999;4:133–52. [44] Willam K et al. Localization in micropolar continua, continuum models for materials with microstructure. Wiley; 1995. p. 297–339. [45] Iordache MM, Willam K. Localized failure analysis in elastoplastic Cosserat continua. Comput Meth Appl Mech Eng 1998;151: 559–86. [46] Cosserat EF. Theorie des corps deformable. Paris: Hermann; 1909. [47] Muhlhaus HB, Vardoulakis I. The thickness of shear bands in granular materials. Geotechnique 1987;37(3):71–283. [48] Muhlhaus HB. Application of Cosserat theory in numerical solutions of limit load problems. Ing Arch 1989;59:124–37. [49] de Borst R. Simulation of strain localization: a reappraisal of the Cosserat continuum. Eng Comput 1991;8:317–32. [50] Eringen AC, Suhubi ES. Nonlinear theory of simple micro-elastic solids – I. Int J Eng Sci 1964;2:189–203. [51] Eringen AC. A unified theory of thermomechnical materials. Int J Eng Sci 1966;4:179. [52] Eringen CA. Microcontinuum field theories I: foundations and solids. New York: Springer-Verlag, Inc; 1999. [53] Forest S, Barbe F, Cailletaud G. Cosserat modelling of size effects in the mechanical behaviour of polycrystals and multiphase materials. Int J Solids Struct 2000;37:7105–26. [54] Walsh SDC, Tordesillas A. A thermomechanical approach to the development of micropolar constitutive models of granular media. Acta Mech 2004;167:145–69. [55] Arslan H, Sture S, Willam K. Analytical and geometrical representation of localization in granular materials. Acta Mech 2007;194(1–4): 159–73. [56] Meyerhof GG. The ultimate bearing capacity of foundations. Geotechnique 1951;2(4):301–32. [57] ABAQUS Version 6.5. Manuals, Hibbitt, Karlsson and Sorensen Inc; 2004. [58] Arslan H, Sture S. Evaluation of a physical length scale for granular materials. Comput Mater Sci. doi:10.1016/j.commatsci.2007.08.016. [59] Arslan H, Sture S. Finite element analysis of localization and micromacro structure relation in granular materials: Part I: formulation. Acta Mech [in press]. [60] Arslan H, Sture S. Finite element analysis of localization and micromacro structure relation in granular materials: Part II: implementation and simulations. Acta Mech [in press].
Computers and Geotechnics 35 (2008) 548–562 www.elsevier.com/locate/compgeo
Finite element simulation of localization in granular materials by micropolar continuum approach Haydar Arslan a
a,*
, Stein Sture
b
Exxon Mobil Upstream Research Company, P.O. Box 2189, Greenway Plaza-3, Room 704-A, Houston, TX, United States b University of Colorado, 26 UCB, Boulder, Colorado 80309-0026, United States Received 23 March 2007; received in revised form 13 October 2007; accepted 24 October 2007 Available online 8 January 2008
Abstract There are two different approaches in continuum mechanics: a continuum considered as solid state – classical continuum mechanics – or a continuum with the microstructure. In order to explain the influences of microstructures on the behavior of materials, higher order effects such as moments that reflect the grain geometry and kinematics (e.g. rotations) must be incorporated, which in turn lead to a basic asymmetry of shear stress and micropolar (or Cosserat) effects. Cosserat continuum incorporates the characteristic length that relates the couple stresses and micro-curvatures at the constitutive model. Finite Element solutions of localization analysis have been carried out in Cosserat continuum to investigate length scale-size effect on the overall behavior of granular materials. The formulations of Cosserat elastoplasticity have been implemented into the nonlinear commercial finite element code using the user defined element utility. The mesh dependency effect in FE solution is reduced. The numerical model was used to predict the load–displacement behavior of a shallow footing and failure behavior of a slope. A proper treatment of the localization phenomenon in granular materials has been observed and simulated in Finite Element analysis. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Localization; Length scale; Size effect; Micropolar; Granular materials
1. Introduction Strain localization in metals, composite materials and geomaterials has been studied since the early 1900s. The strain localization phenomenon has received considerable research attention lately, and much of recent activities have focused on developing mathematical and numerical models for capturing or simulating the experimentally observed strain localization phenomena. The theory of strain localization is an analytical framework to analyze shear bands in materials. It was first proposed by Hadamard [1] and later developed by Thomas [2], Hill [3], Mandel [4], and Rice [5]. It was applied to predict the existence and orientation of shear bands within various *
Corresponding author. Tel.: +1 720 231 9239. E-mail addresses: [email protected] (H. Arslan), stein.sture@ colorado.edu (S. Sture). 0266-352X/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2007.10.006
types of materials, including elastoplastic soils and rocks [6–14]. Rudnicki and Rice [6] established that the elastoplastic Drucker–Prager model predicts localization during triaxial compression of pressure-sensitive materials. Vardoulakis [10] extended the bifurcation analysis of Hill [15], Hill and Hutchinson [16] to sands. By using a constitutive model of the deformation theory of plasticity, he analytically derived the empirical relation that Arthur et al. [17] proposed for defining the inclination of shear bands. Desrues [18] extensively and accurately measured the direction of shear bands during plane-strain compression of sand specimens. However, he could not describe or explain his experimental results by using the relationships developed by Arthur et al. [17] and Vardoulakis [10]. Finite Element analysis studies related to the strain localization process using conventional continuum mechanics principles have shown that shear band formation is mesh
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
dependent. The mesh dependent solutions might lead to a vanishing shear band thickness as the mesh size is refined in the finite element solution. Tejchman [19], Tejchman et al. [20], Maier [21], Manzari [22], Voyiadis et al. [23], Alsaleh et al. [24,25], Alshibli et al. [26], Arslan [27] confirmed that the shear band formation is mesh independent in Cosserat continuum analysis due to the regularization of the deformations by additional (rotational) degrees of freedom. 2. Background on localization analysis solution techniques Coulomb [28] first introduced the concept of slip planes in frictional materials. The original concept of Coulomb has since evolved and reached substantial maturity after the development and implementation of the theory of plasticity as it applies to soils. Roscoe [29] developed a slip plane theory for pressure sensitive materials such as soils, where the subsequently termed Roscoe criterion is based on kinematics rather than static or failure conditions. The Mohr–Coulomb criterion, which is only based on statics would explain the orientation of the failure plane but does not provide any information related to geometry or kinematics. The classical Mohr–Coulomb criterion as developed for two dimensions, calculates the orientation of shear band as: hc ¼ ð45 þ /=2Þ
ð1Þ
where / is the internal friction angle and hc is the angle between the major principal direction and shear band (Fig. 1). It should be noted that the equation does not include any variables related to failure, except failure plane orientation based on static concepts. Roscoe [29] introduced a kinematical assumption for the orientation of the shear band as: hR ¼ ð45 þ w=2Þ
Fig. 1. Mohr–Coulomb solution to shear band orientation.
ð2Þ
549
where hR is the angle between the shear band and the direction of the major principle stress and w is the angle of dilatancy, which is calculated by: w¼
e_ v e_ q
ð3Þ
where e_ v is volumetric strain rate and e_ q is the shear strain rate. Arthur et al. [17] proposed an empirical equation for shear band orientation based on plane strain experimental results: 1 1 hA ¼ 45 þ ð/ þ wÞ ¼ ðhC þ hR Þ ð4Þ 4 2 Alshibli [30] have conducted plane strain experiments on different sands and the experimental results verified the Eq. (4) proposed by Arthur et al. [17]. The equation was theoretically confirmed by Vermer [31]. By using plasticity concepts, Mandel [4] derived the critical hardening parameter for a granular material subjected to loading in plane strain. The critical hardening parameter is associated with the possible existence of discontinuous bifurcation modes and possible loss of uniqueness. The orientation of possible band of localized deformation is determined by the values of a critical hardening parameter Hcr. Mandel did not derive expressions for the bifurcation directions. Vermeer [31] rederived Mandel’s equation that includes the orientation of the shear band as follows: tan2 hV ¼
2 þ ðsin / þ sin wÞ 2 ðsin / þ sin wÞ
ð5Þ
As can bee seen, Arthur’s and Vermeer’s expressions are in good agreement. Vermeer derived the critical hardening parameter as: H cr 1 2 ðsin / sin wÞ ¼ 2G 16ð1 tÞ
ð6Þ
Vardoulakis [8] performed theoretical and experimental investigations for the shear band orientation by using the Mohr–Coulomb failure criterion. He found an expression for the shear band orientation similar to those developed by Arthur and Vermeer. Drucker–Prager and Mohr–Coulomb plasticity are fundamentally different failure theories. Formulations and the terms that used for the yield functions are different. However, the terms are based on friction and cohesion and they both predict similar results. Ottosen and Runesson [32] showed that the localization angle coincides with the failure angle of the yield condition for the Mohr–Coulomb yield criterion based on associatedplastic flow concepts. The main disadvantage of the many different classical Mohr–Coulomb criteria is that the failure plane is independent of the intermediate principal stress. Furthermore, it does not explain the material behavior after the material failure initiation. The Drucker–Prager yield condition, which takes into account all three principal stresses generalizes the
550
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Mohr–Coulomb failure condition. Thus, the stress dependent localization analysis may be investigated by this model. However, the Drucker–Prager model does not provide failure mode information. Hill [15], introduced material bifurcation concepts and processes to plasticity. Later on, Rudnicki and Rice [6], Ranecki and Bruhns [33], Runesson and Mroz [35], Rynesson et al. [34], Ottosen and Runesson [36] considered the localization as a discontinuous bifurcation problem. The bifurcation analysis considers a necessary condition for discontinuous bifurcation, which is expressed as a singularity of the localization tensor which is similar to the elastic acoustic tensor [55]. The localization condition signals the potential formation of a strain discontinuity. From a mathematical point of view, these complicating features are related to the socalled loss of ellipticity of the governing differential equation. The boundary value problem becomes ill-posed, i.e., it does not have a unique solution with continuous dependence on the given data. Thus, a discontinuity approach is necessary to develop a continuum formulation that follows
the failure behavior correctly. The readers are referred to chapter three for detail explanations. The initial thickness of the localization band depends on the material microstructure. Mathematically, narrow zones of highly concentrated evolving strain can be represented in many different ways. There are mainly three types of discontinuity in continuum mechanics as shown in Figs. 2 and 3. The first type of discontinuity incorporates strong discontinuities, i.e., jumps in displacements across a discontinuity curve (in two dimensions) or discontinuity surface (in three dimensions). In the simplest case, the strong discontinuity can be considered as stress-free, which necessarily leads to a stress singularity at the discontinuity line. In physical terms, the strong discontinuity corresponds to a sharp crack. The strong discontinuity approach has been used as a solution technique for quasi-static solid mechanics problems. It exhibits jumps in the displacement field across a line in two dimensional problems. The line can be named as slip line or discontinuity line. Chakrabarty [38] introduced the strong discontinuity approach to a von
Fig. 2. Types of discontinuity [37].
Fig. 3. Classification of discontinuities for localization.
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Mises-based plasticity model to explain the observed shear bands in metals. Later on, Simo et al. [39], Oliver [40,41], Larson et al. [42], and Steinmann [43] continued to use the concept of a strong discontinuity as a tool for capturing localization in solid mechanics. Another possible kinematic description represents the region of localized deformation by a band of a small but of finite thickness, separated from the remaining part of the body by two weak discontinuities. Since the displacement is continuous, the strain components in the plane tangential to the discontinuity surface must remain continuous as well, and only the out-of-plane components can have a jump. In physical terms, the band between the weak discontinuities corresponds to a damage process zone with an almost constant density of microdefects. Instead of lumping all the inelastic effects into a surface, it is possible to distribute them uniformly across the width of a band of a finite thickness. As Fig. 3 illustrates, a weak discontinuity is characterized by discontinuous but bounded strain fields. Since the displacement is continuous, the strain components in the plane tangential to the discontinuity surface must remain continuous as well, and only the out-of-plane components can exhibit a jump. Runesson et al. [34], Peric [14], Willam et al. [44], Iordache and Willam [45], and Arslan et al. [55] have used the weak discontinuity concept for bifurcation solutions in elasto-plastic materials. Finally, the most regular description uses a continuously differentiable displacement field, and the strain field remains continuous. Strain localization is manifested by high strains in a narrow band, with a continuous transition to much lower strain levels in the surrounding parts of the body. In physical terms, this corresponds to a damage process zone with a continuously varying concentration of defects. As the experimental and numerical studies have shown, the shear band has a finite thickness inside granular materials. The formation of weak discontinuities characterized by continuous displacements but discontinuous strains, which concentrate into a band of finite width, thereby this approach is more suitable to capture and calculate the thickness of the shear band. The weak discontinuity has been used in this study to model and capture localization. It must be recognized that the predictions of localization conditions is not only dependent on the chosen discontinuity approach but also to the constitutive description. The underlying physical definition of the behavior of granular materials must be considered for a well representative formulation. When the applied stress reaches the reduced strength, softening starts and the stress decreases. Consequently, the material outside the weaker region must unload elastically, because its strength limit has not been reached. This leads to the conclusion that the size of the softening region cannot exceed the size of the region having minimum strength. In the present one-dimensional setting, loss of ellipticity occurs, when the tangent modulus ceases to be
551
positive. From the numerical point of view, ill-posedness is manifested by pathological sensitivity of the results to the size of finite elements. The limit solution represents a displacement jump at the center, with zero strains everywhere else. Plasticity models have been used to capture discontinuities. As Jirasek [37] explained in detail, the strong discontinuity can be formulated as a traction-separation law, which complements the stress–strain law valid for the continuous part of the body. Instead, models with localization bands bounded by weak discontinuities can be considered as simple regularizations of models with strong discontinuities. Full regularization of the localization problem can be achieved by a proper generalization of the continuum theory. Various techniques have been implemented to provide a physically acceptable solution. Some of these techniques impose restrictions on the constitutive moduli in the postlocalization regime, while others artificially restrict the size of the finite element relative to the localization zone. The former is based on enriching the continuum with non-conventional constitutive relations in such a way that an internal or characteristic length scale is introduced. Non-local theories of this kind include the Cosserat continuum theory [46], higher gradient theories, and the integral or gradient theory. Localization analysis in Cosserat continua has typically focused on analytical and finite element simulation [20– 22,26,47–49,55]. Main research interests have been to illustrate the applicability of the higher order micro-continuum theory to elastic–plastic solid mechanics analysis. The progress in implementation of Cosserat continuum in material modeling has led to more accurate and realistic constitutive relations for granular materials. Instead of going into detailed formulation, the nature and progress in Cosserat continuum approach will be summarized. Additional material parameters in Micropolar continuum and different solution approaches by different researchers will be reviewed to clarify the main underlying concepts of the Cosserat continuum. 3. Method of analysis The Cosserat continuum belongs to the larger class of generalized continua which introduce intrinsic length scales in continuum mechanics via higher order gradients, additional degrees of freedom in fully non-local constitutive equations [50–53,58]. Eringen [52] defined a Cosserat medium as a continuous collection of particles that behave like rigid bodies. Accordingly, translational and rotational degrees of freedom describe displacement and rotation of an underlying microstructure. Cosserat theory incorporates a local rotation of points (Fig. 4) as well as the translation assumed in classical elasticity. N and n are the socalled directors that show the orientations of particles in the reference and the actual configuration, respectively. The rotation of the joint from
552
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
moments created by couple stresses. The infinitesimal strain tensor encompasses the displacement gradient and micro-rotation. Thus, the strain tensor can be written as: e ¼ ru e x
ð9Þ
The micro-curvature strain tensor that is conjugate to the couple stress can be derived from the gradient of the independent micro-rotation field as: j ¼ rx Fig. 4. Definition of the quantities used in the Cosserat theory.
particle i to j depends on the movement of the system and Additionally therefore equals the continuum rotation x. the particles can rotate freely and thus add an independent rotation x* to the total rotation x. The stress in Cosserat continua incorporates a couple stress l (a torque per unit area) as well as the traditional force stress (force per unit area) (Fig. 5). As can be seen in Fig. 5, stress components in Cosserat continuum is similar to free body diagram in civil engineering. There are normal and shear forces and couple stresses can be thought as flexure. The kinematic of Cosserat point can be explained by FBD solution. It must be noted that body couples and body moments are ignored in this research. While classical continuum is appropriate for most materials, there are a number of elements that are better represented with rotational degree of freedom in addition to translational. If a material is composed of relatively rigid micro-scale structures that can rotate independently, a micropolar or Cosserat continuum might be more appropriate. The motion of each point in a micropolar continuum is described both by displacement and a rotation. The stress and strain components of the Cosserat continuum can be written as: e ¼ ½exx ; eyy ; ezz ; exy ; eyx ; jxz lc ; jyz lc T r ¼ ½rxx ; ryy ; rzz ; rxy ; ryx ; mxz lc ; myz lc
ð7Þ T
ð8Þ
Note that the stress tensor is no longer symmetric. The nonsymmetric component is needed to balance additional
ð10Þ
As explained in detail in the previous sections, Cosserat continuum is a generalization of classical continuum approach. Thus, similar to classical continua the differential operator matrix L for two dimensions must be determined. De Borst [49] explained in detail that the L operator as: 3 2 o 0 0 ox 7 6 o 6 0 oy 0 7 7 6 60 0 0 7 7 6 7 6 o 7 6 0 1 ð11Þ L¼6 ox 7 6 o 0 þ1 7 7 6 oy 7 6 6 0 0 lc o 7 4 ox 5 o 0 0 lc oy As can be seen in Eq. 11, the differential operator includes the length scale. Stress–strain relations for Cosserat continuum is developed by using a thermomechanical approach [54] r_ ¼ D_e
ð12Þ
where D is stiffness matrix and for the Cosserat continuum can be defined as: 3 2 2Ga1 2Ga2 2Ga2 0 0 0 0 7 6 0 0 0 0 7 6 2Ga2 2Ga1 2Ga2 7 6 6 2Ga2 2Ga2 2Ga1 0 0 0 0 7 7 6 7 6 D¼6 0 0 0 G þ Gc G Gc 0 0 7 7 6 6 0 0 0 G Gc G þ Gc 0 0 7 7 6 7 6 0 0 0 0 2G 0 5 4 0 0 0 0 0 0 0 2G
Fig. 5. Stress components in Cosserat continuum.
ð13Þ
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
where G is shear modulus, a1 = (1 v)/(1 2v), and a2v/ (1 2v). The elastic stiffness matrix includes additional material parameter Gc that is defined as Cosserat shear modulus. The Cosserat shear modulus has a great effect on the non-symmetric stress state. The parameter lc has the dimension of length and therefore a length scale will be introduced in continuum description. The details of FE formulation can be found in Arslan [27] and Arslan et al. [59,60]. The stress rate, plastic flow, and consistency condition are defined as: uu D0 Duw e_ e_ p r_ 0 _ s_ ¼ ¼ D0 : ð_c c_ p Þ ¼ : ð14Þ j_ j_ p Dwu Dww l_ 0 0 u e_ p _ ¼ k_ m ¼ k_ oQ=or ¼ km ð15Þ c_ p ¼ j_ p oQ=ol mw t oF oF F_ ¼ ð16Þ : s_ þ k_ ¼ nu : r_ þ nw : l_ H p k_ os ok where nu is the gradient of the yield function with respect to stress, nw is the gradient of the yield function with respect to couple stress; mu and mw are the gradient of the plastic potential with regard to s and l, k_ is the plastic multiplier. n : D0 : e_ k_ ¼ hp
ð17Þ
where hp = Hp + nt:D0:m The component of the elastoplastic tangent operator Dep can be written as: Duu ep
¼
Duu 0
1 Duu : mu nu : Duu 0 hp 0
1 uu D : mu nu : Duu 0 hp 0 1 uu Dwu D : mu nu : Duu ep ¼ 0 hp 0 1 uu ww Dww D : mu nu : Duu ep ¼ D0 0 hp 0 Duw ep ¼
ð18Þ
where hp reduces to: ww u w w hp ¼ H p þ nu : Duu 0 : m þ n : D0 : m
ð19Þ
where the hardening parameter Hp is defined as: oF ð20Þ ok There are many yield criteria for different types of materials. As the response of granular materials is pressure sensitive, hydrostatic state stress must be taken into account for a realistic yield criterion. The Drucker–Prager yield criterion [15] is a modification of the von Mises yield criterion that accounts for the influence of the hydrostatic stress component. Since the early developments in the 1950s and 1960s by Drucker and Prager, this yield criterion has become an established framework for modeling the inelastic behavior of soils. Nowadays, it has earned a broad interest and acceptance in geomechanics. The yield surface
Hp ¼
553
is a circular cone that overcomes numerical difficulties of the Mohr–Coulomb yield criteria. As it is powerful and capable in simulating the constitutive behavior of granular materials, a Drucker–Prager yield criteria is used for elastoplastic implementation of Cosserat continuum. A Drucker–Prager type yield function is used for the elastoplasticity of Cosserat continua. The yield function can be defined as: pffiffiffiffiffi F ¼ J 2 þ aJ 1 b J 1 ¼ trðrÞ; ! 1 s s 1 rs þ rt s : s þ 2l:l J2 ¼ ð21Þ with ss ¼ 2 2 lc The calibration parameters a and b may be expressed as: 6c cos / b ¼ pffiffiffi 3ð3 sin /Þ
2 sin / a ¼ pffiffiffi 3ð3 sin /Þ
ð22Þ
The plastic potential Q may be expressed the same as the yield function F with different calibration parameters replaces a and b for A and B. Thus, the gradients of the yield and plastic potential function become: u ( p1ffiffiffiffi ss þ aI ) 2 oF=or n oF 2 J2 ¼ ¼ n¼ ¼ ð23Þ 1ffiffiffiffi l w p os oF=ol n 2 J 2 l2c u ( p1ffiffiffiffi ss þ AI ) 2 oQ=or m oQ 2 J2 ¼ ¼ ð24Þ ¼ m¼ l 1 w pffiffiffiffi os oQ=ol m 2 J 2 2
lc
The equivalent plastic strain for isotropic hardening can be defined as: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c_ ¼ e_ p : e_ p þ l2c j_ p : j_ p ð25Þ Thus, the hardening–softening parameter may be calculated as: oF oa ob ocp ¼ I 1 ðrÞ ð26Þ Hp ¼ ok ocp ocp ok here cp is the equivalent plastic strain. The parameters are calibrated by Alshibli’s [30] plane strain experiments on Ottowa sand. The finite element implementation for the discretized equations was done using the nonlinear commercial finite element program ABAQUS. This program does not provide an element with material rotation; therefore a User Element Subroutine (UEL) is needed to solve the system of the finite element equations within the Micropolar framework. A 4-noded isoparametric element with four integration points was used. However, a selective reduced integration technique was used to avoid any possible volumetric locking during the softening regime. In this sense, full integration was used for all the state variables and only reduced integration technique was used for the volumetric strains. The finite element program ABAQUS [57] uses Newton Raphson iteration technique to fulfill the static equilibrium
554
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
equations and the load–displacement increments are updated using an implicit integration scheme within the standard version of ABAQUS. The problem in hand is a mixed control problem (load–displacement control) and all the internal state variables (such as stresses, plastic work, void ratio, etc.) are updated within the UEL using the explicit forward Euler integration scheme. Thereafter the ABAQUS postprocessor is used to show the analysis results. 4. Geotechnical engineering applications The purpose of this section is to evaluate the performance of the proposed Cosserat finite element implementation. In that study the effect of particle size and confinement on the orientation and thickness of the shear band was investigated. The plane strain numerical model is shown in Fig. 6. Instead of applying the displacement from the bottom of specimen as in plane strain tests, the vertical displacement is applied to the top of the specimen to avoid any numerical problems. The elastic material properties that are used in the Finite element analysis are summarized in Table 1. The length scale is in the range of 5–15 times the medium soil grain size which is related with the shear band thickness obtained from experimental studies. The length scales are ranging from 40% to 80% of the maximum size of a typical finite element in the finite element mesh. As Fig. 6 illustrates, the finite element results are in good agreement with the experiments. The model captures both thickness and orientation of localization quite well. The prediction of the stress–strain curves showed that the model fits the experimental results of granular materials under plane strain. The summary of FE and experimental results are given in Table 2. The finite element model is used for two common geotechnical engineering problems which are both plane strain.
Table 1 Material Parameters for the granular material Young’s modulus, E (kPa)
72,000
Poisson’s ratio, m Length scale (mm) Cosserat shear modulus, Gc Friction angle Cohesion (kPa) Dilatancy angle
0.26 Changes 0.5 * G 40° 5.0 15°
The first example will be the simulation of behavior of the granular material beneath a shallow foundation. A strip footing was assumed to be on the ground surface. The footing was assumed to be very stiff (steel material was assumed) and the soil stratum modeled as the foundation soil was assumed to be 10.0 m width and 5.0 m thick to avoid boundary effects. The loading was performed through pushing the footing towards the foundation soil a total of 25.0 mm. The second example is the stability of a slope with a footing resting on its crest. The footing is elastic and 100 times stiffer than the soil. The footing is subjected to a prescribed displacement at its center. The example assumes plane strain and Poisson’s ratio m = 0.3. The same material properties have been used in this study as given in Table 1. 5. Slope stability There are many situations where footings are constructed on sloping surfaces or adjacent to a slope crest, such as footings for bridge abutments on sloping embankments. It is well known that many parameters affect the stability of soil. The impact of length scale (or particle size as a physical meaning) on the failure of a slope will be tested with the finite element implementation of Cosserat continuum. Two different length scales have been used to illustrate the effect of the length scale. Mesh sizes and other
Fig. 6. Finite Element Experiment comparison of deviatoric stress–strain behavior of granular materials under plane strain condition.
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
555
Table 2 Summary of the FE-experiment comparison Soil type
Confining pressure (kPa)
d50 (mm)
F-75 Ottawa F-75 Ottawa Coarse Slica Sand Coarse Slica Sand
15 100 15 100
0.22 0.22 1.60 1.60
Shear band thickness (mm)
Shear band inclination
Experimental
FE
Experimental
FE
2.97 2.91 17.33 17.00
3.00 2.90 17.50 17.50
51.6 53.7 51.4 53.2
52.0 54.0 52.0 52.0
Fig. 7. Finite element mesh.
parameters are fixed and the length scale effect on the failure mechanism of slope will be compared. Fig. 7 illustrates the finite element mesh used in the analysis of slope failure model. As Figs. 8 and 9 illustrate, the thickness of localized zone is increasing with the increase of length scale. Fig. 10 shows the displacement vectors obtained by the Cosserat finite element solution. It can be seen that the largest vertical and horizontal displacements take place on the top surface and failure surface. The micro-rotation distribution is shown in Fig. 11a. As the figure shows, the Cosserat rotations are noticeable only in the shear zone. Fig. 11b illustrates the load–displacement plots for three
different length scales. As the figure shows, increasing length scale gives higher load responses. It is the key objective to relate the discrete level parameters to continuum parameters (stresses and strains). The finite element simulations in elastoplastic Cosserat elements enriched the classical continuum solutions by accommodating the microstructure to macro structure. 6. Footing punching problem The bearing capacity of a surface footing on granular soils has been a challenging problem approached experimentally, analytically and numerically by numerous
Fig. 8. Equivalent plastic strain distribution of slope failure (lc = 1.0 mm).
Fig. 9. Equivalent plastic strain distribution of slope failure (lc = 2.0 mm).
556
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Fig. 10. Velocity field.
Fig. 11a. Micro-rotation distribution of slope failure (lc = 2.0 mm).
Fig. 11b. Reaction force–displacement relation.
researchers. In the present example the behavior of granular soil under a rectangular footing a well-centered vertical displacement is investigated. Dimensions of the model are chosen large enough to eliminate the boundary condition effects. Bearing resistance is mobilized by applying equal displacement increments to the nodes beneath one half of the footing. As discussed by Arslan [27], displacement control provides numerical convenience, better stability, fewer iterations and represents physical reality (Fig. 12). The footing problem is solved for two different length scales. The first example is the finite element solution of the footing punching problem for the length scale, lc = 1.0 mm. Contours of resultant vertical stresses, shear
stresses and plastic strain viewed in plan are shown in Figs. 13–15, 16a, 16b, 17–19. The vertical stress distribution is uniform over most of the footing (Fig. 13). However, due to singularity at the footing edge, a high value of shear stress was predicted near the edge. The shear stress distribution is given in Fig. 14. As the figure shows, shear stress is localized at the corner of the footing. The granular foundation failure surface is not mobilized simultaneously, but in a progressive manner, so called ‘progressive failure’ [56]. The wedge Hill-type mechanism accompanying failure of the smooth footing is seen in Fig. 15. The deepest part of the failure mechanism observed for the smooth square footing takes place at 1.25 times the width of the footing.
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
557
Fig. 12. Footing punching problem.
Fig. 13. Vertical stress distributions inside the granular soil under footing.
Fig. 14. Shear stress distributions inside the granular soil under footing.
The kinematic mechanism at failure observed in the finite element analysis for rectangular footings is illustrated in Fig. 15, indicating a double wedge Hill-type mechanism (Fig. 20). Considering the geometry of the kinematic mechanisms accompanying failure in the finite element analyses
of rectangular footings various observations were made: plane strain conditions prevailed at approximately 1.25B from the ends of the footings (B is the width of the footing). Micro-rotation fields for this analysis are shown in Fig. 16a. The micro-rotational profile along the shear band
558
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Fig. 15. Equivalent plastic strain distributions inside the granular soil under footing.
Fig. 16a. Micro-rotation distributions inside the granular soil under footing.
Fig. 16b. Micro-rotational profile along the shear band.
is shown in Fig. 16b. As the figure shows, maximum microrotation occurs at the center of the shear band. Displacement vectors from the finite element analyses for a smooth square footing is shown in Fig. 17 and clearly indicates the different direction of soil movement directly
below the footing. As shown in the figures, localized plastic zones initiated at footing edges, and then further penetrated vertically in a narrow zone under the edge. The results of the numerical model showed the foundation soil failed with a localized shear failure mode. Unlike classical continuum, the incremental boundary value problem remains well posed in the micropolar solution and analysis can be continued until a collapse load analysis is attained or the vertical displacement of the footing exceeds a prescribed limit. Figs. 17 and 18 illustrate the finite element results of the footing problem for lc = 1.5 mm. The deepest part of the failure mechanism observed for the smooth square footing increased from 1.25 to 1.50 times the width of the footing. The different length scales provide a comparison to assess the particle size effect on the bearing capacity of the soil. The finite element results illustrates that the larger grain size leads to a deeper failure mechanism. 7. Length scale-size effect study on footing punching problem The effect of length scale on the bearing capacity of shallow foundations will be studied next. To demonstrate
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
559
Fig. 17. Velocity field.
Fig. 18. Vertical stress distributions inside the granular soil under footing.
Fig. 19. Equivalent plastic strain distributions inside the granular soil under footing.
significance of the length scale, three different cases are considered. The length scale effect on the analysis of different footing will be examined for three different scales and
three different footing dimensions (Fig. 21). In order to construct a length scale-bearing capacity relation, other parameters are fixed where only the length scales and the
560
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
Fig. 20. Double wedge Hill-type failure mechanism [56].
Fig. 21. Footing and domain dimensions for three different cases.
size of the domain are varied. Footing dimensions and the size of domain changed at the same ratio as shown in Fig. 21. The example can be explained as: A. Analyzing the size effect for a footing by changing the dimensions of footing having the same material properties and loading condition. The footing sizes studied are 1.0, 10.0, and 100.0 m. The ratio of soil media dimension to footing dimension remains constant at 20:1 for all three examples. B. Length scale effect for a footing with the same dimensions but different length scales. The results are presented as plots of reaction force-vertical displacement curves. Fig. 22 shows load displacement curves of the foundations due to vertical displacement controlled loading of the footing for three different length scales and footing width (L) equals to 1. A 12 mm punching displacement has been applied. As Fig. 22 illustrates, an increase of the length scale from 0.1 mm to 10 mm leads to a clear increase of stiffness. While peak reaction forces of the footing for the length scale of 0.1 is 520 kN/m, it is 680 and 800 kN/m for the length scales 1.0 and 10.0 mm respectively. Virtually the same effect of length scale is shown by the Figs. 23 and 24. As can be seen, length scale has a great effect on the stiffness and peak reaction forces of the foundation. Similar to the plane strain example, the increase of the size of the footing leads to a decrease of peak reaction
Fig. 22. Footing reaction force–displacement curve (L = 1.0 m.).
forces. This can be due to higher stress rotations and higher micro-rotations effect for the larger dimensions. The resultant reaction force of the footing increases first, shows a pronounced peak, drops later and reaches slowly a residual state. The maximum reaction force of the footing increases with increasing length scales. The material becomes more ductile after the peak with increasing length scale. A footing on soil produces a strain field, which varies from relatively high strains immediately beneath the footing to infinitesimal values at large distances beneath and laterally away from the footing. Thus, there is a higher strain gradients due to higher stress concentrations at the corner of footing. Micropolar continua enables to take into account the
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
561
The absence of intrinsic length scale leaves the width of localized zone unspecified. Thus, the finite element solution of the classical continuum is mesh sensitive. To solve the mesh sensitivity and to investigate the effect of microstructure, Cosserat continuum theory is implemented into Drucker Prager criterion and subsequently used in finite element analysis. The finite element analyses were carried out using the commercial nonlinear finite element code ABAQUS. The Cosserat element, with the additional degrees of freedom was implemented in ABAQUS using the UEL interface. This formulation was used in order to analyze selected geotechnical problems or configuration. Some of the major findings of this research are: Fig. 23. Footing reaction force–displacement curve (L = 10.0 m.).
The numerical model was used to predict the load–displacement behavior of a shallow footing and failure behavior of a slope. A Hill-type failure mechanism accompanying failure of the smooth footing was observed. The change of length scale affects the failure mechanism. The bearing capacity (maximum reaction force) of the footing increases with increasing length scales. The material becomes more ductile after the peak with increasing length scale. The thickness of the localized zone in the slope failure increases with increasing length scale. FE calculations show that localization in granular materials can realistically be determined with Cosserat elastoplastic constitutive law. Acknowledgements
Fig. 24. Footing reaction force–displacement curve (L = 100 m.).
higher strain gradient and effect of length scale would be observed clearly under such a higher strain gradient. The bearing capacity of soil to support a footing depends on the properties of soil such as density shear strength also mechanical properties such as particle sizes. Length scale has significantly affected the calculated bearing capacity of a foundation. This analysis clearly demonstrates the effect of the micropolar approach in the numerical simulation of bearing capacity and the significance of length scale in the analysis of bearing capacity. It must be noted that the mesh independent results can be obtained within enhanced continuum if the mesh size length scale ratio is small. For large problems, however, it is not possible without remeshing to achieve this due to the limited calculation time. On the other side, if too large length scale is artificially assumed in calculations, the FE results would not be realistic since both strength and ductility change with length scale. 8. Conclusions The basic deficiency of a classical continuum with regard to microstructure phenomena is the lack of a length scale.
The authors would like to thank Prof. Kaspar Willam for his guidance and help during the study. The comments and suggestions of anonymous reviewers that helped to improve the quality of the paper are also greatly appreciated. The authors gratefully acknowledge the financial support provided by NASA under Contract No. NCC8-242. References [1] Hadamard J. Propagation des Ondes et les Equations d’Hydrodynamique. New York: Chelsea; 1949 [reprinted 1903]. [2] Thomas TY. Plastic flow and fracture in solids. New York (NY): Academic Press; 1961. [3] Hill R. Acceleration waves in solids. J Mech Phys Solids 1962;10: 1–16. [4] Mandel J. Conditions de Stabilite et Postulat de Drucker. In: Kravtchenko, Siireys, editors. Proceedings of the IUTAM symposium on rheology and soil mechanics. Berlin: Springer; 1964. p. 58–68. [5] Rice J. The localization of deformation. In: Koiter W, editor. Theoretical and applied mechanics. North Publishing Company; 1976. [6] Rudnicki JW, Rice JR. Conditions for the localization of deformation in pressure sensitive dilatant materials. J Mech Phys Solids 1975;23:371–94. [7] Rice JR, Rudnicki JW. A note on some features of the theory of localization of deformation. Int J Solids Struct 1980;16:597–605. [8] Vardoulakis I. Equilibrium bifurcation of granular earth bodies. Adv Anal Geotech Instabilities 1978;13:65–120.
562
H. Arslan, S. Sture / Computers and Geotechnics 35 (2008) 548–562
[9] Vardoulakis I. Bifurcation analysis of the triaxial test on sand samples. Acta Mech 1979;32:35–54. [10] Vardoulakis I. Shear band inclination and shear modulus of sand in biaxial tests. Int J Numer Anal Meth Geomech 1980;4:103–19. [11] Vardoulakis I. Constitutive properties of dry sand observable in the triaxial test. Acta Mech 1981;38:219–39. [12] Vardoulakis I. Shear-band and liquefaction in granular materials on the basis of Cosserat continuum theory. Ingenieur-Archiv 1989;59: 106–13. [13] Molenkamp F. Comparison of frictional material models with respect to shear band initiation. Geotechnique 1985;35:127–43. [14] Peric D. Localized deformation and failure analysis of pressure sensitive granular materials. Ph.D. dissertation, University of Colorado-Boulder; 1990. [15] Hill R. A general theory of uniquness and stability in elastic–plastic solids. J Mech Phys Solids 1958;6:36–249. [16] Hill R, Hutcinson KW. Bifurcation phenomena in the plane strain test. J Mech Phys Solids 1975;23:239–64. [17] Arthur J, Dunstan T, Al-Ani Q, Asadi A. Plastic deformation and failure in granular media. Geotechnique 1977;27:53–74. [18] Desrues J. Localization de la deformation plastique dans les materiaux granulaires. These d’etat, Universite’ de Grenoble; 1984. p. 180–85. [19] Tejchman J. Modelling of shear localisation and autogeneous dynamic effects in granular bodies. Publication series of the Institute of Soil and Rock Mechanics, vol. 140. University Karlsruhe; 1997. [20] Tejchman J, Herle I, Wehr J. FE-Studies on the influence of initial void ratio, pressure level, and mean grain diameter on shear localization. Int J Numer Anal Meth Geomech 1999;23:2045–74. [21] Maier T. Nonlocal modeling of softening in hypoplasticity. Comput Geotech 2003;30:599–610. [22] Manzari MT. Application of micropolar plasticity to post failure analysis in geomechanics. Int J Numer Anal Meth Geomech 2004;28:1011–32. [23] Voyiadjis JZ, Alsaleh MI, Alshibli KA. Evolving internal length scales in plastic strain localization for granular materials. Int J Plasticity 2005;21:2000–24. [24] Alsaleh MI, Voyiadjis GZ, Alshibli KA. Modeling strain localization in granular materials using micropolar theory: mathematical formulations. Int J Numer Anal Meth Geomech 2006;30(15):1501–24. [25] Alsaleh MI, Alshibli KA, Voyiadjis GZ. Influence of micro-material heterogeneity on strain localization in granular materials. ASCE: Int J Geomech 2006;6(4):248–59. [26] Alshibli KA, Alsaleh MI, Voyiadjis GZ. Modeling strain localization in granular materials using micropolar theory: numerical implementation and verification. Int J Numer Anal Meth Geomech 2006;30(15):1525–44. [27] Arslan H. Localization analysis of granular materials in Cosserat Elastoplasticity –Formulation and finite element implementation. PhD dissertation, University of Colorado-Boulder; 2006. [28] Coulomb C.A. Essai sur une application des re`gles des maximis et minimis a` quelques proble`mes de statique relatifs a` l’architecture. Me`m. Acad. Roy. Pre`s. Divers Savants 7, Paris; 1776. [29] Roscoe KH. Tenth rankine lecture: the influence of strains in soil mechanics. Geotechnique 1970;20:129–70. [30] Alshibli K.A. Localized deformations in granular materials. Ph.D. dissertation, University of Colorado at Boulder, CO; 1995. [31] Vermeer P.A. A simple shear-band analysis using compliances. In: IUTAM conference on deformation and failure of granular materials, Delft; 1982. p. 493–9. [32] Ottosen N, Runesson K. Acceleration waves in elastoplasticity. Int J Solids Struct 1991;28:135–59. [33] Raniecki B, Bruhns OT. Bounds to bifurcation stresses in solids with non-associated flow law at finite strain. J Mech Phys Solids 1981;29:153–72.
[34] Ottosen NS, Runesson K, Peric D. Discontinuous bifurcations of elastic–plastic solutions at plane stress and plane strain. Int J Plasticity 1991;7:99–121. [35] Runesson K, Mroz Z. A note on non associated plastic flow rules. Int J Plasticity 1989;5:639–58. [36] Ottosen N, Runesson K. Properties of discontinuous bifurcation solutions in elasto-plasticity. Int J Solids Struct 1991;27:401–21. [37] Jira´sek M. Modeling of localized inelastic deformation, Short Course LID, Prague; 19–23 September 2005. [38] Chakrabarty J. Theory of plasticity. New York: McGraw-Hill; 1987. [39] Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by softening relations in rate-independent solids. Comp Mech 1993;12:277–96. [40] Oliver J. Continuum modelling of strong discontinuities in solid mechanics. In: Owen DRJ, Onate E, editors. Computational plasticity. Fundamentals and applications, vol. 1. Pineridge Press; 1995. p. 455–79. [41] Oliver J. Continuum modelling of strong discontinuities in solid mechanics using damage models. Comput Mech 1995;17(1/2): 49–61. [42] Larsson R, Runesson K, Sture S. Embedded localization band in undrained soil based on regularized strong discontinuity theory and finite element analysis. Int J Solids Struct 1996;33(20–22):3081–101. [43] Steinmann P. A finite element formulation for strong discontinuities in fluid saturated porous media. Mech Cohes-Frict Mater 1999;4:133–52. [44] Willam K et al. Localization in micropolar continua, continuum models for materials with microstructure. Wiley; 1995. p. 297–339. [45] Iordache MM, Willam K. Localized failure analysis in elastoplastic Cosserat continua. Comput Meth Appl Mech Eng 1998;151: 559–86. [46] Cosserat EF. Theorie des corps deformable. Paris: Hermann; 1909. [47] Muhlhaus HB, Vardoulakis I. The thickness of shear bands in granular materials. Geotechnique 1987;37(3):71–283. [48] Muhlhaus HB. Application of Cosserat theory in numerical solutions of limit load problems. Ing Arch 1989;59:124–37. [49] de Borst R. Simulation of strain localization: a reappraisal of the Cosserat continuum. Eng Comput 1991;8:317–32. [50] Eringen AC, Suhubi ES. Nonlinear theory of simple micro-elastic solids – I. Int J Eng Sci 1964;2:189–203. [51] Eringen AC. A unified theory of thermomechnical materials. Int J Eng Sci 1966;4:179. [52] Eringen CA. Microcontinuum field theories I: foundations and solids. New York: Springer-Verlag, Inc; 1999. [53] Forest S, Barbe F, Cailletaud G. Cosserat modelling of size effects in the mechanical behaviour of polycrystals and multiphase materials. Int J Solids Struct 2000;37:7105–26. [54] Walsh SDC, Tordesillas A. A thermomechanical approach to the development of micropolar constitutive models of granular media. Acta Mech 2004;167:145–69. [55] Arslan H, Sture S, Willam K. Analytical and geometrical representation of localization in granular materials. Acta Mech 2007;194(1–4): 159–73. [56] Meyerhof GG. The ultimate bearing capacity of foundations. Geotechnique 1951;2(4):301–32. [57] ABAQUS Version 6.5. Manuals, Hibbitt, Karlsson and Sorensen Inc; 2004. [58] Arslan H, Sture S. Evaluation of a physical length scale for granular materials. Comput Mater Sci. doi:10.1016/j.commatsci.2007.08.016. [59] Arslan H, Sture S. Finite element analysis of localization and micromacro structure relation in granular materials: Part I: formulation. Acta Mech [in press]. [60] Arslan H, Sture S. Finite element analysis of localization and micromacro structure relation in granular materials: Part II: implementation and simulations. Acta Mech [in press].