Finite element methods for the analysis of strong discontinuities in coupled poro-plastic media

Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400 www.elsevier.com/locate/cma

Finite element methods for the analysis of strong discontinuities in coupled poro-plastic media C. Callari a

a,1

, F. Armero

b,*,2

Dipartimento di Ingegneria Civile, Universit
a di Roma Tor Vergata, Via del Politecnico 1, 00133 Roma, Italy b Department of Civil and Environmental Engineering, University of California, Berkeley, CA 94720, USA Received 1 October 2001; received in revised form 15 May 2002

Abstract This paper presents the formulation of finite element methods for the numerical resolution of strong discontinuities in poro-plastic solids. Fully coupled infinitesimal conditions are considered. These solutions are characterized by a discontinuous displacement field, with the associated singular strains, and a singular distribution of the fluid content. Here, singular distributions refer to Dirac delta functions. The singular component of the fluid content distribution models the fluid accumulated per unit area of the discontinuity surface, and it is directly related with the dilatancy characterizing singular inelastic strains localized along such a surface. It further accounts for a discontinuous fluid flow vector, given by Darcy’s law in terms of a continuous pore pressure field. All these considerations are incorporated in the proposed finite element methods through a local enhancement of the finite element interpolations as these discontinuities appear. The local character of these interpolations lead after the static condensation of the enhanced fields to a large-scale problem exhibiting the same structure as common finite element models of the global poro-plastic problem, but incorporating now crucially the localized dissipative effects characteristic of the localized failures. Several numerical simulations are presented to evaluate the performance of the proposed methods. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Porous media; Elastoplasticity; Strain localization; Strong discontinuity; Enhanced finite element methods; Coupled problems

1. Introduction The deformation of solids is often characterized by the concentration of the strains in narrow bands, usually referred to as shear bands due to the nature of the dominant strains. The propagation of these bands can lead to the failure of the solid or structure. Therefore, their understanding and physical modeling

*

Corresponding author. E-mail addresses: [email protected] (C. Callari), [email protected] (F. Armero). 1 Supported by University ‘‘Tor Vergata’’: ‘‘Progetto giovani ricercatori’’ 2000, 2001. 2 Supported by ONR nos. N00014-96-1-0818, N00014-00-1-0306 and NSF no. CMS-9703000 with UC Berkeley.

0045-7825/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 3 7 4 - 2

4372

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

are of significant importance in practical applications. Looking in particular at civil engineering problems, the study of strain localization has many important applications in the analysis of structural and geotechnical engineering systems. Characteristic examples include the collapse of foundations, the stability of excavations, slopes and tunnels [10,31]. The problem of interest in this paper involves the case of a porous solid. In this case, the localized failure of the system is strongly influenced by the coupling between the mechanical deformation and the interstitial fluid flow. As shown by experimental and numerical studies, the localization of dilatancy along narrow bands results in a decrease of pore pressures and consequent fluid flow toward these shear bands [14,17,20,21,24–26,34,35,37]. The correct and efficient modeling of these effects is then of crucial importance. The approach considered herein is based on the so-called strong discontinuities. This approach has been developed for purely mechanical problems in [4,5,9,11,22,23,29], among others. The consideration of coupled thermomechanical effects can be found in [6]. As presented in [1], this approach is based on a local modeling of the localized dissipative effects observed in the response of the material through the inclusion of a surface of discontinuity of the displacement with the corresponding singular strains locally at the material points. A cohesive-frictional law along the surface of discontinuity accounts for the required localized character of the dissipation. The local character of these considerations, referred to happen in the small scale, allows to efficiently introduce these dissipative effects in a classical continuum problem (referred to as the large-scale problem). This efficiency is translated directly to the finite element methods developed to solve the problems under consideration, involving then a local enhancement of the finite elements by these singular strain modes through a number of enhanced strain parameters. These parameters are condensed out leading to a large-scale problem possessing the same structure and degrees of freedom in classical solutions in solid mechanics, namely nodal displacements in mechanical uncoupled problems. We have presented in [3] an analysis of strong discontinuities in a coupled poro-plastic solid. In particular, we observed the consistency between a continuum poro-plastic model and a localized cohesivefrictional law, relating the traction with the displacement jump and the fluid content along the discontinuity. Besides the characteristic fluid content per unit volume, a singular part of the fluid content (that is, per unit area of the discontinuity surface) appears as the direct consequence of the dilatancy in the poroplastic model. Notably, these considerations involve a continuous pore pressure field with discontinuous gradients reflecting the localized fluid flow produced by the discontinuity surface. The simple one-dimensional context of a poro-plastic dilatant shear layer was considered in this reference for an illustration of these effects in the numerical solution of the problem. In contrast, we consider the much more general case of plane problems in the current paper. The goal of this paper is then the development of enhanced strain finite element methods that incorporate for general plane problems the singular strain associated to a discontinuous displacement field, the singular part of the fluid content and the general localized poro-plastic cohesive/frictional relating these fields with the traction and the pore pressure along the discontinuity. In this situation, the discontinuities are general curves in the domain of interest, with the need to resolve their propagation through the finite element mesh. The requirements on the finite element formulations are then much more stringent than the simple one-dimensional numerical simulations reported in [3], thus deserving this separate treatment. One of the main concerns in this setting is to accommodate the directionality of the discontinuity in the elements and the aforementioned localized poro-plastic relation along it, all in combination with the mixed finite elements required for a locking-free numerical formulation. The infinitesimal range of the deformation is considered. In particular, the finite element strains are enhanced through the Dirac delta distribution resulting from a discontinuous displacement field, with the enhanced parameters corresponding to a piecewise constant approximation of the displacement jumps. The introduction of these added unknowns requires the consideration of the local equilibrium equation along the discontinuity. The fluid flow problem is then en-

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4373

hanced through the consideration of the singular part of the fluid content, given also as a piecewise constant approximation by the local fluid mass relation in terms of the opening of the discontinuity. The local character of all these considerations, including the ones related to the fluid seepage, leads again to the static condensation of these element enhanced parameters. The final numerical implementation requires then only the solution of a system of coupled nonlinear equations reflecting the global equilibrium of the solid and the global fluid mass balance in terms of the nodal displacements and the nodal pore pressures. The proposed approach proves then to be a computationally very efficient way to incorporate the localized dissipative mechanisms characteristic of these limit solutions. Furthermore, the incorporation of the discontinuities inside the element interpolations (in contrast to being defined by element boundaries only) make the proposed methods especially suitable for the use of general unstructured meshes, as they are commonly found in the modeling of the complex practical applications of interest. More importantly, the proposed approach is able to model these observed localized effects in an objective manner, that is, avoiding the pathological mesh-size dependence observed in classical implementations of continuum poro-plastic models with strain softening. This pathological mesh size dependence is nowadays well understood in the context of uncoupled mechanical problems, being directly linked to deficiencies of the local model. More specifically, the absence of a mechanism able to model the localized nature of the energy dissipation has been identified as the cause for these deficiencies in inviscid elastoplastic models with strain softening. This shortfall is a direct consequence of the lack of a material characteristic length in classical local constitutive models. This situation has been shown to lead to mathematically ill-posed problems. Several regularization techniques have been devised to overcome this limitation, including viscous, nonlocal and high-gradient models, as well as the aforementioned strong discontinuity approach incorporating directly the limit localized dissipative mechanisms. The presence in the coupled problem of diffusion effects (translated mathematically in the presence of the Laplacian operator in the fluid flow equation) has been argued to introduce a length scale in the problem and, hence, to lead to its regularization. Representative recent analyses have been presented in [37,38]. Despite these comments, it has been observed in [6] that this type of diffusive effects are not able to regularize the underlying ill-posed problem by themselves for the analog problem of coupled thermoplasticity, that is, in the absence of a viscosity in the material model with strain softening. A linearized spectral analysis of the coupled thermomechanical problem is presented in this reference. This analysis confirms the observation going back to [24] that the high-frequency disturbances (corresponding to the short wavelength limit) in the coupled inviscid poro-plastic problem with strain softening exhibit an infinite rate of growth. Hence, we can conclude the ill-posedness of such problem in the classical sense of Hadamard [36]. Similar conclusions are reported in [7] for the time-discretized form of the problem. The consideration of the limit discontinuous solutions becomes then essential [3,13]. This situation is to be contrasted with the recent results presented in [38], where a frequency dependent length scale has been identified in the context of a spectral analysis of the linearized poro-plastic dynamic problem. The difficulty in identifying the characteristic frequency in a general nonlinear problem strikes first. In contrast, our analyses support the understanding that despite the presence of length scales associated to the diffusive part of the problem, leading in fact to the rate dependence of the overall response of the solid, the lack of regularization of the high frequency results eventually in the observed deficiencies of inviscid coupled poroplastic models with strain softening. These observations are confirmed for the general plane problem by the results presented in this paper. An outline of the rest of the paper is as follows. Section 2 presents the definition of the problem at hand. The description follows the aforementioned approach by which the localized effects associated to the strong discontinuities are introduced in a local neighborhood of the material points (the small-scale problem) in the context of the global equilibrium of the solid together with the fluid mass balance in it (the largescale problem). The complete characterization of the strong discontinuities in a poro-plastic solid is presented in Section 3. More specifically, Section 3.3 considers the model problem defined by an associated

4374

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Drucker–Prager poro-plastic solid, as it is considered in the numerical simulations presented in this paper. The enhanced finite element methods proposed in this work are developed in Section 4, including complete details of the numerical implementation. Section 5 presents several numerical simulations to evaluate the performance of these methods. Finally, Section 6 concludes with some final remarks.

2. Discontinuous solutions in porous media We present in this section the equations governing the deformation of a porous solid saturated with a fluid, up to the constitutive relations of the material which are presented separately in Section 3 below for a poro-plastic medium. The interest herein is the characterization of the limiting states of the deformation of such a solid at failure. To this purpose, we consider the approach presented in [1,2] for the uncoupled mechanical problem by which the localized dissipative mechanisms appearing in the deformation of the solid are incorporated in the large scale problem (understood as defined by the classical governing equations under the usual regularity conditions) via the use of the so-called strong discontinuities. In this context, we summarize in Section 2.1 the equations describing this large-scale problem in the infinitesimal range, leaving for Section 2.2 the characterization of the strong discontinuities in the small scale of a porous media. Section 2.3 concludes by connecting these two problems. 2.1. The large-scale mechanical and fluid flow problems Under the assumption of infinitesimal strains, the reference and the current configurations of the saturated porous solid are identified with a domain X
Rndim , for the space dimension ndim ¼ 1, 2 or 3 considered in the analysis. In the large scale, the solid deformation is characterized by the displacement field of the solid’s skeleton u : X  ½0; T  ! Rndim for the time interval T of interest. The mechanical response of a local solid is then determined by the infinitesimal strains T

e ¼ rs u ¼ 12½ru þ ðruÞ 

ð1Þ

with rs ðÞ being the symmetric part of the gradient operator rðÞ and ðÞT the transpose of a tensor. The large-scale displacement field u is assumed to satisfy the standard regularity conditions, that is, the first derivatives in (1) are regular distributions in the Lebesgue measure defined in X
Rndim (see e.g., [32]). Discontinuous fields not fitting this regularity assumption are introduced in the Section 2.2 below to characterize the failure of the material in the small scale. Denoting the (total) stress field by r 2 Symðndim Þ, with Symðndim Þ being the space of symmetric tensors in Rndim , the quasi-static equilibrium of the solid can be expressed by the principle of virtual work as Z Z Z t  g dA r : rs g dX ¼ f  g dX þ ð2Þ X

X

ot X

for all admissible variations g 2 Vu where Vu ¼ fg : X ! Rndim : g ¼ 0 on ou Xg

ð3Þ

for the part of the boundary ou X
oX with imposed displacements u ¼ u. The admissible variations g 2 Vu satisfy also the aforementioned regularity conditions. In Eq. (2), we introduced the volume dX and area dA elements (i.e., the Lebesgue measures in the domain X
Rndim and its boundary oX
Rndim 1 ), the volumetric body force f and the imposed tractions t on ot X
oX. The usual assumptions ou X \ ot X ¼ ; and ou X [ ot X ¼ oX are considered for a well-posed problem. We use the symbol ‘‘:’’ to denote the standard double contraction between symmetric rank two tensors and ‘‘’’ to denote the standard Euclidean inner

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4375

product between vectors in Rndim . The total stress field r in the weak Eq. (2) is to be determined from the large-scale strain field e in (1) and other small-scale effects described in Section 2.2. Similarly, the fluid flow in the saturated porous space of the porous solid is characterized by the fluid mass velocity (relative to solid’s skeleton) qw : X  ½0; T  ! Rndim . The same regularity conditions are assumed in this case, thus describing the fluid flow through the solid in the large scale and leaving again the localized effects of the small scale for further consideration in Section 2.2. In this way, and under the assumption that no volumetric distributed sources of fluid exist, we can write _ :¼ div qw M

ð4Þ

for the large-scale fluid content M : X  ½0; T  ! R, i.e., the fluid mass increment (over the initial value) per unit initial volume of the porous solid [8]. In Eq. (4), we have denoted by ‘‘ð_Þ’’ the time derivative and by ‘‘div’’ the standard divergence operator in Rndim . We focus our attention in the following developments on the case of a barotropic fluid characterized in the large scale by a pore pressure field p : X  ½0; T  ! R. The large-scale fluid flow qw is then defined by the classical Darcy’s law qw ¼ qw kðrp  qw gÞ

ð5Þ

with qw being the density of the fluid, k the permeability tensor of the porous solid and g the gravity acceleration vector. We note that tensile stresses r and compressive pore pressures p are considered as positive in the convention taken in this paper. Finally, the weak form of mass-balance equation (4) takes the form Z Z Z _ w dX ¼  M qw  rw dX  ð6Þ qwn w dA X

oq X

X

for all admissible variations w 2 Vp of the pore pressure field, with Vp ¼ fw : X ! R : w ¼ 0 on op Xg;

ð7Þ

where op X
oX is the part of the boundary with an imposed pore pressure (i.e., p ¼ p on op X) and oq X
oX is the part of the boundary with imposed fluid flow (i.e., qw  noX ¼ qwn on oq X for the outward unit normal noX to the domain boundary oX). The standard conditions op X \ oq X ¼ ; and op X [ oq X ¼ oX are assumed again for a well-posed problem. 2.2. The small-scale displacement and fluid flow fields The coupled large-scale problem presented in the previous section considered only solutions satisfying standard regularity conditions. In this way, the nonsmoothness of the solutions arising from the localized dissipative mechanisms appearing in the deformation of the solid was not considered. These localized effects can be incorporated efficiently in the previous large-scale problem by the local characterization of the singular fields arising from discontinuous solutions with the corresponding singular derivatives (singular in the sense of distributions), as we present in this section. In this way, we consider for a given material point x 2 X a local neighborhood Xx
X; see Fig. 1 for an illustration. Locally in this neighborhood the small-scale displacement field is denoted by ul and it is different from the displacement observed in the large scale X. We write ul ¼ u þ fWuCx

in Xx ;

ð8Þ

that is, we parameterize the difference in terms of the local fields f : Xx  ½0; T  ! Rndim for a given fixed function WuCx : Xx ! R modeling the local distribution of the displacements in Xx . The interest herein lies in the consideration of solutions exhibiting a discontinuity across a surface Cx with unit normal n at x. To this purpose, we consider

4376

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Fig. 1. Kinematics of strong discontinuities in plane strain conditions.

WuCx ¼ HCx  NCux

ð9Þ

in Xx

for the Heaviside step function HCx across Cx and a general smooth function NCux , so the function WuCx exhibits a unit jump sWuCx t ¼ 1 across the surface Cx . With this parametrization, the local field f corresponds to the displacement jump across Cx . The infinitesimal strains corresponding to the local displacement field ul are given in the distributional sense by s

s

el :¼ rs ul ¼ eðuÞ þ WuCx rs f  ðf  rNCux Þ þ ðf  nÞ dCx |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} ¼: e regular distribution

ð10Þ

singular distribution

after noting the result rHCx ¼ ndCx for the Dirac delta distribution dCx across Cx [32], defined by the relation Z Z ðÞdCx dX ¼ ðÞ dCx ð11Þ Xx

Cx

for the Lebesgue measure dCx on Cx . We can identify the so-called unresolved strains eunr ¼ el  eðuÞ

ð12Þ

which are to be characterized locally in Xx by the constitutive response of the material, as described in Section 3. The small-scale fluid flow in the local neighborhood Xx is affected by the added kinematics described by the previous considerations, the discontinuity in particular. Hence, we consider the small-scale flow field qwl ¼ qw þ tWqCx

ð13Þ

in terms of the large-scale field qw and the local field t : Xx  ½0; T  ! Rndim for a given function WqCx , defined as before WqCx ¼ HCx  NCqx

ð14Þ

in Xx

for a general smooth function NCqx : Xx ! R. We recover then the jump of the fluid flow sqwl t ¼ t across Cx in Xx . The small-scale fluid content corresponding to the local fluid flow field qwl can be obtained again using the relation (4) in the distributional sense as _  WqC rt : 1 þ t  rNCq _ l :¼ div qw ¼ M M l x x |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} regular distribution

t  ndCx ; |fflfflfflfflffl{zfflfflfflfflffl}

singular distribution

ð15Þ

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4377

_ ¼M _ ðqw Þ is given by (4). Here, the symbol 1 denotes the second-order where the large-scale fluid content M ndim symmetric identity tensor in R . The fluid content decomposition (15) can be written more compactly as e dCx ; Ml ¼ M þ M

ð16Þ

e is the localized fluid content on Cx , representing the fluid mass per unit area of Cx . The continuum where M and the localized mass-balance equations are given by the regular and singular parts of (15) as q q _ ¼ div½q j M wl Xx nCx ¼ divqw  WCx rt : 1 þ t  rNCx

in Xx n Cx

ð17Þ

and e_ ¼ sqw t  n ¼ t  n on Cx M l

ð18Þ

respectively. For future use, we introduce also the unresolved fluid content Munr :¼ Ml  M

ð19Þ

for the large-scale fluid content M ¼ Mðqw Þ. 2.3. The connection of the large and small-scale problems As shown in [1,2], a well-posed local continuum formulation of the mechanical problem is obtained by imposing the local weak equation between the traction vector tCx on Cx and the stresses r in Xx Z Z 1 1  c  rn dX þ c  tCx dCx ¼ 0 ð20Þ measureðXx Þ Xx measureðCx Þ Cx for all variations c of the local displacement jumps f. The traction field tCx is determined by the constitutive relation governing the response of the strong discontinuity Cx in terms of the displacement jumps f, in particular. This is developed in Section 3 for a poro-plastic solid. In the limit hx :¼ measureðXx Þ=measureðCx Þ ! 0, Eq. (20) leads to the local equilibrium equation tCx ¼ rnjCx . The whole formulation is then to be understood in this local limit (that is, in the large-scale), thus avoiding the precise definition of the neighborhood Xx . We refer to [1,2] for complete details. As shown in Section 4 in the context of the finite element methods developed for the solution of the problem at hand, the relation (20) allows to solve for the local displacement jumps f in terms of the large-scale displacement u, hence arriving to a large-scale in this displacement u only. Similarly, the Eqs. (17) and (18) relate the jump in the fluid flow t with the large-scale fluid flow el qw , which is given in terms of the pore pressure p by the Eq. (5) through the regular M l and the singular M parts of the small-scale fluid content. This last quantity is defined again by the particular constitutive response of the strong discontinuity, as developed next in the following section.

3. Strong discontinuities in poro-elastoplastic solids The equations governing the quasi-static equilibrium and the flow of the fluid in a porous solid presented in the previous section need to be supplemented with the appropriate constitutive relations. We consider in this work the case of a poro-elastoplastic solid characterized in the infinitesimal range of interest by an additive decomposition of the strains and of the fluid content, in the context of the general thermodynamic framework first proposed by Biot in [8]. Section 3.1 presents a brief summary of these classical models; we refer to the complete account in [13] for further details.

4378

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

We presented recently in [3] a characterization of the strong discontinuities in this poro-plastic context. Section 3.2 includes a brief summary of these results to apply them in Section 3.3 to the particular case of an associated Drucker–Prager poro-elastoplastic model as it used in the newly developed finite element methods proposed in Section 4. 3.1. The continuum poro-elastoplastic model The coupled fluid-mechanical response of a poro-elastoplastic solid can be characterized by the additive elastoplastic decompositions of the strains and the fluid content el ¼ eel þ epl ;

Ml ¼ Mle þ Mlp :

ð21Þ

The sub-index ðÞl refers to the small-scale nature of these quantities in the context discussed in Section 2.2. In the developments that follow, we study the case of an isotropic response of the porous solid saturated with a barotropic fluid. Consistently with the infinitesimal strain assumption, a constant fluid density qw ¼ qwo is considered. In this context, the rate forms of the constitutive relations are given by r_ ¼ r_ 0  bp_ 1 with r_ 0 ¼ Csk e_ el

ð22Þ

0

for the Biot’s effective stress tensor r , and i Q h_e Ml  qwo b_eel : 1 p_ ¼ qwo

ð23Þ

for the pore pressure p. Here, we have introduced the notation of Csk for the solid’s skeleton elasticities (the so-called ‘‘drained’’ tangent), Q for the Biot’s modulus, and b for the Biot’s coefficient. For simplicity, we consider the case of isotropic hardening/softening determined by a single scalar strainlike internal variable al with conjugated stress-like variable q. Similarly, isothermal conditions are assumed. In the associated case, the evolution of the plastic internal variables epl , Mlp and al is given by the rate equations 9 > e_ pl ¼ kor f > = a_ l ¼ koq f ð24Þ > _ p ¼ kqwo op f > ; M l in terms of the yield function f ðr; p; qÞ, where the plastic consistency parameter k is determined by the Kuhn–Tucker loading/unloading conditions and the consistency conditions f 6 0;

k P 0;

kf ¼ 0;

and

kf_ ¼ 0:

ð25Þ 0

The yield condition is commonly defined in terms of effective stresses f ðr; p; qÞ ¼ fsk ðr ; qÞ, that is, characterizing the drained response of the porous solid. This choice leads to the relations or f ¼ or0 fsk

and

op f ¼ bor0 fsk : 1:

ð26Þ

In this case, Eq. (24) (third term) gives the relation _ p ¼ qwo be_p M l lv

ð27Þ

between the rates of the irreversible part of the fluid content Mlp and of the volumetric plastic strain eplv defined as eplv :¼ epl : 1: Combining (23) with (27), we obtain the relation

ð28Þ

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

_ l ¼ qwo p_ þ qwo be_lv M Q

4379

ð29Þ

for the volumetric strain elv :¼ el : 1. 3.2. Strong discontinuities in a poro-elastoplastic solid The constitutive relations of a poro-elastoplastic solid presented in the previous section are consistent with the presence of strong discontinuities, that is, with the associated singular strain and fluid content fields identified in Section 2.2. In particular, the consideration of a localized poro-plastic flow of the form k ¼ k~dCx

ð30Þ

for the plastic multiplier leads to the following considerations, summarizing briefly the results presented in [3]. We first observe that the softening law (24) (second term) must be understood in the distributional sense H1 q_ :¼ a_ ¼ k~oq f dCx

)

~ 1 dCx H1 ¼ H

ð31Þ

given the regularity of the stress-like internal variable defining locally the yield stress of the material through the yield function f ðr; p; qÞ. The combination of this result with the consistency condition f_ ¼ 0 for its regular and singular parts leads, after some algebraic manipulations, to the relation kf_ k ¼

1 ½or f : r_ þ op f p_  ; ~ HNsk |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} or0 fsk

ð32Þ

:r_ 0

where 2

Nsk ¼

ðoq f Þ Csk or f : ðm  nÞ or f : Csk or f

ð33Þ

for the magnitude of the rate of displacement jump kf_ k. The direction of the displacement jump rate m (i.e., f_ ¼ kf_ km with kmk ¼ 1) is obtained by imposing the regular character of the traction rate t_ ¼ r_ n, leading to the necessary condition   1 ep Qep C m ¼ 0 where Q ¼ n  C o f  C o f  n: ð34Þ  sk sk r sk r sk sk or f : Csk or f We can identify the appearance of the drained elastic perfectly plastic acoustic tensor Qep sk in (34). The introduction of the small-scale strain field el defined by (10) in the constitutive (29) leads to the following expression for the rate of small-scale fluid content _l ¼M _ þ qwo bWu rf_ : 1  qwo bf_  rN u þ qwo bf_  ndCx ; M Cx C |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}x |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} regular distribution

ð35Þ

singular distribution

_ ¼M _ ðuÞ is given by (29) for the large-scale volumetric strain ev ðuÞ. where the large-scale fluid content rate M Comparing Eq. (35) with (15), we obtain the following expression for the singular part of the unresolved fluid content e_ l ¼ t  n ¼ qwo bf_  n M

ð36Þ

and its regular part q q q q _ _ _ M unr ¼ WCx rt : 1 þ t  rNCx ¼ qwo bWCx rf : 1  qwo bf  rNCx

ð37Þ

4380

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

_ M e_ v;unr unr ¼ qwo b

ð38Þ

after denoting by ev;unr the regular part of volumetric unresolved strain ev;unr ¼ eunr : 1. Furthermore, the combination of the relations (36) and (37) motivates the choice NCux ¼ NCqx ) WuCx ¼ WqCx

in Xx

ð39Þ

as originally assumed in [3] and as considered in the finite element methods developed in Section 4. In summary, the developments in this section show the consistency of a classical continuum poro-plastic model with the discontinuous small-scale fields of displacements and fluid flow introduced in Section 2.2. In particular, we have a complete identification of the localized stress–displacement–fluid flow relations. A stress–displacement relation is implied by Eq. (32) on the discontinuity surface Cx . Eqs. (36) and (37) represent the relations between displacement jump rate and fluid flow jump on Cx and in Xx n Cx , respectively. To illustrate the significance of these results, a model example is considered in the next section. Remark 1. The pore pressure field p remains continuous in the current approach, in contrast with the analyses presented in [17,18,33] for the coupled case, extending the results presented in [19,25] for the undrained case where a discontinuity of the pressure field was assumed. The nature of the governing equations (involving the divergence of the fluid flow and consequently the second derivative of the pore pressure) do preclude, however, this discontinuity in the pore pressure. In addition, no ‘‘pressure shocks’’ have been observed in experimental tests, as indicated in [34]. Remark 2. We note the appearance of the drained response in the localization condition for the fluidcoupled case, as obtained in [20,24] for the classic case of ‘‘weak’’ discontinuities (i.e., discontinuous strains with continuous displacements). A different localization condition is reported in [17,18] due to the assumed discontinuity of the pore pressure field mentioned in the previous remark. Remark 3. Eqs. (32) and (36) define a rate-independent localized response on the discontinuity surface. One should note again that the final model is to be understood in the large scale. Looking at the material response in the small scales (i.e., small-scale models incorporating shear bands of finite thickness), a more general transient response could be postulated. The framework considered here allows the incorporation of these small-scale effects in the large-scale problem through a refinement of the localized constitutive relations (32) and (36) (e.g., a general rate-dependent viscoplastic relation). We do not explore further these refinements in this paper, since they do not affect its goals, namely, the development of the actual finite element techniques. Remark 4. Similarly, a main effect of the presence of the dilatant discontinuity is a drainage on the fluid flow pattern. Negative excess pore pressures are then observed; see the results reported in Section 5. This situation clearly indicates that another refinement of the models considered herein involves the consideration of a cavitation model for the fluid response in the bulk, a model that we have not considered in the present work. However, we note that such a model can again be formulated in the large-scale framework of interest here. 3.3. Model problem: associated Drucker–Prager As a model example of the considerations presented in the previous section, we examine in this section a Drucker–Prager model defined by the following yield function in the effective stresses:

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

fsk ðr0 ; qÞ ¼ ksk þ b13r0 : 1 þ

qffiffi

2 q; 3

4381

ð40Þ

where s ¼ devðr0 Þ is the deviatoric part of stress tensor and the pressure coefficient b is a material parameter. The continuum plastic evolution equation (24) read in this case   qffiffi s p 1 _ p ¼ kqwo bb: þ 3b1 ; a_ l ¼ 23k and M ð41Þ e_ l ¼ k l ks k The loading/unloading and consistency conditions (25) are added to these equations. As implied by Eq. (34), the inception of localization is determined by the solid’s skeleton (‘‘drained’’) response of the material. Similarly, the post-localization relation (32) can be expressed only in terms of effective stresses. Hence, to particularize these equations for the considered model example, we can employ the main results obtained in [11] for the uncoupled case following the work presented in [29]. Denoting by t the unit tangent to Cx , the localized dilatancy is defined under the assumption of plane strain as the ratio U :¼

f_n jf_t j

where fn :¼ f  n

and

ft :¼ f  t

ð42Þ

are the normal and the tangential components of displacement jump, respectively. See Fig. 1 for an illustration. For the considered associated Drucker–Prager model, the appearance of the localized plastic flow (30) leads to the expression for the localized dilatancy pffiffiffi 2b U ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð43Þ 2 1  2b2 =3 and to the following expression for the angle h between n and the major (in-plane) principal stress direction e1 (see Fig. 1) " # _t Þ signð f h ¼ 12 arctan : ð44Þ U Denoting by r01 , r02 the in-plane principal effective stresses (r01 P r02 ) and by s3 the out-of-plane principal deviatoric stress, the localization condition (34) is satisfied for sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 r0  r02 ks k ¼ ð45Þ with r ¼ 1 2 r 2 1  b =6 or, equivalently s3 þ 13b ¼ 0:

ð46Þ

In particular, the equivalence between conditions (45) and (46) can be shown by means of the elastic stress– strain relations holding for general plane strain elastoplasticity. The localized relation (32) between the stress and the displacement jump rates reads in this case s_ Cx signðf_t Þ þ Ur_ 0Cx ¼

~ H jf_t j 3  2b2

where sCx :¼ t0Cx  t

r0Cx :¼ t0Cx  n

ð47Þ

4382

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

are the tangential and the normal components of effective traction vector t0Cx ¼ r0 njCx on Cx . Finally, the localized fluid content, given by (36), reads in this case e_ l :¼ t  n ¼ qwo bUjf_t j: M

ð48Þ

Eq. (48) clearly shows how the appearance of a localized dilatancy on surface Cx is the origin of the singular e of fluid content. part M Remark 5. The continuum dilatancy is commonly defined as the ratio e_pv =kdevð_ep Þk. It is well known that in plane strain conditions, the following expression of the Drucker–Prager pressure coefficient b in terms of the drained friction angle u sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 sin2 u b¼ ð49Þ 3 þ sin2 u leads to the same continuum dilatancy prediction of Mohr–Coulomb model (e.g., see Ref. [12]). The substitution of expression (49) in Eqs. (43) and (44) gives p u U ¼ tan u and h ¼   ð50Þ 4 2 as expected. Note that, as a consequence of (50), the relation (47) assumes the same form of the localized response law obtained in [3] for the Mohr–Coulomb model.

4. The enhanced finite element formulation We develop in this section finite element methods in the context of the enhanced strain finite element formulation [30] to solve the localized coupled poro-plastic equations developed in the previous sections. In this way, these developments extend the results presented in [4,5] for the mechanical uncoupled case. 4.1. The finite element interpolations The large-scale displacement field u is approximated by a standard isoparametric interpolation of the general form ue ðxÞ ¼ Nue ðxÞde

ð51Þ

for the shape functions Nue ðxÞ and nodal displacements de in the typical element Xe of a regular partition Sthe elem of the domain of interest X ¼ ne¼1 Xe in nelem finite elements. The simulations presented in Section 5 consider the case given by a quadratic six-noded triangle. The multi-scale framework developed in the previous sections is implemented in the proposed finite element method by identifying the local neighborhood Xx with the element Xe;loc
X where localization has been detected. In particular, the function WuCx in (9) is approximated by WCe ðxÞ ¼ HC ðxÞ  N ðiÞ ðxÞ

where

N ðiÞ ðxÞ ¼ 1 

ðxðiÞ  xÞ  nðiÞ hðiÞ

ð52Þ

is the linear shape function associated to node ðiÞ sustaining the discontinuity line Ce , assumed piecewise straight at the different elements. The element size hðiÞ in (52) is defined in Fig. 2, depicting also the unit normal nðiÞ to the side opposite to node ðiÞ. Noting that rN ðiÞ ¼ nðiÞ =hðiÞ , the assumed enhanced strain field is then constructed as

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4383

Fig. 2. Displacement and fluid flow discontinuous interpolation function WCe .

ee ¼ Bue de Ge ze þ Pe ze dC |ffl{zffl} |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}

ð53Þ

enhanced

conforming

for the standard strain matrix Bue :¼ rNue ðxÞ associated with the conforming part of the displacements (the large-scale strain field in the notation of Section 2), and for the enhanced strain operators   1 1 s ðiÞ s ðiÞ s Ge ¼ ðiÞ ðn  n Þ ðt  n Þ ; ð54Þ Pe ¼ ½n  n ðt  nÞ : h hðiÞ The vector of the normal and tangential displacement jump components is written as   f ze ¼ ne : fte

ð55Þ

Expression (53) corresponds to a constant approximation over the element of the displacement jump vector fe ¼ ½n tze . A similar treatment is adopted for the fields related to the fluid in the porous space. The pore pressure is approximated in a typical element Xe using the isoparametric interpolation pe ðxÞ ¼ Npe ðxÞpe

ð56Þ Npe ðxÞ

for the shape functions and the nodal pore pressures pe . In the notation of Section 2, the large-scale fluid flow field is then approximated as qw ¼ qwo kðBpe pe  qwo gÞ

ð57Þ

for the conforming gradient operator Bpe :¼ rNpe ðxÞ. In a localized element Xe;loc , the small-scale flow field is approximated through the enhanced interpolation qwl ðxÞ ¼ qwo kðBpe pe  qwo gÞ þ te WCe ðxÞ |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} conforming

ð58Þ

enhanced

for the constant approximation over the element of the fluid flow jump vector te . Note that WCe is the same discontinuous interpolation function (52) adopted for the displacement field after the result (39). Remark 6. The expression (53) can be generalized to mixed assumed strain finite element methods through u the consideration of a general strain operator Be . For example, the numerical simulations presented in Section 5 consider a quadratic interpolation of the displacements based on a six-noded triangle, with a

4384

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

discontinuous linear interpolation of the volumetric strain and stress. Even though a mixed formulation may be explicitly considered in the actual numerical implementation, the forthcoming developments apply to this general case, noting the equivalent assumed strain formulation based on the assumed strain operator u Be [16]. The same considerations apply to the interpolation of the pressure field and the fluid mass-balance equation. 4.2. The finite element equations The introduction of the finite element interpolations into the weak Eq. (2) leads, after standard algebraic manipulations, to the residual equation nelem

rd ¼ f

ext



A

Z

e¼1

BuT e r dX ¼ 0

ð59Þ

Xe

for the nodal external force vector f ext , including imposed boundary displacements. A standard quadrature rule is employed to evaluate the integral in (59). The symbol A has been used in (59) to denote the assembly of all the nelem element contributions. The total stresses r are obtained through the relation (22) (first term) in terms of effective stresses r0 and pore pressure pe given by (56). The effective stresses are calculated by (22) (second term) in terms of the regular part ee :¼ Bue de  Ge ze of the enhanced strain field (53). The enhanced strain parameters fe are obtained in terms of the large-scale displacement ue by considering the local equilibrium Eq. (20). Denoting the area of the localized element Xe;loc by Ae , the particular consideration of the piecewise constant approximation of the enhanced and traction fields assumed in the present developments leads to the relation Z 1 ^renh ¼  PT r0 dX þ t0C ¼ 0 ð60Þ Ae Xe;loc e after writing (20) in terms of the effective stresses (and corresponding effective tractions t0 ) for the poroplastic model considered in Section 3.3. As a matter of fact, the direct known relation between the two components of the displacement jump (42) by the dilatancy present in the current model allows the reduction of the problem to a single enhanced strain parameter fte , with the corresponding projection of the residual Eq. (60), that is, Z 1 renh ¼  D PTe r0 dX þ s0C ¼ 0 ð61Þ Ae Xe;loc for the projection operator [11] D ¼ ½signðf_te ÞU

1

ð62Þ

with the localized dilatancy U defined by (42). Here, we have introduced the notation s0C :¼ Dt0C

so s0C ¼ sC þ signðf_te ÞUr0C :

ð63Þ

The projected traction s0C ¼ ^s0C ðfte Þ is obtained by integrating the inelastic relation (47), that is, the given localized softening law on the discontinuity Ce ^s_ 0C ðfte Þ ¼

~ H f_te : 3  2b2

A standard quadrature rule is employed again in the evaluation of the integral in (61).

ð64Þ

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4385

The coupled fluid problem is determined by the weak Eq. (6) in the large scale. The consideration of the finite element interpolations considered in the previous section leads to the finite element equations in residual form nelem

rp ¼ sext 

A ½Se p_ e þ QTe d_ e þ He pe  ¼ 0;

e¼1

ð65Þ

where we have introduced the fluid compressibility matrix Se , the permeability matrix He and the coupling matrix Qe , given elementwise by Z Z qwo pT p p Se ¼ Ne Ne dX; He ¼ qwo BpT ð66Þ e kBe dX Xe Q Xe and Qe ¼

Z

p qwo bBuT e 1Ne dX

ð67Þ

Xe

as well as the vector of imposed external flows sext , including the corresponding contributions due to the imposed boundary pore pressures. A standard quadrature rule is again used in the evaluation of the integrals in (66) and (67). Eq. (65) imposes weakly the global (large-scale) balance of fluid mass as given by the conforming contributions of the fluid flow fields, whereas the localized balance of fluid mass on the strong discontinuity is handled locally at Xe;loc element level. In particular, the appearance of additional fluid-source terms in the localized element is a consequence of the volumetric part of strain enhancement (53) associated to the displacement jump fe . These additional terms of the fluid content can then be obtained as a post-processing in terms of the enhanced parameter fe as e_ e ¼ te  n ¼ qwo b f_ e  n M |ffl{zffl}

ð68Þ

Uf_te

e_ e , and for the singular part of fluid content enhancement M 1 1 _ ðiÞ M ¼ qwo b ðiÞ f_ e  nðiÞ eunr ¼ ðiÞ te  n h h

ð69Þ

_ . These very same relations (68) and (69) allow also for the regular part of fluid content enhancement M eunr to obtain the actual enhanced fluid flow field (58) as a post-processing in Xe;loc , if desired. 4.3. Some remarks on the numerical implementation The final set of nonlinear equations, that is, Eqs. (59), (61), (65) together with the constitutive relations (the localized stress–displacement relation (64) and the bulk poro-elastoplastic model described in Section 3.1) is solved through a standard Newton–Raphson incremental/iterative procedure. In a typical time step (say, ½tn ; tnþ1  with a time step Dt ¼ tnþ1  tn ), the rate Eq. (65) is approximated through a backwardEuler scheme, with the rate terms given by ð_Þ ¼ ½ðÞnþ1  ðÞn =Dt and all other terms evaluated at tnþ1 . A standard closest-point projection return mapping algorithm is employed in the integration of the bulk elastoplastic model at each quadrature point to obtain the effective stresses r0nþ1 leading to the total stresses rnþ1 appearing in the equilibrium Eq. (59) with the added contributions of the pore pressure pnþ1 . We refer to [27] for complete details on these integration algorithms including their linearization. A similar procedure is used for the integration of the localized softening model. Table 1 summarizes the step in this procedure leading to the local update of the enhanced parameter and the calculation of the resolved stress s0Cnþ1 . At the beginning of localization, s0C0 is calculated by imposing fC ¼ 0, that is, s0C0 ¼ jsC j þ Ur0C .

4386

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Table 1 ~ =ð3  2b2 Þf_t Integration algorithm in a localized element for the softening law on the discontinuity: s_ 0C ¼ ½H e Let ften be the initial enhanced parameter at tn and deðkÞ the given displacements at iteration k in ½tn ; tnþ1 . Fix this iteration and compute fkte , for each localized element Xe;loc , by means of the following sub-iteration at the element level: (1) Initialize: fðiÞ te ¼ ften and i ¼ 1 (2) Compute the normal component of displacement jump:   fðiÞ ¼ fn þ UfðiÞ  ft  ne

te

en

en

ðiÞ

u ðkÞ ðiÞ 0 and (3) Outside the discontinuity (Xe;loc nCe ), compute strains eðiÞ e ¼ Be ðkÞde  Ge ze , corresponding effective stresses r ðiÞ drained tangent Csk (elastic response) ~ ðiÞ and (4) On the discontinuity (Ce ), compute H ðiÞ

s0C :¼ s0C signðf_te Þ ¼ s0Cn þ

~ ðiÞ  ðiÞ  H f  ft  t en 3  2b2 e

(5) Compute " # Z ðiÞ 1 ðiÞ r0C PT r0 dX; ¼ ðiÞ Ae Xe;loc e sC

HðiÞ e ¼

1 Ae

Z

ðiÞ

PTe Csk Ge dX þ

Xe;loc

~ ðiÞ H 3  2b2

and ðiÞ

ðiÞ

ðiÞ

ðiÞ

fC :¼ jsC j þ Ur0C  s0C

(6) Check sliding condition ðiÞ

ðkÞ

ðiÞ

ðkÞ ðiÞ 0 0 IF fC < TOL THEN fkte ¼ fðiÞ te , fne ¼ fne , sC ¼ sC EXIT

(7) Update enhanced parameter 1

ðiÞ

ðiÞ ¼ fðiÞ fC signðsC Þ ftðiþ1Þ te þ He e

Set i

i þ 1 and GO TO 2

The final implementation requires the linearization of the governing Eqs. (59) and (65), leading to the linear system of algebraic equations  nelem  1 ðkÞ rd  A Ke Dde  Qe Dpe  LGe Dfte ¼ 0; ð70Þ e¼1 qwo rðkÞ p

nelem 



A

e¼1

   1 T 1 Q Dde þ He þ Se Dpe ¼ 0 Dt e Dt ðkþ1Þ

ð71Þ

ðkÞ

in the increments DðÞ ¼ ðÞnþ1  ðÞnþ1 , where the index ðkÞ denotes the previous iteration in the given interval ½tn ; tnþ1  where all the values are known, with ðk þ 1Þ being the newly obtained iterative values. All ðkÞ the matrices in (70) and (71) are evaluated at the known values ðÞnþ1 (indices omitted to simplify the notation). The element stiffness matrix Ke is given by Z u BuT ð72Þ Ke ¼ e Csk Be dX Xe

with the additional contribution of the localized fields defined as Z LGe ¼ LGe DT with LGe ¼ BuT e Csk Ge dX: Xe;loc

The matrices Se , He and Qe have been defined in (66) and (67).

ð73Þ

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4387

The enhanced parameter increments Dfte are obtained by the linearization of the local Eq. (61), which reads ðkÞ

renh  LTPe Dde þ He Dfte ¼ 0 for the operators Z 1 BuT Csk Pe dX; LPe ¼ Ae Xe;loc e

ð74Þ

1 He ¼ Ae

Z Xe;loc

PTe Csk Ge dX þ

ds0C dfte

ð75Þ

with Pe ¼ Pe DT and Ge ¼ Ge DT (note that Ge Dze ¼ Ge DT Dfte ). For the considered model, the second term of the enhanced tangent operator (75) (second term) is obtained from (64) as ~ ds0C H ¼ : dfte 3  2b2

ð76Þ

The combination of the linearized equations (70) and (74) leads readily to the local elimination of the enhanced strain parameter increments Dfte through a static condensation at the element level. We adopted though the staggered static condensation strategy proposed in [28] and considered in [3,11,18] to solve for these parameters. This strategy consists of solving the nonlinear Eq. (61) (that is, driving renh to zero) for given displacements dðkÞ and pore pressures pðkÞ e e . This local Newton iteration, summarized in Table 1, is based on the return-mapping algorithm presented in [4,5] and involves the scalar Eq. (74) with Dde ¼ 0. The final statically condensed equation reads then  nelem  1 ðkÞ 1 T rd  A ðKe  LGe He LPe Þ Dde  Q Dp ¼ 0 ð77Þ e¼1 qwo e e replacing Eq. (70) in the global solution process. The lack of symmetry of the final statically condensed stiffness matrix (77) in the case of LGe 6¼ LPe , that is, when the discontinuity is not aligned with the element (i.e., n 6¼ nðiÞ in (52)), is apparent. This combination, however, has been observed in practice to lead to a sharper resolution of the nonsmooth deformation patterns when compared with symmetrized versions of this same expression. We observe then that the final formulation involves the global solution of the system of equations (71) and (77) on the increments of the nodal displacements de and the nodal pressures pe only. That is, we need to solve only the large-scale problem, in its defining variables, while the effects of the localized dissipative mechanisms are efficiently introduced locally at the element level through the static condensation. Furthermore, the local fields, including the displacement jumps and localized fluid flow and fluid content, are readily available for post-processing. Remarkably also, these localized effects corresponding to a surface of discontinuity are introduced without the need of additional considerations on the quadrature rules employed in the evaluation of the different matrices and arrays involved in the numerical implementation. The efficiency of the proposed numerical treatment, based on unregularized discontinuous fields, is to be contrasted with alternative numerical approaches based on regularized discontinuities [17– 19,23]. Remark 7. The strong discontinuity is propagated through the mesh, without the spatial discretization knowing a priori its location. The procedure described in [5] is employed, based on the acoustic tensor condition (45) (or, equivalently, (46)) and (44) for the determination of the normal in the Drucker–Prager model considered in Section 3.3. In particular, the choice of a particular orientation of the normal is obtained by identifying the propagating discontinuity as an a or a b slip line in the classical terminology of slip line theory (see e.g., [15]). This propagation is done at the convergence of each time step, with the

4388

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

discontinuity propagating for as many consecutive elements as this localization condition has been satisfied and, hence, involving only geometric arguments. This is efficiently accomplished through the consideration of the connectivity graph of the finite element mesh.

5. Representative numerical simulations The purpose of this section is to evaluate the performance of the finite element formulation proposed in Section 4. In particular, we present the results obtained for two different problems: a plane strain compression test, studied in Section 5.1, and a foundation on a vertical cut, discussed in Section 5.2. As noted in Remark 6, the enhanced finite element method developed in the previous section is implemented over the mixed triangle consisting of a quadratic interpolation of the displacements with a linear interpolation of the volumetric strain and stress, that is, the so-called P2/P1 element based on a six-noded triangle [16]. The same quadratic interpolation as well as a linear interpolation of the pore pressure field have been considered, with no qualitative differences in the different aspects of the solutions discussed below. As noted in Remark 7, the initial spatial discretizations have no information of the discontinuities to be resolved. The enhanced fields are added after propagation of the discontinuities through the use of the localization condition (45). A tolerance of 1:0  103 is considered in the evaluation of this condition. After localization, the linear localized softening law ~ al q~ðal Þ ¼ ry  H

ð78Þ

~ . We is considered in all the examples, for the initial yield stress ry and the localized softening modulus H adopt a linear isotropic law characterized by the Young modulus Esk and by the Poisson coefficient msk for the solid’s skeleton elastic response. Similarly, isotropic conditions are assumed for the permeability, that is, k ¼ k1. Of particular interest is the study of the pathological mesh size dependence of the solutions and its elimination through the incorporation of the localized dissipative mechanisms, as discussed in Section 1. To further investigate this important issue, the numerical tests described in this section are also performed by means of the standard finite element coupled formulation. We adopt again a mixed triangular finite element with a quadratic interpolation of displacement field and linear discontinuous interpolation of the volumetric strain and stress part [16]. We consider this improved element to avoid the eventual locking of standard isoparametric elements. In conditions close to the undrained limit (low permeability, high load rate) and for a quasi-incompressible fluid (high Biot modulus), such a locking is due to the quasi-isochoric response of the saturated porous solid. 5.1. ‘‘Drained’’ plane strain compression tests We consider in this section the numerical simulation of a plane strain compression test. The geometry of the sample and the boundary conditions imposed on the displacements and the pore pressure fields are illustrated in Fig. 3. In particular, the so-called ‘‘drained’’ test is considered, that is, with pervious sample basis. A vertical displacement is imposed to the top base with the constant rate of v ¼ 5:00  106 m/s. The assumed material parameters, reported in Table 2, can be realistically considered as characterizing an overconsolidated clay saturated with water. According to relation (49), the assumed b parameter corresponds to u ¼ 24°; the hydraulic permeability obtained from parameters reported in Table 2 is kh ¼ qwo jgjk ¼ 1:0  107 m/s. The test is performed with two different discretizations of the sample: 2  4  14 (Fig. 4a and b) and 2  8  28 (Fig. 4c and d) elements. The solutions obtained with the two meshes are practically coincident,

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4389

Fig. 3. Plane strain compression test. Configuration of the problem with the assumed boundary conditions.

Table 2 Material parameters considered in the numerical simulations of plane strain compression tests Drained Young modulus Drained Poisson coefficient Pressure coefficient Initial yield stress Localized softening modulus Continuum softening modulus Biot coefficient Biot modulus Permeability Fluid density

Esk msk b ry ~ H H b Q k qwo

20 000 0.25 0.56 63.41 50 000 0 1.0 3:33  106 1:0  108 1:0  103

kPa

kPa kPa/m kPa kPa m2 /(kPa s) kg/m3

both in terms of the vertical reaction (Fig. 5a) and the pore pressures (Figs. 4b,d and 5b), showing the objectivity of the proposed strong discontinuity method. The horizontal displacement distributions shown in Fig. 4a and c point out the localized elements and demonstrate the ability of the enhanced method in capturing localization. In this figure and in the other deformed meshes reported in the following, a displacement amplification factor equal to 3 is adopted to visualize better the different deformation patterns in this infinitesimal setting. We have marked the times corresponding to the first and last localization registered in the elements on the curves plotted in Fig. 5. After the last localization, i.e., when the discontinuity line is fully formed, a super-convergence is attained (e.g.: kRk k=kR0 k ¼ 1:0  100 ; 3:2  101 ; 1:6  1012 ). No imperfection is introduced in the sample to trigger the localization onset in this numerical test. On the contrary, localization is naturally triggered by the heterogeneous effective stress state induced by consolidation. To investigate the coupling effects on the calculated response, the compression test is repeated for other two different rates of imposed vertical displacement: v ¼ 7:34  109 and 1:00  106 m/s. Note that the

4390

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Fig. 4. Plane strain compression test: enhanced finite elements with strong discontinuities. Deformed 2  4  14 (a,b) and 2  8  28 (c,d) meshes with horizontal displacement (a,c) and pore pressure (b,d) distributions (top base vertical displacement: 1:7  103 m).

lowest considered rate is in the range normally adopted for ‘‘drained’’ plane strain compression tests on clay. For comparison, the same test is simulated also with the strong-discontinuity uncoupled finite element formulation implemented in [11]. Before localization, results in terms of the vertical reaction (Fig. 6a) are not significantly influenced by the different applied displacement rates, in spite of the different values of excess pore pressures induced at the sample center (Fig. 6b). On the contrary, after localization, significant coupling effects on the calculated reaction are observed for rates v P 1:00  106 m/s (Fig. 6a). We observe that for all the considered displacement rates and for the uncoupled case, the solutions obtained with the two different meshes (2  4  14 and 2  8  28) are coincident (Fig. 6a and b). To trigger localization in the uncoupled computation, a 1% yield stress imperfection is introduced in the lateral element, on the right side, distant 3.0 cm from top base (see Fig. 7). The same imperfection is also introduced for the two considered lower rates (v ¼ 7:34  109 and 1:00  106 m/s). In fact, the pre-localization excess

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4391

Fig. 5. Plane strain compression test: enhanced finite elements with strong discontinuities. Results corresponding to 2  4  14, 2  8  28 regular meshes and unstructured meshes with 110 and 272 elements. Vertical reaction (a) and pore pressure at the sample center (b) vs. vertical imposed displacement.

Fig. 6. Plane strain compression test: enhanced finite elements with strong discontinuities. Results obtained for different rates of imposed vertical displacement and corresponding to 2  4  14 and 2  8  28 meshes. Vertical reaction (a) and pore pressure at the sample center (b) vs. vertical imposed displacement.

pore pressures in these coupled analyses are very small (see Fig. 6b) and the corresponding effective stress state is not ‘‘heterogeneous enough’’ to avoid the numerical sensitivity of the localization onset. Only minor coupling effects are observed affecting the orientation of the discontinuity. In particular, the same set of localized elements characterizes the uncoupled case (Fig. 7) and the coupled simulations with v ¼ 7:34  109 and 1:00  106 m/s. Some slight differences can be pointed out when comparing it with the set of localized elements obtained for v ¼ 5:00  106 m/s (Fig. 4a and c). The same ‘‘drained’’ plane strain compression test is also simulated by means of standard finite elements with continuum strain softening (H ¼ 400 kPa) with no enhanced fields, for v ¼ 5:00  106 m/s. The deformed 2  4  14 and 2  8  28 meshes in Fig. 8 show that the standard formulation has difficulties in

4392

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Fig. 7. Plane strain compression test: enhanced finite elements with strong discontinuities in the uncoupled case. Deformed 2  4  14 (a) and 2  8  28 (b) meshes with horizontal displacement distributions (top base vertical displacement: 1:1  103 m). Note that the same set of localized elements is obtained from coupled analyses with v ¼ 7:34  109 and 1:00  106 m/s.

Fig. 8. Plane strain compression test: standard coupled finite element formulation. Deformed 2  4  14 (a) and 2  8  28 (b) meshes with pore pressure distributions (top base vertical displacement: 1:5  103 m).

capturing localization. In spite of the fluid flow coupling, a significant mesh-size dependence of the standard finite element solution is observed, both in terms of the vertical reaction (Fig. 9a) and the pore pressures (Figs. 8 and 9b). The finer the mesh the less energy dissipation can be observed globally in the test. This situation is to be contrasted with the results reported above for the enhanced solution. To study the coupling effects on the mesh-size dependence of the standard finite element solution, the test is also performed for v ¼ 1:00  106 m/s. It is observed that the displacement rate increase (from v ¼ 1:00  106 to 5:00  106 m/s), i.e., the increase of fluid flow coupling, has no reducing effect on the mesh dependence of the solution, both in terms of the vertical reaction (Fig. 10a) and the pore pressures

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4393

Fig. 9. Plane strain compression test: standard finite elements. Results corresponding to 2  4  14 and 2  8  28 meshes. Vertical reaction (a) and pore pressure at the sample center (b) vs. imposed vertical displacement.

Fig. 10. Plane strain compression test: standard finite elements. Results obtained for different rates of imposed vertical displacement and corresponding to 2  4  14 and 2  8  28 meshes. Vertical reaction (a) and pore pressure at the sample center (b) vs. vertical imposed displacement.

(Fig. 10b). Due to the convergence problems shown by the standard finite element formulation, it is very difficult to perform the considered tests after the attainment of the peak value of vertical reaction. To investigate the effects of the mesh alignment on the proposed solution method, the compression test is simulated also by means of two unstructured meshes consisting of 110 and 272 finite elements (Fig. 11). Coincident results are obtained, as shown in Fig. 5 (v ¼ 5:00  106 m/s). In terms of vertical reaction, a small difference is observed between the two couples of curves corresponding to the regular and the unstructured meshes, respectively (Fig. 5a). We refer also to [4,5,10,11] for a discussion of this issue in the uncoupled framework. In particular, the objectivity of the strong-discontinuity method is shown by comparisons between results obtained from structured and unstructured meshes in [4,5,10] and by comparisons between results obtained from regular meshes with different orientations in [11].

4394

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Fig. 11. Plane strain compression test: enhanced finite elements with strong discontinuities. Deformed unstructured meshes of 110 elements (a,b) and 272 elements (c,d) with horizontal displacement (a,c) and pore pressure (b,d) distributions (top base vertical displacement: 1:7  103 m).

5.2. Foundation on a vertical cut We analyze the plane strain problem of a foundation on a vertical cut in this section. The geometry of the considered domain and the imposed boundary conditions on the displacements and the pore pressures are depicted in Fig. 12. A vertical displacement is imposed to the foundation center with the constant rate v ¼ 5:00  107 m/s. The adopted soil parameters are reported in Table 3. According to Eq. (49), the considered b is obtained for u ¼ 12°. A stiff and free-draining behaviour is assumed for the foundation material (Young modulus 2:0  107 kPa; Poisson coefficient 0.3; permeability 1.0 m2 /(kPa s)). At the contact between foundation and soil, a perfect adhesion is imposed. We consider two different discretizations of the domain in this numerical example: a 2  15  15 elements mesh (Fig. 13a and b) and a 2  25  25 elements mesh (Fig. 13c and d). The response obtained with

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4395

Fig. 12. Foundation on a vertical cut. Configuration of the problem with the assumed boundary conditions.

Table 3 Material parameters considered in the analysis of the foundation on vertical cut Drained Young modulus Drained Poisson coefficient Pressure coefficient Initial yield stress Localized softening modulus Continuum softening modulus Biot coefficient Biot modulus Permeability Fluid density

Esk msk b ry ~ H H b Q k qwo

20 000 0.25 0.30 63.41 100 0 1.0 3:33  106 1:0  108 1:0  103

kPa

kPa kPa/m kPa kPa m2 /(kPa s) kg/m3

the two meshes are practically coincident, both in terms of the vertical reaction (Fig. 14a) and the pore pressures (Fig. 13b and d), verifying again the lack of pathological mesh-size dependence of the proposed formulation. The enhanced method is also able to capture effectively the localization pattern in this more severe case, as shown by the vertical displacement distributions shown in Fig. 13a and c. This situation is to be contrasted with the significant mesh dependence exhibited by standard finite elements with continuum strain softening (H ¼ 200 kPa), as depicted in Fig. 14b (vertical reaction) and in Fig. 15b and d (pore pressures). Furthermore, Fig. 15a and c show that the standard finite element solution is also strongly dependent from the mesh orientation. In fact, the predicted position of the shear band is practically coincident with the ‘‘softer’’ p=4 direction. The proposed strong-discontinuity method is used to explore the effects of the fluid flow coupling on the foundation response. Results obtained for v ¼ 1:00  108 , 1:00  107 , 5:00  107 and 1:00  106 m/s are plotted in Fig. 16. The response obtained for the lowest considered rate is coincident with the uncoupled formulation results. Coupling effects become significant for v > 1:00  107 m/s, especially in the postlocalization phase. The pore pressures reported in Fig. 16b are calculated in point A under the foundation center (see Fig. 12). Finally, the comparison of Figs. 17 a,b and 13a shows that the discontinuity orientation is significantly influenced by the fluid flow coupling.

4396

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Fig. 13. Foundation on a vertical cut. Strong-discontinuity formulation. Deformed 2  15  15 (a,b) and 2  25  25 (c,d) meshes with vertical displacement (a,c) and pore pressure (b,d) distributions (foundation center vertical displacement: 1:54  101 m).

Fig. 14. Foundation on a vertical cut. Vertical reaction vs. imposed vertical displacement obtained by 2  15  15 and 2  25  25 meshes. Strong-discontinuity formulation (a) and standard finite element method (b) results.

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4397

Fig. 15. Foundation on a vertical cut. Standard finite element formulation. Deformed 2  15  15 (a,b) and 2  25  25 (c,d) meshes with vertical displacement (a,c) and pore pressure (b,d) distributions (foundation center vertical displacement: 9:4  102 m).

Fig. 16. Foundation on a vertical cut. Strong-discontinuity method. Results obtained for different displacement rates v of the foundation center. (a) Vertical reaction, (b) pore pressure in point A (see Fig. 12).

4398

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

Fig. 17. Foundation on a vertical cut. Strong-discontinuity formulation. Deformed 2  15  15 meshes with vertical displacement distributions. (a) Uncoupled formulation (the same set of localized elements is obtained from coupled analyses with v ¼ 1:00  108 and 1:00  107 m/s). (b) Coupled finite elements with v ¼ 1:00  106 m/s (top base vertical displacement: 1:54  101 m).

6. Concluding remarks We have presented in this paper a general procedure for the incorporation of the singular localized fields at failure in standard finite element formulations of coupled poroplastic media. In particular, the discontinuous displacement fields with the corresponding singular strain fields and the singular distribution of the fluid content are introduced at the element level as a local enhancement of the interpolation fields. The appearance of new enhanced parameters requires the introduction of the local equilibrium and the local mass balance. This local structure allows the static condensation of the enhanced parameters, leading altogether to an efficient solution of these problems. The numerical simulations presented herein indicate the lack of any pathological mesh size dependence as well as an efficient solution of the highly oriented deformation patterns, independent of the mesh alignment. This situation has been contrasted with standard finite element treatments of the problem at hand. We are currently considering several extensions of the methods presented herein, starting with the consideration of general models in the finite strain range. Furthermore, the appearance of negative pore pressure associated to the dilatant discontinuity, as observed in some of the above simulations, indicate the appropriateness to consider more complex fluid flow models in the bulk response, accounting e.g., for the cavitation of the fluid, as noted in Remark 4. The refinement outlined in Remark 3 involving a rate-dependent localized law is also of our main interest. Similarly, higher order interpolation of the localized fields (e.g., not only piecewise-constant jumps) are currently under further investigation.

Acknowledgements The finite element methods proposed in this work have been implemented in the general code FEAP, courtesy of Professor R.L. Taylor at UC Berkeley. C. Callari is grateful for the invaluable guidance and continuous support of professors F. Maceri, E. Sacco and F. Auricchio. F. Armero acknowledges the financial support of the ONR under grant numbers N00014-96-1-0818 and N00014-00-1-0306 with UC Berkeley, and the NSF under grant number CMS-9703000 with UC Berkeley. C. Callari was partially supported by a CNR (Italian National Research Council) scholarship (Short Term Mobility Program, 1999) and by the ‘‘Progetto giovani ricercatori’’ 2000, 2001 of the University of Rome ‘‘Tor Vergata’’.

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

4399

References [1] F. Armero, Large-scale modeling of localized dissipative mechanisms in a local continuum: applications to the numerical simulation of strain localization in rate-dependent inelastic solids, Mech. Cohes. Frict. Mater. 4 (1999) 101–131. [2] F. Armero, On the characterization of localized solutions in inelastic solids: an analysis of wave propagation in a softening bar, Comp. Meth. Appl. Mech. Engrg. 191 (2001) 181–213. [3] F. Armero, C. Callari, An analysis of strong discontinuities in a saturated poro-plastic solid, Int. J. Numer. Meth. Engrg. 46 (10) (1999) 1673–1698. [4] F. Armero, K. Garikipati, Recent advances in the analysis and numerical simulation of strain localization in inelastic solids, in: D.R.J. Owen, E. O~ nate, E. Hinton (Eds.), Proc. Comp. Plast. COMPLAS IV, CIMNE, Barcelona, 1995. [5] F. Armero, K. Garikipati, Analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids, Int. J. Solids Struct. 33 (1996) 2863–2885. [6] F. Armero, J.G. Park, An analysis of strain localization in a thermoplastic shear layer under thermally coupled dynamic conditions. I: Continuum thermoplastic models. II: Localized thermoplastic models, Int. J. Numer. Meth. Engrg., to appear. [7] A. Benallal, C. Comi, Properties of finite-step problem in numerical analysis of unstable saturated porous continua, in: D.R.J. Owen, E. O~ nate, E. Hinton (Eds.), Proc. Comp. Plast. COMPLAS V, CIMNE, Barcelona, 1997, pp. 1611–1616. [8] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155–164. [9] R.I. Borja, R.A. Regueiro, Strain localization of frictional materials exhibiting displacement jumps, Comp. Meth. Appl. Mech. Engrg. 190 (20–21) (2001) 2555–2580. [10] C. Callari, The application of a strong-discontinuity FEM to the analysis of strain localization induced by underground openings, in: G.N. Pande, S. Pietruszczak (Eds.), Proc. Eighth Int. Symp. Numer. Models in Geomech. NUMOG VIII, Balkema, Rotterdam, 2002, pp. 163–170. [11] C. Callari, A. Lupoi, Localization analysis in dilatant elasto-plastic solids by a strong-discontinuity method, in: F. Maceri, M. Fremond (Eds.), Lagrange Symposium 2000, Springer, Berlin, in press. Available at http://alpha2.civ.uniroma2.it/callari/Callari & Lupoi4.pdf. [12] W.F. Chen, Plasticity in Reinforced Concrete, McGraw-Hill, New York, 1982. [13] O. Coussy, Mechanics of Porous Continua, Wiley, Chichester, 1995. [14] C. Han, I. Vardoulakis, Plane-strain compression experiments on water-saturated fine-grained sand, G
eotechnique 41 (1) (1991) 49–78. [15] R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, Oxford, 1950. [16] T.J.R. Hughes, The Finite Element Method, Prentice Hall, Englewood Cliffs, NJ, 1987. [17] J. Larsson, R. Larsson, Localization analysis of a fluid-saturated elastoplastic porous medium using regularized discontinuities, Mech. Cohes. Frict. Mater. 5 (7) (2000) 565–582. [18] J. Larsson, R. Larsson, Finite-element analysis of localization of deformation and fluid pressure in an elastoplastic porous medium, Int. J. Solids Struct. 37 (48–50) (2000) 7231–7257. [19] R. Larsson, K. Runesson, S. Sture, Embedded localization band in undrained soil based on regularized strong discontinuity: theory and FE analysis, Int. J. Solids Struct. 33 (1996) 3081–3101. [20] B. Loret, J.H. Prevost, Dynamic strain localization in fluid-saturated porous media, J. Engrg. Mech. 117 (1993) 907–922. [21] M. Mokni, J. Desrues, Strain localization measurements in undrained plane strain biaxial tests on Hostun RF sand, Mech. Cohes. Frict. Mater. 4 (4) (1999) 419–441. [22] J. Mosler, G. Meschke, 3D FE analysis of cracks by means of the strong discontinuity approach, in: E. O~ nate et al. (Eds.), Europ. Congr. Comp. Methods Appl. Sci. Engrg., ECCOMAS 2000, Barcelona, 2000. [23] J. Oliver, Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals. Part 2: Numerical simulation, Int. J. Numer. Meth. Engrg. 39 (1996) 3575–3623. [24] J.R. Rice, On the stability of dilatant hardening for saturated rock masses, J. Geophys. Res. 80 (1975) 1531–1536. [25] K. Runesson, D. Peric, S. Sture, Effect of pore fluid compressibility on localization in elastic–plastic porous solids under undrained conditions, Int. J. Solids Struct. 33 (1996) 1501–1518. [26] B. Schrefler, C.E. Majorana, L. Sanavia, Shear band localization in saturated porous media, Arch. Mech. 47 (1995) 577–599. [27] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer, New York, 1998. [28] J.C. Simo, F. Armero, R.L. Taylor, Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems, Comp. Meth. Appl. Mech. Engrg. 110 (1993) 359–386. [29] J.C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comp. Mech. 12 (1993) 277–296. [30] J.C. Simo, S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Meth. Engrg. 29 (1990) 1595–1638. [31] D. Sterpi, An analysis of geotechnical problems involving strain softening effects, Int. J. Numer. Anal. Meth. Geomech. 23 (13) (1999) 1427–1454.

4400

C. Callari, F. Armero / Comput. Methods Appl. Mech. Engrg. 191 (2002) 4371–4400

[32] I. Stakgold, Green’s Functions and Boundary Value Problems, Wiley, New York, 1979. [33] P. Steinmann, A finite element formulation for strong discontinuities in fluid-saturated porous media, Mech. Cohes. Frict. Mater. 4 (2) (1999) 133–152. [34] I. Vardoulakis, Deformation of water-saturated sand: I. Uniform undrained deformation and shear banding; II. Effect of pore water flow and shear banding, Geotechnique 46 (1996) 441–472. [35] G. Viggiani, R.J. Finno, W.W. Harris, Experimental observations of strain localization in plane strain compression of a stiff clay, in: R. Chambon et al. (Eds.), Localization and Bifurcation Theory for Soils and Rocks, Balkema, Rotterdam, 1994, pp. 189–198. [36] E. Zauderer, Partial Differential Equations of Applied Mathematics, Wiley, New York, 1989. [37] H.W. Zhang, L. Sanavia, B.A. Schrefler, An internal length scale in dynamic strain localization of multiphase porous media, Mech. Cohes. Frict. Mater. 4 (5) (1999) 443–460. [38] H.W. Zhang, B.A. Schrefler, Uniqueness and localization analysis of elastic–plastic saturated porous media, Int. J. Numer. Anal. Meth. Geomech. 25 (1) (2001) 29–48.