- Email: [email protected]
Modelling of diffuse failure mechanisms of catastrophic landslides
Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 www.elsevier.com/locate/cma
Modelling of diffuse failure mechanisms of catastrophic landslides J.A. Fernandez Merodo a,b,*, M. Pastor a,b, P. Mira a,b, L. Tonni c, M.I. Herreros a,b, E. Gonzalez a,b, R. Tamagnini d a Centro de Estudios y Experimentacion de Obras Publicas, Alfonso XII 3, 28014 Madrid, Spain M2 i (Math. Model. Eng. Group), Department of Applied Mathematics, ETS de Ingenieros de Caminos, UPM, Ciudad Universitaria, 28040 Madrid, Spain c Dept. DISTART, Universita di Bologna, Viale Risorgimento 2, 40136 Bologna, Italy Dept. Structural and Geotechnical Engineering, Universita di Roma ‘La Sapienza’, Via Montedoro 28, 00100 Roma, Italy b
d
Received 15 April 2003; received in revised form 15 July 2003; accepted 1 September 2003
Abstract This paper presents a theoretical and numerical framework to model the initiation mechanisms of catastrophic landslides. Attention is paid to the mechanism of failure referred to as diffuse, characteristic of soils presenting very loose or metastable structures, as in the case of the liquefaction of a sand layer induced by an earthquake, where the effective stress approaches zero as the pore pressure increases. The equations describing the coupling between the solid skeleton and the pore fluids are presented following an Eulerian approach based on mixture theory that can provide a unified formulation for both initiation and propagation phases. The system of partial differential equations is then discretized using the classical Galerkin finite element method and neglecting the convective terms. A simple improvement of the Generalized Plasticity model is introduced to deal with bonded geomaterials and to represent the mechanical collapse of such materials. The use of an implicit integration scheme guarantees the efficiency of the numerical convergence. 2004 Elsevier B.V. All rights reserved. Keywords: Landslides; Diffuse failure; Liquefaction; Bonded geomaterials; Finite element analysis; Implicit integration
1. Introduction Landslides are one of natural catastrophes causing important losses of human lives and extreme damage to property. There is a wide variety of types of landslides, depending on the materials involved and the triggering mechanism. The time scale involved can range from minutes to years. For a classification, it is worth consulting Ref. [1] where a detailed description is provided.
*
Corresponding author. Address: Centro de Estudios y Experimentaci on de Obras Publicas, Alfonso XII 3, 28014 Madrid, Spain. E-mail address: [email protected] (J.A. Fernandez Merodo).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.09.016
2912
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Landslides and failure of slopes are caused by changes in the effective stresses, variation of material properties or changes in the geometry. Changes of effective stresses can be induced either directly, as a consequence of variation of the external forces (earthquakes, human action), or indirectly through pore pressures (rainfall effects). Variations in material properties can be caused by processes of degradation (weathering and chemical attack). Finally, geometry can change because of natural causes (erosion) or human action (excavation, construction, reshaping. . .). In some cases, the failure mechanism consists of a clearly defined surface where shear strain concentrates. From a mathematical point of view, the inception of a shear band is characterized by a discontinuity in the strain field, which can evolve towards a discontinuity in the displacements at a later stage. Much effort has been devoted during the past years to better understand this phenomenon. The problem is ill posed for elastoplastic materials, and the results obtained in numerical models depend on the mesh size and alignment. The interested reader is addressed to Ref. [2] or [3] and [4] where a detailed explanation is provided. If the soil has a low relative density, and tends to compact when sheared, failure mechanism can be different from the localized mode. It affects not only the material on a narrow band or failure surface, but a much larger mass of soil. Failure can be described as of a flow type. Pore pressures developed may cause liquefaction of the soil, and liquefaction––or conditions close to it––may have played an important role on failure of slopes such as those of Aberfan in South Wales, Anhui in China, Jupille in Belgium and Rocky Mountains coal mine dumps. This type of mechanism has been described as diffuse [5–7]. The key to properly model this type of landslide is the ability of the constitutive model to reproduce liquefaction. After presenting the mathematical and the numerical model in a general mixture framework, we present an enhanced Generalized Plasticity model that can reproduce the mechanical behaviour of bonded geomaterials. This model is not only able to characterize the basic features of sand behaviour under monotonic and cyclic loading but also to represent the mechanical collapse of bonded materials. An implicit integration of such kind of Generalized Plasticity model is presented to assure the effectiveness of the numerical convergence. Finally, diffuse mechanisms of failure are illustrated with some examples.
2. Mathematical model for soil skeleton–pore fluid coupling 2.1. Introduction Soils and rocks are geomaterials with voids which can be filled with water, air, and other fluids. They are, therefore, multiphase materials, exhibiting a mechanical behaviour governed by the coupling between all the phases. Pore pressures of fluids filling the voids play a paramount role in the behaviour of a soil structure, and indeed, their variations can induce failure. The first mathematical model describing the coupling between solid and fluid phases was proposed by Biot [8,9] for linear elastic materials. This work was followed by further development at Swansea University, where Zienkiewicz and co-workers [10–14] extended the theory to non-linear materials and large deformation problems. It is also worth mentioning the work of Lewis and Schrefler [15], Coussy [16] and de Boer [17]. Among the different alternative ways which can be used to describe the coupling between solid skeleton and pore fluids, we have chosen an approach closer to mixture theories than to the more classical approach used in computational geotechnics, and it provides a more general description which can be used not only for initiation of failure but also for propagation of catastrophic landslides. The soil skeleton consists of particles of density qs having a porosity n (volume percent of voids in the mixture). The voids can be filled with air and water, but in some cases there will be mixtures of water and
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2913
very fine particles which can be considered as a fluid phase. We will assume that fluid phases ðaÞ are not miscible, defining Sa as the degree of saturation for the phase a, i.e., the volume fraction of voids occupied by it. In the case of air and water, we will define Sw and Sa as Sw ¼
volume of water ; volume of voids
Sa ¼
volume of air volume of voids
ð1Þ
with Sw þ Ss ¼ 1: We will introduce qðaÞ to denote the density of the phase ðaÞ: If the fluid density is qa ; qðaÞ ¼ Sa nqa ;
ð2Þ
while for the solid we can write qðsÞ ¼ ð1 nÞqs :
ð3Þ
Density of the mixture is therefore obtained by adding the densities of all constituents in the mixture. Concerning the kinematics, we will use u for displacements and v for velocities, with superindexes referring to the specific constituents. In this way, uðsÞ will be the displacement field of the solid phase, and vðaÞ the velocity of pore fluid ðaÞ: In most geotechnical applications, the movement of the fluid phases is described by the velocity relative to the soil skeleton. The so-called Darcy’s velocity of phase ðaÞ is defined as wðaÞ ¼ nSa ðvðaÞ vðsÞ Þ:
ð4Þ
2.2. Effective and partial stresses If we denote by rs and ra the Cauchy stress tensors acting on the solid particles and on the phase ðaÞ, we can define the partial stresses rðsÞ and rðaÞ as rðsÞ ¼ ð1 nÞrs ;
ð5Þ
rðaÞ ¼ nSa ra : The total stress tensor acting on the mixture is obtained as r ¼ rðsÞ þ
nX phases
rðaÞ :
ð6Þ
a¼1
The partial stresses rðaÞ can be decomposed into hydrostatic and deviatoric components as rðaÞ ¼ nSa pa I þ nSa sa , where sa ¼ devðra Þ is the deviatoric part of ra , and I is the identity tensor of second order. It is important to note that if the saturating fluid is water, the deviatoric stresses can be neglected, but in other cases where viscous contributions are important, this term has to be taken into account. It is convenient to introduce an averaged pore pressure p as p¼
nX phases a¼1
pðaÞ ¼
nX phases a¼1
Sa pa :
ð7Þ
2914
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Two important particular cases are those of (i) non-saturated soils, and (ii) saturated soils, where pore fluids are (i) water and air, and (ii) water. The above relations specialize to r ¼ rðsÞ þ rðwÞ þ rðaÞ ; rðsÞ ¼ ð1 nÞrs ; rðwÞ ¼ nSw pw I;
ð8Þ
rðaÞ ¼ nSa pa I; p ¼ Sw pw þ Sa pa and r ¼ rðsÞ þ rðwÞ ; rðsÞ ¼ ð1 nÞrs ;
ð9Þ
rðwÞ ¼ npw I; p ¼ pw ; respectively. Another important case is that of non-saturated soils having pa ¼ 0; for which r ¼ rðsÞ þ rðwÞ ; rðsÞ ¼ ð1 nÞrs ;
ð10Þ
rðwÞ ¼ nSw pw I; p ¼ Sw pw :
This approach is valid for many engineering cases of interest. A more complete description can be found in Ref. [15]. The effective stress r0 is defined as p: r0 ¼ r þ I
ð11Þ
In the case of viscous pore fluids where deviatoric stresses are important, the effective stress is given by nX phases r0 ¼ r þ I pn sðaÞ : ð12Þ a¼1
Another point worth mentioning is that the effective stress may be modified to account for compressibility of soil grains relative to that of solid skeleton. Here we will make the assumption that solid skeleton is incompressible. 2.3. Balance of mass We will introduce the material derivative of a magnitude / following phase a as DðaÞ / o/ ¼ þ vðaÞ $/: Dt ot
ð13Þ
One important relation which will be used later is DðaÞ / DðsÞ / ¼ þ ðvðaÞ vðsÞ Þ $/ Dt Dt DðsÞ / wðaÞ þ ¼ $/; Dt nSa where we have introduced the Darcy’s velocity of phase ðaÞ relative to soil skeleton defined in (4).
ð14Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2915
The balance of mass equation for the solid phase is given by DðsÞ qðsÞ þ qðsÞ $ vðsÞ ¼ 0 Dt
ð15Þ
and for the fluid phases ðaÞ DðaÞ qðaÞ þ qðaÞ $ vðaÞ ¼ 0; Dt
ð16Þ
where qðsÞ ¼ ð1 nÞqs and qðaÞ ¼ nSa qa : In the case of non-saturated soils, we can write DðsÞ fð1 nÞqs g þ ð1 nÞqs $ vðsÞ ¼ 0; Dt DðaÞ ðnSa qa Þ þ nSa qa $ vðaÞ ¼ 0; Dt
ð17Þ
with a ¼ fa; wg. From here, we obtain DðsÞ n 1 n DðsÞ qs ¼ þ ð1 nÞ$ vðsÞ ; Dt qs Dt DðaÞ n n DðaÞ ¼ ðSa qa Þ þ n$ vðaÞ ; Dt Sa qa Dt If we now add both equations and use (14) we arrive at ðaÞ w 1 n DðsÞ qs n DðaÞ $ þ ðSa qa Þ þ $ vðsÞ ¼ 0 þ qs Sa qa Dt Sa Dt
ð18Þ
ð19Þ
for a ¼ fa; wg. If we neglect now thermal effects, we can elaborate further the second terms as follows: p 1 n DðsÞ qs 1 n DðsÞ : ¼ Dt qs Ks Dt In the case of a non-saturated soil with pa ¼ 0 and p ¼ Sw pw we obtain ( ) 1 n DðsÞ DðsÞ pw DðsÞ Sw p 1n ¼ Sw þ pw : Ks Ks Dt Dt Dt The third term (water phase) can be obtained for non-saturated soils with pa ¼ 0 as ( ) n DðwÞ n DðwÞ pw n DðwÞ Sw ðSw qw Þ ¼ þ : Sw qw Dt Kw Dt Sw Dt Therefore, the balance of mass for the water (non-saturated, pa ¼ 0) can be written as w n DðwÞ pw ð1 nÞ DðsÞ pw n DðwÞ Sw 1 n DðsÞ Sw þ þ þ þ $ vðsÞ ¼ 0: pw $ þ Sw Kw Dt Ks Sw Dt Ks Dt Dt
ð20Þ
ð21Þ
ð22Þ
ð23Þ
Sometimes, it is introduced for convenience the specific storage coefficient of the soil Cs as Cs ¼ n
oSw : opw
ð24Þ
2916
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
From here we can obtain some particular cases of interest: • Saturated soil $wþ
n DðwÞ pw ð1 nÞ DðsÞ pw þ þ $ vðsÞ ¼ 0: Kw Dt Ks Dt
ð25Þ
This equation can be further simplificated by taking into account the relationship between material derivatives following the fluid and the solid phases, Eq. (14) $wþ
n DðsÞ pw ð1 nÞ DðsÞ pw w þ þ $ vðsÞ þ $pw ¼ 0: Kw Dt Ks Kw Dt
The last term is usually neglected in geotechnical analysis, and above equation reduces to $wþ where 1 ¼ Q
1 DðsÞ pw þ $ vðsÞ ¼ 0; Q Dt
n ð1 nÞ þ : Kw Ks
ð26Þ
ð27Þ
• Flow of dry geomaterials In some circumstances, there may exist an important coupling between the pore air and the soil skeleton. Should this be the case, the balance of mass equation for the pore gas can be written in a similar way as $wþ with 1 ¼ Qa
1 DðsÞ pa þ $ vðsÞ ¼ 0 Qa Dt
ð28Þ
n ð1 nÞ þ ; Ka Ks
ð29Þ
where have introduced the volumetric stiffness of the air Ka . The balance of mass equation of the mixture can be obtained by adding the equations for all constituents: nX phases nX phases DðsÞ qðsÞ DðaÞ qðaÞ þ qðsÞ $ vðsÞ þ þ qðaÞ $ vðaÞ ¼ 0 Dt Dt a¼1 a¼1
from where, using (14) we get nX phases ðaÞ DðsÞ q w ðsÞ þ q$ v þ $qðaÞ ¼ 0: Dt nS a a¼1
ð30Þ
If all wðaÞ are negligible, above equation reduces to DðsÞ q þ q$ vðsÞ ¼ 0: Dt
ð31Þ
A particular case of interest is that of very fluid slurries, for which the mixture density can be assumed to be constant. The balance of mass reduces to $ vðsÞ ¼ 0
ð32Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2917
2.4. Balance of momentum of pore fluid The balance of momentum for a fluid phase ðaÞ reads qðaÞ
DðaÞ vðaÞ ¼ qðaÞ b þ div rðaÞ þ Ra ; Dt
ð33Þ
where b is the body forces vector, and Ra is a force representing the interaction between the solid skeleton and the pore fluid ðaÞ; which can be assumed to follow Darcy’s law Ra ¼ ka1 wðaÞ :
ð34Þ
In above, ka is the permeability matrix for phase ðaÞ, ka ¼
karel Kaintr la
ð35Þ
where karel is a relative permeability which depends on Sa , la is the dynamic viscosity of fluid ðaÞ, and Kaintr is the intrinsic permeability matrix which has dimensions of L2 . If the material is isotropic, then the matrix reduces to a constant. The permeability ka has dimensions of L3 T =M. It should be noted that in geotechnical analysis the permeability is often defined in a slightly different way, a ¼ qa bka k with dimensions of LT 1 : If we take into account the relation (14), we can write the balance of momentum as ! DðsÞ vðsÞ ðaÞ q þ c:t: ¼ qðaÞ b þ div rðaÞ ka1 wðaÞ ; Dt where the correcting terms c.t. are given by ðaÞ DðsÞ wðaÞ wðaÞ ðsÞ wðaÞ w $v þ $ þ : c:t: ¼ Dt nSa nSa nSa nSa In many geotechnical problems the correcting terms can be neglected , and wðaÞ is obtained as ! ðsÞ ðsÞ ðaÞ ðaÞ D v ðaÞ ðaÞ w ¼ ka q þ q b þ div r : Dt
ð36Þ
ð37Þ
ð38Þ
ð39Þ
This situation arises when (i) acceleration of fluid phases are small, and (ii) when cross products of wðaÞ and vðsÞ can be neglected. The interested reader is referred to the work of Zienkiewicz et al. [10] where the reasons are detailed. If the soil is saturated by water, and no viscous effects are present in the pore water, we can write ! DðsÞ vðsÞ w ¼ k w qw þ qw b $pw ð40Þ Dt and in the case of dry soils we will have ! DðsÞ vðsÞ þ qa b $pa : w ¼ k a qa Dt
ð41Þ
2918
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2.5. Balance of momentum of the mixture We obtain the balance of linear momentum for the mixture by combining the balance equations of the solid and fluid phases. Concerning the solid, we can write qðsÞ
DðsÞ vðsÞ ¼ qðsÞ b þ div rðsÞ þ Rs ; Dt
ð42Þ
where Rs stands for the drag forces due to all fluid constituents Rs ¼
nX phases
Ra :
a¼1
Eq. (42) is now added to (37), obtaining the balance of momentum for the mixture as ( ðaÞ ) nX phases DðsÞ vðsÞ DðsÞ wðaÞ wðaÞ ðsÞ wðaÞ w þ nSa qa $v þ $ q þ ¼ qb þ div r: Dt Dt nSa nSa nSa nSa a¼1
ð43Þ
This equation can be simplified as q
DðsÞ vðsÞ ¼ qb þ div r; Dt
ð44Þ
which is the form usually encountered in geotechnics.
3. Some useful forms in computational soil mechanics: the u–pw formulation In the preceding section, we have presented the balance of mass and momentum equations for a mixture of solid particles and nphases immiscible fluid phases occupying the pores, with some possible simplifications. The equations have been cast in a Eulerian formulation, from which simple Lagrangian forms involving small deformations can be immediately derived by neglecting the convective components of the material derivatives. The approach we propose in this work consists on modelling catastrophic landslides using a small deformation approach for the initiation phase. Therefore, we will not consider convective terms, and we will further assume that the accelerations of the pore fluids relative to the soil skeleton can be neglected. The validity of this assumption has been studied by Bettess and Zienkiewicz [10,14]. The two cases we will study next are those of (i) non-saturated soils with air at atmospheric pressure, and (ii) saturated soils. We will use displacements of the soil skeleton u instead of velocities. In the former case, the system of equations describing the coupling are the following: • Balance of mass for pore water $ w þ Sw
oev opw 1 opw þ Cs þ ¼ 0; Q ot ot ot
ð45Þ
where ev is the volumetric strain $ vðsÞ ¼ and
oev ot
ð46Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
1 ¼ Q
n
Sw 1 n Cs Sw Sw þ pw þ : Ks Kw n
• Balance of linear momentum (pore water) 2 ou o w $pw þ qw b kw1 w qw þ ¼ 0: ot2 ot nSw • Balance of linear momentum (mixture) o2 u o w div r þ qb q 2 nqw Sw ¼ 0: ot ot nSw
2919
ð47Þ
ð48Þ
ð49Þ
• Constitutive equation dr0 ¼ D : de: • Kinematics 1 oui ouj deij ¼ þ S du: 2 oxj oxi
ð50Þ
ð51Þ
Under certain conditions [10], we can neglect the relative accelerations of the fluid phase, and eliminate the relative displacements of the fluid phase w. This can be done by substituting the value of $ w obtained from (48) o2 u $ w ¼ $ kw $pw þ qw b qw 2 ð52Þ ot into the balance of mass of the fluid phase (45) from where we obtain o2 u 1 opw oev þ Sw ¼ 0: þ Cs þ $ kw $pw þ qw b qw 2 ot Q ot ot
ð53Þ
The main advantage of this formulation is than the field variables reduce to two, displacements u and pressures pw . Therefore, the basic equations of this u–pw formulation are (49) and (53), together with the constitutive equation (50) and the kinematic relations between strains and displacements. The simplification for saturated soils is immediate. Eq. (53) reduces to o2 u 1 opw oev þ ¼ 0: ð54Þ $ kw $pw þ qw b qw 2 þ ot Q ot ot The equations presented so far have to be complemented by suitable boundary and initial conditions.
4. Numerical model: discretization of the u–pw model The system of partial differential equations can be discretized using standard Galerkin techniques, as described in [12]. After approximating the fields u and pw as: u ¼ Nu u, pw ¼ Np pw , it results in two ordinary differential equations:
2920
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
d2 u M 2þ dt
Z
BT r0 dX Q pw f u ¼ 0;
ð55Þ
X
where B ¼ S Nu , and QT
d u d p þ H pw þ C w f p ¼ 0: dt dt
The matrices given above are defined as Z M¼ qNTu Nu dX; ZX 1 T C¼ Np Np dX; X Q Z Q¼ Sw aBT mNp dX; X Z H ¼ $NTp kw $Np dX
ð56Þ
ð57Þ
and Z
NTu b dX þ
Z
NTu t dC; Z Z Z Z op dC þ $NTp kw qw b dX fp ¼ NTp kw $Np kw Nu dX €u NTp s0 dX: on Cq X X X
fu ¼
X
Ct
ð58Þ
The term including accelerations in f p is usually disregarded (Ref. Thesis Chan [18]). The time derivatives of u and pw are approximated in a typical step of computation using the Generalized Newmark GN22 scheme for displacements and a GN11 for the water pressure [19,20]. If we introduce the following notation [19]: unþ1 € un ; D€ un ¼ € Dp_ wn ¼ p_ wnþ1 p_ wn ; u_ nþ1 ¼ u_ p;nþ1 þ b1 DtD€ un ; unþ1 ¼ up;n þ 12b2 Dt2 D€ un ; pwnþ1
¼
pwp;nþ1
þ
ð59Þ
hDtDp_ wn ;
where un ; u_ p;nþ1 ¼ u_ n þ Dt€ up;nþ1 ¼ un þ Dtu_ n þ 12Dt2 € un ; pwp;nþ1
¼
pwn
þ
we obtain the discretized system of equations valid in each time step: Z 0 p_ nw Funþ1 Uu ¼ 0; MD€ un þ BT r nþ1 dX hDtQD X
ð60Þ
Dtp_ wn ;
ð61Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
p_ nw Fpnþ1 Up ¼ 0; b1 DtQT D€ un þ ðDthH þ CÞD
2921
ð62Þ
€ n and D p_ nw . where the unknown values are Du If this system is non-linear, it can be solved by using a Newton Raphson method with a suitable Jacobian matrix: 2
oUu 6 oD€ u 6 6 4 oUp oD€ u
3 oUu ðiÞ " #ðiþ1Þ p_ 7 oD Uu ðiÞ uÞ 7 dðD€ ¼ : 7 Up p_ Þ oUp 5 dðD p_ oD
ð63Þ
Using Eqs. (61) and (62) we can write the above step as "
M þ 12 Dt2 b2 KT
hDtQ
b1 DtQT
DthH þ S
#ðiÞ "
dðD€ uÞ p_ Þ dðD
#ðiþ1Þ
where KT is the tangent stiffness matrix. KT ¼
¼
R
Uu
ðiÞ
Up
;
ð64Þ
BT Dep B dX.
4.1. A note on element type The proper choice of the element type is of paramount importance as some elements introduce errors resulting on (i) unrealistic limit loads and (ii) spurious failure elements. These phenomena have been the subject of an extensive research during the past years [21–23]. Based on the authors experience, the use of enhanced strain elements together with a suitable stabilization technique provides optimal results [24]. Alternatively, reduced integration with control of spurious modes may be used.
5. An enhanced generalized plasticity model for bonded geomaterials 5.1. Basic theory The elastoplastic behaviour of the material is described by the general incremental relationship, relating the strain increment de to the stress increment dr: dr ¼ D de;
ð65Þ
where the matrix D, elastoplastic stiffness matrix, depends not only on the present state conditions but on the stress–strain history and the loading direction. Taking into account this last dependency two different matrices are introduced DL and DU ; that correspond to the loading and unloading situations. A unit vector n is introduced to define the loading/unloading condition: drTe n > 0 loading; drTe n < 0 unloading; drTe n ¼ 0 neutral loading; where dre ¼ De de and De is an elastic matrix dependent only on the present state of the material.
ð66Þ
2922
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
1 The continuity between loading and unloading implies the general form of D1 L and DU :
ngL nT ; HL ngU nT ¼ D1 ; e þ HU
1 D1 L ¼ De þ
D1 U
ð67Þ
where ngL and ngU are arbitrary directions vectors and HL and HU are scalars defined as the loading and unloading plastic moduli respectively. The strain increment de can be decomposed into elastic and plastic parts as de ¼ dee þ dep
ð68Þ
with dee ¼ De dre ; ngL nT dr HL ngU nT dr dep ¼ HU
dep ¼
for loading;
ð69Þ
for unloading:
The values of the DL and DU matrices can be obtained by suitable manipulation, giving: DL ¼ De þ
De ngL nT De ; HL þ nTgL De n
DU ¼ De þ
De ngU nT De : HU þ nTgU De n
ð70Þ
A complete elastoplastic behaviour can thus be specified by prescribing at each point of the stress space (and for a given history) three directions (the loading direction n and the plastic flow directions ngL and ngU ), two scalars (the plastic moduli HL and HU ) and the elastic stiffness matrix De . It has to be remarked here that in the framework of the Generalized Plasticity, neither a yield surface nor a plastic potential have been introduced. It is also noted that plastic strains during unloading are permissible. Furthermore, it must be noted that the classical plasticity, the bounding surface models and the kinematic hardening models are particular cases of such a generalized plasticity. 5.2. The Pastor–Zienkiewicz model for sands In what follows, a sand model developed in the framework of the Generalized Plasticity by Pastor et al. [25] will be described. Using the conventional notation of soil mechanics, three stress invariants are defined: 1 p ¼ I1 ; p3ffiffiffiffiffiffiffi q ¼ 3J2 ; 1 h ¼ arcsin 3
! pffiffiffi 3 3 J3 ; 2 J23=2
ð71Þ
where I1 is the first invariant of the stress tensor and J2 and J3 are the second and third invariants of the deviatoric stress tensor (s ¼ r 13 I1 I).
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2923
Soil dilatancy dg has been expressed as linear function of the stress ratio g ¼ q=p0 : dg ¼ ð1 þ ag ÞðMg gÞ;
ð72Þ
where ag is a material constant and Mg is the slope of the critical state line which is related to the residual angle /0 by Mg ¼
6 sin /0 : 3 sin /0 sin 3h
ð73Þ
As a result, for loading stress increments the plastic flow direction vector is given by T
ngL ¼ ðnpgL ; nqgL ; nhgL Þ ;
ð74Þ
whose components are expressed as dg npgL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ dg2 1 nqgL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ dg2
ð75Þ
qMg cos 3h nhgL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 1 þ dg2 In the case of unloading, it is noted that irreversible strains are of contractive or densifying nature. Therefore the unloading plastic flow direction vector ngU is defined as ngU ¼ ðnpgU ; nqgU ; nhgU ÞT with
ð76Þ
npgU ¼ npgL ; ð77Þ
nqgU ¼ nqgL ; nhgU ¼ nhgL : T
Assuming a non-associated flow rule, the loading direction n ¼ ðnp ; nq ; nh Þ must be different from ngL , but with similar expressions for components: df np ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2 ; 1 þ df 1 nq ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2 ; 1 þ df qM g cos 3h nh ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2 ; 2 1 þ df
ð78Þ
df ¼ ð1 þ af ÞðMf gÞ
ð79Þ
where
and af and Mf are constitutive parameters.
2924
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
In order to take into account the main features of sand response––i.e. the existence of a critical state condition, dilative response after peak, liquefaction in loose sands, memory of previous stress path––Pastor et al. [25] proposed for the plastic modulus HL the following relationship: HL ¼ H0 p0 Hf ðHv þ Hs Þ HDM
ð80Þ
with 4 g 1 Hf ¼ 1 and gf ¼ 1 þ Mf ; gf af g Hv ¼ 1 ; Mg Z Hs ¼ b0 expðb1 nÞ and n ¼ jdnps j accumulated deviatoric plastic strain; HDM ¼
fmax f
cDM
0
and f ¼ p 1
af 1 þ af
g Mf
1=a mobilized stress function;
where H0 ; cDM ; b0 ; b1 are constitutive parameters. For the plastic modulus in unloading condition HU , according to the experimental results during undrained loading-unloading triaxial tests, Pastor et al. [25] proposed the following simple expression: c Mg Mg u > 1; HU ¼ Hu0 for gu gu ð81Þ Mg 6 1; HU ¼ Hu0 for gu where Hu0 , cu are material parameters. gu is the stress ratio from which unloading takes place. The material has a non-linear elastic response, according to the following relationships: dp ¼ Kt deev ;
dq ¼ 3Gt dees ;
ð82Þ
where the tangent bulk and shear modulus Kt and Gt are assumed to be dependent only on the hydrostatic part of the stress tensor: K ¼ K0
p0 ; p00
G ¼ G0
p0 p00
ð83Þ
with K0 , G0 and p00 references values. For the present model there are 12 material parameters requiring definition. Generally, all parameters are identified from monotonic and cyclic triaxial tests, though in certain cases some parameters are adopted from previous experiences if full test records are unavailable. This model is able to characterize the basic features of the sand behaviour under monotonic and cyclic loading. 5.3. An improvement of the model for bonded geomaterials An improvement of the model described in the preceding section can be introduced to reproduce the mechanical behaviour of bonded soils, weak rocks and other materials of a similar kind. Following the framework introduced by Gens and Nova [26] and Lagioia and Nova [27] two basic concepts lie the representation of this mechanical behaviour: the fundamental role played by yield pheno-
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
M
q
A pt
2925
B pc
C pc0
p’
Fig. 1. Successive yield surfaces for increasing degrees of bounding. Surface A corresponds to unbounded material.
mena and the need for considering the observed behaviour of the bonded material in relation with the behaviour of the equivalent unstructured one. A typical yield surface valid for the behaviour of the soil in an unbonded state, in the p0–q plane, is represented by curve A in Fig. 1. As the amount of bonding increases the field surface must grow towards the right to account for the fact that higher mean stresses can be applied to the material without causing it to yield. Furthermore, bonding also provides the sample with real cohesion and tensile strength that is reflected in the fact that yield surfaces are enlarged also towards the left of the stress diagram. Fig. 1 shows the successive yield surfaces for increasing degrees of bonding. Two parameters define the new enlarged yield locus: pc0 that controls the yielding of the bonded soil in isotropic compression and pt which is related to the cohesion and tensile strength of the material. Both pc0 and pt increase with the magnitude of bonding. We can assume that the degradation of the material (decrease in bonding) is related to some kind of damage measure, that will in turn depend on plastic strains. Lagioia and Nova [26] proposed simple laws to describe the debonding effect on a calcarenite material. The evolution of pt is governed by pt ¼ pt0 expðqt ed Þ;
ð84Þ
where pt0 and qt are two constitutive parameters and ed is the accumulated plastic volumetric strain. It appears reasonable to assume that changes of the yield locus will be controlled by two different phenomena: conventional plastic hardening (or softening) for an unbonded material and bond degradation. In that case, the plastic modulus introduced in 80 can be improved introducing Hb such as : HL ¼ ½H0 p Hb Hf ðHv þ Hs Þ HDM
New factors p
; Hf ; Hv
and
HDM
ð85Þ
must be defined by introducing the parameter pt :
p ¼ p0 þ pt ; 4 g Hf ¼ 1 ; gf g Hv ¼ 1 ; Mg cDM fmax ¼ ; HDM f g
1=a
af where g ¼ q=ðp0 þ pt Þ is the new stress ratio, f ¼ ðp0 þ pt Þ ½1 ð1þa Þ is the new mobilized stress f Mf function and fmax is the maximum value of the new mobilized stress function that can be replaced by the distance pc0 in Fig. 1. The term Hb is also controlled by the evolution of bonding. We propose a simple definition that depends only on the plastic volumetric strain:
2926
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Hb ¼ b1 ed expðb2 ed Þ:
ð86Þ
It can be seen that value of Hb decreases when the volumetric plastic strain increases (i.e. when debonding occurs) and in the limit case, when destructuration is complete, Hb becomes zero. In this case, the new plastic modulus defined in (85) coincides with the plastic modulus introduced in (80). To check the proposed model performance we will use the results of the laboratory tests of Lagioia and Nova on the Gravina calcarenite. Tables 1 and 2 list the 15 parameters necessary to characterize the constitutive model used in the following back prediction tests. The 4 constants in Table 2 correspond to the new parameters needed to represent the debonding phenomena. Fig. 2 compares the experimental data from Lagioa and Nova [26] with the model predictions for a specimen in an isotropic compression test. In this case the expression of the plastic modulus reduces to: cDM fmax HL ¼ ½H0 p Hb : ð87Þ f Figs. 3 and 4 compare experimental data and model predictions for specimens sheared in drained conditions under initial isotropic confining pressures of 1300 and 2000 kPa. The agreement is in general
Table 1 Pastor–Zienkiewicz model parameters for the Gravina calcarenite K (kPa)
G (kPa)
Mg
ag
Mf
af
H0
b0
b1
pc0 (kPa)
c
76,923
75,757
1.657
0.2
0.9
0.2
20
20
3000
2300
10
Table 2 New model parameters for the Gravina calcarenite b1
b2
pt0 (kPa)
qt
5e)6
36
142
1000
Experimental Calculated
1.15 1.1
void ratio
1.05 1 0.95 0.9 0.85 0.8 0
1000 2000 3000 mean effective stress p’ [kPa]
4000
Fig. 2. Isotropic compression test: experimental data from Lagioa and Nova [26] and calculated curves.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 4000
0
2927
Experimental Calculated
-0.05 volume strain
q [kPa]
3000
2000
-0.1 -0.15 -0.2
1000 -0.25 0
0
0.1 0.2 axial strain
0.3
0
0.1 0.2 axial strain
0.3
Fig. 3. Drained constant cell pressure test (r03 ¼ 1300 kPa): experimental data from Lagioa and Nova [26] and calculated curves.
5000
0
volume strain
q [kPa]
Calculated
-0.05
4000
3000
2000
1000
0
Experimental
-0.1 -0.15 -0.2 -0.25
0
0.1 0.2 axial strain
0.3
0
0.1 0.2 axial strain
0.3
Fig. 4. Drained constant cell pressure test (r03 ¼ 2000 kPa): experimental data from Lagioa and Nova [27] and calculated curves.
Experimental
1000
1000
Calculated
q [kPa]
1500
q [kPa]
1500
500
0
500
0
1000
2000 p’ [kPa]
3000
0
0
0.05 0.1 axial strain
0.15
Fig. 5. Oedometric test: experimental data from Lagioa and Nova [26] and calculated curves.
2928
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
good for both the stress–strain and the volumetric strain–axial strain curves. Fig. 5 shows the results for an oedometric test.
6. Implicit integration of generalized plasticity models In this section, a method of performing the fully implicit integration of a Generalized Plasticity model for sands is outlined. The algorithm presented here is general, since it has been devised for the integration of constitutive equations in the whole stress space [28,29]. The formulation involves an elastic predictor to give a first estimate of the stress change, coupled with an iterative sub-algorithm whose objective is to ensure that, on convergence, the constitutive behaviour is satisfied. In the following, the basic steps of the algorithm will be explained. Let us consider the finite element step n þ 1 and the incremental strain Denþ1 . The material is initially assumed to behave elastically and a trial stress state renþ1 is computed by adding the elastic stress increment associated with the incremental strain to the accumulated stress at the beginning of the increment: Z renþ1 ¼ rn þ Drenþ1 ¼ rn þ Des;nþ1 : Denþ1 ¼ rn þ Det;nþ1 : de; ð88Þ Denþ1
where Des;nþ1 is the elastic constitutive matrix. According to the non-linear elastic law proposed by Pastor et al. [25], the procedure assumes that the elastic moduli are not constant and accounts for their variation over the increment. Hence in Eq. (88) a secant elastic matrix has been used, which can be expressed just as a function of the secant bulk modulus. On the basis of the projection of the trial stress increment Drenþ1 onto the direction vector nn , it is possible to check whether loading, unloading or neutral loading is occurring. Under neutral loading ðnn : Drenþ1 ¼ 0Þ the material is assumed to have remained elastic at the Gauss point level. In that case the trial stress state computed in Eq. (88) is accepted as the final solution and there are neither any incremental plastic strains nor any incremental change in the hardening modulus. Under loading ðnn : Drenþ1 > 0Þ or unloading ðnn : Drenþ1 < 0Þ conditions, plastic deformations occur and an integration procedure must be adopted. Using an implicit scheme, the integration of the constitutive law results in a set of non-linear equations for the unknowns rnþ1 (stresses), HL=U ;nþ1 (plastic modulus), Dknþ1 (plastic multiplier) and Ks;nþ1 (bulk modulus): rnþ1 ¼ renþ1 Dknþ1 Des;nþ1 : ngL;nþ1 ;
ð89Þ
HL;nþ1 ¼ HL ðrnþ1 ; ep ðDknþ1 ÞÞ;
ð90Þ
Dknþ1 ¼
nnþ1 : Des;nþ1 : Denþ1 ; HL;nþ1 þ nnþ1 : Des;nþ1 : ngL;nþ1
Ks;nþ1 ¼ Ks ðpnþ1 ; Deenþ1 Þ;
ð91Þ ð92Þ
where pnþ1 is hydrostatic part of the stress tensor. To solve this non-linear local problem, a numerical method must be used following Crisfield [30]. It is advantageous to define a vector of residuals r:
ð93Þ rTnþ1 ¼ rT1;nþ1 r2;nþ1 r3;nþ1 r4;nþ1
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2929
in which r1;nþ1 , r2;nþ1 , r3;nþ1 and r4;nþ1 represent the difference between the current values of the abovementioned primary variables and the values obtained from the integration scheme. Iterations are then introduced to reduce residuals to zero, while the primary variables must satisfy the constitutive equations. After having collected the unknowns of the problem in vector v: vTnþ1 ¼ rTnþ1
HL;nþ1
Dknþ1
Ks;nþ1
ð94Þ
the iterative scheme can be expressed as follows: iþ1 iþ1 i riþ1 nþ1 ¼ rnþ1 þ Jr;v dvnþ1 ¼ 0;
ð95Þ
where the superscript i refers to the iteration number at the Gauss point level and Jr;v represents the Jacobian of the residual rnþ1 with respect to the unknowns. Solving the linear system (95) we obtain, at each iteration, the update of the primary variables, given by !1 orinþ1 iþ1 dvnþ1 ¼ rinþ1 : ð96Þ ovinþ1 Hence the unknowns of the problem can be simultaneously updated: i iþ1 riþ1 nþ1 ¼ rnþ1 þ drnþ1 ; iþ1 i iþ1 ¼ HL;nþ1 þ dHL;nþ1 ; HL;nþ1 iþ1 iþ1 Dknþ1 ¼ Dkinþ1 þ dDknþ1 ; iþ1 i iþ1 ¼ Ks;nþ1 þ dKs;nþ1 : Ks;nþ1
The iteration is considered to have converged if riþ1 nþ1 ¼ 0 within some prescribed tolerance. 6.1. Consistent tangent operator An important issue in computational plasticity is the derivation of the tangent elastoplastic matrix. In what follows we shall derive the so-called consistent elastoplastic matrix [31,32] that results in a quadratic rate of convergence when the Newton–Raphson method is used at the global equilibrium level, provided the current solution lies within the convergence region. To compute this matrix, the consistent tangent modulus at each Gauss point is needed. Differentiation of Eq. (89) with respect to the strains, taking into account all the variables appearing in the integration algorithm, gives: drnþ1 ornþ1 ornþ1 ovnþ1 ¼ þ : denþ1 oenþ1 ovnþ1 oenþ1
ð97Þ
Furthermore, during the iterative procedure we have drnþ1 ornþ1 ornþ1 ovnþ1 ¼ þ ¼ 0; denþ1 oenþ1 ovnþ1 oenþ1
ð98Þ
which can be rewritten as 1 ovnþ1 ornþ1 ornþ1 ¼ : oenþ1 ovnþ1 oenþ1
ð99Þ
2930
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
By substituting Eq. (99) in Eq. (97) we obtain the elastoplastic consistent tangent stiffness tensor: drnþ1 ornþ1 ornþ1 1 ornþ1 ¼ ðJr;v Þ ¼ Dcon t;nþ1 ; denþ1 oenþ1 ovnþ1 oenþ1
ð100Þ
where Jr;v stands for the Jacobian of the residual rnþ1 with respect to the primary variables, as specified above.
7. Diffuse mechanisms of failure This section includes a series of numerical examples to illustrate the problem of diffuse failure mechanisms under different circumstances. 7.1. Liquefaction of a sand layer under earthquake loading Fig. 6 shows a fully saturated sand layer subjected to a horizontal earthquake. The seismic input is the accelerogram given in Fig. 7 that correspond to the NS component accelerogram of the El Salvador earthquake 13/01/2001 in Santa Tecla. The sand layer is modelled by a sand column with both sides and the bottom assumed impermeable. Pore pressures are assumed to be zero at the surface of the sand layer.
Fig. 6. Soil layer problem: finite element mesh.
5
Comp. NS
4 3 acc. [m/s-2]
2 1 0 -1 -2 -3 -4 -5 0
10
20
30
40 50 time [s]
60
70
80
Fig. 7. Seismic input: NS component accelerogram of the El Salvador earthquake 13/01/2001 in Santa Tecla.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2931
Table 3 Data used in sand layer analysis qs (kg/m3 )
qw (kg/m3 )
Ks (Pa)
Kw (Pa)
n
k (m/s)
2720
980
1e)20
1.092e)9
0.363
2.1e)3
Table 4 Pastor–Zienkiewicz soil model parameters K0 (kPa)
G0 (kPa)
Mg
Mf
ag ¼ af
H0
b0
b1
HU (kPa)
cDM ¼ cU
45,000
22,500
1.5
0.4
0.45
350
4.2
0.2
6e)3
2
Repeated boundary conditions are assumed for the lateral nodes, which implies that the displacement of a right hand side node equals that of the corresponding left hand side node. The finite element mesh is shown in Fig. 6 and consists of stabilized four node quadrilateral elements, where bilinear shape functions are used both for displacements and pressures. The material is assumed to be a very loose sand. Constitutive behaviour has been modelled using the Pastor–Zienkiewicz model described in the previous section. Material properties and other relevant data used in the analysis are listed in Tables 3 and 4. The results can be seen in Fig. 8, which depicts the evolution of the pore pressure and the mean effective stress at stations A–E. It is interesting to note that the pore pressure increases and reaches a plateau between 7 and 17 s depending on the depth, and the mean effective stress decreases to zero. Fig. 9 provides further insight into the phenomenon of liquefaction. We have plotted profiles of pore pressure at several instants together with the vertical stress. In this way, when the excess pore pressure coincides with the vertical effective stress, it is possible to see how liquefaction extends from the surface. At 17 s all the column is liquefied. 7.2. Liquefaction failure of an embankment under earthquake action The case we will be considered next is that of an earthquake induced flowslide in very loose saturated sand. Here it is important to notice that proposed model is intended to describe the initiation of failure. In x 10
5
5
x 10 2 E
E
D
1.5 D
2 C
p’ [Pa]
excess pore-pressure [Pa]
3
C
1
B
B
1
0.5 A
A 0
0
10 time [s]
20
0
0
10
20
time [s]
Fig. 8. Evolution of excess pore pressure and mean effective stress at several stations.
2932
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 30 t=4s t=8s t=12s t=16s t=20s
25
height [m]
20
15
10
5
0
0
1 2 3 excess pore-pressure [Pa]
4 x 10
5
Fig. 9. Vertical profiles of excess pore pressure at different times showing extent of liquefaction.
order to cope with the propagation of the flowslide, the authors have recently proposed a fluid-like, depth integrated model formulated in an Eulerian framework. Flowing material is described by a suitable rheological law [5,33]. The problem consists of an embankment 10 m in height with slopes 2:1, founded on a sand layer which extends 10 m in depth and lies on a rigid rockbed. The material of both the embankment and the foundation is a very loose saturated sand. Initial conditions correspond to geostatic equilibrium under gravity forces. Suction at the surface has been assumed to be equal to )20 kPa. The finite element mesh can be seen in Fig. 10 and consists on 500 quadrilaterals with eight nodes for displacements and 4 for pore pressures. A reduced integration rule has been used in the solid part to avoid locking. The number of nodes is 1611, with 3535 degrees of freedom (Fig. 10). Loading is applied by prescribing horizontal accelerations at the base. We have used the input motion defined in Fig. 7. A simplified absorbing boundary condition has been applied at lateral boundaries. Concerning pore pressures, it has been assumed that no flux occurs at artificial boundaries, and the constant value of )20 kPa has been kept at the surface. The behaviour of the loose material is represented using the Pastor–Zienkiewicz model for sand, Tables 3 and 4. To better understand this behaviour we simulate a loading-unloading-reloading undrained triaxial test. The material loses progressively its resistant capacity, the effective stress decreases and the pore pressure increases, Fig. 11. The results can be seen in Figs. 12–14, where the contours of pore pressure, deviatoric plastic strain 0 0 invariant epl D , p =p0 and displacements are given at different times.
Fig. 10. Finite element mesh of the embankment.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 20
B
20
10 5
F
0 -5
A
E
0
20
40 60 p' (en kPa.)
80
C
10 F
5 0 -5
D
-10
B
15
C q (en kPa.)
q (en kPa.)
15
A
E
-10 -6.E-06
100
-3.E-06
100 E pw (en kPa.)
2933
D 0.E+00 3.E-06 E22 (en %)
6.E-06
C
80 60 40 20 0 -6.E-06
A -3.E-06
0.E+00
3.E-06
6.E-06
E22 (en %)
Fig. 11. Non-drained triaxial simulation of the loose sand: loading, unloading, reloading path.
Fig. 12. Evolution of excess pore pressure contour [Pa].
The deviatoric plastic strain invariant epl D is defined as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i pl pl 2 pl pl 2 pl pl 2 2 ¼ ðe e Þ þ ðe e Þ þ ðe e Þ epl : D 1 2 2 3 3 1 3
2934
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Fig. 13. Evolution of p0 =p00 contour.
Fig. 14. Deformation contour [m].
The ratio between the mean effective confining pressure and its initial value p0 =p00 has been used as an indicator of the extent of the liquefied zones. From these results, it can be concluded that failure of the embankment is caused by liquefaction of the outer liquefied zones. 7.3. Failure of an embankment on a collapsible soil under earthquake loading The following example illustrates the collapse of a bonded material under earthquake loading. The problem is similar to that presented in the previous section but drained conditions are now assumed and the bonded material is now represented using the enhanced Generalized Plasticity model with the new parameters listed in Table 5. The evolution of plastic strain during the earthquake is plotted in Fig. 15. The plastic strain is concentrated in the outer zones of the slopes. Fig. 16 shows the evolution of the pt parameter. This parameter decreases as the volumetric plastic strain increases and becomes zero. This indicates the evolution of the debonding and the progressive reduction of the apparent cohesion of the material.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2935
Table 5 New model parameters for the enhanced Generalized Plasticity model b1
b2
pt0 (kPa)
qt
pc0 (kPa)
5e)6
36
100
1000
2000
Fig. 15. Evolution of deviatoric plastic strain invariant epl D contour.
Another interesting result shows evidence of the mechanical collapse in Figs. 17 and 18, where we have plotted the evolution of the void ratio throughout the embankment and at a critical point A. Finally, the deformation contour at t ¼ 0:88 s is plotted in Fig. 19. 7.4. A note on air liquefaction failure If the soil is partially saturated, its behaviour depends on the coupling between the solid skeleton and the pore air and water. In the limiting case of a dry soil, the air has to flow out from the pores for the material to consolidate, but typical air consolidation times are much shorter than those of the water. Therefore, in practical cases, the role of air pressures is neglected, as the characteristic time of loading is much longer than consolidation time. However, it is possible to imagine situations with much shorter loading times, where coupling between pore air and soil skeleton plays a paramount role. This is the case of fluidized granular beds, just to mention a particular example of industrial interest. Bishop [34] describes the case of Jupille flowslide, which happened in Belgium in February 1961. A tip of uncompacted fly ash located in the upper part of a narrow valley collapsed, and the subsequent flowslide travelled for about 600 m at very high speeds (130 km/h) until it stopped. A triggering mechanism was suggested to be ‘‘collapse due to undermining of a steep, partly saturated and slightly cohesive face’’. Bishop referred to Calembert and Dantinne [35], who pointed out the role of the entrapped air. The mechanism of pore air consolidation can explain the ‘‘sort of fog’’ which formed above the flowing material, as warmer air met the colder winter air in the exterior. The ‘‘honeycomb’’ like structure of soil in the tip was responsible for its sudden collapse when movement started.
2936
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Fig. 16. Evolution of pt parameter [kPa].
Fig. 17. Evolution of void ratio contour.
Dry liquefaction can be exhibited also by soils of volcanic origin. Some of the catastrophic landslides during the El Salvador earthquake of 13th of January 2001 can be explained by this mechanism, the collapse of the material under the dynamic loading could make the pore air pressure increase.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2937
collapse
0.57
void ratio
0.56
0.55
0.54
0.53 0
0.2
0.4 0.6 time [s]
0.8
1
Fig. 18. Evolution of void ratio at point A.
Fig. 19. Deformation contour [m].
8. Conclusions We have proposed in this paper a theoretical framework to model the initiation mechanisms of catastrophic landslides. First, the equations describing the coupling have been presented following an Eulerian approach based on mixture theory. The main advantage is that they provide a unified formulation that can also be used for the study of a subsequent propagation phase [5,6]. The system of partial differential equations is then discretized using the classical Galerkin Finite Element Method. An enhanced model in the framework of the Generalized Plasticity has been presented to describe the mechanical behaviour of bonded soils. This model is able to reproduce the mechanical collapse phenomenon and to characterize the basic features of a sand behaviour under monotonic and cyclic loading when debonding occurred. An implicit integration scheme for generalized plasticity has been introduced to guarantee the efficiency of the numerical convergence. Finally, applications to diffuse failure have been presented. Acknowledgements The authors gratefully acknowledge the financial support provided by the Spanish Agency for International Cooperation (AECI), by the Spanish Ministry of Science and Technology (Project Andes and Ram on y Cajal contract) and by the European Union (Project Diga).
2938
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
References [1] R. Dikau, D. Brundsen, L. Schrott, M.L. Ibsen, Landslide Recognition, John Wiley and Sons, 1996. [2] I. Vardoulakis, J. Sulem, Bifurcation Analysis in Geomechanics, Blakie Academic and Professional, 1995. [3] M. Pastor, C. Tamagnini (Eds.), Numerical Modelling in Geomechanics, Revue Francßaise du genie civil, Hermes-Lavoisier, Paris, vol. 6, 2002. [4] P. Mira, An alisis por elementos finitos de problemas de rotura de geomateriales, Ph.D. thesis, Universidad Politecnica de Madrid, Madrid, 2002. [5] M. Pastor, M. Quecedo, J.A. Fernandez Merodo, M.I. Herreros, E. Gonzalez, P. Mira, Modelling tailing dams and mine waste dumps failures, Geotechnique 52 (8) (2002) 579–591. [6] J.A. Fern andez Merodo, Une Approche a la modelisation des glissements et des effonfrements de terrains: Initiation et Propagation, Ph.D. thesis, Ecole Centrale de Paris, December 2001. [7] F. Darve, F. Laouafa, Constitutive equations and instabilities of granular materials, in: G. Capriz, V.N. Ghionna, P. Giovine (Eds.), Modeling and Mechanics of Granular and Porous Materials, Birkh€auser, 2002, pp. 3–44. [8] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155–164. [9] M.A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26 (1955) 182–185. [10] O.C. Zienkiewicz, C.T. Chang, P. Bettess, Drained, undrained, consolidating dynamic behaviour assumptions in soils, Geotechnique 30 (1980) 385–395. [11] O.C. Zienkiewicz, T. Shiomi, Dynamic behaviour of saturated porous media: the generalised Biot formulation and its numerical solution, Int. J. Num. Anal. Methods Geomech. 8 (1984) 71–96. [12] O.C. Zienkiewicz, A.H.C. Chan, M. Pastor, D.K. Paul, T. Shiomi, Static and dynamic behaviour of soils: a rational approach to quantitative solutions. I. Fully saturated problems, Proc. R. Soc. Lond. A 429 (1990) 285–309. [13] O.C. Zienkiewicz, Y.M. Xie, B.A. Schrefler, A. Ledesma, N. Bicanic, Static and dynamic behaviour of soils: a rational approach to quantitative solutions. II. Semi-saturated problems, Proc. R. Soc. Lond. A 429 (1990) 311–321. [14] O.C. Zienkiewicz, A.H.C. Chan, M. Pastor, B. Schrefler, T. Shiomi, Computational Geomechanics, John Wiley and Sons, 2000. [15] R.L. Lewis, B.A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, John Wiley and Sons, 1998. [16] O. Coussy, Mechanics of Porous Media, John Wiley and Sons, Chichester, 1995. [17] R. de Boer, Theory of Porous Media, Springer-Verlag, Berlin, 2000. [18] A.H.C. Chan, A unified finite element solution to static and dynamic geomechanics problems, Ph.D. dissertation, University college of Swansea, Wales, 1988. [19] M.G. Katona, O.C. Zienkiewicz, A unified set of single-step algorithms, Part 3: the beta-m method, a generalisation of the Newmark scheme, Int. J. Num. Methods Engrg. 21 (1985) 1345–1359. [20] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method (4th ed.), vol. 2, McGraw-Hill, 1991. [21] S.W. Sloan, M.F. Randolph, Numerical prediction of collapse loads using finite element methods, Int. J. Numer. Anal. Methods Geomech. 6 (1982) 47–76. [22] M. Pastor, M. Quecedo, O.C. Zienkiewicz, A patch test for mesh alignment effects in localized failure, Commun. Numer. Methods Engrg. 11 (1995) 1015–1024. [23] M. Pastor, T. Li, X. Liu, O.C. Zienkiewicz, Stabilized low order finite elements for failure and localization problems in undrained soils and foundations, Comput. Methods Appl. Mech. Engrg. 174 (1999) 219–234. [24] P. Mira, M. Pastor, T. Li, X. Liu, A new stabilized enhanced strain element with equal order of interpolation for soil consolidation problems, Comput. Methods Appl. Mech. Engrg. 192 (2003) 37–38, and 4257–4277. [25] M. Pastor, O.C. Zienkiewicz, A.H.C. Chan, Generalized plasticity and the modelling of soil behaviour, Int. J. Numer. Anal. Methods Geomech. 14 (1990) 151–190. [26] A. Gens, R. Nova, Conceptual bases for a constitutive model for bonded soils and weak rocks, in: Proceedings of International Symposium on Geotechnical Engineering of Hard Soils–Soft Rocks, Balkema, Rotterdam, 1993, pp. 485–494. [27] R. Lagioia, R. Nova, An experimental and theoretical study of the behaviour of a calcarenite in triaxial compression, Geotechnique 45 (4) (1995) 633–648. [28] L. Tonni, Modellazione numerica del comportamento di terreni granulari con la Plasticita Generalizzata. Ph.D. thesis, Politecnico di Torino, 2002. [29] O.M. Heeres, R. de Borst, Implicit integration of standard and non-standard elasto-plastic models, ECCOMAS 2000, September 2000, Barcelona (Spain), pp. 11–14. [30] M.A. Crisfield, Non-Linear Finite Element Analysis of Solids and Structures, vol. 1–2: Essentials, John Wiley and Sons, 1991. [31] K. Runesson, A. Samuelsson, Aspects on numerical techniques in small deformation plasticity, in: J. Middleton et al. (Ed.), NUMETA 85, Numerical Methods in Engineering: Theory and Applications, vol. 1, A.A. Balkema, Rotterdam, 1985, pp. 337– 348.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2939
[32] J.C. Simo, R.L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Methods Appl. Mech. and Engrg. 48 (1985) 101–118. [33] M. Quecedo, M. Pastor, M.I. Herreros, J.A. Fernandez Merodo, Numerical modelling of the propagation of fast landslides using the finite element method, Int. J. Num. Meth. Eng. 59 (2004) 755–794. [34] A.W. Bishop, The stability of tips and spoil heaps, Quart. J. Engrg. Geol. 6 (1973) 335–376. [35] L. Calembert, R. Dantinne, The avalanche of ash at Jupille (Liege) on February 3rd, 1961. From: The commemorative volume dedicated to Professeur F. Campus, Liege, Belgium, 1964, pp. 41–57.
Modelling of diffuse failure mechanisms of catastrophic landslides J.A. Fernandez Merodo a,b,*, M. Pastor a,b, P. Mira a,b, L. Tonni c, M.I. Herreros a,b, E. Gonzalez a,b, R. Tamagnini d a Centro de Estudios y Experimentacion de Obras Publicas, Alfonso XII 3, 28014 Madrid, Spain M2 i (Math. Model. Eng. Group), Department of Applied Mathematics, ETS de Ingenieros de Caminos, UPM, Ciudad Universitaria, 28040 Madrid, Spain c Dept. DISTART, Universita di Bologna, Viale Risorgimento 2, 40136 Bologna, Italy Dept. Structural and Geotechnical Engineering, Universita di Roma ‘La Sapienza’, Via Montedoro 28, 00100 Roma, Italy b
d
Received 15 April 2003; received in revised form 15 July 2003; accepted 1 September 2003
Abstract This paper presents a theoretical and numerical framework to model the initiation mechanisms of catastrophic landslides. Attention is paid to the mechanism of failure referred to as diffuse, characteristic of soils presenting very loose or metastable structures, as in the case of the liquefaction of a sand layer induced by an earthquake, where the effective stress approaches zero as the pore pressure increases. The equations describing the coupling between the solid skeleton and the pore fluids are presented following an Eulerian approach based on mixture theory that can provide a unified formulation for both initiation and propagation phases. The system of partial differential equations is then discretized using the classical Galerkin finite element method and neglecting the convective terms. A simple improvement of the Generalized Plasticity model is introduced to deal with bonded geomaterials and to represent the mechanical collapse of such materials. The use of an implicit integration scheme guarantees the efficiency of the numerical convergence. 2004 Elsevier B.V. All rights reserved. Keywords: Landslides; Diffuse failure; Liquefaction; Bonded geomaterials; Finite element analysis; Implicit integration
1. Introduction Landslides are one of natural catastrophes causing important losses of human lives and extreme damage to property. There is a wide variety of types of landslides, depending on the materials involved and the triggering mechanism. The time scale involved can range from minutes to years. For a classification, it is worth consulting Ref. [1] where a detailed description is provided.
*
Corresponding author. Address: Centro de Estudios y Experimentaci on de Obras Publicas, Alfonso XII 3, 28014 Madrid, Spain. E-mail address: [email protected] (J.A. Fernandez Merodo).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.09.016
2912
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Landslides and failure of slopes are caused by changes in the effective stresses, variation of material properties or changes in the geometry. Changes of effective stresses can be induced either directly, as a consequence of variation of the external forces (earthquakes, human action), or indirectly through pore pressures (rainfall effects). Variations in material properties can be caused by processes of degradation (weathering and chemical attack). Finally, geometry can change because of natural causes (erosion) or human action (excavation, construction, reshaping. . .). In some cases, the failure mechanism consists of a clearly defined surface where shear strain concentrates. From a mathematical point of view, the inception of a shear band is characterized by a discontinuity in the strain field, which can evolve towards a discontinuity in the displacements at a later stage. Much effort has been devoted during the past years to better understand this phenomenon. The problem is ill posed for elastoplastic materials, and the results obtained in numerical models depend on the mesh size and alignment. The interested reader is addressed to Ref. [2] or [3] and [4] where a detailed explanation is provided. If the soil has a low relative density, and tends to compact when sheared, failure mechanism can be different from the localized mode. It affects not only the material on a narrow band or failure surface, but a much larger mass of soil. Failure can be described as of a flow type. Pore pressures developed may cause liquefaction of the soil, and liquefaction––or conditions close to it––may have played an important role on failure of slopes such as those of Aberfan in South Wales, Anhui in China, Jupille in Belgium and Rocky Mountains coal mine dumps. This type of mechanism has been described as diffuse [5–7]. The key to properly model this type of landslide is the ability of the constitutive model to reproduce liquefaction. After presenting the mathematical and the numerical model in a general mixture framework, we present an enhanced Generalized Plasticity model that can reproduce the mechanical behaviour of bonded geomaterials. This model is not only able to characterize the basic features of sand behaviour under monotonic and cyclic loading but also to represent the mechanical collapse of bonded materials. An implicit integration of such kind of Generalized Plasticity model is presented to assure the effectiveness of the numerical convergence. Finally, diffuse mechanisms of failure are illustrated with some examples.
2. Mathematical model for soil skeleton–pore fluid coupling 2.1. Introduction Soils and rocks are geomaterials with voids which can be filled with water, air, and other fluids. They are, therefore, multiphase materials, exhibiting a mechanical behaviour governed by the coupling between all the phases. Pore pressures of fluids filling the voids play a paramount role in the behaviour of a soil structure, and indeed, their variations can induce failure. The first mathematical model describing the coupling between solid and fluid phases was proposed by Biot [8,9] for linear elastic materials. This work was followed by further development at Swansea University, where Zienkiewicz and co-workers [10–14] extended the theory to non-linear materials and large deformation problems. It is also worth mentioning the work of Lewis and Schrefler [15], Coussy [16] and de Boer [17]. Among the different alternative ways which can be used to describe the coupling between solid skeleton and pore fluids, we have chosen an approach closer to mixture theories than to the more classical approach used in computational geotechnics, and it provides a more general description which can be used not only for initiation of failure but also for propagation of catastrophic landslides. The soil skeleton consists of particles of density qs having a porosity n (volume percent of voids in the mixture). The voids can be filled with air and water, but in some cases there will be mixtures of water and
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2913
very fine particles which can be considered as a fluid phase. We will assume that fluid phases ðaÞ are not miscible, defining Sa as the degree of saturation for the phase a, i.e., the volume fraction of voids occupied by it. In the case of air and water, we will define Sw and Sa as Sw ¼
volume of water ; volume of voids
Sa ¼
volume of air volume of voids
ð1Þ
with Sw þ Ss ¼ 1: We will introduce qðaÞ to denote the density of the phase ðaÞ: If the fluid density is qa ; qðaÞ ¼ Sa nqa ;
ð2Þ
while for the solid we can write qðsÞ ¼ ð1 nÞqs :
ð3Þ
Density of the mixture is therefore obtained by adding the densities of all constituents in the mixture. Concerning the kinematics, we will use u for displacements and v for velocities, with superindexes referring to the specific constituents. In this way, uðsÞ will be the displacement field of the solid phase, and vðaÞ the velocity of pore fluid ðaÞ: In most geotechnical applications, the movement of the fluid phases is described by the velocity relative to the soil skeleton. The so-called Darcy’s velocity of phase ðaÞ is defined as wðaÞ ¼ nSa ðvðaÞ vðsÞ Þ:
ð4Þ
2.2. Effective and partial stresses If we denote by rs and ra the Cauchy stress tensors acting on the solid particles and on the phase ðaÞ, we can define the partial stresses rðsÞ and rðaÞ as rðsÞ ¼ ð1 nÞrs ;
ð5Þ
rðaÞ ¼ nSa ra : The total stress tensor acting on the mixture is obtained as r ¼ rðsÞ þ
nX phases
rðaÞ :
ð6Þ
a¼1
The partial stresses rðaÞ can be decomposed into hydrostatic and deviatoric components as rðaÞ ¼ nSa pa I þ nSa sa , where sa ¼ devðra Þ is the deviatoric part of ra , and I is the identity tensor of second order. It is important to note that if the saturating fluid is water, the deviatoric stresses can be neglected, but in other cases where viscous contributions are important, this term has to be taken into account. It is convenient to introduce an averaged pore pressure p as p¼
nX phases a¼1
pðaÞ ¼
nX phases a¼1
Sa pa :
ð7Þ
2914
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Two important particular cases are those of (i) non-saturated soils, and (ii) saturated soils, where pore fluids are (i) water and air, and (ii) water. The above relations specialize to r ¼ rðsÞ þ rðwÞ þ rðaÞ ; rðsÞ ¼ ð1 nÞrs ; rðwÞ ¼ nSw pw I;
ð8Þ
rðaÞ ¼ nSa pa I; p ¼ Sw pw þ Sa pa and r ¼ rðsÞ þ rðwÞ ; rðsÞ ¼ ð1 nÞrs ;
ð9Þ
rðwÞ ¼ npw I; p ¼ pw ; respectively. Another important case is that of non-saturated soils having pa ¼ 0; for which r ¼ rðsÞ þ rðwÞ ; rðsÞ ¼ ð1 nÞrs ;
ð10Þ
rðwÞ ¼ nSw pw I; p ¼ Sw pw :
This approach is valid for many engineering cases of interest. A more complete description can be found in Ref. [15]. The effective stress r0 is defined as p: r0 ¼ r þ I
ð11Þ
In the case of viscous pore fluids where deviatoric stresses are important, the effective stress is given by nX phases r0 ¼ r þ I pn sðaÞ : ð12Þ a¼1
Another point worth mentioning is that the effective stress may be modified to account for compressibility of soil grains relative to that of solid skeleton. Here we will make the assumption that solid skeleton is incompressible. 2.3. Balance of mass We will introduce the material derivative of a magnitude / following phase a as DðaÞ / o/ ¼ þ vðaÞ $/: Dt ot
ð13Þ
One important relation which will be used later is DðaÞ / DðsÞ / ¼ þ ðvðaÞ vðsÞ Þ $/ Dt Dt DðsÞ / wðaÞ þ ¼ $/; Dt nSa where we have introduced the Darcy’s velocity of phase ðaÞ relative to soil skeleton defined in (4).
ð14Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2915
The balance of mass equation for the solid phase is given by DðsÞ qðsÞ þ qðsÞ $ vðsÞ ¼ 0 Dt
ð15Þ
and for the fluid phases ðaÞ DðaÞ qðaÞ þ qðaÞ $ vðaÞ ¼ 0; Dt
ð16Þ
where qðsÞ ¼ ð1 nÞqs and qðaÞ ¼ nSa qa : In the case of non-saturated soils, we can write DðsÞ fð1 nÞqs g þ ð1 nÞqs $ vðsÞ ¼ 0; Dt DðaÞ ðnSa qa Þ þ nSa qa $ vðaÞ ¼ 0; Dt
ð17Þ
with a ¼ fa; wg. From here, we obtain DðsÞ n 1 n DðsÞ qs ¼ þ ð1 nÞ$ vðsÞ ; Dt qs Dt DðaÞ n n DðaÞ ¼ ðSa qa Þ þ n$ vðaÞ ; Dt Sa qa Dt If we now add both equations and use (14) we arrive at ðaÞ w 1 n DðsÞ qs n DðaÞ $ þ ðSa qa Þ þ $ vðsÞ ¼ 0 þ qs Sa qa Dt Sa Dt
ð18Þ
ð19Þ
for a ¼ fa; wg. If we neglect now thermal effects, we can elaborate further the second terms as follows: p 1 n DðsÞ qs 1 n DðsÞ : ¼ Dt qs Ks Dt In the case of a non-saturated soil with pa ¼ 0 and p ¼ Sw pw we obtain ( ) 1 n DðsÞ DðsÞ pw DðsÞ Sw p 1n ¼ Sw þ pw : Ks Ks Dt Dt Dt The third term (water phase) can be obtained for non-saturated soils with pa ¼ 0 as ( ) n DðwÞ n DðwÞ pw n DðwÞ Sw ðSw qw Þ ¼ þ : Sw qw Dt Kw Dt Sw Dt Therefore, the balance of mass for the water (non-saturated, pa ¼ 0) can be written as w n DðwÞ pw ð1 nÞ DðsÞ pw n DðwÞ Sw 1 n DðsÞ Sw þ þ þ þ $ vðsÞ ¼ 0: pw $ þ Sw Kw Dt Ks Sw Dt Ks Dt Dt
ð20Þ
ð21Þ
ð22Þ
ð23Þ
Sometimes, it is introduced for convenience the specific storage coefficient of the soil Cs as Cs ¼ n
oSw : opw
ð24Þ
2916
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
From here we can obtain some particular cases of interest: • Saturated soil $wþ
n DðwÞ pw ð1 nÞ DðsÞ pw þ þ $ vðsÞ ¼ 0: Kw Dt Ks Dt
ð25Þ
This equation can be further simplificated by taking into account the relationship between material derivatives following the fluid and the solid phases, Eq. (14) $wþ
n DðsÞ pw ð1 nÞ DðsÞ pw w þ þ $ vðsÞ þ $pw ¼ 0: Kw Dt Ks Kw Dt
The last term is usually neglected in geotechnical analysis, and above equation reduces to $wþ where 1 ¼ Q
1 DðsÞ pw þ $ vðsÞ ¼ 0; Q Dt
n ð1 nÞ þ : Kw Ks
ð26Þ
ð27Þ
• Flow of dry geomaterials In some circumstances, there may exist an important coupling between the pore air and the soil skeleton. Should this be the case, the balance of mass equation for the pore gas can be written in a similar way as $wþ with 1 ¼ Qa
1 DðsÞ pa þ $ vðsÞ ¼ 0 Qa Dt
ð28Þ
n ð1 nÞ þ ; Ka Ks
ð29Þ
where have introduced the volumetric stiffness of the air Ka . The balance of mass equation of the mixture can be obtained by adding the equations for all constituents: nX phases nX phases DðsÞ qðsÞ DðaÞ qðaÞ þ qðsÞ $ vðsÞ þ þ qðaÞ $ vðaÞ ¼ 0 Dt Dt a¼1 a¼1
from where, using (14) we get nX phases ðaÞ DðsÞ q w ðsÞ þ q$ v þ $qðaÞ ¼ 0: Dt nS a a¼1
ð30Þ
If all wðaÞ are negligible, above equation reduces to DðsÞ q þ q$ vðsÞ ¼ 0: Dt
ð31Þ
A particular case of interest is that of very fluid slurries, for which the mixture density can be assumed to be constant. The balance of mass reduces to $ vðsÞ ¼ 0
ð32Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2917
2.4. Balance of momentum of pore fluid The balance of momentum for a fluid phase ðaÞ reads qðaÞ
DðaÞ vðaÞ ¼ qðaÞ b þ div rðaÞ þ Ra ; Dt
ð33Þ
where b is the body forces vector, and Ra is a force representing the interaction between the solid skeleton and the pore fluid ðaÞ; which can be assumed to follow Darcy’s law Ra ¼ ka1 wðaÞ :
ð34Þ
In above, ka is the permeability matrix for phase ðaÞ, ka ¼
karel Kaintr la
ð35Þ
where karel is a relative permeability which depends on Sa , la is the dynamic viscosity of fluid ðaÞ, and Kaintr is the intrinsic permeability matrix which has dimensions of L2 . If the material is isotropic, then the matrix reduces to a constant. The permeability ka has dimensions of L3 T =M. It should be noted that in geotechnical analysis the permeability is often defined in a slightly different way, a ¼ qa bka k with dimensions of LT 1 : If we take into account the relation (14), we can write the balance of momentum as ! DðsÞ vðsÞ ðaÞ q þ c:t: ¼ qðaÞ b þ div rðaÞ ka1 wðaÞ ; Dt where the correcting terms c.t. are given by ðaÞ DðsÞ wðaÞ wðaÞ ðsÞ wðaÞ w $v þ $ þ : c:t: ¼ Dt nSa nSa nSa nSa In many geotechnical problems the correcting terms can be neglected , and wðaÞ is obtained as ! ðsÞ ðsÞ ðaÞ ðaÞ D v ðaÞ ðaÞ w ¼ ka q þ q b þ div r : Dt
ð36Þ
ð37Þ
ð38Þ
ð39Þ
This situation arises when (i) acceleration of fluid phases are small, and (ii) when cross products of wðaÞ and vðsÞ can be neglected. The interested reader is referred to the work of Zienkiewicz et al. [10] where the reasons are detailed. If the soil is saturated by water, and no viscous effects are present in the pore water, we can write ! DðsÞ vðsÞ w ¼ k w qw þ qw b $pw ð40Þ Dt and in the case of dry soils we will have ! DðsÞ vðsÞ þ qa b $pa : w ¼ k a qa Dt
ð41Þ
2918
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2.5. Balance of momentum of the mixture We obtain the balance of linear momentum for the mixture by combining the balance equations of the solid and fluid phases. Concerning the solid, we can write qðsÞ
DðsÞ vðsÞ ¼ qðsÞ b þ div rðsÞ þ Rs ; Dt
ð42Þ
where Rs stands for the drag forces due to all fluid constituents Rs ¼
nX phases
Ra :
a¼1
Eq. (42) is now added to (37), obtaining the balance of momentum for the mixture as ( ðaÞ ) nX phases DðsÞ vðsÞ DðsÞ wðaÞ wðaÞ ðsÞ wðaÞ w þ nSa qa $v þ $ q þ ¼ qb þ div r: Dt Dt nSa nSa nSa nSa a¼1
ð43Þ
This equation can be simplified as q
DðsÞ vðsÞ ¼ qb þ div r; Dt
ð44Þ
which is the form usually encountered in geotechnics.
3. Some useful forms in computational soil mechanics: the u–pw formulation In the preceding section, we have presented the balance of mass and momentum equations for a mixture of solid particles and nphases immiscible fluid phases occupying the pores, with some possible simplifications. The equations have been cast in a Eulerian formulation, from which simple Lagrangian forms involving small deformations can be immediately derived by neglecting the convective components of the material derivatives. The approach we propose in this work consists on modelling catastrophic landslides using a small deformation approach for the initiation phase. Therefore, we will not consider convective terms, and we will further assume that the accelerations of the pore fluids relative to the soil skeleton can be neglected. The validity of this assumption has been studied by Bettess and Zienkiewicz [10,14]. The two cases we will study next are those of (i) non-saturated soils with air at atmospheric pressure, and (ii) saturated soils. We will use displacements of the soil skeleton u instead of velocities. In the former case, the system of equations describing the coupling are the following: • Balance of mass for pore water $ w þ Sw
oev opw 1 opw þ Cs þ ¼ 0; Q ot ot ot
ð45Þ
where ev is the volumetric strain $ vðsÞ ¼ and
oev ot
ð46Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
1 ¼ Q
n
Sw 1 n Cs Sw Sw þ pw þ : Ks Kw n
• Balance of linear momentum (pore water) 2 ou o w $pw þ qw b kw1 w qw þ ¼ 0: ot2 ot nSw • Balance of linear momentum (mixture) o2 u o w div r þ qb q 2 nqw Sw ¼ 0: ot ot nSw
2919
ð47Þ
ð48Þ
ð49Þ
• Constitutive equation dr0 ¼ D : de: • Kinematics 1 oui ouj deij ¼ þ S du: 2 oxj oxi
ð50Þ
ð51Þ
Under certain conditions [10], we can neglect the relative accelerations of the fluid phase, and eliminate the relative displacements of the fluid phase w. This can be done by substituting the value of $ w obtained from (48) o2 u $ w ¼ $ kw $pw þ qw b qw 2 ð52Þ ot into the balance of mass of the fluid phase (45) from where we obtain o2 u 1 opw oev þ Sw ¼ 0: þ Cs þ $ kw $pw þ qw b qw 2 ot Q ot ot
ð53Þ
The main advantage of this formulation is than the field variables reduce to two, displacements u and pressures pw . Therefore, the basic equations of this u–pw formulation are (49) and (53), together with the constitutive equation (50) and the kinematic relations between strains and displacements. The simplification for saturated soils is immediate. Eq. (53) reduces to o2 u 1 opw oev þ ¼ 0: ð54Þ $ kw $pw þ qw b qw 2 þ ot Q ot ot The equations presented so far have to be complemented by suitable boundary and initial conditions.
4. Numerical model: discretization of the u–pw model The system of partial differential equations can be discretized using standard Galerkin techniques, as described in [12]. After approximating the fields u and pw as: u ¼ Nu u, pw ¼ Np pw , it results in two ordinary differential equations:
2920
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
d2 u M 2þ dt
Z
BT r0 dX Q pw f u ¼ 0;
ð55Þ
X
where B ¼ S Nu , and QT
d u d p þ H pw þ C w f p ¼ 0: dt dt
The matrices given above are defined as Z M¼ qNTu Nu dX; ZX 1 T C¼ Np Np dX; X Q Z Q¼ Sw aBT mNp dX; X Z H ¼ $NTp kw $Np dX
ð56Þ
ð57Þ
and Z
NTu b dX þ
Z
NTu t dC; Z Z Z Z op dC þ $NTp kw qw b dX fp ¼ NTp kw $Np kw Nu dX €u NTp s0 dX: on Cq X X X
fu ¼
X
Ct
ð58Þ
The term including accelerations in f p is usually disregarded (Ref. Thesis Chan [18]). The time derivatives of u and pw are approximated in a typical step of computation using the Generalized Newmark GN22 scheme for displacements and a GN11 for the water pressure [19,20]. If we introduce the following notation [19]: unþ1 € un ; D€ un ¼ € Dp_ wn ¼ p_ wnþ1 p_ wn ; u_ nþ1 ¼ u_ p;nþ1 þ b1 DtD€ un ; unþ1 ¼ up;n þ 12b2 Dt2 D€ un ; pwnþ1
¼
pwp;nþ1
þ
ð59Þ
hDtDp_ wn ;
where un ; u_ p;nþ1 ¼ u_ n þ Dt€ up;nþ1 ¼ un þ Dtu_ n þ 12Dt2 € un ; pwp;nþ1
¼
pwn
þ
we obtain the discretized system of equations valid in each time step: Z 0 p_ nw Funþ1 Uu ¼ 0; MD€ un þ BT r nþ1 dX hDtQD X
ð60Þ
Dtp_ wn ;
ð61Þ
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
p_ nw Fpnþ1 Up ¼ 0; b1 DtQT D€ un þ ðDthH þ CÞD
2921
ð62Þ
€ n and D p_ nw . where the unknown values are Du If this system is non-linear, it can be solved by using a Newton Raphson method with a suitable Jacobian matrix: 2
oUu 6 oD€ u 6 6 4 oUp oD€ u
3 oUu ðiÞ " #ðiþ1Þ p_ 7 oD Uu ðiÞ uÞ 7 dðD€ ¼ : 7 Up p_ Þ oUp 5 dðD p_ oD
ð63Þ
Using Eqs. (61) and (62) we can write the above step as "
M þ 12 Dt2 b2 KT
hDtQ
b1 DtQT
DthH þ S
#ðiÞ "
dðD€ uÞ p_ Þ dðD
#ðiþ1Þ
where KT is the tangent stiffness matrix. KT ¼
¼
R
Uu
ðiÞ
Up
;
ð64Þ
BT Dep B dX.
4.1. A note on element type The proper choice of the element type is of paramount importance as some elements introduce errors resulting on (i) unrealistic limit loads and (ii) spurious failure elements. These phenomena have been the subject of an extensive research during the past years [21–23]. Based on the authors experience, the use of enhanced strain elements together with a suitable stabilization technique provides optimal results [24]. Alternatively, reduced integration with control of spurious modes may be used.
5. An enhanced generalized plasticity model for bonded geomaterials 5.1. Basic theory The elastoplastic behaviour of the material is described by the general incremental relationship, relating the strain increment de to the stress increment dr: dr ¼ D de;
ð65Þ
where the matrix D, elastoplastic stiffness matrix, depends not only on the present state conditions but on the stress–strain history and the loading direction. Taking into account this last dependency two different matrices are introduced DL and DU ; that correspond to the loading and unloading situations. A unit vector n is introduced to define the loading/unloading condition: drTe n > 0 loading; drTe n < 0 unloading; drTe n ¼ 0 neutral loading; where dre ¼ De de and De is an elastic matrix dependent only on the present state of the material.
ð66Þ
2922
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
1 The continuity between loading and unloading implies the general form of D1 L and DU :
ngL nT ; HL ngU nT ¼ D1 ; e þ HU
1 D1 L ¼ De þ
D1 U
ð67Þ
where ngL and ngU are arbitrary directions vectors and HL and HU are scalars defined as the loading and unloading plastic moduli respectively. The strain increment de can be decomposed into elastic and plastic parts as de ¼ dee þ dep
ð68Þ
with dee ¼ De dre ; ngL nT dr HL ngU nT dr dep ¼ HU
dep ¼
for loading;
ð69Þ
for unloading:
The values of the DL and DU matrices can be obtained by suitable manipulation, giving: DL ¼ De þ
De ngL nT De ; HL þ nTgL De n
DU ¼ De þ
De ngU nT De : HU þ nTgU De n
ð70Þ
A complete elastoplastic behaviour can thus be specified by prescribing at each point of the stress space (and for a given history) three directions (the loading direction n and the plastic flow directions ngL and ngU ), two scalars (the plastic moduli HL and HU ) and the elastic stiffness matrix De . It has to be remarked here that in the framework of the Generalized Plasticity, neither a yield surface nor a plastic potential have been introduced. It is also noted that plastic strains during unloading are permissible. Furthermore, it must be noted that the classical plasticity, the bounding surface models and the kinematic hardening models are particular cases of such a generalized plasticity. 5.2. The Pastor–Zienkiewicz model for sands In what follows, a sand model developed in the framework of the Generalized Plasticity by Pastor et al. [25] will be described. Using the conventional notation of soil mechanics, three stress invariants are defined: 1 p ¼ I1 ; p3ffiffiffiffiffiffiffi q ¼ 3J2 ; 1 h ¼ arcsin 3
! pffiffiffi 3 3 J3 ; 2 J23=2
ð71Þ
where I1 is the first invariant of the stress tensor and J2 and J3 are the second and third invariants of the deviatoric stress tensor (s ¼ r 13 I1 I).
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2923
Soil dilatancy dg has been expressed as linear function of the stress ratio g ¼ q=p0 : dg ¼ ð1 þ ag ÞðMg gÞ;
ð72Þ
where ag is a material constant and Mg is the slope of the critical state line which is related to the residual angle /0 by Mg ¼
6 sin /0 : 3 sin /0 sin 3h
ð73Þ
As a result, for loading stress increments the plastic flow direction vector is given by T
ngL ¼ ðnpgL ; nqgL ; nhgL Þ ;
ð74Þ
whose components are expressed as dg npgL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ dg2 1 nqgL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ dg2
ð75Þ
qMg cos 3h nhgL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 1 þ dg2 In the case of unloading, it is noted that irreversible strains are of contractive or densifying nature. Therefore the unloading plastic flow direction vector ngU is defined as ngU ¼ ðnpgU ; nqgU ; nhgU ÞT with
ð76Þ
npgU ¼ npgL ; ð77Þ
nqgU ¼ nqgL ; nhgU ¼ nhgL : T
Assuming a non-associated flow rule, the loading direction n ¼ ðnp ; nq ; nh Þ must be different from ngL , but with similar expressions for components: df np ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2 ; 1 þ df 1 nq ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2 ; 1 þ df qM g cos 3h nh ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi2 ; 2 1 þ df
ð78Þ
df ¼ ð1 þ af ÞðMf gÞ
ð79Þ
where
and af and Mf are constitutive parameters.
2924
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
In order to take into account the main features of sand response––i.e. the existence of a critical state condition, dilative response after peak, liquefaction in loose sands, memory of previous stress path––Pastor et al. [25] proposed for the plastic modulus HL the following relationship: HL ¼ H0 p0 Hf ðHv þ Hs Þ HDM
ð80Þ
with 4 g 1 Hf ¼ 1 and gf ¼ 1 þ Mf ; gf af g Hv ¼ 1 ; Mg Z Hs ¼ b0 expðb1 nÞ and n ¼ jdnps j accumulated deviatoric plastic strain; HDM ¼
fmax f
cDM
0
and f ¼ p 1
af 1 þ af
g Mf
1=a mobilized stress function;
where H0 ; cDM ; b0 ; b1 are constitutive parameters. For the plastic modulus in unloading condition HU , according to the experimental results during undrained loading-unloading triaxial tests, Pastor et al. [25] proposed the following simple expression: c Mg Mg u > 1; HU ¼ Hu0 for gu gu ð81Þ Mg 6 1; HU ¼ Hu0 for gu where Hu0 , cu are material parameters. gu is the stress ratio from which unloading takes place. The material has a non-linear elastic response, according to the following relationships: dp ¼ Kt deev ;
dq ¼ 3Gt dees ;
ð82Þ
where the tangent bulk and shear modulus Kt and Gt are assumed to be dependent only on the hydrostatic part of the stress tensor: K ¼ K0
p0 ; p00
G ¼ G0
p0 p00
ð83Þ
with K0 , G0 and p00 references values. For the present model there are 12 material parameters requiring definition. Generally, all parameters are identified from monotonic and cyclic triaxial tests, though in certain cases some parameters are adopted from previous experiences if full test records are unavailable. This model is able to characterize the basic features of the sand behaviour under monotonic and cyclic loading. 5.3. An improvement of the model for bonded geomaterials An improvement of the model described in the preceding section can be introduced to reproduce the mechanical behaviour of bonded soils, weak rocks and other materials of a similar kind. Following the framework introduced by Gens and Nova [26] and Lagioia and Nova [27] two basic concepts lie the representation of this mechanical behaviour: the fundamental role played by yield pheno-
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
M
q
A pt
2925
B pc
C pc0
p’
Fig. 1. Successive yield surfaces for increasing degrees of bounding. Surface A corresponds to unbounded material.
mena and the need for considering the observed behaviour of the bonded material in relation with the behaviour of the equivalent unstructured one. A typical yield surface valid for the behaviour of the soil in an unbonded state, in the p0–q plane, is represented by curve A in Fig. 1. As the amount of bonding increases the field surface must grow towards the right to account for the fact that higher mean stresses can be applied to the material without causing it to yield. Furthermore, bonding also provides the sample with real cohesion and tensile strength that is reflected in the fact that yield surfaces are enlarged also towards the left of the stress diagram. Fig. 1 shows the successive yield surfaces for increasing degrees of bonding. Two parameters define the new enlarged yield locus: pc0 that controls the yielding of the bonded soil in isotropic compression and pt which is related to the cohesion and tensile strength of the material. Both pc0 and pt increase with the magnitude of bonding. We can assume that the degradation of the material (decrease in bonding) is related to some kind of damage measure, that will in turn depend on plastic strains. Lagioia and Nova [26] proposed simple laws to describe the debonding effect on a calcarenite material. The evolution of pt is governed by pt ¼ pt0 expðqt ed Þ;
ð84Þ
where pt0 and qt are two constitutive parameters and ed is the accumulated plastic volumetric strain. It appears reasonable to assume that changes of the yield locus will be controlled by two different phenomena: conventional plastic hardening (or softening) for an unbonded material and bond degradation. In that case, the plastic modulus introduced in 80 can be improved introducing Hb such as : HL ¼ ½H0 p Hb Hf ðHv þ Hs Þ HDM
New factors p
; Hf ; Hv
and
HDM
ð85Þ
must be defined by introducing the parameter pt :
p ¼ p0 þ pt ; 4 g Hf ¼ 1 ; gf g Hv ¼ 1 ; Mg cDM fmax ¼ ; HDM f g
1=a
af where g ¼ q=ðp0 þ pt Þ is the new stress ratio, f ¼ ðp0 þ pt Þ ½1 ð1þa Þ is the new mobilized stress f Mf function and fmax is the maximum value of the new mobilized stress function that can be replaced by the distance pc0 in Fig. 1. The term Hb is also controlled by the evolution of bonding. We propose a simple definition that depends only on the plastic volumetric strain:
2926
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Hb ¼ b1 ed expðb2 ed Þ:
ð86Þ
It can be seen that value of Hb decreases when the volumetric plastic strain increases (i.e. when debonding occurs) and in the limit case, when destructuration is complete, Hb becomes zero. In this case, the new plastic modulus defined in (85) coincides with the plastic modulus introduced in (80). To check the proposed model performance we will use the results of the laboratory tests of Lagioia and Nova on the Gravina calcarenite. Tables 1 and 2 list the 15 parameters necessary to characterize the constitutive model used in the following back prediction tests. The 4 constants in Table 2 correspond to the new parameters needed to represent the debonding phenomena. Fig. 2 compares the experimental data from Lagioa and Nova [26] with the model predictions for a specimen in an isotropic compression test. In this case the expression of the plastic modulus reduces to: cDM fmax HL ¼ ½H0 p Hb : ð87Þ f Figs. 3 and 4 compare experimental data and model predictions for specimens sheared in drained conditions under initial isotropic confining pressures of 1300 and 2000 kPa. The agreement is in general
Table 1 Pastor–Zienkiewicz model parameters for the Gravina calcarenite K (kPa)
G (kPa)
Mg
ag
Mf
af
H0
b0
b1
pc0 (kPa)
c
76,923
75,757
1.657
0.2
0.9
0.2
20
20
3000
2300
10
Table 2 New model parameters for the Gravina calcarenite b1
b2
pt0 (kPa)
qt
5e)6
36
142
1000
Experimental Calculated
1.15 1.1
void ratio
1.05 1 0.95 0.9 0.85 0.8 0
1000 2000 3000 mean effective stress p’ [kPa]
4000
Fig. 2. Isotropic compression test: experimental data from Lagioa and Nova [26] and calculated curves.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 4000
0
2927
Experimental Calculated
-0.05 volume strain
q [kPa]
3000
2000
-0.1 -0.15 -0.2
1000 -0.25 0
0
0.1 0.2 axial strain
0.3
0
0.1 0.2 axial strain
0.3
Fig. 3. Drained constant cell pressure test (r03 ¼ 1300 kPa): experimental data from Lagioa and Nova [26] and calculated curves.
5000
0
volume strain
q [kPa]
Calculated
-0.05
4000
3000
2000
1000
0
Experimental
-0.1 -0.15 -0.2 -0.25
0
0.1 0.2 axial strain
0.3
0
0.1 0.2 axial strain
0.3
Fig. 4. Drained constant cell pressure test (r03 ¼ 2000 kPa): experimental data from Lagioa and Nova [27] and calculated curves.
Experimental
1000
1000
Calculated
q [kPa]
1500
q [kPa]
1500
500
0
500
0
1000
2000 p’ [kPa]
3000
0
0
0.05 0.1 axial strain
0.15
Fig. 5. Oedometric test: experimental data from Lagioa and Nova [26] and calculated curves.
2928
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
good for both the stress–strain and the volumetric strain–axial strain curves. Fig. 5 shows the results for an oedometric test.
6. Implicit integration of generalized plasticity models In this section, a method of performing the fully implicit integration of a Generalized Plasticity model for sands is outlined. The algorithm presented here is general, since it has been devised for the integration of constitutive equations in the whole stress space [28,29]. The formulation involves an elastic predictor to give a first estimate of the stress change, coupled with an iterative sub-algorithm whose objective is to ensure that, on convergence, the constitutive behaviour is satisfied. In the following, the basic steps of the algorithm will be explained. Let us consider the finite element step n þ 1 and the incremental strain Denþ1 . The material is initially assumed to behave elastically and a trial stress state renþ1 is computed by adding the elastic stress increment associated with the incremental strain to the accumulated stress at the beginning of the increment: Z renþ1 ¼ rn þ Drenþ1 ¼ rn þ Des;nþ1 : Denþ1 ¼ rn þ Det;nþ1 : de; ð88Þ Denþ1
where Des;nþ1 is the elastic constitutive matrix. According to the non-linear elastic law proposed by Pastor et al. [25], the procedure assumes that the elastic moduli are not constant and accounts for their variation over the increment. Hence in Eq. (88) a secant elastic matrix has been used, which can be expressed just as a function of the secant bulk modulus. On the basis of the projection of the trial stress increment Drenþ1 onto the direction vector nn , it is possible to check whether loading, unloading or neutral loading is occurring. Under neutral loading ðnn : Drenþ1 ¼ 0Þ the material is assumed to have remained elastic at the Gauss point level. In that case the trial stress state computed in Eq. (88) is accepted as the final solution and there are neither any incremental plastic strains nor any incremental change in the hardening modulus. Under loading ðnn : Drenþ1 > 0Þ or unloading ðnn : Drenþ1 < 0Þ conditions, plastic deformations occur and an integration procedure must be adopted. Using an implicit scheme, the integration of the constitutive law results in a set of non-linear equations for the unknowns rnþ1 (stresses), HL=U ;nþ1 (plastic modulus), Dknþ1 (plastic multiplier) and Ks;nþ1 (bulk modulus): rnþ1 ¼ renþ1 Dknþ1 Des;nþ1 : ngL;nþ1 ;
ð89Þ
HL;nþ1 ¼ HL ðrnþ1 ; ep ðDknþ1 ÞÞ;
ð90Þ
Dknþ1 ¼
nnþ1 : Des;nþ1 : Denþ1 ; HL;nþ1 þ nnþ1 : Des;nþ1 : ngL;nþ1
Ks;nþ1 ¼ Ks ðpnþ1 ; Deenþ1 Þ;
ð91Þ ð92Þ
where pnþ1 is hydrostatic part of the stress tensor. To solve this non-linear local problem, a numerical method must be used following Crisfield [30]. It is advantageous to define a vector of residuals r:
ð93Þ rTnþ1 ¼ rT1;nþ1 r2;nþ1 r3;nþ1 r4;nþ1
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2929
in which r1;nþ1 , r2;nþ1 , r3;nþ1 and r4;nþ1 represent the difference between the current values of the abovementioned primary variables and the values obtained from the integration scheme. Iterations are then introduced to reduce residuals to zero, while the primary variables must satisfy the constitutive equations. After having collected the unknowns of the problem in vector v: vTnþ1 ¼ rTnþ1
HL;nþ1
Dknþ1
Ks;nþ1
ð94Þ
the iterative scheme can be expressed as follows: iþ1 iþ1 i riþ1 nþ1 ¼ rnþ1 þ Jr;v dvnþ1 ¼ 0;
ð95Þ
where the superscript i refers to the iteration number at the Gauss point level and Jr;v represents the Jacobian of the residual rnþ1 with respect to the unknowns. Solving the linear system (95) we obtain, at each iteration, the update of the primary variables, given by !1 orinþ1 iþ1 dvnþ1 ¼ rinþ1 : ð96Þ ovinþ1 Hence the unknowns of the problem can be simultaneously updated: i iþ1 riþ1 nþ1 ¼ rnþ1 þ drnþ1 ; iþ1 i iþ1 ¼ HL;nþ1 þ dHL;nþ1 ; HL;nþ1 iþ1 iþ1 Dknþ1 ¼ Dkinþ1 þ dDknþ1 ; iþ1 i iþ1 ¼ Ks;nþ1 þ dKs;nþ1 : Ks;nþ1
The iteration is considered to have converged if riþ1 nþ1 ¼ 0 within some prescribed tolerance. 6.1. Consistent tangent operator An important issue in computational plasticity is the derivation of the tangent elastoplastic matrix. In what follows we shall derive the so-called consistent elastoplastic matrix [31,32] that results in a quadratic rate of convergence when the Newton–Raphson method is used at the global equilibrium level, provided the current solution lies within the convergence region. To compute this matrix, the consistent tangent modulus at each Gauss point is needed. Differentiation of Eq. (89) with respect to the strains, taking into account all the variables appearing in the integration algorithm, gives: drnþ1 ornþ1 ornþ1 ovnþ1 ¼ þ : denþ1 oenþ1 ovnþ1 oenþ1
ð97Þ
Furthermore, during the iterative procedure we have drnþ1 ornþ1 ornþ1 ovnþ1 ¼ þ ¼ 0; denþ1 oenþ1 ovnþ1 oenþ1
ð98Þ
which can be rewritten as 1 ovnþ1 ornþ1 ornþ1 ¼ : oenþ1 ovnþ1 oenþ1
ð99Þ
2930
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
By substituting Eq. (99) in Eq. (97) we obtain the elastoplastic consistent tangent stiffness tensor: drnþ1 ornþ1 ornþ1 1 ornþ1 ¼ ðJr;v Þ ¼ Dcon t;nþ1 ; denþ1 oenþ1 ovnþ1 oenþ1
ð100Þ
where Jr;v stands for the Jacobian of the residual rnþ1 with respect to the primary variables, as specified above.
7. Diffuse mechanisms of failure This section includes a series of numerical examples to illustrate the problem of diffuse failure mechanisms under different circumstances. 7.1. Liquefaction of a sand layer under earthquake loading Fig. 6 shows a fully saturated sand layer subjected to a horizontal earthquake. The seismic input is the accelerogram given in Fig. 7 that correspond to the NS component accelerogram of the El Salvador earthquake 13/01/2001 in Santa Tecla. The sand layer is modelled by a sand column with both sides and the bottom assumed impermeable. Pore pressures are assumed to be zero at the surface of the sand layer.
Fig. 6. Soil layer problem: finite element mesh.
5
Comp. NS
4 3 acc. [m/s-2]
2 1 0 -1 -2 -3 -4 -5 0
10
20
30
40 50 time [s]
60
70
80
Fig. 7. Seismic input: NS component accelerogram of the El Salvador earthquake 13/01/2001 in Santa Tecla.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2931
Table 3 Data used in sand layer analysis qs (kg/m3 )
qw (kg/m3 )
Ks (Pa)
Kw (Pa)
n
k (m/s)
2720
980
1e)20
1.092e)9
0.363
2.1e)3
Table 4 Pastor–Zienkiewicz soil model parameters K0 (kPa)
G0 (kPa)
Mg
Mf
ag ¼ af
H0
b0
b1
HU (kPa)
cDM ¼ cU
45,000
22,500
1.5
0.4
0.45
350
4.2
0.2
6e)3
2
Repeated boundary conditions are assumed for the lateral nodes, which implies that the displacement of a right hand side node equals that of the corresponding left hand side node. The finite element mesh is shown in Fig. 6 and consists of stabilized four node quadrilateral elements, where bilinear shape functions are used both for displacements and pressures. The material is assumed to be a very loose sand. Constitutive behaviour has been modelled using the Pastor–Zienkiewicz model described in the previous section. Material properties and other relevant data used in the analysis are listed in Tables 3 and 4. The results can be seen in Fig. 8, which depicts the evolution of the pore pressure and the mean effective stress at stations A–E. It is interesting to note that the pore pressure increases and reaches a plateau between 7 and 17 s depending on the depth, and the mean effective stress decreases to zero. Fig. 9 provides further insight into the phenomenon of liquefaction. We have plotted profiles of pore pressure at several instants together with the vertical stress. In this way, when the excess pore pressure coincides with the vertical effective stress, it is possible to see how liquefaction extends from the surface. At 17 s all the column is liquefied. 7.2. Liquefaction failure of an embankment under earthquake action The case we will be considered next is that of an earthquake induced flowslide in very loose saturated sand. Here it is important to notice that proposed model is intended to describe the initiation of failure. In x 10
5
5
x 10 2 E
E
D
1.5 D
2 C
p’ [Pa]
excess pore-pressure [Pa]
3
C
1
B
B
1
0.5 A
A 0
0
10 time [s]
20
0
0
10
20
time [s]
Fig. 8. Evolution of excess pore pressure and mean effective stress at several stations.
2932
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 30 t=4s t=8s t=12s t=16s t=20s
25
height [m]
20
15
10
5
0
0
1 2 3 excess pore-pressure [Pa]
4 x 10
5
Fig. 9. Vertical profiles of excess pore pressure at different times showing extent of liquefaction.
order to cope with the propagation of the flowslide, the authors have recently proposed a fluid-like, depth integrated model formulated in an Eulerian framework. Flowing material is described by a suitable rheological law [5,33]. The problem consists of an embankment 10 m in height with slopes 2:1, founded on a sand layer which extends 10 m in depth and lies on a rigid rockbed. The material of both the embankment and the foundation is a very loose saturated sand. Initial conditions correspond to geostatic equilibrium under gravity forces. Suction at the surface has been assumed to be equal to )20 kPa. The finite element mesh can be seen in Fig. 10 and consists on 500 quadrilaterals with eight nodes for displacements and 4 for pore pressures. A reduced integration rule has been used in the solid part to avoid locking. The number of nodes is 1611, with 3535 degrees of freedom (Fig. 10). Loading is applied by prescribing horizontal accelerations at the base. We have used the input motion defined in Fig. 7. A simplified absorbing boundary condition has been applied at lateral boundaries. Concerning pore pressures, it has been assumed that no flux occurs at artificial boundaries, and the constant value of )20 kPa has been kept at the surface. The behaviour of the loose material is represented using the Pastor–Zienkiewicz model for sand, Tables 3 and 4. To better understand this behaviour we simulate a loading-unloading-reloading undrained triaxial test. The material loses progressively its resistant capacity, the effective stress decreases and the pore pressure increases, Fig. 11. The results can be seen in Figs. 12–14, where the contours of pore pressure, deviatoric plastic strain 0 0 invariant epl D , p =p0 and displacements are given at different times.
Fig. 10. Finite element mesh of the embankment.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939 20
B
20
10 5
F
0 -5
A
E
0
20
40 60 p' (en kPa.)
80
C
10 F
5 0 -5
D
-10
B
15
C q (en kPa.)
q (en kPa.)
15
A
E
-10 -6.E-06
100
-3.E-06
100 E pw (en kPa.)
2933
D 0.E+00 3.E-06 E22 (en %)
6.E-06
C
80 60 40 20 0 -6.E-06
A -3.E-06
0.E+00
3.E-06
6.E-06
E22 (en %)
Fig. 11. Non-drained triaxial simulation of the loose sand: loading, unloading, reloading path.
Fig. 12. Evolution of excess pore pressure contour [Pa].
The deviatoric plastic strain invariant epl D is defined as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i pl pl 2 pl pl 2 pl pl 2 2 ¼ ðe e Þ þ ðe e Þ þ ðe e Þ epl : D 1 2 2 3 3 1 3
2934
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Fig. 13. Evolution of p0 =p00 contour.
Fig. 14. Deformation contour [m].
The ratio between the mean effective confining pressure and its initial value p0 =p00 has been used as an indicator of the extent of the liquefied zones. From these results, it can be concluded that failure of the embankment is caused by liquefaction of the outer liquefied zones. 7.3. Failure of an embankment on a collapsible soil under earthquake loading The following example illustrates the collapse of a bonded material under earthquake loading. The problem is similar to that presented in the previous section but drained conditions are now assumed and the bonded material is now represented using the enhanced Generalized Plasticity model with the new parameters listed in Table 5. The evolution of plastic strain during the earthquake is plotted in Fig. 15. The plastic strain is concentrated in the outer zones of the slopes. Fig. 16 shows the evolution of the pt parameter. This parameter decreases as the volumetric plastic strain increases and becomes zero. This indicates the evolution of the debonding and the progressive reduction of the apparent cohesion of the material.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2935
Table 5 New model parameters for the enhanced Generalized Plasticity model b1
b2
pt0 (kPa)
qt
pc0 (kPa)
5e)6
36
100
1000
2000
Fig. 15. Evolution of deviatoric plastic strain invariant epl D contour.
Another interesting result shows evidence of the mechanical collapse in Figs. 17 and 18, where we have plotted the evolution of the void ratio throughout the embankment and at a critical point A. Finally, the deformation contour at t ¼ 0:88 s is plotted in Fig. 19. 7.4. A note on air liquefaction failure If the soil is partially saturated, its behaviour depends on the coupling between the solid skeleton and the pore air and water. In the limiting case of a dry soil, the air has to flow out from the pores for the material to consolidate, but typical air consolidation times are much shorter than those of the water. Therefore, in practical cases, the role of air pressures is neglected, as the characteristic time of loading is much longer than consolidation time. However, it is possible to imagine situations with much shorter loading times, where coupling between pore air and soil skeleton plays a paramount role. This is the case of fluidized granular beds, just to mention a particular example of industrial interest. Bishop [34] describes the case of Jupille flowslide, which happened in Belgium in February 1961. A tip of uncompacted fly ash located in the upper part of a narrow valley collapsed, and the subsequent flowslide travelled for about 600 m at very high speeds (130 km/h) until it stopped. A triggering mechanism was suggested to be ‘‘collapse due to undermining of a steep, partly saturated and slightly cohesive face’’. Bishop referred to Calembert and Dantinne [35], who pointed out the role of the entrapped air. The mechanism of pore air consolidation can explain the ‘‘sort of fog’’ which formed above the flowing material, as warmer air met the colder winter air in the exterior. The ‘‘honeycomb’’ like structure of soil in the tip was responsible for its sudden collapse when movement started.
2936
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
Fig. 16. Evolution of pt parameter [kPa].
Fig. 17. Evolution of void ratio contour.
Dry liquefaction can be exhibited also by soils of volcanic origin. Some of the catastrophic landslides during the El Salvador earthquake of 13th of January 2001 can be explained by this mechanism, the collapse of the material under the dynamic loading could make the pore air pressure increase.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2937
collapse
0.57
void ratio
0.56
0.55
0.54
0.53 0
0.2
0.4 0.6 time [s]
0.8
1
Fig. 18. Evolution of void ratio at point A.
Fig. 19. Deformation contour [m].
8. Conclusions We have proposed in this paper a theoretical framework to model the initiation mechanisms of catastrophic landslides. First, the equations describing the coupling have been presented following an Eulerian approach based on mixture theory. The main advantage is that they provide a unified formulation that can also be used for the study of a subsequent propagation phase [5,6]. The system of partial differential equations is then discretized using the classical Galerkin Finite Element Method. An enhanced model in the framework of the Generalized Plasticity has been presented to describe the mechanical behaviour of bonded soils. This model is able to reproduce the mechanical collapse phenomenon and to characterize the basic features of a sand behaviour under monotonic and cyclic loading when debonding occurred. An implicit integration scheme for generalized plasticity has been introduced to guarantee the efficiency of the numerical convergence. Finally, applications to diffuse failure have been presented. Acknowledgements The authors gratefully acknowledge the financial support provided by the Spanish Agency for International Cooperation (AECI), by the Spanish Ministry of Science and Technology (Project Andes and Ram on y Cajal contract) and by the European Union (Project Diga).
2938
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
References [1] R. Dikau, D. Brundsen, L. Schrott, M.L. Ibsen, Landslide Recognition, John Wiley and Sons, 1996. [2] I. Vardoulakis, J. Sulem, Bifurcation Analysis in Geomechanics, Blakie Academic and Professional, 1995. [3] M. Pastor, C. Tamagnini (Eds.), Numerical Modelling in Geomechanics, Revue Francßaise du genie civil, Hermes-Lavoisier, Paris, vol. 6, 2002. [4] P. Mira, An alisis por elementos finitos de problemas de rotura de geomateriales, Ph.D. thesis, Universidad Politecnica de Madrid, Madrid, 2002. [5] M. Pastor, M. Quecedo, J.A. Fernandez Merodo, M.I. Herreros, E. Gonzalez, P. Mira, Modelling tailing dams and mine waste dumps failures, Geotechnique 52 (8) (2002) 579–591. [6] J.A. Fern andez Merodo, Une Approche a la modelisation des glissements et des effonfrements de terrains: Initiation et Propagation, Ph.D. thesis, Ecole Centrale de Paris, December 2001. [7] F. Darve, F. Laouafa, Constitutive equations and instabilities of granular materials, in: G. Capriz, V.N. Ghionna, P. Giovine (Eds.), Modeling and Mechanics of Granular and Porous Materials, Birkh€auser, 2002, pp. 3–44. [8] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155–164. [9] M.A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26 (1955) 182–185. [10] O.C. Zienkiewicz, C.T. Chang, P. Bettess, Drained, undrained, consolidating dynamic behaviour assumptions in soils, Geotechnique 30 (1980) 385–395. [11] O.C. Zienkiewicz, T. Shiomi, Dynamic behaviour of saturated porous media: the generalised Biot formulation and its numerical solution, Int. J. Num. Anal. Methods Geomech. 8 (1984) 71–96. [12] O.C. Zienkiewicz, A.H.C. Chan, M. Pastor, D.K. Paul, T. Shiomi, Static and dynamic behaviour of soils: a rational approach to quantitative solutions. I. Fully saturated problems, Proc. R. Soc. Lond. A 429 (1990) 285–309. [13] O.C. Zienkiewicz, Y.M. Xie, B.A. Schrefler, A. Ledesma, N. Bicanic, Static and dynamic behaviour of soils: a rational approach to quantitative solutions. II. Semi-saturated problems, Proc. R. Soc. Lond. A 429 (1990) 311–321. [14] O.C. Zienkiewicz, A.H.C. Chan, M. Pastor, B. Schrefler, T. Shiomi, Computational Geomechanics, John Wiley and Sons, 2000. [15] R.L. Lewis, B.A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, John Wiley and Sons, 1998. [16] O. Coussy, Mechanics of Porous Media, John Wiley and Sons, Chichester, 1995. [17] R. de Boer, Theory of Porous Media, Springer-Verlag, Berlin, 2000. [18] A.H.C. Chan, A unified finite element solution to static and dynamic geomechanics problems, Ph.D. dissertation, University college of Swansea, Wales, 1988. [19] M.G. Katona, O.C. Zienkiewicz, A unified set of single-step algorithms, Part 3: the beta-m method, a generalisation of the Newmark scheme, Int. J. Num. Methods Engrg. 21 (1985) 1345–1359. [20] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method (4th ed.), vol. 2, McGraw-Hill, 1991. [21] S.W. Sloan, M.F. Randolph, Numerical prediction of collapse loads using finite element methods, Int. J. Numer. Anal. Methods Geomech. 6 (1982) 47–76. [22] M. Pastor, M. Quecedo, O.C. Zienkiewicz, A patch test for mesh alignment effects in localized failure, Commun. Numer. Methods Engrg. 11 (1995) 1015–1024. [23] M. Pastor, T. Li, X. Liu, O.C. Zienkiewicz, Stabilized low order finite elements for failure and localization problems in undrained soils and foundations, Comput. Methods Appl. Mech. Engrg. 174 (1999) 219–234. [24] P. Mira, M. Pastor, T. Li, X. Liu, A new stabilized enhanced strain element with equal order of interpolation for soil consolidation problems, Comput. Methods Appl. Mech. Engrg. 192 (2003) 37–38, and 4257–4277. [25] M. Pastor, O.C. Zienkiewicz, A.H.C. Chan, Generalized plasticity and the modelling of soil behaviour, Int. J. Numer. Anal. Methods Geomech. 14 (1990) 151–190. [26] A. Gens, R. Nova, Conceptual bases for a constitutive model for bonded soils and weak rocks, in: Proceedings of International Symposium on Geotechnical Engineering of Hard Soils–Soft Rocks, Balkema, Rotterdam, 1993, pp. 485–494. [27] R. Lagioia, R. Nova, An experimental and theoretical study of the behaviour of a calcarenite in triaxial compression, Geotechnique 45 (4) (1995) 633–648. [28] L. Tonni, Modellazione numerica del comportamento di terreni granulari con la Plasticita Generalizzata. Ph.D. thesis, Politecnico di Torino, 2002. [29] O.M. Heeres, R. de Borst, Implicit integration of standard and non-standard elasto-plastic models, ECCOMAS 2000, September 2000, Barcelona (Spain), pp. 11–14. [30] M.A. Crisfield, Non-Linear Finite Element Analysis of Solids and Structures, vol. 1–2: Essentials, John Wiley and Sons, 1991. [31] K. Runesson, A. Samuelsson, Aspects on numerical techniques in small deformation plasticity, in: J. Middleton et al. (Ed.), NUMETA 85, Numerical Methods in Engineering: Theory and Applications, vol. 1, A.A. Balkema, Rotterdam, 1985, pp. 337– 348.
J.A. Fernandez Merodo et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2911–2939
2939
[32] J.C. Simo, R.L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Methods Appl. Mech. and Engrg. 48 (1985) 101–118. [33] M. Quecedo, M. Pastor, M.I. Herreros, J.A. Fernandez Merodo, Numerical modelling of the propagation of fast landslides using the finite element method, Int. J. Num. Meth. Eng. 59 (2004) 755–794. [34] A.W. Bishop, The stability of tips and spoil heaps, Quart. J. Engrg. Geol. 6 (1973) 335–376. [35] L. Calembert, R. Dantinne, The avalanche of ash at Jupille (Liege) on February 3rd, 1961. From: The commemorative volume dedicated to Professeur F. Campus, Liege, Belgium, 1964, pp. 41–57.