Elasticplastic consolidation analysis

of saturated soil
by a time discontinuous method
H. Cramer
Universitat Rostock, Fachgebiet f
ur Baustatik und Baudynamik
PF 1210, 23952 Wismar, Germany
[Link]@[Link]

R. Findei, W. Wunderlich
Technische Universitat M
unchen, Lehrstuhl f
ur Statik
Arcisstr. 21, 80333 M
unchen, Germany
findeiss@[Link]
ww@[Link]

Keywords: TimeDiscontinuousGalerkin, spacetime finite elements, consolidation,

saturated soil, time integration, localization, Cosserat continuum, regularization
Abstract. The subject of this paper is the numerical solution of initialboundaryvalue
problems that arise in the elasticplastic consolidation analysis of saturated soil. A consistent implementation of a time discontinuous Galerkin finite element procedure is discussed.
The basic idea of the presented method is to permit the unknown fields to be discontinuous at the discrete time levels. The integration over discrete time intervals is performed
with appropriate variational formulations of the governing equations using Galerkintype
weighting functions. Concerning the integration of the rate equations of elastoplasticity
a numerical procedure is proposed using Gaussintegration in time for the discretization
of the variational equations. For the global NewtonRaphson iterative method a tangent
operator is derived by a consistent linearization of the nonlinear equations. The general applicability of the proposed integration scheme to general consolidation problems of
saturated soil is demonstrated by numerical examples.

Introduction

The coupled finite element analysis of elasticplastic consolidation of saturated soil requires the solution of discretized equations in space and time. Most of the applied time
integration procedures for quasistatic problems are based on a space discretization by
finite elements and a type of generalized trapezoidal rule in the time domain as described
by the authors [1]. These result in time stepping schemes where the unkown variables
at the end of each interval are calculated with a continuous progression of displacements
and pressures in time. A superior method has been presented by Ehlers [3] where the
resulting system of differential and algebaic equations is solved using a diagonally implicit
RungeKutta scheme.
In the paper an alternative approach using finite elements in both space and time is
discussed. In this approach the displacements as well as the pore water pressures are
approximated as piecewise linear functions within each time step. In contradiction to
continuous time stepping schemes, the basic idea of the presented method is to permit the unknown fields to be discontinuous at the discrete time levels. The integration
over discrete time intervals is performed with appropriate variational formulations of the
governing equations using Galerkintype weighting functions. The continuity of displacements and pore water pressures at the end of the time intervals is satisfied only in a
weak sense. This so called time discontinuous Galerkin (TDG) method originally was
proposed by Hulbert and Hughes [5] and was applied in structural dynamic analysis
by Wiberg and X.D. Li [9]. The presented paper deals with the application of the
TDG to the quasistatic analysis of consolidation and considers a consistent extension to
elasticplastic material behaviour.
Concerning the integration of the rate equations of elastoplasticity the generalized trapezoidal rule as well as the generalized midpoint rule have been suggested by Stein [7]. In
combination with a continuous time stepping algorithm the backward Euler scheme has
gained a wide acceptance (see Simo and Hughes [6]). However, this scheme is based on
a continuous approximation of the displacements and thus not applicable to the discontinuous Galerkin method. Therefore, in this paper a numerical procedure is proposed
using Gauss integration in time for the discretization of the variational equations. This
results in an implicit scheme in which the elasticplastic stressstrain relations are enforced at two Gauss points in each time interval. In order to ensure consistency with the
TDG method, the strain displacement relation as well as the rate equations are used in
appropriate variational forms.
In the first part of the following presentation the governing equations for saturated soil and
the formulation of a space time finite element is presented. The second part is concerned
with the integration of the rateequations of plasticity and the derivation of a consistent
tangent operator. The paper is completed with numerical examples that demonstrate the
applicability of the proposed method.

Governing equations

2.1

Formulation of the initial-boundary-value problem

The governing equations of single phase flow in a deforming porous medium can be formulated on the basis of the theory of porous media and the theory of mixtures. For
a detailed representation of this theory the interested reader is referred to the work of
Truesdell and Toupin [8], de Boer and Ehlers [2] and Ehlers [4].
The complete initialboundaryvalue problem in a domain bounded by and in a
time intervall I = [0, T ] is represented by a system of differential equations defining the
balance of linear momentum and the conservation of mass for the mixture as a whole:

div

= 0,
div ( + pI) + b

(1)




u k
F
grad p + b = 0 .
+
t

(2)

with the appropriate initial and boundary conditions:

u(t = t0 ) = u0 in
p(t = t0 ) = p0
= 0
uu

in
on u

(3)

( + pI) n t = 0

on t

u t = p Q =

p p = 0

on p

u t = p Q =

= 0
QnQ

on Q

Herein u and represent the displacements of the soil matrix and the effective stress
tensor actiong on the soil skeleton, p the pore water pressure, I the fundamental tensor
the volume force density. The dynamic viscosity of the pore fluid
of second order and b
and the permeability tensor are denoted by and k, respectively.

2.2

Equations of rate-independent plasticity

The material behaviour of soil is modelled by a nonlinear elasticplastic constitutive law.

The basic equations of rateindependent plasticity include the following aspects:
elastic stressstrain relationship with additive decomposition of the strainrates


(4)
= C : pl
kinematic equations

(t)
=


S
1

u=
u + (u)T =
Duu
t
t 2
t
3

(5)

 non associative flow rule and isotropic hardening law

pl =
yield condition

q = D

g
q

(6)

f (, q) 0

(7)

KuhnTucker complementarity conditions

0

f (, q) 0

consistency condition

f
(, q) = 0

(8)

f(, q) = 0

(9)

A more detailed discussion of the plasticity equations can be found in e.g. Simo and
Hughes [6]. The following presentation will be restricted to perfect plasticity, however,
the extension to hardening material is possible in a straightforward manner.

A space-time formulation for the consolidation problem

In order to solve the set of differential equations (1) and (2) of the quasistatic consolidation problem, a consistent approach using finiteelements in space and in time is discussed
in this chapter. The displacement and the pore water pressure fields are interpolated by
standard finiteelement shape functions in the spatial domain. In time they are approximated as piecewise linear functions over each time interval In = [tn , tn+1 ]. In contradiction
to continuous time stepping schemes, the basic idea of the presented method is to permit
jumps in the unknown fields at the discrete timelevels (see figure 1).
PSfrag replacements

FEsolution
(discontinuous)

displacement

exact solution
(continuous)

[un ]

u
n

u+
n

time
tn1

tn

tn+1
t

Figure 1: Illustration of the time discontinuous Galerkin method

t
P

discrete structure

PSfrag replacements X 2

timeslab
tn+1
t

tn
X

X2

X1
Figure 2: Linear 8noded spacetime finite element
Therefore it is appropriate to introduce the following notation in order to distinguish
between the two values at a discrete time tn :
u+
n =
u
n =

lim u(tn + ) = u(t+

n)

0+

(10)

lim u(tn + ) = u(t

n)

[un ] = u+
n un

The continuity of displacements and pore water pressures is satisfied in a weak sense
between the time intervals. In this way a time slab in spacetime is discretized by the use of
standard quadrilateral finiteelements in space and by linear finite element approximations
in time.
The set of differential equations (1) and (2) can be transformed into an equivalent variational form using Galerkintype weighting functions u and p. Notice that the jumps
in time for the displacements and pressures [un ], [pn ] and in addition for the stresses
[ n ] have to be taken into account. For convenience the vector notation is used in the
following:
Z
Z
Z
Z Z
T
T
T
T
T
u D u d + u D u m p d u b d uT t d dt+
In

T
T
(u+
n ) Du

[ n ] d +

Z Z

In

T
T
(u+
n ) Du

(11)

m [pn ] d = 0 ,

p m D u u d

pD Tp

k
D p p d +

dQ dt +
p Q

T
p+
n m D u [un ] d = 0 ,

(12)
where m = [1, 1, 1, 0]T and D u and D p are the standard gradient operators. In these
equations the spatial coordinates x as well as the time coordinate t are independent
5

variables. For a discretization in spacetime shape functions are used in an appropriate

product form:
,
u(x, t) = N u N t u

(13)

,
p(x, t) = N p N t p

(14)

and p
represent the generalized vector of discrete values at the element nodes
where u
for time tn and tn+1 , respectively.

+
+
n
n
u
p
,
.
=
=
u
p
(15)

n+1
n+1
u
p
In the following biquadratic shape functions for the displacements and bilinear shape
functions for the pore water pressures are used. This is equivalent to the 9/4node
quadrilateral element proposed in [1]:
N u = N u (x) ,

(16)

N p = N p (x) .

(17)

The shape functions in time N t represent a linear variation of u and p within each time
step and therefore takes the form:


tn+1 t t tn
,
N t = N t (t) =
t In .
(18)
t
t
After applying these discretizations to the variational formulation (11) and (12) a system
of algebraic equations is obtained. For implementational aspects especially in the context
of plasticity this system is transformed to an incremental form. With the definitions of
incremental values of displacements and pore water pressures:

un = u
n+1 un ,

(19)

pn = p
n+1 pn ,

(20)

the following system of equations can finally be obtained for a Newton scheme by a
consistent linearization:
k
k

k
ext +
int +
F (tn ) F (tn )
[un ]
K
K 12
L
0

11

ext

int
F (tn+1 ) F (tn+1 )
un
K 21 K 22

0
L
=

1
1

[p ]
0 6 tH 6 tH
0

p
T
1
1
pn
F
0
L
2 tH 2 tH
(21)
with the update:

k+1

[un ]

un

[pn ]

pn

[un ]

un

+
=

[pn ]

pn
6

[un ]

un
.

[pn ]

pn

(22)

Herein the following definitions are applied:

nel

K =

e=1

(D u N u )T C D u N u d e ,

(23)

(D u N u )T m N p d e ,

(24)

(D p N p )T

e
nel

L =

e=1

e
nel

H =

e=1

k
D p N p d e ,

(25)

where
denotes the standard finiteelement assembly operator, nel is the number of
elements and F ext , F int and F p are the vectors of external forces, internal forces and
equivalent fluidpressure forces, respectively, and are defined as follows:
nel

ext

(t) =

e=1

N Tu

N Tu D Tu (t) d e ,

d e +
b(t)

e
nel

int

(t) =

e=1

N Tu t(t) de ,

(26)

(27)


k

k 1

.
F p = L T u
+ tH p+
n + pn+1
n un+1
2

(28)

Since the first two lines in equation (21) include the nonlinearity of the material characteristics in the form of the path dependent matrices C obtained as consistent linearizations,
this system is solved iteratively by a NewtonRaphson algorithm. Therefore in each
iterationstep k the systemmatrix as well as the vector of internal forces have to be reassembled. In order to obtain a quadratic convergence rate a consistent tangent operator
according to the applied variational equations is derived.
In order to ensure consistency with the TDG method the straindisplacement relation and
the stressstrain relation have to be transformed into appropriate variational forms. With
the permission of an unsteady displacement field in time the jumps can be considered in
a variational sense by the integral form:
t
n+1

t+
n

T
S
S u dt + ( +
n ) [un ] = 0.

(29)

From the equations of rateindependent plasticity (4) and (6) the local form
C + C

g
= 0

t In

(30)

can be obtained. A variational form which takes unsteady fields at the discrete time levels
into account can be formulated by integration over the time interval In = [tn , tn+1 ] and
7

PSfrag replacements

Figure 3: Position of integrationpoints in spacetime. (2x2x2) scheme

using as a weighting function:
t
n+1

C + C

T
dt + (+
[ n ] = 0.
n)

(31)

t+
n

For the approximation of , , and an interpolation represented by discrete values

at Gauss points in space is used. By a similar approach , and are also represented
replacements
atPSfrag
two Gauss
points in time. Their positions in space are defined by the local coordinates
and and depend on the integration scheme (2x2 or 3x3) and their positions in time
are:
1
3
(32)
1/2 =
3
with the local time coordinate [1, 1]
tn
1

tn+1
1

time t

Figure 4: Position of Gauss points in time

Under the assumption of a linear variation within each time step can be expressed as



1
1
(33)
= (1 3 ), (1 + 3 ) ,
2
2
where represents the generalized vector of the Gauss point values:

=1
|
.
=
=2
|

(34)

Equivalent representations are obtained for and . This representation has the advantage that a corresponding Gauss integration of equations (29) and (31) results in
decoupled equations for the Gauss point values in space and for the variational form of
the kinematic relations in time.
8

Thus, when the displacement gradient u is given in the form:

1
+
u
u =
n+1 un ,
t

(35)

the values of the strain rates at the Gauss points in time can be calculated directly
from the displacement increments [un ] and un by:

3 1
D N
0
[un ]
1
u u

=
(36)

t
3 1
0
DuN u
un
With the assumption of a linear variation of the fields within In the stress values at the

Gauss points can be calculated from the values +

n and n+1 at the ends of a time interval
by

1 1
1
1
+
+

3
3

| =1
2
6
2
6

n
.
=
=
(37)

1
1
1
1
n+1
| =2
3
3

+
2

Setting
=

1
+

n
t n+1

(38)

and performing the numerical Gauss integration over the actual time intervall under the
assumption of linear elastic behaviour the following equations are obtained relating the
and
Gauss point values ,
:

n
C
G
+C
=

(39)

n
G and
where C,
are given by:

1
C
3

1+ 3
C
6

= t
C

G =

g 
1

|
1
|
2

0

g 

1 3
6
1
3

(40)

(41)

(42)

As the stress relation is formulated in terms of the discrete values at the spatial and
temporal Gauss points it is natural to ensure the yield condition also at the Gauss
points:

f ()|1
=
f ()
= 0.
f ()|2
9

(43)

Iterative calculation of the path dependent variables

4.1

Formulation of the implicit algorithm

This chapter is concerned with the determination of the correct stressfield (x, t) as
a consequence of the strainfield (x, t) which is the driving variable. Regarding the
restriction of the yieldcondition (7) has to be determined by the solution of a nonlinear
initialvalue problem. The fields D u u, and are evaluated at a distinct temporal
and plastic
Gauss point. The determination of the Gauss point values of stresses

multipliers
is performed by a local NewtonRaphson scheme based on a consistent
linearization of the basic equations. The algorithm includes the following steps:
trial from (36) and (39) by first
Initialize i = 0 and compute elastic trial stresses
assuming elastic behaviour and setting:

(0) = 0

(0)

(44)
trial

(45)

Check flow rule and determine set of active (plastic) Gauss points:


(i)
(i)
f
= f
n
o


(i)
Jact = {1, 2}| f (i)  > 0

(46)

= 0 where f () < 0

(48)

(47)

It has to be recognized that the set of active (plastic) Gauss points within a time
interval can change during the iteration. Therefore, an algorithm has to be applied
which is usually used in the context of non smooth multisurface plasticity. For all
Gauss points satisfying the yield condition is set to zero:

The further steps are:

Calculate residual in stressspace:
G(i)
trial
(i) C
R(i) =
(i)

(49)

Check for convergence at all active Gauss points and finish iteration when convergence is archieved:


(i)
(50)
f (i)  < tol1 for all Jact
||R(i) || < tol2

Calculate consistent tangent moduli:

A(i) = I + C

2 g (i) 
2 

10

(51)

|
1

2
(i)

g
|


2
2

(52)

 (i) and d
Calculate increments d
(i) :
h
i1  (i)

G(i)
d
(i) = F (i) A(i) 1 C
f + F (i) A(i) 1 R(i)


(i) 1
(i)
(i) (i)
(i)

d
= A
R C G d

where

F (i) =

f (i) 


f (i) 


(53)
(54)

(55)

Update variables and perform next iteration step ( i i + 1):

(i+1) =
(i) + d
(i)

(56)

(i+1) =
(i) + d
(i)

(57)

After the iteration has converged at every spatial sampling point, the vector of internal
forces can be assembled by equation (27) for each time t.

4.2

Consistent tangent operator

With known internal forces the next step in the global iteration can be performed. However, to archieve a quadratic rate of asymptotic convergence that is characteristic for the
NewtonRaphson method the stiffness matrix K has to be calculated by use of the
consistent tangent operator C ct which can be obtained as a byproduct in the last step of
the iteration:

C ct =
.
(58)

The exact linearization of the given algorithm yields the following formulation:


C
G B 1 F A1 C
,
C ct = A1 C
where

G.
B = F A1 C

(59)

(60)

For use in the global iteration procedure with the system of equations (21) a transformation has to be performed from the Gauss points in time to the time interval limits tn and
tn+1 . Together with equation (36) this leads to:

ct
ct
1
1

(1
+
3)
(1

3)
3
1
2
2
ct
C 11 C 12
ct =
C
(61)

C
=
ct
ct

1
1
C
C
(1 3) 2 (1 + 3)
3 1
21
22
2
Finally the stiffness matrices K can be calculated in accordance to equation (23):
nel

K =

e=1

ct D u N u d e ,
(D u N u )T C

11

(62)

where the time interval size t has been cancelled out. It can be seen that the influence
of both Gauss points is strongly coupled and in the general case the stiffness matrix K
is nonsymmetric. In the special case of an elastic behaviour at both temporal Gauss
points the coupling terms become zero and the consistent tangent operator reduces to a
blockdiagonal form:

ct =
C

(63)

The stiffness matrices K 11 and K 22 then take the standard form of linear finiteelement
calculations and K 12 as well as K 21 reduce to zero. This shows that the proposed formulation is a consistent extension to elastoplasticity which of course also includes the
linearelastic case.

5
5.1

PSfrag replacements

Numerical Examples
One-dimensional consolidation

The first example deals with the application of the proposed method in a onedimensional
consolidation problem. A soil layer of 10m height is compressed by a ramptype loading
process. The system characteristics and material parameters are given in figure 5. A
Coulombtype yield surface with a smooth shape in the plane and a smooth tension
cutoff together with an nonassociated flow rule was used to descibe the elastoplastic
behaviour of the soil.
P
E = 15000 kN
m2
= 0.25

u(t)
P

loadingfunction

elementtype

ky

ky = 1.0 1010 ms
= 1.5 103

kg
ms

H = 10 m
P = 20 kN
m2

Figure 5: System for onedimensional consolidation analysis

As a result the surface displacements and the porewater pressures at the bottom are
plotted in a logarithmic time scale (figures 6 and 7). An exact solution can be given in
the elastic case which is also shown in figure 6
Two calculations with time steps of different size are performed. The first with very small
time intervals shall serve as a reference solution. Within the second calculation larger
time steps are used, to show the effect of jumps that arise from the timediscontinuous
formulation. It can be seen that a significant increase in the resulting jumps is associated
with larger time steps is when the evolution of the curve is not linear. Therefore they
could be used as an errorindicator in the time domain. A second advantage of the time
discontinuous method lies in the damping of temporal oscillations. This gets obvious in
12

surface displacements u(t) [cm]

-0.5
-1

small time steps

elastic
discontinuous

-1.5
-2

-2.5
-3

-3.5
-4

10

100

1000

time [h]

Figure 6: Displacements

excess pore pressure p(t) [kPa]

0
-5
-10
small time steps
large time steps
-15
-20
-25

10

100

time [h]

1000

Figure 7: Pore water pressures

13

0.002

excess pore pressure p(t) [kPa]

0
-0.002
-0.004
continuous
discontinuous
-0.006
-0.008
-0.01
time [h]

Figure 8: Oscillation at the and of the consolidation process

the later stage of consolidation when the porewater pressures should have been reduced to
zero. Figure 8 demonstrates the difference to a calculation with continuous timestepping
where the trapezoidal rule is applied. The latter shows oscillations which are damped out
only very slowly.

5.2

Two-dimensional localization problem

This example shows the applicability of the timediscontinuous method to a more complicated strain localization problem in the analysis of twodimensional consolidation. A
stripfooting on halfspace is again loaded by the same type of loading function. The
geometric characteristics are shown in figure 9. The material parameters and the elasto
plastic properties remain unchanged to the onedimensional problem.
5m
p(t)

rigid, undrained

drained

loading function
p(t)
undrained

undrained

8m

rough, undrained

10m

Figure 9: Twodimensional consolidation

14

Figure 10: Strain localization in shearbands

Figure 10 shows the accumulation of plastic strain increments in distinct shearbands that
start from the edge of the foundation and procede in the diagonal direction throughout
the halfspace.
The introduction of a time dependency by considering the viscous fluid flow represents
a regularization of the otherwise illposed problem. Therefore, the width of the shear
band is now is influenced by the permeability of the soil. Future investigations will have to
demonstrate the correlation between the soil permeability and the calculated bandwidth.

5.3

Convergence properties of the proposed integration scheme

The determination of the state of stress at the end of the time interval is accomplished
by the integration of the rateequations of plasticity. The resulting nonlinear problem is
solved by a local NewtonRaphson iterative scheme that enforces the yield condition
at both Gauss points in time. Because the trial elastic predictor does not necessarily
lie inside of the elastic regime, several iteration passes have to be performed to assure
f () 0. The following table shows the progress within the iterative process:
Iteration

Gauss point 1 (t = 1 )
||Res||

trialstate

Gauss point 2 (t = 2 )
||Res||

1.6861

f
4.2688

5.3593 102

1.5885 102

3.7441 101

8.9990 102

3.9152 104

1.5821 106

6.4567 103

6.9415 105

3.0690 105

7.5648 1011

1.7868 105

1.2914 109

3.4748 107

< 1015

2.2194 107

< 1015

Table 1: Iteration progress at both Gauss points

It can be seen that the state of stress at both active Gauss points is iteratively returned
to the yieldsurface f () = 0. The subsequent values of the yieldfunction demonstrate
15

the quadratic convergence rate of the applied NewtonRaphson scheme until the convergence criteria (50) and (51) are satisfied.

Conclusions

In this paper a complete formulation of the timediscontinuous Galerkin procedure for

the analysis of saturated soil is developed. The use of spacetime finite elements combines the discretization in space and time to a consistent approach. The extension to
elastoplastic problems has been presented. Within a global NewtonRaphson iterative scheme the importance of the formulation of a consistent algorithmic tangent operator
has been pointed out. Since the standard Eulerbackward method for the integration
of the rateequations of plasticity is incompatible with the underlying formulation, a new
method is presented using variational forms for the kinematic equations and the elastic
plastic stressstrain relation. The restriction by the yieldcondition is satisfied at Gauss
points in time and in space. Within the presented algorithm this is done by an implicit
local iteration scheme.
The numerical examples demonstrate the applicability in one- and twodimensional consolidation analysis. For further improvements the reliability and accuracy of the finite
element calculations it may be necessary to apply adaptive procedures. On the basis of
the presented formulation it is possible on one hand to use standard errorestimators like
residualtype or postprocessing methods in the spatial domain. On the other hand the
presented method enables the control of the time discretization error by the jumps at the
discrete time levels which will be the object of further investigations.

References
[1] H. Cramer, R. Findei, G. Steinl, and W. Wunderlich. An adaptive finite element approach in associated and nonassociated plasticity considering localization phenomena. In P. Ladeveze and J.T. Oden, editors, Advances in Adaptive Computational
Methods in Mechanics. Elsevier Science, Amsterdam, 1998.
[2] R. de Boer and W. Ehlers. Theorie der Mehrkomponentenkontinua mit Anwendung auf
bodenmechanische Probleme. Forschungsbericht aus dem Fachbereich Bauwesen 40,
UniversitatGHEssen, 1986.
[3] S. Diebels, P. Ellsiepen, and W. Ehlers. A twophase model for viscoplastic geomaterials. Bericht aus dem Institut f
ur Mechanik (Bauwesen) 96-II-6, Universitat Stuttgart,
1996.
[4] W. Ehlers. Porose Medien ein kontinuumsmechanisches Modell auf der Basis der
Mischungstheorie. Forschungsbericht aus dem Fachbereich Bauwesen 47, Universitat
GHEssen, 1989.
[5] G. Hulbert and T.J.R. Hughes. Spacetime finite element methods for second order
hyperbolic equations. Comput. Meth. Appl. Mech. Eng., 84:327348, 1990.

16

[6] J.C. Simo and T.J.R. Hughes. Computational Inelasticity. Springer Verlag, New York,
1998.
[7] E. Stein and R. Mahnken. Adaptive timestep control in creep analysis. Int. J. Numer.
Meth. Engng., 28:16191633, 1989.
[8] C. Truesdell and R.A. Toupin. The Classical Field Theories. In S. Fl
ugge, editor,
Handbuch der Physik, volume III/1, pages 226902. Springer Verlag, Berlin, 1960.
[9] N.-E. Wiberg and X.D. Li. Structural dynamic analysis by a timediscontinuous
Galerkin finite element method. Int. J. Numer. Meth. Engng., 39:21312152, 1996.

17