Deformation and localization analysis of partially saturated soil

Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910 www.elsevier.com/locate/cma

Deformation and localization analysis of partially saturated soil W. Ehlers *, T. Graf, M. Ammann Institute of Applied Mechanics (Civil Engineering), University of Stuttgart, Pfaffenwaldring 7, 70550 Stuttgart, Germany Received 21 March 2003; received in revised form 2 September 2003; accepted 2 September 2003

Abstract Deformation and localization analysis is a crucial issue and has thus been intensively investigated in the last decades. However, in contrast to solid mechanical problems, geotechnical applications do not only concern a single solid material, the soil, but they also affect the pore-fluids, water and air, and, consequently, the coupling of the solid deformation with the pore-fluid flow. As a result, both the deformation and the localization analysis must be applied to a triphasic material consisting of the soil skeleton, the pore-water and the pore-gas, which, in geotechnical engineering, is known as unsaturated or partially saturated soil. Based on a continuum mechanical approach, unsaturated soil can be described within the well-founded framework of the Theory of Porous Media (TPM), thus including saturated soil (solid matrix and pore-water) as well as empty soil (solid matrix and pore-gas) as special cases. It is the goal of the present contribution to investigate the deformation and the localization behavior of unsaturated soil and to exhibit the influence of the solid–fluid coupling on the localization analysis. In the framework of a triphasic formulation, unsaturated soil is considered as a materially incompressible elasto-plastic or elasto-viscoplastic skeleton saturated by two viscous pore-fluids, a materially incompressible pore-liquid and a materially compressible pore-gas. Assuming quasi-static situations, the numerical computations proceed from weak formulations of the momentum balance of the overall triphasic material together with the mass balance equations of the pore-fluids and Darcy-like relations for the seepage velocities. As a result, a system of strongly coupled differential-algebraic equations (DAE) occurs, which is solved by use of the FE tool PANDAS. In particular, various initial boundary-value problems are treated on the basis of time- and space-adaptive methods, thus demonstrating the efficiency of the overall formulation. Furthermore, the influence of the pore-gas constituent on the material behavior of partially saturated soil is studied with respect to fluid-flow simulations or embankment and slope failure problems.
2004 Elsevier B.V. All rights reserved. Keywords: Theory of Porous Media (TPM); Triphasic material; Solid elasto-plasticity; Viscous pore-fluids; Deformation and localization analysis; Adaptive computations

*

Corresponding author. E-mail address: [email protected] (W. Ehlers). URL: www.mechbau.uni-stuttgart.de/ls2.

0045-7825/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.09.026

2886

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

1. Introduction In geotechnical engineering, there is a rapidly growing interest in the coupled analysis of the soil deformation and the pore-fluid flow. To set an example, it is evident that dykes or embankments, which are built to protect the environment from the elements of water, are generally loaded by gravitation and by a water table on one side. As a consequence, one observes both the deformation of the porous embankment structure and a flow process of the pore-content, water and air, cf. Fig. 1. Furthermore, if the water table rises or decreases rapidly driven, e.g., by natural hazards, stability problems occur, which are initiated by the localization of plastic solid deformations in narrow bands (shear bands). In order to correctly predict the overall behavior of such constructions, it is necessary to carefully investigate the coupled solid–fluid behavior of a biphasic material (solid matrix and pore-liquid) or of a triphasic material (solid matrix, poreliquid and pore-gas). In the triphasic case, which is of course more general and contains the ‘‘saturated’’ material (solid matrix and pore-liquid) or the ‘‘empty’’ material (solid matrix and pore-gas) as special cases, a partially saturated porous solid material is considered and described within the well-founded Theory of Porous Media (TPM). Concerning the general TPM approach, the reader is referred, e.g., to the work by Truesdell and Toupin [1], Bowen [2], de Boer [3] or Ehlers [4,5]. In the present contribution, the triphasic material under consideration is assumed to consist of a materially incompressible, elasto-plastic or elastoviscoplastic, cohesive-frictional solid skeleton saturated by two viscous fluids, a materially incompressible pore-liquid (water) and a materially compressible pore-gas (air). Concerning the localization analysis of solids as well as of multiphasic materials, it is well known that the computation of shear band localizations generally leads to a mathematically ill-posed problem that has to be regularized by means of additional assumptions as, e.g., the consideration of additional kinematical degrees of freedom in the sense of a micropolar continuum [6–9], by the assumption of viscoplastic skeleton behavior [10–12] or by proceeding from further techniques like the consideration of gradient plasticity models [13] or of non-local approaches [14]. With regard to the organization of the following contribution, there is firstly the triphasic material for the description of partially saturated soil presented in Section 2. Therein, proceeding from a geometrically linear framework of the solid deformation, the soil matrix is considered as a cohesive-frictional material governed by a general elasto-viscoplastic description of the solid stresses based on a Hookean type elasticity law and a single-surface yield criterion [15] together with an additional plastic potential function in order to catch the non-associativity of the plastic behavior of geomaterials. Viscoplasticity is not only a convenient tool for the regularization of the shear band problem, but it also represents the basic matrix behavior of typical embankment materials like clayey silt. Furthermore, it additionally allows for a simple transfer to geomaterials plasticity. Concerning the pore-fluids, the effective pressure of the materially incompressible liquid acts as a Lagrangean, whereas the effective gas pressure is assumed to be governed by the ideal gas law. The interaction between the constituents, solid matrix, liquid and gas, is taken into consideration by additional constitutive relations for the momentum production terms, thus leading to Darcy-like relations for the seepage velocities. Finally, the fluids interact by a capillary-pressure-saturation relation based on relative permeabilities to consider the mutually interacting pore-fluid mobilities. In Section 3, the numerical

Fig. 1. Pore-liquid distribution in an embankment.

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2887

treatment of the strongly coupled solid–fluid problem is described. Based on quasi-static initial boundaryvalue problems, the numerical computations proceed from weak formulations of the momentum balance of the overall triphasic material together with the volume and mass balance equations of the pore-fluids. As a result, a system of strongly coupled differential-algebraic equations (DAE) occurs, which can be solved by use of the FE tool PANDAS. 1 This tool [16] is designed for the solution of volumetrically coupled initial boundary-value problems, where, if necessary, time- and space-adaptive methods [12,17,18] can be taken into consideration. Numerical examples are presented in Section 4, thus demonstrating the efficiency of the overall formulation. In particular, the draining of a soil column shows the differences of the triphasic in comparison to a biphasic formulation based on the Reynolds assumption, whereas the two-dimensional computation of diverse embankment problems together with the fully three-dimensional investigation of a slope failure situation demonstrate the full capability of the model by exhibiting the coupled behavior of the water and gas flow with the solid deformation ranging from very small strains to shear band localizations.

2. A triphasic continuum embedded in the TPM 2.1. Preliminaries and basic concepts Based on the fundamental concepts of the Theory of Porous Media, one proceeds from the hypothesis that the individual constituents, solid, liquid and gas, composing the triphasic material under consideration are in a state of ideal disarrangement, i.e., they are statistically distributed over the control space. This assumption together with the prescription of a real or of a virtual averaging process leads to a model u of superimposed and interacting continua ua : [ u ¼ ua : ð1Þ a

Therein, a is the constituent index with a ¼ S (solid), a ¼ L (liquid) and a ¼ G (gas). Furthermore, as a result of the assumed averaging process, the mathematical functions for the description of the geometrical and physical properties of the individual materials are field functions defined all over the control space. Thus, the volume V of the triphasic aggregate B is obtained from the sum of the partial volumes of the constituents ua in B: Z Z Z X V ¼ dv ¼ V a ; where V a ¼ dva ¼: na dv: ð2Þ B

a

B

B

Following this, the quantity na is defined as the local ratio of the volume element dva of a given constituent ua with respect to the volume element dv of the overall medium u: na ¼

dva : dv

ð3Þ

The relations (2) and (3) represent the concept of volume fractions. Since, in general, there is no vacant space in the overall medium, this concept directly leads to the saturation condition X na ¼ 1: ð4Þ a

1 Porous media Adaptive Nonlinear finite element solver based on Differential Algebraic Systems (http://www.getpandas.com).

2888

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

Proceeding from the volume fractions na , two different density functions of each constituent ua can be introduced: qaR ¼

dma dva

and

qa ¼

dma : dv

ð5Þ

Therein, the material (realistic or effective) density qaR relates the local mass dma to the volume element dva , whereas the partial (global or bulk) density qa relates the same mass to the volume element dv. Following this, the density functions are related to each other by
aR a  q dv a dm ¼ ð6Þ ! qa ¼ na qaR : qa dv Based on the above relation, it is immediately evident that the property of material incompressibility (qaR ¼ const.) of a constituent ua is not equivalent to the property of bulk incompressibility, since the partial density qa can still vary through changes in the volume fraction na . In addition to the volume fractions, the triphasic model requires the introduction of saturation functions sb ¼

nb ; nF

where nF ¼ nL þ nG

ð7Þ

and b 2 fL; Gg. Note in passing that the overall fluid volume fraction nF also represents the porosity variable. Furthermore, one easily concludes to the saturation constraint sL þ sG ¼ 1:

ð8Þ

Following this, it is seen that (8) relates the saturation condition to the pore content, whereas (4) relates it to the overall material. 2.2. Kinematical relations Proceeding from the concept of superimposed continua with internal interactions and individual states of motion, the Theory of Porous Media combines the volume fraction concept with the classical results of mixture theories [1,19]. In this framework, each spatial point x of the current configuration is at any time t simultaneously occupied by material particles (material points) P a of all constituents ua . These particles proceed from different reference positions at time t0 , cf. Fig. 2. Thus, each constituent is assigned its own motion function x ¼ va ðXa ; tÞ:

ð9Þ

Fig. 2. Motion of a triphasic continuum.

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2889

As a result, each spatial point x can only be occupied by one single material point P a of each constituent ua . The assumption of unique motion functions, where each material point P a of the current configuration has a unique reference position Xa at time t0 , requires the existence of unique inverse motion functions v
1 a based on non-singular Jacobian determinants Ja : Xa ¼ v
1 a ðx; tÞ;

Ja ¼ det

ova 6¼ 0: oXa

ð10Þ

It follows from (10) that each constituent has its own velocity and acceleration fields. In the basic Lagrangean setting, these fields are given by x0a ¼

ova ðXa ; tÞ ; ot

x00a ¼

o2 va ðXa ; tÞ : ot2

ð11Þ

With the aid of the inverse motion function (10)1 , an alternative formulation of (11) leads to the Eulerian description x0a ¼ x0a ðx; tÞ;

x00a ¼ x00a ðx; tÞ:

ð12Þ

Suppose that C is an arbitrary, steady and sufficiently often steadily differentiable scalar function of ðx; tÞ. Then, the material time derivative of C following the motion of ua reads 0

ðCÞa ¼

da oC þ grad C  x0a : C¼ ot dt

ð13Þ

Therein, the operator ‘‘gradðÞ’’ denotes the partial derivative of (Æ) with respect to the local position x. Describing coupled solid–fluid problems, it is convenient to proceed from a Lagrangean description of the solid matrix by the solid displacement vector uS as the primary kinematic variable, whereas the porefluids are better described in a modified Eulerian setting by use of the seepage velocities wb describing the fluid motions with respect to the deforming skeleton material: uS ¼ x
X S ;

wb ¼ x0b
x0S :

ð14Þ

From (9) and (10)1 , one obtains the material deformation gradient Fa and its inverse F
1 a by Fa ¼ Grada x;

F
1 a ¼ grad Xa :

ð15Þ

Therein, the operator ‘‘Grada ðÞ’’ denotes the partial derivative of (Æ) with respect to the reference position Xa of ua . 2.3. Balance relations The discussion of balance relations for multiphasic materials is based on Truesdell’s ‘‘metaphysical principles’’ of mixture theories [20, p. 221]. These principles proceed from the idea that both the balance relations of the constituents ua and the balance relations of the overall medium u can be given in analogy to the balance relations of classical continuum mechanics of single-phase materials, extended by so-called production terms to describe the interaction mechanisms between the constituents. These productions are constrained by sum relations. The reader, who is interested in an extended discussion on the development of balance equations in frame of thermodynamics of porous and multiphasic materials, is referred, e.g., to [3–5,19]. ^a vanish and the partial Excluding mass exchanges between the constituents, the density productions q balance equations of mass and of linear momentum yield ðqa Þ0a þ qa div x0a ¼ 0; qa x00a ¼ div Ta þ qa ba þ ^ pa :

ð16Þ

2890

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

Therein, ‘‘divðÞ’’ is the divergence operator corresponding to ‘‘gradðÞ’’, and ba is the mass specific body force, which is usually identified with the overall gravitation g. Furthermore, ^pa is the direct term of the linear momentum production. In the above setting, the momentum productions are constrained by ^ pS þ ^ pF ¼ 0;

^F ¼ ^ where p pL þ ^ pG :

a

If u is assumed to be materially incompressible (q the volume balance

ð17Þ aR

¼ const:), the mass balance equation (16)1 reduces to

0

ðna Þa þ na div x0a ¼ 0:

ð18Þ

In case of the skeleton material, an integration of (18) yields the relation
1

nS ¼ nS0S ðdetFS Þ ;

ð19Þ

which, in the framework of a small strain approach, can be formally linearized around the natural state of uS . The result is nS ¼ nS0S ð1
DivS uS Þ;

ð20Þ nS0S

where DivS ðÞ is the divergence operator corresponding to GradS ðÞ. Furthermore, is the volume fraction of uS in the solid reference configuration. Note in passing that, in case of a small strain approach, it is generally not necessary to distinguish between GradS ðÞ and gradðÞ or between DivS ðÞ and divðÞ. Consequently, the remainder of this article is based on spatial gradient and divergence operators. Given (19) or (20), the overall fluid volume fraction nF (porosity) is coupled with the solid deformation by the solid volume fraction nS via the saturation condition (4). However, to determine the portions of the liquid and the gas constituents, an additional constitutive equation for the liquid saturation function sL or for the gas saturation function sG , respectively, is required such that nL and nG can be determined through (7) and (8). It is typical for geotechnical applications that most of the problems under study proceed from isothermal and quasi-static situations described with respect to the skeleton motion. Following this firstly means 0 0 that the material time derivative ðÞb of the fluid motion has to be related to the skeleton time derivative ðÞS by a modification of the convective part. Considering an arbitrary field function C, this transformation implies 0

0

ðCÞb ¼ ðCÞS þ grad C  wb :

ð21Þ

Furthermore, an isothermal triphasic formulation of partially saturated soil is governed by five primary variables given by the solid displacement uS , the seepage velocities wL and wG and the effective pore-fluid pressures pLR and pGR . However, proceeding from initial boundary-value problems under quasi-static conditions given by x00a  0, one obtains a coupling of the seepage velocities and the effective fluid pressures resulting from the individual fluid momentum balances in the shape of Darcy-like relations. Following this reduces the set of primary variables from five to three: the solid displacement uS and the fluid pressures pLR and pGR . The corresponding set of governing equations is then given by the vector-valued overall momentum balance corresponding to uS and the scalar-valued fluid volume or mass balance equations corresponding to pLR and pGR . In particular, taking the sum of the momentum balances (16)2 of uS , uL and uG results in the overall momentum balance 0 ¼ div T þ qg;

ð22Þ

where the constraint (17) has been taken into consideration. In addition T ¼ TS þ TL þ TG ; q ¼ nS qSR þ nL qLR þ nG qGR

ð23Þ

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2891

are defined as the overall Cauchy stress and the overall density of the triphasic material. In addition to (22), one proceeds from the liquid volume balance 0

0

ðnL ÞS þ nL divðuS ÞS þ divðnL wL Þ ¼ 0

ð24Þ

and from the gas mass balance nG ðqGR Þ0S þ ðnG Þ0S qGR þ nG qGR divðuS Þ0S þ divðnG qGR wG Þ ¼ 0:

ð25Þ

Note in passing that, by use of (6)2 , (18) and (21), the last two balances are written with respect to the skeleton motion. 2.4. Constitutive equations 2.4.1. General setting To close the triphasic model under consideration, constitutive equations are required for the partial Cauchy stress tensors Ta , the linear momentum productions ^pb of the pore-fluids, the liquid saturation sL and the effective gas pressure pGR . However, since pGR is assumed as a primary variable, the constitutive equation for pGR will be given in an inverse form as an equation for the effective density qGR . Proceeding from standard arguments of Rational Thermodynamics, an evaluation of the overall entropy inequality of the triphasic model yields the solid and the fluid stresses as well as the linear momentum productions to consist of two terms [2,4,5], where the first one is governed by the pore-pressure variables, while the second one, the so-called ‘‘extra term’’, results from the solid deformation (effective stress) or from the pore-fluid flow (frictional stress): TS ¼
nS pFR I þ TSE ; Tb ¼
nb pbR I þ TbE ; ^ pb ¼ pbR grad nb þ

ð26Þ

^ pbE :

Therein, the effective pore-pressure p :¼ pFR is obtained by the well-known Dalton’s law via p ¼ sL pLR þ sG pGR :

ð27Þ

Note in passing that Dalton’s law is more than only a convenient expression, since it can be recovered from thermodynamical considerations to satisfy the entropy principle. Furthermore, as a result of the materially incompressible pore-liquid, the effective liquid pressure pLR acts as a Lagrangean multiplier and is thus determined by the boundary conditions of the problem under study. Finally, in the remainder of this article, the pressures pLR and pGR are understood as ‘‘effective excess pressures’’ exceeding a typical surrounding pressure like, e.g., the atmospheric pressure p0 . In geotechnical applications, it can be shown by arguments of dimensional analysis [21] that the frictional fluid forces f bE ¼ div TbE can be neglected in comparison to the viscous interaction terms ^pbE . Thus, proceeding from TbE  0, the overall Cauchy stress T yields the well-known concept of effective stress [22,23]: T ¼
pI þ TSE :

ð28Þ

2.4.2. The fluid constituents Concerning the fluid materials, liquid and gas, the extra momentum productions ^pbE , the liquid saturation L s and the effective gas density qGR have to be specified by constitutive equations. In the present setting, the extra momentum productions are related to the seepage velocities wb through ^ pbE ¼
ðnb Þ2 cbR ðKbr Þ
1 wb ;

ð29Þ

2892

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

where cbR is the specific weight of ub and Kbr the corresponding relative permeability tensor. The relative permeabilities included in Kbr are related to the Darcy permeabilities included in Kb by Kbr ¼ jbr Kb :

ð30Þ

Therein, jbr is the so-called relative permeability factor depending on the saturation of ub , whereas Kb is understood as the permeability tensor of ub specified under fully saturated conditions (sb ¼ 1). Finally, Kb is related to the intrinsic permeability KS of the porous skeleton material through Kb ¼

cbR S K : lbR

ð31Þ

Note that this equation can also be used to determine the Darcy permeability tensor Kb of various fluids from the intrinsic permeability KS through the specific weights cbR and the effective shear viscosities lbR . In order to describe the deformation dependence of KS , it is assumed following Eipper [24] that  p 1
nS KS0S ; ð32Þ KS ¼ 1
nS0S where KS0S is the intrinsic permeability tensor of the undeformed skeleton, and p is a material parameter governing the exponential function (32). If an initially isotropic solid is concerned, KS0S reduces to S KS0S ¼ K0S I

ð33Þ

S governed by a single initial intrinsic permeability coefficient K0S . In analogy to (31), this coefficient is related b to an initial Darcy permeability coefficient k0S by b k0S ¼

cbR S K : lbR 0S

ð34Þ

Concerning partially saturated soil, three zones have to be distinguished, cf. Fig. 3. In the zone beneath the ground-water table (saturated domain), most of the pore-space is filled with the pore-liquid, which is mobile governed by the Darcy permeability measured under fully saturated conditions. Nevertheless, there is a small amount of trapped pore-gas in this zone with a residual saturation sG res . In a certain height above the ground-water table (empty domain), the pore-gas is mobile, whereas a small amount of the pore-liquid is trapped with the residual saturation sLres . In between these zones, there is the unsaturated or the partially

Fig. 3. Zones of a partially saturated soil.

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2893

saturated domain, respectively, where both fluids are mobile. The particular height of this domain depends on the suction properties of the soil material under study. In order to obtain a relation between the porefluid mobilities through a relation between the liquid saturation and the pore-pressures, one introduces the capillary pressure pC ¼ pGR
pLR :

ð35Þ

Following this, a relation between the liquid saturation sL and the capillary pressure pC can be given following Brooks and Corey [25] or van Genuchten [26]. In the present contribution, use is made of the van Genuchten model given by j

sLeff ðpC Þ ¼ ½1 þ ðagen pC Þ gen 


hgen

:

ð36Þ

Therein, agen , jgen and hgen are material parameters, whereas the effective saturation function describing the area between the two residual saturations is given with respect to Finsterle [27], cf. Fig. 4: sLeff :¼

sL
sLres : 1
sLres
sG res

ð37Þ

In the van Genuchten model, the relative permeability functions are given by
 h  L 1=hgen ihgen 2 L L
gen jr ¼ ðseff Þ 1
1
seff ; jG r

¼ ð1


sLeff Þcgen

h

1




1=hgen i2hgen sLeff ;

ð38Þ

where
gen and cgen are additional parameters governing the hydraulic behavior of the soil in the unsaturated domain, cf. Fig. 5. If the effective saturation vanishes, jLr is zero, and one obtains an immobile pore-liquid. On the other hand, if the effective saturation is one, jLr is one, and the pore-liquid is fully mobile with the Darcy permeability under fully saturated conditions. Concerning the pore-gas, equivalent statements hold. The reader, who is interested in further details of the mutual fluid behavior, is referred, e.g., to the work by Helmig [28]. Inserting (29) into the quasi-static fluid momentum balances (16)2 yields the Darcy-like equations nb wb ¼


Kbr ðgrad pbR
qbR gÞ; cbR

Fig. 4. Capillary-pressure-saturation relation according to [26] with agen ¼ 2  10
4 , jgen ¼ 2:3 and hgen ¼ 1:5.

ð39Þ

2894

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

Fig. 5. Relative permeability factor according to [26] with
gen ¼ 0:5 and cgen ¼ 0:333.

where nb wb is the filter velocity of ub and  1  ib ¼
bR grad pbR
qbR g c

ð40Þ

the corresponding hydraulic gradient. Finally, the effective density function of the pore-gas is assumed to be governed by Boyle’s ideal gas law: qGR ¼

pGR þ p0 G

ð41Þ

:

R h

G

Therein, R denotes the specific gas constant of the pore-gas and h the absolute Kelvin’s temperature, which, in the present investigations, is constant (h ¼ const:) due to the assumption of an overall isothermal problem. 2.4.3. The solid constituent Following the geometrically linear approach to elasto-plasticity, the solid strain tensor i 1h T eS ¼ grad uS þ ðgrad uS Þ 2

ð42Þ

is additively decomposed into an elastic and a plastic part: eS ¼: eSe þ eSp :

ð43Þ

The solid extra stress is governed by the generalized Hookean law rSE ¼ 2lS eSe þ kS ðeSe  IÞI; S

S

ð44Þ rSE

is the solid extra stress under small strain where l and k are the Lam
e constants of the porous soil and conditions (rSE  TSE ), where no difference between the Cauchy, the Kirchhoff, the Piola and the Piola– Kirchhoff stress must be made. In order to describe the plastic or the viscoplastic material properties of the skeleton material, one has to consider a convenient yield function to bound the elastic domain. Following this, the single-surface yield criterion by Ehlers [4,15] is applied, cf. Fig. 6: F ¼ U1=2 þ bI þ
I2
j ¼ 0; m

U ¼ IID ð1 þ c#Þ þ 12aI2 þ d2 I4 ; D

D 3=2

# ¼ III =ðII Þ

:

ð45Þ

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

–σ2

2895

–I

–σ1

I

–σ3

Fig. 6. Single-surface yield criterion for cohesive-frictional materials; r1 , r2 , r3 : principal stresses of rSE (tension positive).

In this yield function, which has been designed for cohesive-frictional materials, I, IID and IIID are the first principal invariant of rSE and the (negative) second and third principal invariants of the effective stress D deviator ðrSE Þ . The material parameter sets Sh ¼ fa; b; d;
; jg;

Sd ¼ fc; mg

ð46Þ

govern the shape of the yield surface in the hydrostatic plane (Sh ) and in the deviatoric plane (Sd ). Proceeding either from the viscoplastic approach or from the perfect plasticity concept, the parameters included in Sh and in Sd are kept constant during the deformation process and can be computed from standard experimental data by use of an optimization procedure [29]. Since the associated plasticity concept cannot be applied to frictional materials [30], an additional plastic potential [31] G ¼ C1=2 þ w2 I þ
I2 ; C ¼ w1 IID þ 12aI2 þ d2 I4

ð47Þ

is considered, where w1 and w2 serve to relate the dilatation angle to experimental data. From the concept of the plastic potential, it is straight forward to obtain the evolution equation (flow rule) for the plastic strain eSp via 0

ðeSp ÞS ¼ K

oG ; orSE

ð48Þ

where K is the plastic multiplier. In the framework of viscoplasticity using the overstress concept of Perzyna-type [32], the plastic multiplier included in (48) is given by   1 F ðrSE Þ r K¼ : ð49Þ g r0 Therein, hi are the Macaulay brackets, g is the relaxation time, r0 the reference stress, and r is the viscoplastic exponent. However, in the framework of elasto-plasticity, where the plastic strains are rateindependent, K has to be computed from the Kuhn-Tucker conditions F 6 0;

K P 0;

rather than from (49).

KF ¼ 0

ð50Þ

2896

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

3. Numerical treatment 3.1. Weak formulations of the governing field equations The numerical treatment of initial boundary-value problems of the triphasic material under study is based on the weak formulations of the governing field equations together with discretization methods in the space and time domains. As was mentioned above, the isothermal problem is basically governed by five independent fields: the solid displacement uS , the seepage velocities wL and wG and the effective fluid pressures pLR and pGR . However, under quasi-static conditions, Darcy-like relations (39) have been found to eliminate the seepage velocities by the effective fluid pressures. Consequently, the problem is finally governed by the variables uS , pLR and pGR corresponding, in the framework of the standard Galerkin procedure (Bubnov–Galerkin) to the balance relations (22), (24) and (25). Following the above remarks, the overall momentum balance (22) is multiplied by a test function duS and integrated over the domain X. Thus, one obtains Z ½divðrSE
pIÞ  duS þ qg  duS dv ¼ 0; ð51Þ X

where (28) and TSE  rSE have been used. Integration by parts of the first term in (51) together with the Gaussian integral theorem finally yields Z Z Z S t  duS da: ðrE
pIÞ  grad duS dv ¼ qg  duS dv þ ð52Þ X

X

Ct

Therein, t is the external load vector acting on the Neumann boundary Ct of the overall medium. Analogously, multiplication of the pore-liquid volume balance (24) with a test function dpLR , integration over the domain X, integration by parts and application of the Gaussian integral theorem yields Z h Z Z i  L 0 0 vL dpLR da; n S þ nL divðuS ÞS dpLR dv
nL wL  grad dpLR dv ¼
ð53Þ X

X

Cv

where vL ¼ nL wL  n is the efflux of liquid volume through the Neumann boundary Cv with the outward oriented unit surface normal n. Applying the same procedure as above to the pore-gas mass balance, a multiplication of (25) with a test function dpGR leads to Z h Z Z i  G GR 0 0 G GR GR G GR LR qG dpGR da: n q q divðu Þ dv
n q w  grad dp dv ¼
ð54Þ þ n dp S G S S X

X

Cq

Therein,  qG ¼ nG qGR wG  n is the efflux of gaseous mass through the Neumann boundary Cq . 3.2. Discretization in space and time In the framework of the finite element method (FEM), the spatial discretization (semi-discretization with respect to the space variable x) of the field equations (52)–(54) proceeds from extended Taylor–Hood elements based on quadratic shape functions for the solid displacement uS and linear shape functions for the effective fluid pressures pLR and pGR . In contrast to the computation of the external variables uS , pLR and pGR , the plastic strains eSp and the plastic multiplier K determined by (48) and (49) act as internal variables and must thus be computed, in the sense of the collocation method, at the integration points of the numerical quadrature. For a mesh of Nu nodes and Nq integration points, the space-discrete variables of the semi-discrete problem are collected in the vectors

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

   T u1S ; p1LR ; p1GR ; . . . ; uNS u ; pNLR ; pNGR ; u u h iT Nq q ¼ ðe1Sp ; K1 Þ; . . . ; ðeNq Sp ; K Þ :

u¼

2897



0

0

ð55Þ T

Using the abbreviation ðÞ :¼ ðÞS and the vector y :¼ ðuT ; qT Þ , one obtains the semi-discrete initial-value problem   F 1 ðt; u; u0 ; qÞ Mu0 þ kðu; qÞ
f ! Fðt; y; y0 Þ  ð56Þ  ¼0 F 2 ðt; q; q0 ; uÞ Aq0
gðq; uÞ of first order in the time variable t, where t P t0 and yðt0 Þ ¼ y0 are the corresponding initial conditions. In (56), the first equation (F 1 ¼ 0) represents the discretization of the governing field equations, where M is the generalized mass matrix, k is the generalized stiffness vector, and f is the vector of the external forces. The second equation (F 2 ¼ 0) exhibits the plastic or the viscoplastic evolution equations together with the constraints resulting from the Kuhn-Tucker conditions in case of the elasto-plastic formulation. The introduction of the matrix A formally allows for a joint formulation of elasto-viscoplastic and elasto-plastic problems. Finally, g represents the right-hand side of the evolution equations and constraints, which are element-wise decoupled as a result of their evaluation at the integration points of the finite elements [12]. As a result of the quasi-static problem under consideration, it may occur that the generalized mass matrix M is not regular. Then, the system (56) turns out to be a system of differential-algebraic equations (DAE) of index one in the time variable. Furthermore, in case of elasto-plastic solid material behavior, the matrix A is also not regular, thus additionally yielding a DAE system on the Gauss point level. Details on the solution of DAE systems can be taken from the literature, cf. the work by Ellsiepen [12] and the quotations therein. The numerical computations presented in the following section proceed from the above scheme together with time- and space-adaptive methods. In particular, the time integration as well as the time-adaptive strategy are based on one step methods together with an embedded time step control, where the solution at time tnþ1 only depends on the solution at time tn . This choice is of essential importance with respect to space-adaptive methods (refinements as well as coarsenings), since the transfer of the numerical solution thus only includes two meshes [11]. Based on the fact that it may occur that the system (56) is a DAE system of index one, it is convenient to apply diagonally implicit Runge–Kutta methods (DIRK) with suitable stability properties. Concerning the application of space-adaptive methods, no mathematically founded methods are known so far to estimate the spatial error [18]. Thus, the following procedure is applied, cf. Fig. 7, where n is a measure for the spatial error, whereas ey represents the total error in the time domain. A time step of the non-stationary problem is treated as a stationary problem, where the initial conditions are taken from the solution of the previous step. In order to estimate the spatial error of the discretized problem, the gradient-based error indicator of Zienkiewicz–Zhu type [33] is extended in such a way that all the driving quantities of saturated and unsaturated elasto-plastic and elasto-viscoplastic materials are included. Apart from the standard consideration of the effective solid stresses representing the elastic part of the problem, the error indicator is extended towards the plastic part of the strain representing the accumulated plasticity and towards the seepage velocities representing the viscosities of the pore-fluids together with the intrinsic permeability of the solid material. It was furthermore found that, in order to find a precise demarcation between the fully saturated, the partially saturated and the empty zones of initial boundaryvalue problems under study, the space-adaptivity concept has to be extended by including the effective saturation as an additional input to the overall error indicator. Further details on numerical methods applied to saturated and partially saturated porous materials can be found in the work by Schrefler and coworkers, cf., e.g., [34] and quotations therein, as well as in [5,17– 18,35–37,41,42].

2898

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

TimeStep

FullStep Mesh Coarsening (if possible)

Embedded Time Step

∆t ← ∆t/2

Embed.Step Adaptive Time Step

Ok?

TimeStep Estimate Space Error

Mesh Refinement

Estimate Time Error

MeshEstim

TimeEstim

ξ <– 1

Set ∆t new

Update State Variables

ey <– 1

Return

Return Fig. 7. Algorithm for a time- and space-adaptive step.

4. Numerical examples The numerical examples presented in this article exhibit typical problems in the framework of the deformation and localization analysis of unsaturated soil. As was described above, the computations proceed from a triphasic TPM formulation together with the numerical tools presented in the preceding section and implemented in the FE code PANDAS. In particular, the examples concern the draining of a soil column following the Liakopoulos test and the deformation and localization behavior of a partially saturated soil embankment under gravitational forces and varying heights of the water table. Finally, the triphasic formulation is applied to a three-dimensional slope failure problem driven by a raising water table. The material parameters summarized in Table 1 describe a clayey silt. These parameters are assumed to govern most of the numerical computations. However, since different values of the intrinsic permeabilities S K0S are used for the various examples, these values are given within the description of the individual problems. Furthermore, the numerical computation of the Liakopoulos test describing the draining of a soil column follows the line of the ALERT 2. Geomaterials network, where different elastic constants and relative permeability parameters have been used.

2

Alliance of Laboratories in Europe for Research and Technology (http://alert.epfl.ch).

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2899

Table 1 Material parameters of clayey silt Parameter

Symbol

Value

Lam
e constants

lS kS

5583 kN/m2 8375 kN/m2

Effective densities

qSR qLR qGR 0

2720 kg/m3 1000 kg/m3 1.23 kg/m3

Gas phase

R h p0

287.17 J/(kg K) 283 K 105 N/m2

Volume fractions

nS0S nF0S

0.54 0.46

Gravitation value

g

10 m/s2

Fluid viscosities and relative permeability parameters

lLR lGR p agen jgen hgen
gen cgen sLres sG res

10
3 Ns/m2 1.8 · 10
5 Ns/m2 1.0 2 · 10
4 2.3 1.5 0.5 0.333 0.1 0.1

Single-surface yield criterion

a b c d
j m

1.0740 · 10
2 0.1195 1.555 1.377 · 10
4 m2 /kN 4.330 · 10
6 m2 /kN 10.27 kN/m2 0.5935

Viscoplasticity

g r0 r

5 · 103 s 10.27 kN/m2 1

Plastic potential

w1 w2

1.33 0.107

G

4.1. Draining of a soil column The draining of a soil column is a typical problem for the validation of a triphasic model describing partially saturated soil. In his Ph.D. thesis of 1964 [38], Liakopoulos experimentally investigated the leaking of a fully saturated soil column of very fine Del Monte sand under the effect of gravitational forces. In the experiment, the soil column was situated in an impermeable cylindrical vessel open at the bottom and at the top such that a one-dimensional problem was obtained. This problem was chosen by the ALERT network as a benchmark for the description of coupled solid–fluid problems. Numerical solutions can be found, e.g., in the book by Lewis and Schrefler [34] or in the Ph.D. thesis by Klubertanz [39] following the line of the geomechanics group in Lausanne [40]. In particular, the benchmark is based on the assumption of a deformable linear elastic soil. However, since Liakopoulos did not measure the elastic constants of the solid

2900

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

Table 2 Material parameters of the ALERT benchmark Parameter

Symbol

Value

Lam
e constants

lS kS

464 kN/m2 1857 kN/m2

Volume fractions

nS0S nF0S

0.63 0.37

Relative permeability parameters

p agen jgen hgen
gen cgen sLres sG res

7.0 2 · 10
5 1.5 1.03 3.5 0.333 0 0.01

material, the Young’s modulus and the Poisson’s ratio have been assumed by members of the ALERT network to yield E ¼ 1300 kN/m2 and m ¼ 0:4, thus leading to the Lam
e constants included in Table 2. Furthermore, the remainder of the material parameters of this benchmark different from the values summarized in Table 1 can also be taken from Table 2. Based either on the fully triphasic formulation or on a pseudo triphasic formulation proceeding from the Reynolds assumption of a so-called static gas-phase, where pGR  0, Fig. 8 shows the time-dependent development of the pore-liquid saturation sL . Note in passing that the Reynolds assumption with pGR  0 obviously results in a biphasic description of the problem, since the mass balance (54) vanishes in the numerical setting. The present computations are based on a soil column with a height of h ¼ 1 m and an S L intrinsic permeability of K0S ¼ 4:5  10
13 m2 implying a Darcy pore-liquid permeability of k0S ¼ 4:5  10
6 L m/s. Finally, the initial condition for the pore-liquid saturation has been set to s ¼ 0:99. Obviously, a 5 min

30 min

2h

5h

7h

20 h

88 h

0.990 0.981 0.972 0.963 0.954 0.945 0.936 0.927 0.918 0.909 0.900

Fig. 8. Saturation of the pore-liquid during the draining process, upper series: biphasic model (pGR ¼ 0), lower series: triphasic model.

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2901

1

column height [m]

0.8

0.6

0.4 5 min 30 min 2h 5h 7h 20 h

0.2

0 0

–1400

–2800

pressure

–4200 –5600 –7000

pGR [N/m2]

Fig. 9. Time dependent development of the effective gas-pressure pGR , triphasic formulation.

distribution of the pore-gas suction can only occur in case of the triphasic formulation, cf. Fig. 9. Furthermore, it is evident that the maximal gas suction of around pGR ¼
6200 N/m2 at t ¼ 30 min significantly influences the liquid-flow situation. This effect, although it vanishes with respect to the proceeding time due to the gas inflow through the specimen’s top surface, exhibits the major difference between the triphasic and the biphasic formulations. The advantage of the fully triphasic formulation can furthermore be seen from Fig. 10, where the water efflux is plotted versus time. In particular, a comparison of the experimental data by Liakopoulos with the bi- and triphasic formulations reveals that the biphasic computation with the static gas phase overestimates the water efflux at the beginning of the draining process, whereas the fully triphasic computation yields a slight underestimation of the primary outflow. Furthermore, in contrast to the experimental evidence, the biphasic formulation predicts the final distribution of the liquid saturation at

filter velocity [10–6 m/s]

5 experiment biphasic computation triphasic computation

4 3 2 1 0 0

2.5

5 time [h]

7.5

10

Fig. 10. Liquid filter velocity nL wL versus time of the Liakopoulos experiment and of the bi- and triphasic computations.

2902

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

t ¼ 5 h, while the triphasic formulation fits the experiments very well. In this formulation, the final saturation distribution is reached after approximately 88 h, cf. Fig. 8. 4.2. The embankment problem The second numerical example concerns the 2-dimensional description of the liquid flow through an embankment, cf. Fig. 11. In particular, the embankment has a height of 10 m and a slope gradient on both sides of s ¼ 1=3. Furthermore, the ground level is assumed to be more or less impermeable governed by an S L intrinsic permeability of K0S ¼ 10
15 m2 implying a liquid Darcy permeability of k0S ¼ 10
8 m/s. The S
12 2 L
5 embankment itself proceeds from the values of K0S ¼ 10 m and k0S ¼ 10 m/s, whereas the filter at the S L right-hand side of the structure is assumed to be governed by K0S ¼ 10
9 m2 and k0S ¼ 10
2 m/s. In order to exhibit the differences between typical embankment constructions, the embankment problem is discussed twicely, namely, with and without a central seal unit governed by the same permeability properties like the ground level. Based on the triphasic formulation, it is assumed that a water table of 2 m at the left side (water side) of the embankment is rapidly increased towards 8 m, cf. Fig. 11. As a result, water streams into the embankment and finally reaches a stationary state. In particular, Fig. 12 shows the pore-fluid flow at different times prior to the stationary state. Whereas the upper image exhibits the pore-liquid saturation

8m

S =10 –12 m 2 K 0S

S =10 –15 m 2 K 0S S =10 –9 m 2 K 0S

2m

S =10 –15 2 K 0S m

Fig. 11. Embankment with a central seal unit.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 12. Saturation of the pore-liquid and vector arrows of the seepage velocity of the pore-gas, embankment with a central seal unit at t ¼ 3 h (top) and at t ¼ 5 d (bottom).

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2903

together with the vector arrows of the outstreaming pore-gas three hours after the water level has been increased, the lower one illustrates the evolution of the water front five days later, when the water reaches the seal unit and rises along the line of different permeabilities at the barrier between the embankment and the seal unit materials. In the stationary state, the pore-liquid saturation reaches its final distribution. In particular, Fig. 13 (top) represents the stationary state of the embankment without a seal unit reached 57 d after the water table increase, whereas, at the same time, the embankment with a seal unit is far away from the final saturation distribution, cf. Fig. 13 (bottom) and Fig. 14. To demonstrate the effect of the space-adaptive computation, it is seen from Fig. 15 that the mesh is strongly refined in the zones with high gradients of the effective saturation. To obtain this result, where the element size only depends on the distribution of the saturation,

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 13. Pore-liquid saturation 57 d after the increase of the water table: embankment without a seal unit in the stationary state (top), embankment with a central seal unit prior to the stationary state (bottom).

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 14. Pore-liquid saturation of the embankment with a central seal unit in the stationary state reached 6.3 years after the increase of the water table.

Fig. 15. Final mesh in the stationary state of the embankment without a seal unit.

2904

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

the computations have been carried out with an error indicator that was only governed by the saturation input. Finally, the above representations realize that the seal unit prevents the right part of the embankment (air side) from an increasing water table. As a result, the stability of the air side of the embankment is significantly improved. The investigation of embankment stability problems is a crucial issue and has to be investigated very carefully. In particular, instabilities are usually initiated by the localization of plastic deformations in narrow bands, the so-called shear bands, driven by gravitational forces and local saturation values. Furthermore, the onset of localizations can appear on both sides of the embankment, the water side as well as the air side. However, both a moderate embankment slope angle and the existence of a filter unit at the air side help to prevent the embankment construction from stability problems. Concerning the construction and the geometry of the examples computed above, no shear band development could be found. However, there are geometries and technical layouts of embankments which are very sensitive to stability phenomena. To support this statement, two further boundary-value problems have been computed for an embankment without a central seal unit. The first of the following examples concerns a sudden decrease of the water table at the water side from a top level of 8 m to zero, whereas the second example describes an embankment without a filter unit such that the water with a level of 8 m at the water side can leak from the slope at the air side of the embankment. Concerning the first of these examples, a slope angle on both sides of the embankment of 30 is assumed. Furthermore, the initial water table is at a height of 8 m until the stationary state of the saturation distribution is reached. Then, the water table is suddenly lowered towards the ground level, thus causing the water to leak from the water side of the embankment. During the course of the decreasing saturation in the embankment construction, cf. Fig. 16 (top), the soil beneath the water table is still under the effect of buoyancy, while the excess liquid pressure at the slope is already zero. As a result, plastic strains begin to localize starting at the kink between the slope and the ground level and forming a nearly circular shear band, cf. Fig. 16 (bottom). In the second example, the embankment is again considered without a filter unit. It is furthermore assumed that the slope gradient at the water side keeps its original value of s ¼ 1=3, whereas, at the air side, there is a stronger inclination with an angle of 35. As a result of the water table at the left side of the embankment in combination with the missing filter unit, water is leaking from the embankment at the air

Fig. 16. Saturation of the pore-liquid (top) and localization of the accumulated plastic strains [10
1 ] (bottom).

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2905

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.4

0.8

1.2

1.6

2.0

2.4

2.8

3.2

3.6

4.0

Fig. 17. Saturation of the pore-liquid (top) and localization of the accumulated plastic strains [10
1 ] (bottom).

side through the slope, cf. Fig. 17. As a consequence, one observes the onset of a shear band at the air side, starting again at the kink between the slope and the ground level. Note in passing that this situation is very typical to appear in case of natural hazards initiated, e.g., by extremely heavy rainfall and comparable situations. In contrast to the previous example, where the saturation of the embankment is decreasing, thus yielding a stabilization of the localized structure, the present example needs to be prevented from destruction. For example, this can be done by loading the air side of the embankment by additional weights given, e.g., by sandbags. Finally, it should be mentioned with respect to the well-known ill-posedness of the localization problem that a non-regularized space-adaptive computation would reveal a shear band thickness shrinking with the element size. In order to overcome this unphysical behavior, a regularization of the problem is mandatory and has been applied to the present computations by the consideration of viscoplastic solid behavior. Further regularization strategies are possible and have been commented in the introduction to this article. 4.3. The excavation problem The final example concerns the application of the triphasic model to a 3-dimensional computation of a standard excavation problem carried out at a natural slope, cf. Fig. 18. In particular, the top figure shows the total geometrical situation after the excavation process has been realized, whereas the lower figure gives a perspective image of half of the problem. Due to symmetry conditions, the geometry of half of the total slope is sufficient for the following numerical computations. Assuming the ground-water table to be close beneath the excavation ground, Fig. 19 (left) exhibits the saturation distribution in the stationary state, whereas Fig. 19 (right) shows the distribution of the accumulated plastic strains as a result of the rearrangement of stresses driven by gravitational forces due to the excavation process. For convenience, these figures are given both with and without the underlying FE mesh. In addition to Fig. 19, Fig. 20 represents a situation, where the water table has been raised up to the top level in the back of the slope. Then, the computations have been carried out until the stationary state was reached. In the stationary state, Fig. 20 (left) exhibits the saturation distribution of the pore-water, whereas Fig. 20 (right) shows the distribution of the accumulated plastic strains. As before, both figures are given with and without the underlying FE mesh. Furthermore, it is again seen that a leaking of the pore-water at

2906

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

18.0 m 5.0 m 24.0 m

30.0 m

14.0 m 7.0 m

cutting planes

21.0 m

12.0 m

Fig. 18. Sketch of the whole slope (top) and perspective view on half of the excavated slope (bottom).

Fig. 19. Ground-water table (left) and accumulated plastic strains (right) after the excavation process (scaled 2.5 times), upper series without and lower series with the underlying FE mesh.

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2907

Fig. 20. Ground-water table (left) and accumulated plastic strains (right) after an increase of the water table (scaled 2.5 times), upper series without and lower series with the underlying FE mesh.

0.02 0.04 0.05 0.07 0.09 0.10 0.12 0.14 0.15 0.17 0.18 0.20 0.22 0.23 0.25

Fig. 21. Accumulated plastic strains in the cutting planes after an increase of the water table (scaled 2.5 times).

the slope is sufficient for the onset of a shearing domain with a certain thickness. Separating the upcoming shearing domain into individual shear bands by cutting off slices from the slope at distinct positions reveals that the localizing domain is obtained in a shape comparable to a double cranked shell, cf. Fig. 21.

2908

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

However, due to the huge amount of approximately 110 700 degrees of freedom, it was found inconvenient to furthermore refine the mesh size in order to obtain a shearing domain with a lower thickness. Whereas the present computations have been carried out on a single processor system, comparable studies have also been analyzed on multi-processor systems, however, restricted to a biphasic analysis, cf. the work by Wieners et al. [41,42].

5. Conclusion In the present contribution, a triphasic model for the description of partially saturated soil has been considered on the basis of the Theory of Porous Media. In particular, a thermodynamically consistent framework of an elasto-plastic or an elasto-viscoplastic cohesive-frictional soil skeleton saturated by the pore-fluids water and air has been presented. In order to describe the volumetrically fully coupled triphasic material with individual states of motion driven by external loads, pressures, displacements and saturations in combination with mutual interactions between the individual constituents, it has been assumed that the solid skeleton and the pore-water are materially incompressible, whereas the pore-gas follows the ideal gas law. Furthermore, the interaction between the solid and the fluid materials has been based on the capillarypressure-saturation relation together with the relative-permeabilities-saturation relation by van Genuchten [26] and constitutive equations for the momentum supply terms leading to Darcy-like filter laws. Based on an elasto-plastic or an elasto-viscoplastic deformation behavior of the cohesive-frictional skeleton materials, a realistic model of partially saturated soil was obtained. On the numerical side, the triphasic model was implemented in the FE tool PANDAS such that time- and space-adaptive computations of a variety of initial boundary-value problems could be carried out. In particular, the 1-dimensional computation of a leaking soil column exhibited the difference of the fully triphasic approach in comparison to the biphasic analysis based on Reynolds’ assumption of a static gasphase. This example clearly pointed out that the gas phase in a partially saturated soil is not generally negligible. The following numerical examples additionally demonstrated the wide range of useful applications of the coupled solid–fluid behavior in the framework of the triphasic formulation. In particular, it could be shown how close the fluid saturation is coupled to the solid deformation up to the onset and the development of localization phenomena prior to a possible destruction of embankments and slopes.

References [1] C. Truesdell, R.A. Toupin, The classical field theories, in: S. Fl€ ugge (Ed.), Handbuch der Physik, vol. III/1, Springer-Verlag, Berlin, 1960, pp. 226–902. [2] R.M. Bowen, Incompressible porous media models by use of the theory of mixtures, Int. J. Engrg. Sci. 18 (1980) 1129–1148. [3] R. de Boer, Theory of Porous Media, Springer-Verlag, Berlin, 2000. [4] W. Ehlers, Constitutive equations for granular materials in geomechanical context, in: K. Hutter (Ed.), Continuum Mechanics in Environmental Sciences and Geophysics, CISM Courses and Lectures No. 337, Springer-Verlag, Wien, 1993, pp. 313–402. [5] W. Ehlers, Foundations of multiphasic and porous materials, in: W. Ehlers, J. Bluhm (Eds.), Porous Media: Theory, Experiments and Numerical Applications, Springer-Verlag, Berlin, 2002, pp. 3–86. [6] R. de Borst, Simulation of strain localization: A reappraisal of the Cosserat continuum, Engrg. Comput. 8 (1991) 317–332. [7] P. Steinmann, A micropolar theory of finite deformation and finite rotation multiplicative elasto-plasticity, Int. J. Solids Struct. 31 (1994) 1063–1084. [8] W. Ehlers, W. Volk, On theoretical and numerical methods in the theory of porous media based on polar and non-polar elastoplastic solid materials, Int. J. Solids Struct. 35 (1998) 4597–4617. [9] W. Volk, Untersuchung des Lokalisierungsverhaltens mikropolarer por€ oser Medien mit Hilfe der Cosserat-Theorie, Dissertation, Bericht Nr. II-2 aus dem Institut f€ ur Mechanik (Bauwesen), Universit€at Stuttgart, 1999.

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

2909

[10] A. Needleman, Material rate dependence and mesh sensitivity in localization problems, Comput. Methods Appl. Mech. Engrg. 67 (1988) 69–85. [11] S. Diebels, P. Ellsiepen, W. Ehlers, Error-controlled Runge-Kutta time integration of a viscoplastic hybrid two-phase model, Tech. Mech. 19 (1999) 19–27. [12] P. Ellsiepen, Zeit- und ortsadaptive Verfahren angewandt auf Mehrphasenprobleme por€ oser Medien, Dissertation, Bericht Nr. II3 aus dem Institut f€ ur Mechanik (Bauwesen), Universit€at Stuttgart, 1999. [13] H.-B. M€ uhlhaus, E.C. Aifantis, A variational principle for gradient plasticity, Int. J. Solids Struct. 28 (1991) 845–857. [14] R. de Borst, L.J. Sluys, H.-B. M€ uhlhaus, J. Pamin, Fundamental issues in finite element analysis of localization of deformation, Engrg. Comput. 10 (1993) 99–121. [15] W. Ehlers, A single-surface yield function for geomaterials, Arch. Appl. Mech. 65 (1995) 246–259. [16] W. Ehlers, P. Ellsiepen, PANDAS: Ein FE-System zur Simulation von Sonderproblemen der Bodenmechanik, in: P. Wriggers, U. Meißner, E. Stein, W. Wunderlich (Eds.), Finite Elemente in der Baupraxis: Modellierung, Berechnung und Konstruktion, Ernst & Sohn, Berlin, 1998, pp. 391–400. [17] W. Ehlers, P. Ellsiepen, M. Ammann, Time- and space-adaptive methods applied to localization phenomena in empty and saturated micropolar and standard porous materials, Int. J. Numer. Methods Engrg. 52 (2001) 503–526. [18] W. Ehlers, M. Ammann, S. Diebels, h-Adaptive FE methods applied to single- and multiphase problems, Int. J. Numer. Methods Engrg. 54 (2002) 219–239. [19] R.M. Bowen, Theory of mixtures, in: A.C. Eringen (Ed.), Continuum Physics, vol. III, Academic Press, New York, 1976, pp. 1– 127. [20] C. Truesdell, Thermodynamics of diffusion, in: C. Truesdell (Ed.), Rational Thermodynamics, second ed., Springer-Verlag, New York, 1984, pp. 219–236. [21] W. Ehlers, P. Ellsiepen, P. Blome, D. Mahnkopf, B. Markert, Theoretische und numerische Studien zur L€ osung von Rand- und Anfangswertproblemen in der Theorie Por€ oser Medien, Abschlußbericht zum DFG-Forschungsvorhaben Eh 107/6-2, Bericht aus dem Institut f€ ur Mechanik (Bauwesen), Nr. 99-II-1, Universit€at Stuttgart, 1999. [22] A.W. Bishop, The effective stress principle, Teknisk Ukeblad 39 (1959) 859–863. [23] A.W. Skempton, Significance of Terzaghi’s concept of effective stress (Terzaghi’s discovery of effective stress), in: L. Bjerrum, A. Casagrande, R.B. Peck, A.W. Skempton (Eds.), From Theory to Practice in Soil Mechanics, Wiley, New York, 1960, pp. 42–53. [24] G. Eipper, Theorie und Numerik finiter elastischer Deformationen in fluidges€attigten por€ osen Medien, Dissertation, Bericht Nr. II-1 aus dem Institut f€ ur Mechanik (Bauwesen), Universit€at Stuttgart, 1998. [25] R.N. Brooks, A.T. Corey, Properties of porous media affecting fluid flow, ASCE: J. Irrigat. Drain. Div. 92 (1966) 61–68. [26] M.T. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Amer. J. 44 (1980) 892–898. [27] S. Finsterle, Inverse Modellierung zur Bestimmung hydrogeologischer Parameter eines Zweiphasensystems, Dissertation, Technischer Bericht der Versuchsanstalt f€ ur Wasserbau, Hydrologie und Glaziologie der ETH Z€ urich, 1993. [28] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface, Springer-Verlag, Berlin, 1997. [29] H. M€ ullersch€ on, Spannungs-Verzerrungsverhalten granularer Materialien am Beispiel von Berliner Sand, Dissertation, Bericht Nr. II-6 aus dem Institut f€ ur Mechanik (Bauwesen), Universit€at Stuttgart, 2000. [30] P.V. Lade, J.M. Duncan, Cubical triaxial tests on cohesionless soil, ASCE: J. Soil Mech. Found. Div. 99 (1973) 793–812. [31] D. Mahnkopf, Lokalisierung fluidges€attigter por€ oser Festk€ orper bei finiten elastoplastischen Deformationen, Dissertation, Bericht Nr. II-5 aus dem Institut f€ ur Mechanik (Bauwesen), Universit€at Stuttgart, 2000. [32] P. Perzyna, Fundamental problems in viscoplasticity, Adv. Appl. Mech. 9 (1966) 243–377. [33] O.C. Zienkiewicz, J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Numer. Methods Engrg. 24 (1987) 337–357. [34] R.W. Lewis, B.A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, second ed., Wiley, Chichester, 1998. [35] W. Ehlers, P. Ellsiepen, Theoretical and numerical methods in environmental continuum mechanics based on the theory of porous media, in: B.A. Schrefler (Ed.), Environmental Geomechanics, CISM Courses and Lectures No. 417, Springer-Verlag, Wien, 2001, pp. 1–81. [36] W. Ehlers, B. Markert, O. Klar, Biphasic description of viscoelastic foams using an extended Ogden-type formulation, in: W. Ehlers, J. Bluhm (Eds.), Porous Media: Theory, Experiments and Numerical Applications, Springer-Verlag, Berlin, 2002, pp. 275– 294. [37] W. Ehlers, B. Markert, A. Acart€ urk, A continuum approach for 3-d finite viscoelastic swelling of charged tissues and gels, in: H.A. Mang, F.G. Rammerstorfer, J. Eberhardsteiner (Eds.), Proceedings of the Fifth World Congress on Computational Mechanics (WCCM V), TU Wien, ISBN 3-9501554-0-6, 2002, http://wccm.tuwien.ac.at, Paper-ID 80236. [38] A. Liakopoulos, Transient flow through unsaturated porous media, PhD thesis, University of California at Berkeley, 1964. [39] G. Klubertanz, Zur hydromechanischen Kopplung in dreiphasigen por€ osen Medien: Modellbildung und Anwendung auf die
Ausl€ osung von Murg€angen, Dissertation Nr. 2027, Ecole Polytechnique F
ed
erale de Lausanne, 1999.

2910

W. Ehlers et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2885–2910

[40] G. Klubertanz, L. Laloui, L. Vulliet, Numerical modelling of the hydro-mechanical behaviour of unsaturated porous media, in: A. Creechan, M. Gillot, T. Kenny, D.R. Owen, E. Ramm, J. Wood (Eds.), Proceedings of NAFEMS World Congress on Design, Simulation and Optimisation, Reliability and Applicability of Computational Methods, NAFEMS Ltd, Glasgow, 1997, pp. 1302– 1313. [41] C. Wieners, M. Ammann, S. Diebels, W. Ehlers, Parallel 3-d simulations for porous media models in soil mechanics, Comput. Mech. 29 (2002) 73–87. [42] C. Wieners, M. Ammann, W. Ehlers, Distributed point objects: a new concept for parallel finite elements applied to a geomechanical problem, Future Gener. Comput. Syst., accepted for publication.