Jordan Journal of Civil Engineering, Volume 5, No.

2, 2011

Viscoplastic Finite Element Analysis of Complex Geotechnical Problems

Waddah Salman Abdullah

Associate Professor of Civil Engineering, Department of Civil Engineering, Jordan University of Science and
Technology, Irbid, P.O. Box 3030, Jordan. E-Mail: waddah@just.edu.jo

ABSTRACT
The principles of elasto/viscoplastic finite element analysis were presented and explained, and the results and
the analysis of some geotechnical problems were presented in this work. Plasticity models such as von-Mises,
Tresca, Drucker-Prager and Mohr-Coulomb models with associated and non-associated flow rules were
incorporated in the viscoplastic algorithm. The ultimate bearing capacity of a rigid surface footing on
weightless clayey soil, Tresca material, predicted by the elasto/viscoplasticity approach agrees very well with
that obtained by Prandtl exact solution (only 1% above Prandtl exact solution). Solution was also presented
for complex problems with no available solution such as the problem of an anchor buried in sands (Mohr-
Coulomb materials), where large zones within the soil domain are dominated by tensile stresses as well as
sharp changes in shear stresses.

KEYWORDS: Viscoplasticity algorithm, Plasticity models, Associated flow rule, Non-associated

flow rule, Foundations, Anchors.

INTRODUCTION deficiency in simulating and predicting soil

deformations as well as stresses developed within soil
Soils are highly non-linear materials with domain due to applied loads and stresses. Plasticity
anisotropic and inhomogeneous nature. Therefore, exact modeling, however, represented an enormous step in
solutions for problems with such complex nature are not pursuit of accurate prediction of soil behavior even
achievable. Early developments were to incorporate the under very complex stress systems such as tensile
non-linear elastic models such as the hyperbolic model stresses and zones of very sharp changes in shear
in the finite element method (Kondner, 1963; Duncan stresses. Early works on plasticity modeling with finite
and Chang, 1970). Other types of nonlinear functions element method was mainly applied to metal stress
such as the spline function were also incorporated in analysis (Zienkiewicz and Nayak, 1971; Nayak and
finite element analysis (Desai, 1971). Non-linear Zienkiewicz, 1972). Implementation of viscoplasticity
elasticity provided solutions, but with limited success. with finite elements (Zienkiewicz and Cormeau, 1972,
Nonlinear elasticity could not model or simulate soil 1974; Cormeau, 1974) was intended for solving general
dilation on shearing which is an important phenomenon problems in solid mechanics. Applications of
for cohesionless soils. Accordingly, analysis viscoplasticity with finite elements using diversity of
incorporating such models experience immense plasticity models were advanced, specifically for
solving complex problems in geotechnical engineering
(Abdullah, 1983, 1987a, 1987b, 2008). The various
Accepted for Publication on 15/4/2011.
plasticity models used were; von-Mises, Tresca,

- 302 - © 2011 JUST. All Rights Reserved.

Jordan Journal of Civil Engineering, Volume 5, No. 2, 2011

Drucker-Prager, Mohr-Coulomb, Rowe’s stress- the viscoplastic algorithm coupled with plasticity
dilatancy model, and the critical state model, with models in the finite element method.
associated and non-associated flow rules. 2. To demonstrate the powerful nature of the
elasto/viscoplastic approach for precise analysis of
Scope of the Work complex geotechnical engineering problems.
1. To illustrate the principles and implementation of

σ3
σ1 = σ2 = σ3
Mohr-Coulomb

Tresca

σ2

σ1
Figure 1: Mohr-Coulomb and Tresca plastic yield surfaces

PLASTICITY MODELS yielding is represented by a surface called the plastic

yield surface, which separates the elastic state from the
The use of plasticity is intended to model the plastic state. There are several types of plastic yield
inelastic behavior of soils. This includes the detection of surfaces. The shape of each of these surfaces is
inelastic behavior, finding out the amount of dependent on the type of soil it represents. Broadly
irrecoverable strains and accounting for changes in speaking, the shapes of the plastic yield surfaces are
material strength after plastic yielding had occurred; divided into two categories. The first category is the
usually called hardening/softening (Nayak and hydrostatic stress independent surfaces, such as the von-
Zienkiewicz, 1972; Zienkiewicz and Nayak, 1971; Mises and Tresca surfaces which are appropriate for
Zienkiewicz and Cormeau, 1972; Zienkiewicz and modeling cohesive soils (Fig.1, Fig. 2 and Fig.3). The
Cormeau, 1974; Cormeau, 1974; Abdullah, 1983, second category is the hydrostatic stress dependent
1987a, 1987b, 2008). In a uniaxial stress system, the surfaces, such as Mohr-Coulomb and Drucker-Prager
stress that limits elastic behavior is called the yield surfaces which are appropriate for cohesionless soils or
strength. In three- dimensional stress systems, plastic mixed soils (Fig.1, Fig. 2 and Fig.3). Mathematically,

- 303 -
Viscoplastic Finite Element… Waddah Salman Abdullah

the plastic yield surface (Table 1) is represented by, where: σ represents the stress components, εp
represents the accumulated plastic strains and κ is the
(1) hardening coefficient.

σ3 Drucker-Prager σ1 = σ2 = σ3

C
Von-Mises

σ2

σ1

Figure 2: Drucker-Prager and Von-Mises plastic yield surfaces

Mohr-coulomb and Tresca plastic yield surfaces and (2)

Drucker-Prager and von-Mises plastic yield surfaces are
shown in Fig.1 and Fig.2, respectively. The π-plane of
some commonly used plastic yield functions is shown in and (3)
Fig. 3.
The conditions for F are; where: Q is the plastic potential function and may be
F < 0.0 Elastic state of stress; represented in a fashion similar to the plastic yield
F = 0.0 Plastic flow; surface (Eq. 1). d is a non-negative proportionality
F > 0.0 Not permitted. constant which may vary throughout the loading history.
The flow theory of plasticity establishes a The plastic potential may assume the same function
relationship between increments of stress and strain as that of the plastic yield function, and in this case we
rates, usually known as “normality condition”. The have what is commonly known as “associated flow
latter does not give the magnitude of plastic strain but rule”, such that;
establishes the ratio of the strain components, such that;
(4)

- 304 -
Jordan Journal of Civil Engineering, Volume 5, No. 2, 2011

σ3
π
θ0 = +
Mohr-Coulomb 6 θ0 = 0

External θ0
envelope θ0 = -

Internal
Internal cone envelope
coincident at θ0 =-

σ1
σ2

Figure 3: π-plane for some commonly used plastic yield functions

- y

εvp
Inactive
if

εe

Figure 4: One-dimensional elasto/viscoplastic model

- 305 -
Viscoplastic Finite Element… Waddah Salman Abdullah

- 306 -
Jordan Journal of Civil Engineering, Volume 5, No. 2, 2011

Figure 6: Load settlement relationship for a rigid smooth foundation on weightless soil;
Tresca yield criterion, Plane strain condition, critical time stepping (Abdullah, 1983)

Changes in the state of the yield surface whether in measure of hardening, then the process is known as
shape, size and location compared to its previous state “Strain Hardening”. It establishes a relationship
are loosely called “hardening”. Hardening is called between the accumulated plastic strain and hardening,
isotropic if and only if the size of the plastic yield which is related to parameter κ (Eq. 1), such that;
surface changes, whilst the location and shape does not
change. On the other hand, hardening is called (6)
“kinematic” if the location of the yield surface translates
along the direction of the plastic strain increment. H is always an increasing function of plastic strain.
Isotropic hardening usually involves one simple However, the rate of H may become negative as in the
parameter denoting the state of the material. If plastic case of softening. The incremental plastic stress-strain
work is used as the hardening parameter, the process is law may be given directly from the normality principle
called “Work-Hardening”. It establishes a relationship (Eq. 3), such that;
between the plastic work done by the external agency
and the hardening produced by the plastic deformation, (7)
which is eventually related to the parameter κ (Eq. 1),
such that; and

(5)
or (8)
where: Wp is the amount of plastic work per unit
volume. However, the rate of plastic work may be
positive, zero or even negative depending on the type
and state of the material. Negative rate of plastic work and
signifies softening rather than hardening.
If the accumulated plastic strain is taken as the (9)

- 307 -
Viscoplastic Finite Element… Waddah Salman Abdullah

Equation 8 may be rewritten as; Abdullah, 1983, 1987a, 1987b, 2008).

The concept of the elasto/viscoplastic modeling is
or (10) that the material behaves elastically on loading as long
as the stresses stay within the plastic yield function
where, a1 , a 2 and a 3 are the derivatives of (F<0). Beyond that, there is an additional time
σ m , σ and J 3 with respect to the components of dependent component of strains (called viscoplastic) the
stresses and such that; rate of which is dependent on the amount of excess
stress beyond the plastic yield surface. Therefore, the
total strain rate may be obtained as;

; (12)

where: [D] = elasticity matrix, = fluidity parameter,

F = plastic yield function and Q = plastic potential
function.
Also, the following means;

+ (11) for F < 0 (13a)

for F > 0 (13b)

For a uniaxial model, the viscoplastic concept is

θ0, σm, σ , Sx, Sy, Sz and J3 are given in the shown in Fig. 4. It consists of an elastic spring
Appendix. connected in parallel with a Bingham unit, which
B1, B2 and B3 are constant parameters related to a consists of a plastic unit capable of taking a stress σy
particular plastic yield function (Table 2). and a dashpot. The total stress σ may exceed the yield
value σy by an amount, which is a fact that has a
Viscoplasticity Algorithm physical justification (Perzyna, 1966) and is the
Viscoplasticity is a physically more logical approach fundamental difference between plasticity and
for simulating material inelasticity behavior than the viscoplasticity. The plastic yield surface F, the plastic
plasticity approach. Moreover, numerically the potential surface Q, the normality condition and the
viscoplastic approach eliminates some numerical hardening/softening phenomenon have the same
difficulties associated with the elasto-plastic modeling, meaning as in the context of plasticity. The function F
such as strain softening and non-associated flow rules. may assume a number of forms, linear, exponential,
The time factor involved in the viscoplastic model may logarithmic …etc. The linear form of (φ ( F ) δ = F ) is
be used to obtain real time dependent solutions, given fairly adequate to describe the behavior of many
that the parameters involved in the model are materials (Zienkiewicz and Cormeau, 1972; Abdullah,
experimentally determined. It could also be used as a 1983). The viscoplastic solution involves the integration
pure computational artifice leading to an elasto/plastic of matrix equations in time. There are a large number of
solution (Nayak and Zienkiewicz, 1972; Zienkiewicz these stepping schemes. Euler’s forward time marching
and Nayak, 1971; Zienkiewicz and Cormeau, 1972; scheme has been shown to be efficient (Cormeau, 1974;
Zienkiewicz and Cormeau, 1974; Cormeau, 1974; Abdullah, 1983). In this scheme, it is assumed that

- 308 -
Jordan Journal of Civil Engineering, Volume 5, No. 2, 2011

viscoplastic strain rate ε&ν p remains constant within may occur if the time interval length is too large. Firstly
each time interval after it has been determined from a empirical rules were introduced (Zienkiewicz, and
state of stress assumed to remain constant at the Cormeau, 1972, 1974), such as;
beginning of that time interval. Numerical instability

ε
∆t ≤ τ (14)
ε&νp
and
ε
∆t ≤ τ (15)
ε& νp

where;

(16a)

(16b)

(16c)

(16d)

It was found from numerical experiments that τ Mohr-Coulomb

values ranging from 0.01 to 0.15 gave stable solutions.
Experiments have shown that none of the empirical (17c)
rules was satisfactory as far as stability and economy of
the solution were concerned. For associated flow rules, where: E = modulus of elasticity, ν = Poisson’s ratio,
and based on the viscous properties of the material, a φ = friction angle and γ = fluidity parameter.
critical time interval length was obtained (Zienkiewicz
and Cormeau, 1974). For von-Mises, Tresca and Mohr- VERIFICATION OF
Coulomb plastic yield functions, the critical time THE VISCOPLASTIC APPROACH
interval length is given as;
In order to demonstrate the powerful nature of the
Von-Mises (17a) viscoplastic approach, a rigid strip surface footing
resting on weightless clayey soil was considered.
Tresca (17b) Prandtl (1920 and 1921) obtained an exact analytic
solution for a rigid strip surface footing resting on

- 309 -
Viscoplastic Finite Element… Waddah Salman Abdullah

weightless clayey soil. Prandtl (1920 and 1921) found and the ultimate bearing capacity qu = (2+π) C or qu/C =
that yielding occurs at a bearing pressure equal to πC, (2+π); (C is the soil cohesion).

Foundation

Logarithmic
spiral
Figure 7: Failure mechanism for surface footing elasto-viscoplastic approach and
Prandtl solution (Abdullah, 1983)

Pulling-up force

Mohr-Coulomb material H Rigid rough

E/γH = 333 side
H/B = 3
5B ν = 0.3
B
φ = 30o
γH/C = 120 Rigid anchor

Area dominated by
Rigid rough tensile and very sharp
base changes in shear stresses

Figure 8: Vertically pulled rigid anchor

For the elasto/viscoplastic finite element solution, result (Fig. 6) as compared with Prandtl exact solution
the soil domain was idealized with a carefully designed (only 1% above Prandtl exact solution). The load
finite element mesh (Fig. 5). The mesh provided small- settlement relationship (Fig. 6) for the rigid strip surface
sized elements at areas of load discontinuity and very footing resting on weightless clayey soil demonstrates
sharp shear stress changes in order to account for such that initial yielding started at q/C = π as predicted by
sharp changes. The finite elements used were 8-noded Prandtl exact solution (Abdullah, 1983). The plastic
isoparametric elements, and the material was modeled zone or the failure mechanism predicted by Prandtl
as Tresca type material (Abdullah, 1983). The ultimate exact method and the plastic zone predicted by the
bearing capacity predicted by the viscoplastic finite viscoplastic finite element solution (Fig. 7) are both
element approach (Abdullah, 1983) yielded an excellent identical (Abdullah, 1983).

- 310 -
Jordan Journal of Civil Engineering, Volume 5, No. 2, 2011

Figure 9: Finite element mesh for the anchor problem (Abdullah, 1983)

Figure 10: Load displacement relationships for the anchor problem (Abdullah, 1983)

- 311 -
Viscoplastic Finite Element… Waddah Salman Abdullah

Mohr-Coulomb material Zone dominated by tensile

E/γH = 333 stresses and very sharp
H/B = 3 changes in shear stresses
ν = 0.3
Rigid rough
φ = 30o
side
γH/C = 120

Rigid anchor

Rigid rough
base

Figure 11: Vector displacement field for the anchor problem (Abdullah, 1983)

Table 1: Common plasticity models

Type of plastic
Plastic yield function
model
Von-Mises
Tresca

Drucker-Prager

Mohr-Coulomb

Table 2: Parameters for some common plasticity models

Plastic
B1 B2 B3
model
Von-
0 0
Mises

Tresca 0

Drucker-
1.0 0
Prager

Mohr-
Coulomb

- 312 -
Jordan Journal of Civil Engineering, Volume 5, No. 2, 2011

Undoubtedly, this example demonstrated that the solutions with a great deal of accuracy. Viscoplasticity
elasto/viscoplastic finite element solution yielded very approach is more realistic than plasticity in terms of
accurate results as far as initial yielding, ultimate reproducing actual material behavior and yields highly
bearing capacity and failure mechanism were accurate results even for very complex problems. In this
concerned. Therefore, for any type of geotechnical work, the viscoplasticity modeling and algorithm were
problem where no exact solution is available, the incorporated in the finite element method (8-noded
elasto/viscoplastic finite element approach may be used isoparametric elements). Various plasticity models were
to obtain a high precision solution. incorporated. These were; von-Mises, Tresca, Drucker-
Prager, Mohr-Coulomb, Rowe’s stress-dilatancy model
ANCHOR PROBLEM and the critical state model, with associated and non-
associated flow rules. Evaluation of the performance of
The second problem considered in this work was the the viscoplasticity analysis was conducted on a problem
anchor problem where, due to its complex nature, there with known exact solution. The considered problem was
is no available exact solution. In addition to the plastic a rigid surface footing resting on a cohesive weightless
nature of soils, the anchor problem produces large areas soil (Tresca type material, φ = 0). Prandtl provided exact
dominated by tensile stresses and very sharp changes in solution for the mentioned problem with the ultimate
shear stresses (Fig. 8) which render the problem as bearing capacity qu = (2+π) C. The elasto/viscoplastic
insolvable by any analytical method. The soil domain finite element solution was in very close agreement with
was idealized by an appropriate finite element mesh Prandtl exact solution (only 1% above the exact
(Fig. 9), that takes into account the very sharp changes solution). This kind of result represents significant
in shear stresses as designated by the problem development for finding highly accurate results for very
description (Fig. 8). The finite elements used were 8- complex geotechnical problems which have no known
noded is oparametric elements and the material was exact solutions. The anchor problem represents one such
modeled by Mohr-Coulomb yield criterion (Abdullah, complex problem with zones dominated by tensile
1983). The normalized force versus normalized vertical stresses as well as very sharp changes in shear stresses.
displacement for an associated flow rule and a non-
associated flow rule (ψ = 15°, and ψ = 0°) differs, fairly, APPENDIX
in behavior up till near failure, getting closer at later
stages of loading, where the effect of dilation due to = Lode angle, and is given as;
associated flow rule diminishes more and more at later
stages(Fig. 10). The displacement field for the anchor with
problem shows a large amount of movements around
the anchor as well as at the soil surface right above the
anchor (Fig. 11). (A1)

SUMMARY AND CONCLUSIONS

= (A2)
Soils are highly nonlinear materials, therefore,
assuming soils as linear elastic materials produces a
gross amount of error. Non-linear elasticity was
introduced, but only with limited success. Plasticity = (A3)
models incorporated in finite element analysis yielded

- 313 -
Viscoplastic Finite Element… Waddah Salman Abdullah

2 2 2
S x =σ x −σ m S y = σ y −σ m S z = σ z − σ m and J 3 = S x S y S z + 2τ xy τ yz τ zx − S x τ yz − S y τ xz − S z τ xy
(A4)

REFERENCES and Foundations Division, ASCE, 96 (SM5).

Kondener, R.L. 1963. Hyperbolic stress-strain response:
Abdullah, W.S. 2008. New elastoplastic method for cohesive soils, Journal of the Soil Mechanics and
calculating the contact pressure distribution under rigid Foundations Division, ASCE, 89 (SM10): 115-143.
foundations, Jordan Journal of Civil Engineering, 2 Nayak, G.C. and Zienkiewicz, O.C. 1972. Elasto-plastic
(1), January. stress analysis. A generalization for various constitutive
Abdullah, W.S. 1987a. Shakedown analysis of footings on relations including strain softening, International
elastoplastic continua, Proceedings of the Second Arab Journal for Numerical Methods in Engineering.
Conference on Structural Engineering, Jordan Perzyna, P. 1966. Fundamental problems in visco-
University, Department of Civil Engineering, pp.32/1- plasticity. Recent advances in applied mechanics,
32/25, 19-22 April, Amman. Academic Press, 9, New York.
Abdullah, W.S. 1987b. Vertically pulled anchors buries in Prandtl, L. 1920. Über die Härte plastischer Körper, Nachr.
no-tension materials, Proceedings of the Second Arab K. Ges. Wiss. Gött., Math.-Phys., Kl, 74-85.
Conference on Structural Engineering, Jordan Prandtl, L. 1921. Über die Eindringungsfestigkeit (Härte)
University, Department of Civil Engineering, pp.29/1- plastischer Baustoffe und die Festigkeit von Schneiden,
29/17, 19-22 April, Amman. Z. Angew. Math. Mech., 15-20.
Abdullah, W.S. 1983. Non-linear finite element analysis of Zienkiewicz, O.C. and Nayak, G.C. 1971. A general
shallow footings, Ph.D. thesis submitted to the approach to problems of large deformation and
University of Wales, School of Engineering, Division plasticity using isoparametric elements, 3rd Conference
of Civil Engineering, University College of Swansea, on Matrix Methods, Wright-Patterson Air Force Base,
Wales. Ohio.
Cormeau, I.C. 1974. Viscoplasticity and plasticity in the Zienkiewicz, O.C. and Cormeau, I.C. 1972. Viscoplasticity
finite element method, Ph.D. Thesis, University solution by finite element process, Archives of
College of Swansea, Wales. Mechanics, 24 (5-6): 873-889, Warszawa.
Desai, C.S. 1971. Nonlinear analysis using spline Zienkiewicz, O.C. and Cormeau, I.C. 1974.
functions, Journal of the Soil Mechanics and Viscoplasticity-plasticity and creep in elstic solids – a
Foundations Division, ASCE, 97 (SM10): 1461-1480. unified numerical solution approach, International
Duncan, J.M. and Chang, C.Y. 1970. Nonlinear analysis of Journal for Numerical Methods in Engineering, 8: 821-
stress and strain in soils, Journal of the Soil Mechanics 845.

- 314 -