(5) The problem of element performance in finite element computations of elasto-plastic media exhibiting incompressible or dilatant behavior has been investigated intensively in recent years with the aim of developing low order elements which do not lock. Early approaches to overcome locking due to incompressibility used reduced selective integration, mean dilatation (Nagtegaal & al. 1974), and later on the B-bar approach (Hughes 1980). More recent approaches adopted mixed formulations, derived from the Hu-Washizu variational principle, and enhanced strain approximations; this led to the so-called Enhanced Assumed Strain approach (EAS) (Simo & Rifai 1990, Taylor & al. 1976). Unfortunately, this method does not apply to many well known finite elements (like linear triangles in 2D or linear tetrahedrons in 3D) and, moreover, it requires the specific design of the enhanced part of the strain field for each element separately. Yet another approach is presented in this paper. It is based on a mixed continuous displacement-pressure formulation, and complements the Galerkin scheme by least-squares terms which enhance its stability.
The B-bar approach introduces a modification in the volumetric contribution to the standard B matrix of space derivatives of shape functions Bvol, which is underintegrated or averaged over the element. This formulation has proven to be appropriate to overcome locking due to incompressibility in very general situations, but poor performance was observed in dilatant cases (de Borst & Groen 1995).
Several comprehensive studies on EAS elements, based on the Hu-Washizu functional, have recently been done by different authors (Andelfinger & Ramm 1993) (Groen 1994) (de Borst & Groen 1995), mainly in conjunction with simple isotropic elastoplastic models. In these cases, it has been shown that a simple enhancement of normal strains only, is sufficient to overcome the problem of volumetric locking due to dilatant constitutive behavior. In the case of real rock-like layered materials, the need of enhancement for both normal and shear strains is confirmed by a number of numerical tests (Truty & al. 1997). In the case of an enhanced shear strain field, some spurious mechanisms may however be produced and, so far, it has been impossible to overcome locking when a mixture of elements of different types are to be used.
This new approach follows ideas which were developed and proven successful for fluid mechanics in (Hughes & al., 1986, Franca, 1987) and following papers. A least-squares term is added to the standard Galerkin formulation for stabilization. The following general form of the least square term is discussed in Truty & Zimmermann 1997, Truty 2002, and Commend 2001).
(5)
where the differential operator L takes the form:
(6)
and [σx,σy, σxy, σz, σxz, σyz]T is the corresponding stress vector. uh, ph are the approximations of the solution, and f is the body force.
The proposed formulation is designed for a class of rate independent multisurface elastoplastic models. The constitutive equation for the stress and the pressure p are as follows:
; 
(7)
where ε is the
total strain and εp is the plastic strain. The
general deviatoric projection 
of
the elastic stiffness matrix D and the generalized elastic
bulk modulus 
are
defined as follows :
;
(8)
where symbol 1 is the vector representation of the Kronecker delta. The exact deviatoric-volumetric split is retrieved for ξ = 0 and such a version was initially used (Truty & Zimmermann 1997). Unfortunately this split is not acceptable for plastic surfaces depending exclusively on I1. Thus some perturbation factor ξ has to be introduced but its value should be small enough to be close to the pure deviatoric-volumetric split (10-6 to 10-5). The flow rule and the hardening law take the following form:
(9)
(10)
where the plastic surface index is denoted by α, the flow direction by r, the hardening parameters by q, the hardening functions by h, and summation is performed only over active surfaces (α Î Jact). The loading- unloading conditions read:
(11)
where f designates yield surfaces.
The weak form of the momentum equation, constitutive equation for the mean stress can be expressed in the following form:
(12)
where 
is the imposed boundary traction,
(13)
The above equations can be written in combined form:
(14)
The minus sign for R2 is introduced here in order to preserve the positive definiteness of the resulting Jacobian matrix in Newton's scheme while adding the least-square term. Since the plastic flow rule and the hardening law are formulated in rate form they need to be integrated first and then substituted into the variational equations. A fully implicit integration scheme has been adopted and it results in the following formulae for plastic strains and hardening parameters at time step n+1:
(15)
(16)
The corresponding stress state defined at step n+1
can be expressed as follows:
(17)
The system of variational equations is nonlinear in arguments u,p. It can be solved with the aid of Newton's method but first it has to be linearized. Consistent linearization is an essential step to get an efficient iterative scheme. Also, inserting the integrated formulae for state variables σn+1, εpn+1, qn+1 into variational equations prior to linearization should preserve the quadratic convergence rate of the resulting Newton scheme.
The stabilization term is meant to enhance stability. The proposed expression (5) is the most general form of least-squares term. The general aim of adding this term is to prevent the effect of volumetric locking, to reduce the effect of stress oscillations, to render possible the use of equal order continuous u-p interpolation (including linear triangles), and, moreover, to enable the mixing of different types of elements to fit the desired geometry of the domain. In order to define the form of the stabilization term, let us split the stress state into deviatoric and volumetric parts, and rewrite the stabilization term in more detailed form as follows:
(18)
If the weighting term related to plastic strains
is omitted the following stabilization term can be obtained:
(19)
An even simpler version is obtained by taking as
a weighting operator only the gradient of the pressure and only this type
of stabilization is used in all the following applications. This term added
after the introduction of the finite element approximation has to be linearized
and finally brings a contribution to both the left- and the right-hand side
of the resulting linearized system of equations. It should be noted that
for some linear elements the first weighting term disappears automatically
(eg. for linear triangles). For Navier-Stokes problems the stabilization factor is defined as τ = αh2 / 2G (Franca 1987),
where α, the optimal stabilization factor, should be a constant for a given
element type, and can be obtained by a convergence analysis; G is the shear
modulus. The same τ factor is used here.
Detailed derivations
can be found in Commend (2001), Truty (2002), Commend & al.(2002), and
alternative equivalent approaches are described under the name of Finite
Increment Calculus by Onate(2000), and
Laplacian pressure operator scheme (LPOS) by Pastor
& al.(1997). All three approaches are compared in Commend (2001).
A superficial foundation on an incompressible (or dilatant) medium is subjected to a uniform loading, which is increased until a failure is detected. The problem is run here as a displacement driven on a mesh with 1296 elements, using a Drucker-Prager criterion, with size adjustment to a Mohr-Coulomb criterion, such that ultimate loads coincide, under plane strain conditions (Z_Soil 2002). E = 50000 kN/m2, ν = 0.3, γ = 0, c = 1kN/m2, friction angle φ=20o, dilatancy angle ψ=0o.
Analytical solutions to this problem
can be found in (Terzaghi,1951, Chen & Liu 1990), which yield an ultimate
pressure qu = ~18kN/m and also in Matar
and Salençon (1979) which yields qu = ~15.6 kN/m.
When standard bilinear displacement quad elements are used, results show an overshoot of theoretical solutions (Fig. 4), characteristic of locking; this also holds for constant strain triangles, or when a mixture of the two different types of elements is being used (Zimmermann & Commend 2001) .
Stabilized formulations are shown,
in Figure 4 (here LPOS formulation with 
= 0.5), to provide a distinct ultimate load in the
incompressible case. This is also true for a mesh composed of quads, triangles,
or both. The same formulation performs as
well in the incompressible and in the dilatant case (Commend 2001, Truty
2002) .
Earlier solutions to overcome locking,
like strain-projection 
(Fig.4),
enhanced assumed strains EAS or cross-diagonal triangles all solve the incompressible
case for quads, EAS also overcomes locking in dilatant media when the mesh
is composed of quads, but all fail when quads and triangles are used simultaneously
within the same mesh (Truty & al. 1997).
Stabilized
formulations may, however, require additional mesh refinement as compared
to
,
and a residual locking tendency may show up when coarse meshes are used (Truty
2002).
The use of a stabilized formulation in order to overcome pressure oscillations in 2-phase consolidation problems has also proven successful (Truty 2001, Zimmermann & al. 2001), as illustrated next on an essentially one-dimensional consolidation problem.
A
Figure 5. Pore pressure distributions using different weak forms.
The comparison of the numerical solutions obtained with a standard Galerkin formulation (SG) (which shows oscillations) and a Galerkin/Least-squares formulation (GLS) (which coïncides with the analytical solution), both obtained with a time step of 0.2 days, much smaller than the critical value, proves the effectiveness of the novel approach (Fig.5).