- Email: [email protected]
Stabilized low-order finite elements for failure and localization problems in undrained soils and foundations
Computer methods in applied mechanics and engineering Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
ELSEVIER
Stabilized low-order finite elements for failure and localization problems in undrained soils and foundations M. P a s t o r " ' * , T. Li b, X. L i u b, O.C. Z i e n k i e w i c z c aCentro de Estudios y Experimentacifn de Obras Pt~blicas and ETS de lngenieros de Caminos, Madrid, Spain bUniversity of Hohai Nanjing, PR China ~Universib' College of Swansea, Swansea, UK Received 20 April 1998; revised 10 July 1998
Abstract Geomaterials in general, and soils in particular, are highly nonlinear materials presenting a very strong coupling between solid skeleton and intersticial water. In the limit of zero compressibility of water and soil grains and zero permeability (which correspond to the classical 'undrained' assumption of Soil Mechanics), the functions used to interpolate displacements and pressures must fulfill either the Babuska-Brezzi conditions or the much simpler patch test proposed by Zienkiewicz and Taylor. These requirements exclude the use of elements with equal order interpolation for pressures and displacements, for which spurious oscillations may appear. The simplest elements with continuous pressures which can be used in 2D are the quadratic triangle and quadrilateral with linear and bilinear pressures, respectively. The purpose of this paper is to present a stabilization technique allowing the use of both linear triangles for displacements and pressures (T3P3) and bilinear quadrilaterals (Q4P4). The proposed element will be applied to obtain limit loads and failure surfaces in simple boundary value problems for which analytical solutions exist. © 1999 Elsevier Science S.A. All fights reserved.
1. Introduction
Low-order finite elements, such as linear triangles and bilinear quadrilaterals in 2D, provide some interesting features such as simplicity and low computational cost. In recent years, they have been found to provide an excellent definition for shocks and discontinuities in Fluid Dynamics. An additional advantage of triangles and tetrahedra is their flexibility for use in adaptive remeshing codes [1]. The situation is far less favourable in Solid Mechanics, as these advantages are outweighed by the poor behaviour exhibited by linear triangles and quadrilaterals in bending dominated conditions, or when the material behaviour is quasi-incompressible. Concerning the first problem, Pastor and Quecedo described in [2] how the failure mechanism depended on bending behaviour. Nagtegaal [3] and Sloan [4] described in their early work the second problem, which has been known to introduce important errors in the computation of collapse loads and failure mechanisms. The means to overcome this difficulty abound. In a previous paper [5] it was shown how mixed displacement-pressure formulations can improve dramatically the accuracy of collapse load computations. In the incompressible limit, the function spaces used to interpolate both displacements and pressures must fulfill the Babuska-Brezzi conditions [6,7] or the much simpler Zienkiewicz-Taylor patch test for mixed formulations [8]. Otherwise, spurious 'checker-board' oscillations may appear in the pressure field. These conditions exclude elements with equal order of interpolation of both variables. Therefore, elements such as the T3P3 triangle (linear displacements and pressures) or the Q4P4 quadrilateral (bilinear for * Corresponding author. 0045-7825/99/$ - see front matter PII: S 0 0 4 5 - 7 8 2 5 ( 9 8 ) 0 0 3 1 6-8
© 1999 Elsevier Science S.A. All fights reserved.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
220
displacements and pressures) cannot be safely used. The simplest triangles and quadrilaterals passing the test are the T6P3 (6-noded quadratic triangle for displacements and 3-noded linear triangle for pressure) and the Q8P4 quadrilateral. In fact, the T6P3 element was used in the computations described in [5] showing an excellent behaviour. In Soil Mechanics, the problem is further complicated as there exists a strong coupling between the solid skeleton and pore fluid which has to be taken into account [9-11]. One of the most widespread formulations of the problem is the so-called u-pw formulation, in which both displacements and pore pressures are chosen as field variables to describe the problem [12,13]. When approaching zero permeability and incompressibility of soil grains and pore water, the system of discretized equations is of a similar nature to that found when using mixed formulations for Solid Mechanics problems, and again, either the Babuska-Brezzi conditions, or the patch test have to be fulfilled. Again, the simplest 2D elements which can be used are the T6Pw3 and the Q8P~4 quadrilateral. The restrictions imposed by the Babuska-Brezzi conditions on both types of problems can be circumvented if special stabilization techniques are used. Such methods have been investigated in the past in the context of Fluid Dynamics 114-17], and have been extended to both Solid Mechanics and Geomechanics problems [18-21]. The purpose of this paper is to introduce a simple stabilization technique which can be used to develop suitable elements with equal order of interpolation of pressures and displacements for the computation of collapse loads of foundations and earth structures. The method is based on the fractional step algorithm proposed by the authors [21,29] for dynamic problems, and it will be applied here to problems of failure by localization. Localization of plastic deformation in shear bands is one of the most frequently found mechanisms leading to failure in soils. At a local level, it is possible to know the instant at which a shear band can appear and its orientation. In the case of elasto-plastic materials following the Tresca or Von Mises yield criteria, it corresponds to the instant at which the plastic modulus becomes zero, i.e. at the yield surface, and the inclination in 2D problems is 45 ° relative to the principal stress axes. It should be mentioned that the problem becomes ill-posed from a mathematical point of view, resulting in non-uniqueness and mesh-dependence of the solution. For a review of the research work done in the past years in this area, the reader may consult Vardoulakis and Sulem [22] and the references given therein. This paper will not focus on non-uniqueness mesh-dependence issues, but on how to develop simple Finite Elements suitable for computation of collapse loads and failure mechanisms of soils and foundations in the undrained state in which the pore pressure Pw plays the role of the mean pressure p.
2. The fractional-step algorithm for soil dynamics The fractional step algorithm was introduced by Chorin [23] as a device to allow the use of standard time integration techniques in fluid dynamics problems. Among the several alternative formulations some forms were found later to provide the required stabilization for elements with equal order of interpolation of velocities and pressures. The discovery was first made by Schneider et al. [24], and Kawahara and Ohmiya [25], and later justified by Zienkiewicz [26]. In fluid dynamics, Chorin's projection step method is based on splitting the increment of velocity into two orthogonal parts, consisting on a solenoidal field and the gradient of a scalar function:
v=w+V~
(1)
In the case of incompressible fluid flow the governing equations are 013
Ot
- STs - V p
vTv =
0
where v is the velocity, p the pressure, and S the strain operator:
(2)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
"0
0
X
S=
0
0
0
o
0
0
o
0
0 Ox 0 Oz
Oy 0 0
221
(3) 0 0 Oy 0
0
~x-
Fractional step equations introduce an intermediate velocity v* obtained in a first step: U*
- - 1.3n
_ _ At
_ STS"
(4)
where s is the deviatoric stress tensor written in vector form. The second step is D
n+l
D: ~
--
At
--719 n + ]
with
~T
D n+ ]
=0
(5) n+l
It can be seen that v* has the two orthogonal components mentioned above, the divergence free field v and the gradient field Vp. Although in its original form the scheme consisted of an explicit step followed by an implicit step (pressure was considered at time t "+1 in the second equation), further research showed that better results could be obtained by evaluating pressure at some intermediate time t" +02 and imposing the incompressibility condition at t "+°' [17], where t "+° is defined as (1 - O)t" + Ot"+1, with 0 ~< 0 ~< 1. The soil-fluid problem will be solved here using this approach. First, the balance of linear momentum for the mixture is written as dv STar ' --Vp.. +pm b -- p,. ~ - = 0
(6)
where o" is the effective stress tensor, Pw the pore pressure, p,. the mixture's density and b the body forces. Dirichlet boundary conditions are imposed on F. and tractions are prescribed on Ft. The balance of momentum and conservation of fluid phase equations can be combined, yielding: v T v -- V T k V p w + Q* [J"' = O
(7)
where we have assumed that body forces are divergence-free and that the permeability tensor is isotropic: k~--k/. 1/Q* characterizes the combined compressibility of the fluid and solid phases, and depends on voids ratio n, and bulk moduli of grains and water, Ks and K t 1
n
Q*-Kj
1-n
+
K,.
Pressures and fluxes are assumed to be given on boundaries Fp and Fq, respectively. In the incompressible undrained limit equation 7 reduces to vTv = 0. The first step of the proposed extension for saturated soils consists in discretizing the equations in the time domain using the intermediate velocity v*: V*
P"
- - V '~
At
-
AU* P" At
_
S,To .'~
+ pmb
(8)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
222
V n+l -- O*
P~
At
AV**
- Pm
=
At
~]] "+02
(9)
- p"
where we have introduced: n+O.
n
(10)
Pw " : P w + O 2 A p w
Therefore, the increment of velocity Av" is decomposed into two parts, Av* and Av**. Eq. 7 is written as 1 Ap w T n+Ol__VTvn+O Q* At = V (kVpw) [
I
(11)
where we have introduced O"+°~ = V" + 01 Av" = (1 - 01)v" + rio
n+l
Substituting now the value of v "+°~ in (11) we arrive to 1 Ap~ Q* At
--
VT.kT, .+o,. t p~ ) - (1 - O~)VTv n
~..,v .+1
-fflv
v
(12)
If we obtain velocity at time n + 1 from (9) 1 V.+ 1 = O* -- - - AtVp~ +°~ P..
(13)
and substitute it in (12), it results in 1 AAt p . , _ V T . k( ~ .p~ + o l " ) - (1 _ Q*
01)VTv
n
n+02 ~ - 01VT ( v * - - - 1 At Pw ) P.,
(14)
from which we finally get At 0102V2 ] Apw : V T ( k V p ; ) --1-~At 0 1V 2P~n _ V T v n -- 0~VT(v * - v") 0 j V T k V - - -p,.
Q* At
(15)
The solution of the problem is therefore obtained in three steps (1) Solve Eq. (8) for v* explicitly m tl n+ 1 ~! tt (2) Obtain Pw from (15), and Pw =Pw + APw (3) Use (9) to determine v "+l once Apw and v* are known The above equations can be discretized in space applying standard Galerkin procedures, and the result is the following system of algebraic equations: ij*=On+AtM]
{
1 f~+R~+
f To
N.pwdF-
f BVo-'nd/2 }
with =
T
N . o'.n d / '
with (o'.n); = ~r/jnj, or introducing f * " =f.n + R n + f r N.pwT . d E
O*=O"+AtMfl{f*.- f Bvtr'" dS2} where
(16)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
M. = f pmN uN. T d+Q B =SN.
223
(17)
It is important to note here that, in addition to external forces f : and reactions R", f . n now includes integrals over the boundary of pressures at time n, as proposed in [28] and [29]. The second discrete equation is
/'
~ Mp + Oln + At OlO2H,
)
Aft w =~r f.+o,
- Hp+n - 01 At H , f f nw
- QTt3" + OiGT A0* - fr. +rq NTp AO n d F
(18)
where Mp :
f,- ~ N~N
e dl2
H = f VTNpk.,VNp d n H* = f VTNpVNv d$2 Q = f B TmNp dO
G =f --~1 N.TVNe d~ REMARK 1. In the discretization process of terms involving second-order derivatives with respect to space, we obtain the following boundary integrals:
1 T 0 (Ap~)dF+OIAt O~02AtfFv +rq~mNP~n -01
+Gp1---NT m P OP~W On d F
frv+G 1P"N T P ao* d r
(19)
(20)
which can be written as
f
T[-At Op+ .+o~_ A v * ] d E
If we now write Eq. (13) as
1
Av"-Av*+--At
~jp~, n+02 = 0
Pm
and project it along the normal to the boundary, we obtain At ap+ In+°2
Pm On
- A v * +Av~=O
from which we can write the contribution of the boundary terms as
fFp+G NTp Av.n d F which will be assumed to be zero following Codina et al. [28]. It follows from here that Q T AV- nn :
- G pT A ~ n.
(21)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
224
REMARK 2. The force t e r m f ~ +°~ is given by f.+o, =
ap
f~, + ]'q
c9 n + O~ k-~n Pw - dF
and reduces to zero in the incompressible-undrained limit. Finally, the third equation is
~"+J = 6" - At M~lGpff~, +°2
(22)
The computational scheme requires of course that the stresses are evaluated at the end of the interval, i.e. at time n + 1. This can be done by using Ali" = t~ "+° At and accordingly computing o" '~" +~) from the constitutive relation.
3. An asymptotic form for steady state conditions in the incompressible undrained limit In the incompressible undrained limit matrices Alp and H are zero as k w and 1/Q* are zero. Eq. (16)
- ~ M.(tS* - 15")
+ Jr
.Pw a ~
(23)
reduces to
1 A--tM,(O * - ts") = f " + R" - Kti
(24)
when the material behaviour is elastic. Here, K is the stiffness matrix and ti n the nodal displacements at time t n. In steady state conditions (6* = 6n) Eq. (23) yields
f
BTtT
'
dO-
dF=f. + R
(25)
where we have omitted superscripts referring to time step. Second (18) and third (22) steps are written now as (At O~02H*) A,pw-" = -01 At H*ff~
-
QT o" + oia f A~*
(26)
and 1
--M
At
~lr - n + 8 2 . t"o- n + l __ [ ) , ) = -- pPw
Eqs. (25) and (27) can be combined to yield in the steady state
fBT~r'dO+G,Pw-fr
NT .pwdF=f. +R
or, taking into account that
(27)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
225
Gpfiw = f NTVp~ d~2 = - f VN~N~ d~2 fiw + f r N ~ P w d F = - Q TPw + fl N .Tp w dE f BTo ,
- =L+R d~f2 -Q Tm
which for elastic materials reduces to
K.a-Q
T-
p.=f
+R
(28)
Next, we will eliminate the intermediate velocity field iS* using Eqs. (26) and (27). From the latter, 1
~ M . ( v-~÷l - t3* - t~" + 0 ") = -Gp#-.+o~ w and Ate* = AO" +
A t M . JGFff~.,+°2
(29)
Substituting this value in (26) it results
O, Oz At H* A - "
-01 At H*ff"., _ QTe,, + O,GT AO" + OIAtGTM~ lG -"+°~ppw"
from which, after rearranging, we obtain
QT6" + OIAt(H* - GTM-1GI, ,,
t,)P'-"+°2~, = 0
where we have taken into account that fe = 0 as the permeability is zero. In the case of solids, steady state conditions imply ~" = 0, and above equation reduces to (H*
--
G pTM u
J
Gp)lJw = 0
which can be combined with the incompressibility constraint QTo = 0 or Qrt~ = 0 to give
QTa + a(H* - Gi~M~'Gp)ffw = 0 where a is a suitable constant. The system can be written as
f B . ¢r' - Q.15.,-f. - R = ~ = 0
(30)
n . - = ¢Pp = 0 QTti + ce(H* - GT= , , f i t I op)pw
(31)
If this system is nonlinear, it can be solved by using the Newton-Raphson method
F0"u o4
/
o~.
[ =
(32)
with a suitable Jacobian matrix j~k> which is given by (33) where K r is a consistent tangent stiffness matrix and A is:
A
=
a(H*
-- G pT M
u
J
G.) "
(34)
226
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
It is important to note that A has provided the stabilization required to circumvent Babuska-Brezzi conditions. Otherwise, the resulting system would have adopted the form
This type of system is frequently found in constrained problems, presenting pressure oscillations wherever the number of sc variables n~ ~< n,~, i.e. the number of ~b variables, as occurs when using equal order of interpolation for both variables. Concerning the stabilization matrix A added, it consists of two parts, H * which is a discrete Laplacian T --1 operator acting on the pressure, and G p M , Gp, which requires the inversion of M,. A can be obtained more economically using patches of elements surrounding each node npatch
A = ~
OjO2 At 2
H*-
GeM"
(36)
It is interesting to note that other stabilization terms proposed included the discrete Laplacian multiplied by a scalar factor [14,19,29].
4. Applications: numerical computation of collapse loads and failure mechanisms of foundations The purpose of this section is to show how simple linear triangles can provide accurate predictions of collapse loads and failure mechanisms of foundations. We will assume that both soil grains and pore water are incompressible and permeability is zero, so the behaviour will correspond to undrained conditions. The tests we propose are simple boundary value problems for which there exist analytical solutions of collapse loads against which finite element predictions can be compared [27]. Indeed, they have been used in a previous paper [5], where the objective was to show how the use of mixed elements together with adaptive remeshing techniques improved the results obtained with classical displacement based formulations. Two new aspects have been introduced in this paper: • First, we will use elements with equal order of interpolation for both displacements and pressures • Second, we will use stabilized linear triangles for both displacements and pressures which are simpler than the quadratic triangles used in the previous paper. Of course, we will use adaptive remeshing techniques to obtain an optimal resolution of the failure surface. The interested reader is referred to [1] and [5] where the procedure is explained in some detail. In the two examples that follow the soil has been assumed to be an elastoplastic material with a Tresca yield surface and no hardening (H = 0). The problem can be characterized in terms of the following non-dimensional parameters E / Q = 2000
p = 0.49 where E and ~ are Young's modulus and Poisson's ratio, and Cu is the undrained soil cohesion. Load is applied through rigid and rough footings which have been modelled as elastic, with an elastic moduli one hundred times higher than that of the soil. The problems will be solved controlling the vertical displacement 6 of the point at which load is applied. Load-displacement curves will be presented in terms of a dimensional loads F / ( B C , ) and displacements a = (E6)/BC,, where B is the width of the footing. 4.1. Strip foundation
The first example consists of a rigid and rough strip footing resting on an undrained soil layer. The geometry and boundary conditions are given in Fig. 1. Due to symmetry, only one half of the domain will be discretized, but the results will be given for the whole domain.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
227
F
I
"I
I
~'J
I
Fig. 1. Strip footing foundation: geometry and boundary conditions.
Fig. 2. Initial mesh.
The theoretical collapse load for this problem is given by F BC, = (2 + rr) The analysis begins with the mesh given in Fig. 2. The results at a value of the non-dimensional displacement cr of 25 are given in Figs. 3-5 where the deformed mesh, displacement vectors and contours of effective plastic strain are shown. Load vs. displacement curves can be seen in Fig. 14. From these results two new meshes are obtained using adaptive techniques (Figs. 6 and 10). Results are given in Figs. 7 - 9 for the first mesh and Figs. 11-13 for the second. The process of adaptive remeshing ends here as the two last force-displacement curves coincide.
Fig. 3. Deformed mesh.
...... .
•
~-~
r[hlll i hhli I .,,,~.,,.,,,
, • ,, ~ ' , , ~ J u v
I I.II~4~--,.~/" / / / , ,
. :.
. .
: .
Fig. 4. Displacement vectors.
Fig. 5. Contours of effective plastic strain.
228
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
Fig. 6. Adaptive mesh 1.
Fig. 7. Deformed mesh.
Fig. 8. Displacement vectors.
Fig. 9. Contours of effective plastic strain.
Fig. lO. Adaptive mesh 2.
Fig. 11. Deformed mesh.
Fig. 12. Displacement vectors.
Fig. 13. Contours of effective plastic strain.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 2 1 9 - 2 3 4
229
Limiting loads by different meshes
for foundation problem 6
5.5 5 4.5
,,! 13 al
4
I
l
l
t
l
l
l
l
l
/
l
l
l
l
l
l
..........
l
l
l
l
l
l
l
initialmesh
. . . .
adapevemesh1
-
-
adaptive mesh 2
I
I
3.5 3 2.5 2 I
!
I
3
I
5
I
I
7
I
I
9
I
11
I
l
13
I
15
I
I
i
i
17
|
19
I
21
I
I
23
I
25
dE/BCu Fig. 14, Strip foundations: load vs. displacement.
It is important to notice the good resolution of the failure mechanism, with an error in the computed collapse load below 5%.
4.2. Vertical cut
The last example we will consider is that of the vertical cut shown in Fig. 15, where both the geometry and boundary conditions are shown.
F ~t
=~z =sd
C~
II ,k~,%
II ,vk.~
&
II
l
:3"
u, -- O,t, = 0 ,
J(
IOta
Fig. 15. Vertical cut problem: geometry and boundary conditions.
230
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
Fig. 17 shows the initial mesh used in the analysis, and the corresponding results can be seen in Figs. 18-20. The adaptive mesh refinement process leads to the mesh shown in Fig. 21, with results given in Figs. 22-24. We then perform a last mesh refinement to arrive at the mesh of Fig. 21, for which the results are given in Figs. 26-28. Concerning the limit load, it is found again to agree well with the analytical value given by F
BC,
-
2.0
Limiting loads by different m e s h s
for vertical cut problem 2.5
I
I
I
I
I
t
I
I
/ S; "
3 U
m
..........
'" inil~l mesh
.
~
.
-
.
.
-
m
u
h
l
~p~ve me~ 2
1.5
1
I
I
I
I
I
I
I
l
2
3
4
S
6
7
8
9
t0
dF./BCu Fig. 16. Vertical cut problem: load vs. displacement plot.
1 Fig. 17. Vertical cut problem: initial mesh.
Fig. 18. Deformed mesh.
231
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
~X.s.~.3~%,.
- -
J~ /,,.,.......;... t // I/ /
/
•
,
h Fig. 19. Displacement vectors.
Fig. 20. Contours of effective plastic strain.
i
\
i
J 1
/ Fig. 21. Adaptive mesh
Fig. 22. Deformed mesh.
1
Fig. 23. Displacement vectors.
Fig. 24. Contours of plastic strain.
232
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 2 1 9 - 2 3 4
i
l V Fig. 25, Adaptive mesh 2.
/ /
/
Fig. 26. Deformed mesh.
//./I
,,///f,,#...".
/ /,
Fig. 27. Displacement vectors.
Fig. 28. Contours of plastic strain.
5. Conclusions
The restrictions imposed by Babuska-Brezzi conditions or the Zienkiewicz and Taylor patch test for mixed formulations do not allow the use of equal order of interpolation shape functions for displacements and pressures in mixed formulations. Therefore, the simplest elements with continuous pressures which can be used are the T6P3 triangle with quadratic displacements and the Q8P4 quadrilateral. This restriction can be circumvented if suitable stabilization techniques are used. The purpose of this paper has been to introduce a stabilization technique allowing the use of equal order of interpolation shape functions for pressures and displacements, and, therefore, of elements such as the linear triangle T3P3. These stabilized elements do not present volumetric locking and can be used to obtain limit loads and failure mechanisms of foundations. The T3P3 element which has been used in the proposed boundary value problems, presents an excellent resolution of the failure surfaces, which can be further improved if adaptive remeshing techniques are used.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg, 174 (1999) 219-234
233
Acknowledgments T h e r e s e a r c h p r e s e n t e d in this p a p e r h a s b e e n s u p p o r t e d b y A g e n c i a Espafiola de C o o p e r a c i d n I n t e r n a c i o n a l ( A E C I ) o f the F o r e i g n A f f a i r s M i n i s t r y o f S p a i n a n d the E u r o p e a n U n i o n ( P r o g r a m A L E R T G e o m a t e r i a l s ) . T. Li a n d X. L i u w o u l d like to e x p r e s s t h e i r g r a t i t u d e to H o h a i U n i v e r s i t y ( N a n j i n g , China), for h a v i n g g r a n t e d a leave p e r m i s s i o n for T. Li a n d X. Liu, a n d to C e n t r o de E s t u d i o s y E x p e r i m e n t a c i d n de O b r a s Pfiblicas ( C E D E X ) a n d F u n d a c i d n A g u s t i n de B e t h e n c o u r t for h o s t i n g t h e m in Spain.
References [1] J. Peraire, M. Vahdati, K. Morgan and O.C. Zienkiewicz, Adaptive remeshing for compressible flow computations, J. Comput. Phys. 72 (1987) 449-466. [2] M. Pastor and M. Quecedo, A patch test for mesh alignment effects in localized failure, Comm. Appl. Numer. Methods Engrg. I1 (1995) 1015-1024. [3] J.C. Nagtegaal, D.M. Parks and J.R. Rice, On numerically accurate finite element solutions in the fully plastic range, Comput. Methods Appl. Mech. Engrg. 4 (1974) 153-177. [4] S.W. Sloan and M.F. Randolph, Numerical prediction of collapse loads using finite element methods, Int. J. Numer. Anal. Methods Geomech. 6 (1982) 47-76. [5] O.C. Zienkiewicz, M. Huang and M. Pastor, Localization problems in plasticity using finite elements with adaptive remeshing, Int. J. Numer. Anal. Methods Geomech. 19 (1995) 127-148. [6] 1. Babfiska, The finite element method with Lagrange multipliers, Numer. Math. 20 (1973) 179-192. [7] F. Brezzi, On the existence, uniqueness and approximation of saddle point problems arising from Lagrangian multipliers, RAIRO 8-R2 (1974) 129-151. [8] O.C. Zienkiewicz, S. Qu, R.L. Taylor and S. Nakazawa, The patch test for mixed formulations, Int. J. Numer. Methods Engrg. 23 (1986) 1873-1883. [9] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155-164. [10] M.A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26 (1955) 182-185. [11] O.C. Zienkiewicz, C.T. Chan and P. Bettess, Drained, undrained, consolidating dynamic behaviour assumptions in soils, Geotechnique 30 (1980) 385-395. [12[ O.C. Zienkiewicz and P. Bettess, Soils and other saturated media under transient, dynamic conditions: general formulation and the validity of the various simplifying assumptions, in: G.N. Pande and O,C. Zienkiewicz, eds., Soil mechanics: Transient and Cyclic loads, Chap. 1 (Wiley, New York, 1982). 1113] O.C. Zienkiewicz and T. Shiomi, Dynamic behavinur of saturated porous media: The generalized Blot formulation and its numerical solution, Int. J. Numer. Anal. Methods Geomech. 8 (1984) 71-96. [14] F. Brezzi and J. Pitkaranta, On the stabilization of finite element approximations of the Stokes problem, in: W. Hackbusch, ed., Efficient Solutions of Elliptic Problems, Notes on Numerical Fluid Mechanics 10 (Vieweg, Wiesbaden, 1984) 11-19. [15] T.J.R. Hughes, L.P. Franca and M. Balestra, A new finite element formulation for fluid dynamics. V. Circumventing the BabfiskaBrezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal order interpolation, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85-99. [16] M. Hafez and M. Soliman, Numerical solution of the incompressible Navier-Stokes equations in primitive variables on unstaggered grids, Proc. AIAA Conf. 91 1561 -CP, 368-379, 1991. [17] O.C. Zienkiewicz and J. Wu, Incompressibility without tears! How to avoid restrictions of mixed formulations, Int. J. Numer. Methods Engrg. 32 (1991) 1184-1203. [18] O.C. Zienkiewicz, M. Huang, J. Wu and S. Wu, A new algorithm for the coupled soil - pore fluid problem, Shock Vib. 1 (1993) 3-14. [19] O.C. Zienkiewicz, M. Huang and M. Pastor, Computational soil dynamics--A new algorithm for drained and undrained conditions, Comput. Methods Adv. Geomech, 47-59. H.J. Siriwardne and M.M. Zaman, eds. (Balkema, 1994). [20] M. Pastor, M. Quecedo and O.C. Zienkiewicz, A mixed displacement-pressure formulation for numerical analysis of plastic failure, Comput. Struct. (1996). [21] M. Pastor, O.C. Zienkiewicz, T. Li, L. Xiaoqing and M. Huang, Stabilized finite elements with equal order of interpolation for soil dynamics problems, Arch. Comput. Mech. 6(1) (1999) 3-34. [22] I. Vardoulakis and J. Sulem, Bifurcation Analysis in Geomechanics (Blakie Academic & Professional, 1995). [23] A.J. Chorin, Numerical solution of incompressible flow problems, Studies Numer Anal. 2 (1968) 64-71. [24] G.E. Schneider, G.D. Raithby and M.M. Yovanovich, Finite element analysis of incompressible flow incorporating equal order pressure and velocity interpolation, in: C. Taylor, K. Morgan and C. Brebbia, eds., Numerical Methods for Laminar and Turbulent Flow (Pentech Press, Plymouth, 1978). [25] M. Kawahara and K. Ohmiya, Finite element analysis of density flow using velocity correction method, Int. J. Numer. Methods Fluids 5 (1985) 981-993. [26] O.C. Zienkiewicz, R. Codina, K. Morgan and B.V.P. Satya Sai, A general algorithm for compressible and incompressible flow. Part I: The split characteristic based scheme, Int. J. Numer. Methods Fluids 20 (1995) 869-885.
234
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
[27] W.F. Chen, Limit Analysis and Soil Plasticity (Elsevier, 1975). [28] R. Codina, M. V~zquez and O.C. Zienkiewicz, A fractional step method for compressible flows: boundary conditions and incompressible limit, in: M. Morandi Cecchi, K. Morgan, J. Periaux, B.. Schrefler and O.C. Zienkiewicz, eds., Proc. Int. Conf. on Finite Elements in Fluids--New Trends and Applications, Venezia, 15-21 Oct. 1995, pp. 409-418. [29] O.C. Zienkiewicz, J. Rojek, R.L. Taylor and M. Pastor, Triangles and Tetrahedra in Explicit Dynamic Codes for Solids, Int. J. Numer. Methods Engrg. 43 (1998) 565-583.
ELSEVIER
Stabilized low-order finite elements for failure and localization problems in undrained soils and foundations M. P a s t o r " ' * , T. Li b, X. L i u b, O.C. Z i e n k i e w i c z c aCentro de Estudios y Experimentacifn de Obras Pt~blicas and ETS de lngenieros de Caminos, Madrid, Spain bUniversity of Hohai Nanjing, PR China ~Universib' College of Swansea, Swansea, UK Received 20 April 1998; revised 10 July 1998
Abstract Geomaterials in general, and soils in particular, are highly nonlinear materials presenting a very strong coupling between solid skeleton and intersticial water. In the limit of zero compressibility of water and soil grains and zero permeability (which correspond to the classical 'undrained' assumption of Soil Mechanics), the functions used to interpolate displacements and pressures must fulfill either the Babuska-Brezzi conditions or the much simpler patch test proposed by Zienkiewicz and Taylor. These requirements exclude the use of elements with equal order interpolation for pressures and displacements, for which spurious oscillations may appear. The simplest elements with continuous pressures which can be used in 2D are the quadratic triangle and quadrilateral with linear and bilinear pressures, respectively. The purpose of this paper is to present a stabilization technique allowing the use of both linear triangles for displacements and pressures (T3P3) and bilinear quadrilaterals (Q4P4). The proposed element will be applied to obtain limit loads and failure surfaces in simple boundary value problems for which analytical solutions exist. © 1999 Elsevier Science S.A. All fights reserved.
1. Introduction
Low-order finite elements, such as linear triangles and bilinear quadrilaterals in 2D, provide some interesting features such as simplicity and low computational cost. In recent years, they have been found to provide an excellent definition for shocks and discontinuities in Fluid Dynamics. An additional advantage of triangles and tetrahedra is their flexibility for use in adaptive remeshing codes [1]. The situation is far less favourable in Solid Mechanics, as these advantages are outweighed by the poor behaviour exhibited by linear triangles and quadrilaterals in bending dominated conditions, or when the material behaviour is quasi-incompressible. Concerning the first problem, Pastor and Quecedo described in [2] how the failure mechanism depended on bending behaviour. Nagtegaal [3] and Sloan [4] described in their early work the second problem, which has been known to introduce important errors in the computation of collapse loads and failure mechanisms. The means to overcome this difficulty abound. In a previous paper [5] it was shown how mixed displacement-pressure formulations can improve dramatically the accuracy of collapse load computations. In the incompressible limit, the function spaces used to interpolate both displacements and pressures must fulfill the Babuska-Brezzi conditions [6,7] or the much simpler Zienkiewicz-Taylor patch test for mixed formulations [8]. Otherwise, spurious 'checker-board' oscillations may appear in the pressure field. These conditions exclude elements with equal order of interpolation of both variables. Therefore, elements such as the T3P3 triangle (linear displacements and pressures) or the Q4P4 quadrilateral (bilinear for * Corresponding author. 0045-7825/99/$ - see front matter PII: S 0 0 4 5 - 7 8 2 5 ( 9 8 ) 0 0 3 1 6-8
© 1999 Elsevier Science S.A. All fights reserved.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
220
displacements and pressures) cannot be safely used. The simplest triangles and quadrilaterals passing the test are the T6P3 (6-noded quadratic triangle for displacements and 3-noded linear triangle for pressure) and the Q8P4 quadrilateral. In fact, the T6P3 element was used in the computations described in [5] showing an excellent behaviour. In Soil Mechanics, the problem is further complicated as there exists a strong coupling between the solid skeleton and pore fluid which has to be taken into account [9-11]. One of the most widespread formulations of the problem is the so-called u-pw formulation, in which both displacements and pore pressures are chosen as field variables to describe the problem [12,13]. When approaching zero permeability and incompressibility of soil grains and pore water, the system of discretized equations is of a similar nature to that found when using mixed formulations for Solid Mechanics problems, and again, either the Babuska-Brezzi conditions, or the patch test have to be fulfilled. Again, the simplest 2D elements which can be used are the T6Pw3 and the Q8P~4 quadrilateral. The restrictions imposed by the Babuska-Brezzi conditions on both types of problems can be circumvented if special stabilization techniques are used. Such methods have been investigated in the past in the context of Fluid Dynamics 114-17], and have been extended to both Solid Mechanics and Geomechanics problems [18-21]. The purpose of this paper is to introduce a simple stabilization technique which can be used to develop suitable elements with equal order of interpolation of pressures and displacements for the computation of collapse loads of foundations and earth structures. The method is based on the fractional step algorithm proposed by the authors [21,29] for dynamic problems, and it will be applied here to problems of failure by localization. Localization of plastic deformation in shear bands is one of the most frequently found mechanisms leading to failure in soils. At a local level, it is possible to know the instant at which a shear band can appear and its orientation. In the case of elasto-plastic materials following the Tresca or Von Mises yield criteria, it corresponds to the instant at which the plastic modulus becomes zero, i.e. at the yield surface, and the inclination in 2D problems is 45 ° relative to the principal stress axes. It should be mentioned that the problem becomes ill-posed from a mathematical point of view, resulting in non-uniqueness and mesh-dependence of the solution. For a review of the research work done in the past years in this area, the reader may consult Vardoulakis and Sulem [22] and the references given therein. This paper will not focus on non-uniqueness mesh-dependence issues, but on how to develop simple Finite Elements suitable for computation of collapse loads and failure mechanisms of soils and foundations in the undrained state in which the pore pressure Pw plays the role of the mean pressure p.
2. The fractional-step algorithm for soil dynamics The fractional step algorithm was introduced by Chorin [23] as a device to allow the use of standard time integration techniques in fluid dynamics problems. Among the several alternative formulations some forms were found later to provide the required stabilization for elements with equal order of interpolation of velocities and pressures. The discovery was first made by Schneider et al. [24], and Kawahara and Ohmiya [25], and later justified by Zienkiewicz [26]. In fluid dynamics, Chorin's projection step method is based on splitting the increment of velocity into two orthogonal parts, consisting on a solenoidal field and the gradient of a scalar function:
v=w+V~
(1)
In the case of incompressible fluid flow the governing equations are 013
Ot
- STs - V p
vTv =
0
where v is the velocity, p the pressure, and S the strain operator:
(2)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
"0
0
X
S=
0
0
0
o
0
0
o
0
0 Ox 0 Oz
Oy 0 0
221
(3) 0 0 Oy 0
0
~x-
Fractional step equations introduce an intermediate velocity v* obtained in a first step: U*
- - 1.3n
_ _ At
_ STS"
(4)
where s is the deviatoric stress tensor written in vector form. The second step is D
n+l
D: ~
--
At
--719 n + ]
with
~T
D n+ ]
=0
(5) n+l
It can be seen that v* has the two orthogonal components mentioned above, the divergence free field v and the gradient field Vp. Although in its original form the scheme consisted of an explicit step followed by an implicit step (pressure was considered at time t "+1 in the second equation), further research showed that better results could be obtained by evaluating pressure at some intermediate time t" +02 and imposing the incompressibility condition at t "+°' [17], where t "+° is defined as (1 - O)t" + Ot"+1, with 0 ~< 0 ~< 1. The soil-fluid problem will be solved here using this approach. First, the balance of linear momentum for the mixture is written as dv STar ' --Vp.. +pm b -- p,. ~ - = 0
(6)
where o" is the effective stress tensor, Pw the pore pressure, p,. the mixture's density and b the body forces. Dirichlet boundary conditions are imposed on F. and tractions are prescribed on Ft. The balance of momentum and conservation of fluid phase equations can be combined, yielding: v T v -- V T k V p w + Q* [J"' = O
(7)
where we have assumed that body forces are divergence-free and that the permeability tensor is isotropic: k~--k/. 1/Q* characterizes the combined compressibility of the fluid and solid phases, and depends on voids ratio n, and bulk moduli of grains and water, Ks and K t 1
n
Q*-Kj
1-n
+
K,.
Pressures and fluxes are assumed to be given on boundaries Fp and Fq, respectively. In the incompressible undrained limit equation 7 reduces to vTv = 0. The first step of the proposed extension for saturated soils consists in discretizing the equations in the time domain using the intermediate velocity v*: V*
P"
- - V '~
At
-
AU* P" At
_
S,To .'~
+ pmb
(8)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
222
V n+l -- O*
P~
At
AV**
- Pm
=
At
~]] "+02
(9)
- p"
where we have introduced: n+O.
n
(10)
Pw " : P w + O 2 A p w
Therefore, the increment of velocity Av" is decomposed into two parts, Av* and Av**. Eq. 7 is written as 1 Ap w T n+Ol__VTvn+O Q* At = V (kVpw) [
I
(11)
where we have introduced O"+°~ = V" + 01 Av" = (1 - 01)v" + rio
n+l
Substituting now the value of v "+°~ in (11) we arrive to 1 Ap~ Q* At
--
VT.kT, .+o,. t p~ ) - (1 - O~)VTv n
~..,v .+1
-fflv
v
(12)
If we obtain velocity at time n + 1 from (9) 1 V.+ 1 = O* -- - - AtVp~ +°~ P..
(13)
and substitute it in (12), it results in 1 AAt p . , _ V T . k( ~ .p~ + o l " ) - (1 _ Q*
01)VTv
n
n+02 ~ - 01VT ( v * - - - 1 At Pw ) P.,
(14)
from which we finally get At 0102V2 ] Apw : V T ( k V p ; ) --1-~At 0 1V 2P~n _ V T v n -- 0~VT(v * - v") 0 j V T k V - - -p,.
Q* At
(15)
The solution of the problem is therefore obtained in three steps (1) Solve Eq. (8) for v* explicitly m tl n+ 1 ~! tt (2) Obtain Pw from (15), and Pw =Pw + APw (3) Use (9) to determine v "+l once Apw and v* are known The above equations can be discretized in space applying standard Galerkin procedures, and the result is the following system of algebraic equations: ij*=On+AtM]
{
1 f~+R~+
f To
N.pwdF-
f BVo-'nd/2 }
with =
T
N . o'.n d / '
with (o'.n); = ~r/jnj, or introducing f * " =f.n + R n + f r N.pwT . d E
O*=O"+AtMfl{f*.- f Bvtr'" dS2} where
(16)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
M. = f pmN uN. T d+Q B =SN.
223
(17)
It is important to note here that, in addition to external forces f : and reactions R", f . n now includes integrals over the boundary of pressures at time n, as proposed in [28] and [29]. The second discrete equation is
/'
~ Mp + Oln + At OlO2H,
)
Aft w =~r f.+o,
- Hp+n - 01 At H , f f nw
- QTt3" + OiGT A0* - fr. +rq NTp AO n d F
(18)
where Mp :
f,- ~ N~N
e dl2
H = f VTNpk.,VNp d n H* = f VTNpVNv d$2 Q = f B TmNp dO
G =f --~1 N.TVNe d~ REMARK 1. In the discretization process of terms involving second-order derivatives with respect to space, we obtain the following boundary integrals:
1 T 0 (Ap~)dF+OIAt O~02AtfFv +rq~mNP~n -01
+Gp1---NT m P OP~W On d F
frv+G 1P"N T P ao* d r
(19)
(20)
which can be written as
f
T[-At Op+ .+o~_ A v * ] d E
If we now write Eq. (13) as
1
Av"-Av*+--At
~jp~, n+02 = 0
Pm
and project it along the normal to the boundary, we obtain At ap+ In+°2
Pm On
- A v * +Av~=O
from which we can write the contribution of the boundary terms as
fFp+G NTp Av.n d F which will be assumed to be zero following Codina et al. [28]. It follows from here that Q T AV- nn :
- G pT A ~ n.
(21)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
224
REMARK 2. The force t e r m f ~ +°~ is given by f.+o, =
ap
f~, + ]'q
c9 n + O~ k-~n Pw - dF
and reduces to zero in the incompressible-undrained limit. Finally, the third equation is
~"+J = 6" - At M~lGpff~, +°2
(22)
The computational scheme requires of course that the stresses are evaluated at the end of the interval, i.e. at time n + 1. This can be done by using Ali" = t~ "+° At and accordingly computing o" '~" +~) from the constitutive relation.
3. An asymptotic form for steady state conditions in the incompressible undrained limit In the incompressible undrained limit matrices Alp and H are zero as k w and 1/Q* are zero. Eq. (16)
- ~ M.(tS* - 15")
+ Jr
.Pw a ~
(23)
reduces to
1 A--tM,(O * - ts") = f " + R" - Kti
(24)
when the material behaviour is elastic. Here, K is the stiffness matrix and ti n the nodal displacements at time t n. In steady state conditions (6* = 6n) Eq. (23) yields
f
BTtT
'
dO-
dF=f. + R
(25)
where we have omitted superscripts referring to time step. Second (18) and third (22) steps are written now as (At O~02H*) A,pw-" = -01 At H*ff~
-
QT o" + oia f A~*
(26)
and 1
--M
At
~lr - n + 8 2 . t"o- n + l __ [ ) , ) = -- pPw
Eqs. (25) and (27) can be combined to yield in the steady state
fBT~r'dO+G,Pw-fr
NT .pwdF=f. +R
or, taking into account that
(27)
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
225
Gpfiw = f NTVp~ d~2 = - f VN~N~ d~2 fiw + f r N ~ P w d F = - Q TPw + fl N .Tp w dE f BTo ,
- =L+R d~f2 -Q Tm
which for elastic materials reduces to
K.a-Q
T-
p.=f
+R
(28)
Next, we will eliminate the intermediate velocity field iS* using Eqs. (26) and (27). From the latter, 1
~ M . ( v-~÷l - t3* - t~" + 0 ") = -Gp#-.+o~ w and Ate* = AO" +
A t M . JGFff~.,+°2
(29)
Substituting this value in (26) it results
O, Oz At H* A - "
-01 At H*ff"., _ QTe,, + O,GT AO" + OIAtGTM~ lG -"+°~ppw"
from which, after rearranging, we obtain
QT6" + OIAt(H* - GTM-1GI, ,,
t,)P'-"+°2~, = 0
where we have taken into account that fe = 0 as the permeability is zero. In the case of solids, steady state conditions imply ~" = 0, and above equation reduces to (H*
--
G pTM u
J
Gp)lJw = 0
which can be combined with the incompressibility constraint QTo = 0 or Qrt~ = 0 to give
QTa + a(H* - Gi~M~'Gp)ffw = 0 where a is a suitable constant. The system can be written as
f B . ¢r' - Q.15.,-f. - R = ~ = 0
(30)
n . - = ¢Pp = 0 QTti + ce(H* - GT= , , f i t I op)pw
(31)
If this system is nonlinear, it can be solved by using the Newton-Raphson method
F0"u o4
/
o~.
[ =
(32)
with a suitable Jacobian matrix j~k> which is given by (33) where K r is a consistent tangent stiffness matrix and A is:
A
=
a(H*
-- G pT M
u
J
G.) "
(34)
226
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
It is important to note that A has provided the stabilization required to circumvent Babuska-Brezzi conditions. Otherwise, the resulting system would have adopted the form
This type of system is frequently found in constrained problems, presenting pressure oscillations wherever the number of sc variables n~ ~< n,~, i.e. the number of ~b variables, as occurs when using equal order of interpolation for both variables. Concerning the stabilization matrix A added, it consists of two parts, H * which is a discrete Laplacian T --1 operator acting on the pressure, and G p M , Gp, which requires the inversion of M,. A can be obtained more economically using patches of elements surrounding each node npatch
A = ~
OjO2 At 2
H*-
GeM"
(36)
It is interesting to note that other stabilization terms proposed included the discrete Laplacian multiplied by a scalar factor [14,19,29].
4. Applications: numerical computation of collapse loads and failure mechanisms of foundations The purpose of this section is to show how simple linear triangles can provide accurate predictions of collapse loads and failure mechanisms of foundations. We will assume that both soil grains and pore water are incompressible and permeability is zero, so the behaviour will correspond to undrained conditions. The tests we propose are simple boundary value problems for which there exist analytical solutions of collapse loads against which finite element predictions can be compared [27]. Indeed, they have been used in a previous paper [5], where the objective was to show how the use of mixed elements together with adaptive remeshing techniques improved the results obtained with classical displacement based formulations. Two new aspects have been introduced in this paper: • First, we will use elements with equal order of interpolation for both displacements and pressures • Second, we will use stabilized linear triangles for both displacements and pressures which are simpler than the quadratic triangles used in the previous paper. Of course, we will use adaptive remeshing techniques to obtain an optimal resolution of the failure surface. The interested reader is referred to [1] and [5] where the procedure is explained in some detail. In the two examples that follow the soil has been assumed to be an elastoplastic material with a Tresca yield surface and no hardening (H = 0). The problem can be characterized in terms of the following non-dimensional parameters E / Q = 2000
p = 0.49 where E and ~ are Young's modulus and Poisson's ratio, and Cu is the undrained soil cohesion. Load is applied through rigid and rough footings which have been modelled as elastic, with an elastic moduli one hundred times higher than that of the soil. The problems will be solved controlling the vertical displacement 6 of the point at which load is applied. Load-displacement curves will be presented in terms of a dimensional loads F / ( B C , ) and displacements a = (E6)/BC,, where B is the width of the footing. 4.1. Strip foundation
The first example consists of a rigid and rough strip footing resting on an undrained soil layer. The geometry and boundary conditions are given in Fig. 1. Due to symmetry, only one half of the domain will be discretized, but the results will be given for the whole domain.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
227
F
I
"I
I
~'J
I
Fig. 1. Strip footing foundation: geometry and boundary conditions.
Fig. 2. Initial mesh.
The theoretical collapse load for this problem is given by F BC, = (2 + rr) The analysis begins with the mesh given in Fig. 2. The results at a value of the non-dimensional displacement cr of 25 are given in Figs. 3-5 where the deformed mesh, displacement vectors and contours of effective plastic strain are shown. Load vs. displacement curves can be seen in Fig. 14. From these results two new meshes are obtained using adaptive techniques (Figs. 6 and 10). Results are given in Figs. 7 - 9 for the first mesh and Figs. 11-13 for the second. The process of adaptive remeshing ends here as the two last force-displacement curves coincide.
Fig. 3. Deformed mesh.
...... .
•
~-~
r[hlll i hhli I .,,,~.,,.,,,
, • ,, ~ ' , , ~ J u v
I I.II~4~--,.~/" / / / , ,
. :.
. .
: .
Fig. 4. Displacement vectors.
Fig. 5. Contours of effective plastic strain.
228
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
Fig. 6. Adaptive mesh 1.
Fig. 7. Deformed mesh.
Fig. 8. Displacement vectors.
Fig. 9. Contours of effective plastic strain.
Fig. lO. Adaptive mesh 2.
Fig. 11. Deformed mesh.
Fig. 12. Displacement vectors.
Fig. 13. Contours of effective plastic strain.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 2 1 9 - 2 3 4
229
Limiting loads by different meshes
for foundation problem 6
5.5 5 4.5
,,! 13 al
4
I
l
l
t
l
l
l
l
l
/
l
l
l
l
l
l
..........
l
l
l
l
l
l
l
initialmesh
. . . .
adapevemesh1
-
-
adaptive mesh 2
I
I
3.5 3 2.5 2 I
!
I
3
I
5
I
I
7
I
I
9
I
11
I
l
13
I
15
I
I
i
i
17
|
19
I
21
I
I
23
I
25
dE/BCu Fig. 14, Strip foundations: load vs. displacement.
It is important to notice the good resolution of the failure mechanism, with an error in the computed collapse load below 5%.
4.2. Vertical cut
The last example we will consider is that of the vertical cut shown in Fig. 15, where both the geometry and boundary conditions are shown.
F ~t
=~z =sd
C~
II ,k~,%
II ,vk.~
&
II
l
:3"
u, -- O,t, = 0 ,
J(
IOta
Fig. 15. Vertical cut problem: geometry and boundary conditions.
230
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
Fig. 17 shows the initial mesh used in the analysis, and the corresponding results can be seen in Figs. 18-20. The adaptive mesh refinement process leads to the mesh shown in Fig. 21, with results given in Figs. 22-24. We then perform a last mesh refinement to arrive at the mesh of Fig. 21, for which the results are given in Figs. 26-28. Concerning the limit load, it is found again to agree well with the analytical value given by F
BC,
-
2.0
Limiting loads by different m e s h s
for vertical cut problem 2.5
I
I
I
I
I
t
I
I
/ S; "
3 U
m
..........
'" inil~l mesh
.
~
.
-
.
.
-
m
u
h
l
~p~ve me~ 2
1.5
1
I
I
I
I
I
I
I
l
2
3
4
S
6
7
8
9
t0
dF./BCu Fig. 16. Vertical cut problem: load vs. displacement plot.
1 Fig. 17. Vertical cut problem: initial mesh.
Fig. 18. Deformed mesh.
231
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
~X.s.~.3~%,.
- -
J~ /,,.,.......;... t // I/ /
/
•
,
h Fig. 19. Displacement vectors.
Fig. 20. Contours of effective plastic strain.
i
\
i
J 1
/ Fig. 21. Adaptive mesh
Fig. 22. Deformed mesh.
1
Fig. 23. Displacement vectors.
Fig. 24. Contours of plastic strain.
232
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 2 1 9 - 2 3 4
i
l V Fig. 25, Adaptive mesh 2.
/ /
/
Fig. 26. Deformed mesh.
//./I
,,///f,,#...".
/ /,
Fig. 27. Displacement vectors.
Fig. 28. Contours of plastic strain.
5. Conclusions
The restrictions imposed by Babuska-Brezzi conditions or the Zienkiewicz and Taylor patch test for mixed formulations do not allow the use of equal order of interpolation shape functions for displacements and pressures in mixed formulations. Therefore, the simplest elements with continuous pressures which can be used are the T6P3 triangle with quadratic displacements and the Q8P4 quadrilateral. This restriction can be circumvented if suitable stabilization techniques are used. The purpose of this paper has been to introduce a stabilization technique allowing the use of equal order of interpolation shape functions for pressures and displacements, and, therefore, of elements such as the linear triangle T3P3. These stabilized elements do not present volumetric locking and can be used to obtain limit loads and failure mechanisms of foundations. The T3P3 element which has been used in the proposed boundary value problems, presents an excellent resolution of the failure surfaces, which can be further improved if adaptive remeshing techniques are used.
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg, 174 (1999) 219-234
233
Acknowledgments T h e r e s e a r c h p r e s e n t e d in this p a p e r h a s b e e n s u p p o r t e d b y A g e n c i a Espafiola de C o o p e r a c i d n I n t e r n a c i o n a l ( A E C I ) o f the F o r e i g n A f f a i r s M i n i s t r y o f S p a i n a n d the E u r o p e a n U n i o n ( P r o g r a m A L E R T G e o m a t e r i a l s ) . T. Li a n d X. L i u w o u l d like to e x p r e s s t h e i r g r a t i t u d e to H o h a i U n i v e r s i t y ( N a n j i n g , China), for h a v i n g g r a n t e d a leave p e r m i s s i o n for T. Li a n d X. Liu, a n d to C e n t r o de E s t u d i o s y E x p e r i m e n t a c i d n de O b r a s Pfiblicas ( C E D E X ) a n d F u n d a c i d n A g u s t i n de B e t h e n c o u r t for h o s t i n g t h e m in Spain.
References [1] J. Peraire, M. Vahdati, K. Morgan and O.C. Zienkiewicz, Adaptive remeshing for compressible flow computations, J. Comput. Phys. 72 (1987) 449-466. [2] M. Pastor and M. Quecedo, A patch test for mesh alignment effects in localized failure, Comm. Appl. Numer. Methods Engrg. I1 (1995) 1015-1024. [3] J.C. Nagtegaal, D.M. Parks and J.R. Rice, On numerically accurate finite element solutions in the fully plastic range, Comput. Methods Appl. Mech. Engrg. 4 (1974) 153-177. [4] S.W. Sloan and M.F. Randolph, Numerical prediction of collapse loads using finite element methods, Int. J. Numer. Anal. Methods Geomech. 6 (1982) 47-76. [5] O.C. Zienkiewicz, M. Huang and M. Pastor, Localization problems in plasticity using finite elements with adaptive remeshing, Int. J. Numer. Anal. Methods Geomech. 19 (1995) 127-148. [6] 1. Babfiska, The finite element method with Lagrange multipliers, Numer. Math. 20 (1973) 179-192. [7] F. Brezzi, On the existence, uniqueness and approximation of saddle point problems arising from Lagrangian multipliers, RAIRO 8-R2 (1974) 129-151. [8] O.C. Zienkiewicz, S. Qu, R.L. Taylor and S. Nakazawa, The patch test for mixed formulations, Int. J. Numer. Methods Engrg. 23 (1986) 1873-1883. [9] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155-164. [10] M.A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26 (1955) 182-185. [11] O.C. Zienkiewicz, C.T. Chan and P. Bettess, Drained, undrained, consolidating dynamic behaviour assumptions in soils, Geotechnique 30 (1980) 385-395. [12[ O.C. Zienkiewicz and P. Bettess, Soils and other saturated media under transient, dynamic conditions: general formulation and the validity of the various simplifying assumptions, in: G.N. Pande and O,C. Zienkiewicz, eds., Soil mechanics: Transient and Cyclic loads, Chap. 1 (Wiley, New York, 1982). 1113] O.C. Zienkiewicz and T. Shiomi, Dynamic behavinur of saturated porous media: The generalized Blot formulation and its numerical solution, Int. J. Numer. Anal. Methods Geomech. 8 (1984) 71-96. [14] F. Brezzi and J. Pitkaranta, On the stabilization of finite element approximations of the Stokes problem, in: W. Hackbusch, ed., Efficient Solutions of Elliptic Problems, Notes on Numerical Fluid Mechanics 10 (Vieweg, Wiesbaden, 1984) 11-19. [15] T.J.R. Hughes, L.P. Franca and M. Balestra, A new finite element formulation for fluid dynamics. V. Circumventing the BabfiskaBrezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal order interpolation, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85-99. [16] M. Hafez and M. Soliman, Numerical solution of the incompressible Navier-Stokes equations in primitive variables on unstaggered grids, Proc. AIAA Conf. 91 1561 -CP, 368-379, 1991. [17] O.C. Zienkiewicz and J. Wu, Incompressibility without tears! How to avoid restrictions of mixed formulations, Int. J. Numer. Methods Engrg. 32 (1991) 1184-1203. [18] O.C. Zienkiewicz, M. Huang, J. Wu and S. Wu, A new algorithm for the coupled soil - pore fluid problem, Shock Vib. 1 (1993) 3-14. [19] O.C. Zienkiewicz, M. Huang and M. Pastor, Computational soil dynamics--A new algorithm for drained and undrained conditions, Comput. Methods Adv. Geomech, 47-59. H.J. Siriwardne and M.M. Zaman, eds. (Balkema, 1994). [20] M. Pastor, M. Quecedo and O.C. Zienkiewicz, A mixed displacement-pressure formulation for numerical analysis of plastic failure, Comput. Struct. (1996). [21] M. Pastor, O.C. Zienkiewicz, T. Li, L. Xiaoqing and M. Huang, Stabilized finite elements with equal order of interpolation for soil dynamics problems, Arch. Comput. Mech. 6(1) (1999) 3-34. [22] I. Vardoulakis and J. Sulem, Bifurcation Analysis in Geomechanics (Blakie Academic & Professional, 1995). [23] A.J. Chorin, Numerical solution of incompressible flow problems, Studies Numer Anal. 2 (1968) 64-71. [24] G.E. Schneider, G.D. Raithby and M.M. Yovanovich, Finite element analysis of incompressible flow incorporating equal order pressure and velocity interpolation, in: C. Taylor, K. Morgan and C. Brebbia, eds., Numerical Methods for Laminar and Turbulent Flow (Pentech Press, Plymouth, 1978). [25] M. Kawahara and K. Ohmiya, Finite element analysis of density flow using velocity correction method, Int. J. Numer. Methods Fluids 5 (1985) 981-993. [26] O.C. Zienkiewicz, R. Codina, K. Morgan and B.V.P. Satya Sai, A general algorithm for compressible and incompressible flow. Part I: The split characteristic based scheme, Int. J. Numer. Methods Fluids 20 (1995) 869-885.
234
M. Pastor et al. / Comput. Methods Appl. Mech. Engrg. 174 (1999) 219-234
[27] W.F. Chen, Limit Analysis and Soil Plasticity (Elsevier, 1975). [28] R. Codina, M. V~zquez and O.C. Zienkiewicz, A fractional step method for compressible flows: boundary conditions and incompressible limit, in: M. Morandi Cecchi, K. Morgan, J. Periaux, B.. Schrefler and O.C. Zienkiewicz, eds., Proc. Int. Conf. on Finite Elements in Fluids--New Trends and Applications, Venezia, 15-21 Oct. 1995, pp. 409-418. [29] O.C. Zienkiewicz, J. Rojek, R.L. Taylor and M. Pastor, Triangles and Tetrahedra in Explicit Dynamic Codes for Solids, Int. J. Numer. Methods Engrg. 43 (1998) 565-583.