Geomechanics Module

User´s Guide

VERSION 4.4
Geomechanics Module User’s Guide
© 1998–2013 COMSOL
Protected by U.S. Patents 7,519,518; 7,596,474; 7,623,991; and 8,457,932. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement ([Link]/sla) and may be used or copied only under the terms of the license
agreement.
COMSOL, COMSOL Multiphysics, Capture the Concept, COMSOL Desktop, and LiveLink are either
registered trademarks or trademarks of COMSOL AB. All other trademarks are the property of their
respective owners, and COMSOL AB and its subsidiaries and products are not affiliated with, endorsed by,
sponsored by, or supported by those trademark owners. For a list of such trademark owners, see
[Link]/tm.
Version: November 2013 COMSOL 4.4

Contact Information
Visit the Contact COMSOL page at [Link]/contact to submit general
inquiries, contact Technical Support, or search for an address and phone number. You can
also visit the Worldwide Sales Offices page at [Link]/contact/offices for
address and contact information.

If you need to contact Support, an online request form is located at the COMSOL Access
page at [Link]/support/case.

Other useful links include:

• Support Center: [Link]/support

• Product Download: [Link]/support/download
• Product Updates: [Link]/support/updates
• COMSOL Community: [Link]/community
• Events: [Link]/events
• COMSOL Video Center: [Link]/video
• Support Knowledge Base: [Link]/support/knowledgebase

Part number: CM021801

C o n t e n t s

Chapter 1: Introduction

Geomechanics Module Overview 6

About the Geomechanics Module . . . . . . . . . . . . . . . . 6
Where Do I Access the Documentation and Model Libraries? . . . . . . 6

Chapter 2: Geomechanics Theory

General Geomechanics Theory 12

Sign Conventions for Stress and Strain Analysis . . . . . . . . . . . 12
Defining the Stress Invariants . . . . . . . . . . . . . . . . . . 13
Defining the Yield Surface . . . . . . . . . . . . . . . . . . . 15

Defining Perfectly Plastic Materials 17

von Mises Criterion . . . . . . . . . . . . . . . . . . . . . . 17
Tresca Criterion . . . . . . . . . . . . . . . . . . . . . . . 18

Plasticity Models for Soils 20

Mohr-Coulomb Criterion . . . . . . . . . . . . . . . . . . . 20
Drucker-Prager Criterion . . . . . . . . . . . . . . . . . . . 22
Matsuoka-Nakai Criterion . . . . . . . . . . . . . . . . . . . 25
Lade-Duncan Criterion . . . . . . . . . . . . . . . . . . . . 26
Tension Cut-Off . . . . . . . . . . . . . . . . . . . . . . . 27

Theory for the Cam-Clay Material 28

About the Cam-clay Material . . . . . . . . . . . . . . . . . . 28
Volumetric Elastic Behavior . . . . . . . . . . . . . . . . . . . 30
Hardening and Softening . . . . . . . . . . . . . . . . . . . . 31
Including Pore Fluid Pressure . . . . . . . . . . . . . . . . . . 32

Failure Criteria for Concrete, Rocks, and Other Brittle

CONTENTS |3
 Materials 33
Bresler-Pister Criterion . . . . . . . . . . . . . . . . . . . . 33
Willam-Warnke Criterion . . . . . . . . . . . . . . . . . . . 34
Ottosen Criterion . . . . . . . . . . . . . . . . . . . . . . 35
Original Hoek-Brown Criterion . . . . . . . . . . . . . . . . . 37
Generalized Hoek-Brown Criterion. . . . . . . . . . . . . . . . 38

Elastoplastic Material Theory 40

Introduction to Small and Large Plastic Strains . . . . . . . . . . . . 40
Plastic Flow for Small Strains . . . . . . . . . . . . . . . . . . 42
Isotropic Plasticity . . . . . . . . . . . . . . . . . . . . . . 43
Yield Function . . . . . . . . . . . . . . . . . . . . . . . . 45
Hardening Models . . . . . . . . . . . . . . . . . . . . . . 46
Plastic Flow for Large Strains . . . . . . . . . . . . . . . . . . 49
Numerical Solution of the Elastoplastic Conditions . . . . . . . . . . 50
References for Elastoplastic Materials . . . . . . . . . . . . . . . 51

Creep and Viscoplasticity 52

About Creep . . . . . . . . . . . . . . . . . . . . . . . . 52
Fundamental Creep Material Models . . . . . . . . . . . . . . . 54
Solver Settings for Creep . . . . . . . . . . . . . . . . . . . . 57

References for the Geomechanics Module 58

Chapter 3: The Geomechanics Materials

Working with the Geomechanics Materials 62

Adding a Material to a Solid Mechanics Interface . . . . . . . . . . . 62
Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Soil Plasticity . . . . . . . . . . . . . . . . . . . . . . . . 65
Concrete . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Rocks . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Cam-Clay Material . . . . . . . . . . . . . . . . . . . . . . 69
Creep . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4 | CONTENTS
 1

Introduction

Welcome to the Geomechanics Module User’s Guide. This chapter provides a

Geomechanics Module Overview and the type of modeling you can achieve with
its extensive set of material models for geomechanics and soil mechanics. In
addition, it includes some general information about documentation and models.

The Geomechanics Module extends the capabilities of the Structural Mechanics

Module. For general information about structural analysis, see the Structural
Mechanics Module User’s Guide.

5
 Geomechanics Module Overview
In this section:

• About the Geomechanics Module

• Where Do I Access the Documentation and Model Libraries?

About the Geomechanics Module

The Geomechanics Module is an optional package that extends the Structural
Mechanics Module to the quantitative investigation of geotechnical processes. It is
designed for researchers, engineers, developers, teachers, and students, and suits both
single-physics and interdisciplinary studies within geomechanics and soil mechanics.

The module includes an extensive set of fundamental material models, such as the
Drucker-Prager and Mohr-Coulomb criteria and the Cam-Clay model in soil
mechanics. These material models can also couple to any new equations created, and
to interfaces (for example, heat transfer, fluid flow, and solute transport in porous
media) already built into COMSOL Multiphysics and its other specialized modules.

The Physics Interfaces and Building a COMSOL Model in the COMSOL

Multiphysics Reference Manual

NOTE ABOUT MATERIALS

The material property groups (including all associated properties) can be

added to models from the Material page. See Materials in the COMSOL
Multiphysics Reference Manual.

Where Do I Access the Documentation and Model Libraries?

A number of Internet resources provide more information about COMSOL, including
licensing and technical information. The electronic documentation, topic-based (or

6 | CHAPTER 1: INTRODUCTION
context-based) help, and the Model Libraries are all accessed through the COMSOL
Desktop.

If you are reading the documentation as a PDF file on your computer, the
blue links do not work to open a model or content referenced in a
different guide. However, if you are using the Help system in COMSOL
Multiphysics, these links work to other modules (as long as you have a
license), model examples, and documentation sets.

THE DOCUMENTATION AND ONLINE HELP

The COMSOL Multiphysics Reference Manual describes all core physics interfaces
and functionality included with the COMSOL Multiphysics license. This book also has
instructions about how to use COMSOL and how to access the electronic
Documentation and Help content.

Opening Topic-Based Help

The Help window is useful as it is connected to many of the features on the GUI. To
learn more about a node in the Model Builder, or a window on the Desktop, click to
highlight a node or window, then press F1 to open the Help window, which then
displays information about that feature (or click a node in the Model Builder followed
by the Help button ( ). This is called topic-based (or context) help.

To open the Help window:

• In the Model Builder, click a node or window and then press F1.
• On any toolbar (for example, Home or Geometry), hover the mouse over
a button (for example, Browse Materials or Build All) and then press F1.
• From the File menu, click Help ( ).
• In the upper-right part of the COMSOL Desktop, click the ( )
button.

To open the Help window:

• In the Model Builder, click a node or window and then press F1.
• On the main toolbar, click the Help ( ) button.
• From the main menu, select Help>Help.

GEOMECHANICS MODULE OVERVIEW | 7

 Opening the Documentation Window

To open the Documentation window:

• Press Ctrl+F1.
• From the File menu select Help>Documentation ( ).

To open the Documentation window:

• Press Ctrl+F1.
• On the main toolbar, click the Documentation ( ) button.
• From the main menu, select Help>Documentation.

THE MODEL LIBRARIES WINDOW

Each model includes documentation that has the theoretical background and
step-by-step instructions to create the model. The models are available in COMSOL
as MPH-files that you can open for further investigation. You can use the step-by-step
instructions and the actual models as a template for your own modeling and
applications. In most models, SI units are used to describe the relevant properties,
parameters, and dimensions in most examples, but other unit systems are available.

Once the Model Libraries window is opened, you can search by model name or browse
under a module folder name. Click to highlight any model of interest and a summary
of the model and its properties is displayed, including options to open the model or a
PDF document.

The Model Libraries Window in the COMSOL Multiphysics Reference

Manual.

8 | CHAPTER 1: INTRODUCTION
Opening the Model Libraries Window
To open the Model Libraries window ( ):

• From the Home ribbon, click ( ) Model Libraries.

• From the File menu select Model Libraries.

To include the latest versions of model examples, from the File>Help

menu, select ( ) Update COMSOL Model Library.

• On the main toolbar, click the Model Libraries button.

• From the main menu, select Windows>Model Libraries.

To include the latest versions of model examples, from the Help menu
select ( ) Update COMSOL Model Library.

CONTACTING COMSOL BY EMAIL

For general product information, contact COMSOL at info@[Link].

To receive technical support from COMSOL for the COMSOL products, please
contact your local COMSOL representative or send your questions to
support@[Link]. An automatic notification and case number is sent to you by
email.

COMSOL WEBSITES

COMSOL website [Link]

Contact COMSOL [Link]/contact
Support Center [Link]/support
Product Download [Link]/support/download
Product Updates [Link]/support/updates
COMSOL Community [Link]/community
Events [Link]/events
COMSOL Video Gallery [Link]/video
Support Knowledge Base [Link]/support/knowledgebase

GEOMECHANICS MODULE OVERVIEW | 9

10 | CHAPTER 1: INTRODUCTION
 2

Geomechanics Theory

The Geomechanics Module contains new materials which are used in combination
with the Structural Mechanics Module to account for the plasticity of soils and
failure criteria in rocks, concrete and other brittle materials.

In this chapter:

• General Geomechanics Theory

• Defining Perfectly Plastic Materials
• Plasticity Models for Soils
• Theory for the Cam-Clay Material
• Failure Criteria for Concrete, Rocks, and Other Brittle Materials
• Elastoplastic Material Theory
• Creep and Viscoplasticity
• References for the Geomechanics Module

11
 General Geomechanics Theory
In this section:

• Sign Conventions for Stress and Strain Analysis

• Defining the Stress Invariants
• Defining the Yield Surface

Theory for the Solid Mechanics Interface in the Structural Mechanics

Module User’s Guide

Sign Conventions for Stress and Strain Analysis

The theory of plasticity was initially developed in continuum mechanics where tensile
stresses are usually considered to be positive quantities and compressive stresses are
negative quantities. Under this convention, force and displacement components are
considered positive if directed in the positive directions of the coordinate axes. Tensile
normal strains and tensile normal stresses are treated as positive (Ref. 11).

Engineering analysis and design for soil and rock structures, however, are in most cases
concerned with compressive stresses (Ref. 1). Therefore, in geotechnical applications
the opposite sign convention is usually adopted because compressive normal stresses
are more common than tensile ones (Ref. 5).

The Geomechanics Module follows the sign convention of continuum mechanics as

stated in Ref. 1 (positive tensile stresses). This is consistent with the notation used in
the Structural Mechanics Module. This means the principal stresses often have a
negative sign, and are sorted as 123.

The convention used in Ref. 1 refers to the hydrostatic pressure (trace of the stress
Cauchy tensor) with a positive sign. The use of the first invariant of Cauchy stress
tensor I1() is preferred through this document, in order to avoid misunderstandings
with the convention in the Structural Mechanics Module (where pressure is positive
under compression, or equivalently, it has the opposite sign of the Cauchy stress
tensor’s trace).

12 | CHAPTER 2: GEOMECHANICS THEORY

Defining the Stress Invariants
The starting point for defining the theory of plastic deformation of soils is the
definition of the first, second, and third invariants for the Cauchy stress tensor . These
invariants are defined as follows:

I 1    = trace   
1 2
I 2    = ---  I 1 – :  (2-1)
2
I 3    = det   

The first invariant I1 is the trace of the tensor, the second invariant I2 is the sum of the
principal two-rowed minors of the determinant of , and the third invariant I3 is the
determinant of . (Ref. 1).

When 1, 2, and 3, represent the principal components of the stress tensor, these
invariants can be written as

I1    = 1 + 2 + 3
I2    = 1 2 + 2 3 + 1 3
I3    = 1  2  3

The principal components of the stress tensor are the roots of the characteristic
equation (Cayley–Hamilton theorem)

3 2
 – I1  + I2  – I3 = 0

which comes from calculating the roots for the equation

det  – I = 0

The invariants I1, I2, and I3 can be called in user-defined yield criteria by
referencing the corresponding variables solid.I1s, solid.I2s, and
solid.I3s.

Defining the deviatoric stress as the traceless tensor

1
dev    =  – --- I 1 I
3

GENERAL GEOMECHANICS THEORY | 13

 introduces us to the first, second, and third deviatoric stress invariants

J 1    = trace  dev     = 0
1
J 2    = --- dev   :dev    (2-2)
2
J 3    = det  dev    

As defined above J20. In soil plasticity, the most relevant invariants are I1, J2, and
J3. I1 represents the effect of mean stress, J2 represents the magnitude of shear stress,
and J3 is the direction of the shear stress.

The invariants J2 and J3 can be called in user-defined yield criteria by

referencing the variables solid.II2s and solid.II3s, where solid is
the name of the interface identifier for the physics interface.

OTHER INVARIANTS
It is possible to define other invariants in terms of the primary invariants. One common
auxiliary invariant is the Lode angle 

3 3 J3
cos 3 = -----------  -----------
- (2-3)
2 J3  2
2

The Lode angle is bounded to 03 when the principal stresses are sorted as
123 (Ref. 1).

Following this convention, =corresponds to the tensile meridian, and =3

corresponds to the compressive meridian. The Lode angle is part of a cylindrical
coordinate system (the Haigh–Westergaard coordinates) with height (hydrostatic axis)

 = I 1 / 3 and radius r = 2J 2 .

The Lode angle is undefined at the hydrostatic axis, where all three
principal stresses are equal (1=2=3=I1/3) and J2=0. To avoid
division by zero, the Lode angle is computed from the inverse tangent
function atan2, instead of the inverse cosine, as stated in Equation 2-3.

14 | CHAPTER 2: GEOMECHANICS THEORY

 The Lode angle and the effective (von Mises) stress can be called in
user-defined yield criteria by referencing the variables [Link] and
[Link], where solid is the name of the interface identifier for the
physics interface.

The octahedral plane (also called -plane) is defined perpendicular to the hydrostatic
axis in the Haigh–Westergaard coordinate system. The stress normal to this plane is
oct=I1/3, and the shear stress on that plane is defined by

 oct = 2/3J 2

The functions described in Equation 2-1 and Equation 2-2 enter into expressions that
define various kind of yield and failure surfaces. A yield surface is a surface in the 3D
space of principal stresses which circumscribe an elastic state of stress.

PRINCIPAL STRESSES
The principal stresses 1, 2, and 3) are the eigenvalues of the stress tensor, and when
sorted as 123 they can be written by using the invariants I1 and J2 and the
Lode angle 03 (Ref. 1):

1 4J
 1 = --- I 1 + ---------2- cos 
3 3
4J
 2 = --- I 1 + ---------2- cos   – ------
1 2
3 3  3
4J
 3 = --- I 1 + ---------2- cos   + ------
1 2
3 3 3

here, i0 means tensile stress, and i0 means compressive stress.

The principal stresses can be called in user-defined yield criteria by

referencing the variables solid.sp1, solid.sp2 and solid.sp3.

Defining the Yield Surface

A yield criterion serves to define the stress condition under which plastic deformation
occurs. Stress paths within the yield surface result in purely recoverable deformations

GENERAL GEOMECHANICS THEORY | 15

 (elastic behavior), while paths intersecting the yield surface produces both recoverable
and permanent deformations (plastic strains).

In general, the yield surface can be described as

F = f    – fc = 0

where fc can be a constant value (for perfectly plastic materials), or a variable for
strain-hardening materials. The yield surface F is a surface in the space of principal
stresses, for which the elastic regime (F0) is enclosed.

For brittle materials, the yield surface represents a failure surface, which is a stress
level at which the material collapses instead of deforms plastically.

Some authors define the yield criterion as f ()= fc, while the yield
surface is an isosurface in the space of principal stresses F=0, which can
2 2
be chosen for numerical purposes as F = f    – f c = 0 .

16 | CHAPTER 2: GEOMECHANICS THEORY

Defining Perfectly Plastic Materials
For perfect elastoplastic materials, the yield surface is fixed in the 3D space of principal
stress, and therefore, plastic deformations occur only when the stress path moves on
the yield surface (the regime inside the yield surface is elastic, and stress paths beyond
the yield surface are not allowed).

In general, the yield criterion depends on various parameters. Most of the plasticity
models are based on isotropic assumptions, which require the yield function to be
independent of the chosen coordinate system. This introduces the concept of using the
stress invariants previously defined in Equation 2-1, Equation 2-2, and Equation 2-3.

In this section:

• von Mises Criterion

• Tresca Criterion

von Mises Criterion

The von Mises criterion suggests that the yielding of the material begins when the
second deviatoric stress invariant J2 reaches a critical value. This criterion can be
written in terms of the elements of Cauchy’s stress tensor (Ref. 1)

1 2 2 2 2 2 2
J 2 = ---    11 –  22  +   22 –  33  +   33 –  11   +  12 +  23 +  13 = k 2
6

or equivalently J2 = k .

The von Mises criterion is implemented as

F = 3J 2 –  ys = 0

where ys is the yield stress level (yield stress in uniaxial tension).

The effective or von Mises stress (  mises = 3J 2 ) is available in the

variable [Link], where solid is the name of the interface identifier
for the physics interface.

DEFINING PERFECTLY PLASTIC MATERIALS | 17

 Tresca Criterion
The Tresca yield surface is normally expressed in terms of the principal stress
components

1
--- max   1 –  2   1 –  3   2 –  3  = k
2

The Tresca criterion is a hexagonal prism with its axis equally inclined to the three
principal stress axes. When the principal stresses fulfill 123, this criterion is
written as

1
---   1 –  3  = k
2

By using the representation of principal stresses in term of the invariants J2 and the
Lode angle 03, this criterion can alternatively be written as

4J 
--- ---------2-  cos    – cos   + ------  = J 2 sin   + --- = k
1 2
2 3   3   3

or equivalently


J 2 cos   – --- = k
 6

The maximum shear stress is reached at the meridians (=0 or =). The Tresca
criterion can be circumscribed by setting the Lode angle =0, or equivalently, by a
von Mises criterion

3J 2 = 2k

The minimum shear is reached at=, so the Tresca criterion can be inscribed by
setting a von Mises criterion

J2 = k

When dealing with soils, the parameter k is also called undrained shear strength.

The Tresca effective stress, tresca=13 is implemented in the

variable [Link], where solid is the name of the interface
identifier for the physics interface.

18 | CHAPTER 2: GEOMECHANICS THEORY

 
Upper limit
Lower limit

Figure 2-1: The upper and lower limits of the Tresca criterion.

Figure 2-2: Classical yield criteria for metals. Left: Tresca criterion. Right: von Mises
criterion.

The von Mises and Tresca criteria are independent of the first stress invariant I1 and
are mainly used for the analysis of plastic deformation in metals and ductile materials,
though some researchers also use these criteria for describing fully saturated cohesive
soils (that is, clays) under undrained conditions.

DEFINING PERFECTLY PLASTIC MATERIALS | 19

 Plasticity Models for Soils
In this section:

• Mohr-Coulomb Criterion
• Drucker-Prager Criterion
• Matsuoka-Nakai Criterion
• Lade-Duncan Criterion
• Tension Cut-Off

Elastoplastic Material Theory

Mohr-Coulomb Criterion
The Mohr-Coulomb criterion is the most popular criterion in soil mechanics. It was
developed by Coulomb before the Tresca and von Mises criteria for metals, and it was
the first criterion to account for the hydrostatic pressure. The criterion states that
failure occurs when the shear stress and the normal stress acting on any element in the
material satisfy the equation

 +  tan  – c = 0

here, is the shear stress, c the cohesion, and  denotes the angle of internal friction.

With the help of Mohr’s circle, this criterion can be written as

1
---   1 –  3  + 1
---   1 +  3  sin  – c cos  = 0
2 2

The Mohr-Coulomb criterion defines an irregular hexagonal pyramid in the space of

principal stresses, which generates singularities in the derivatives of the yield function.

20 | CHAPTER 2: GEOMECHANICS THEORY

 ia n
rid
me
an
ile
i
rid
ns

me
Te

si ve
p res
m
Co

Figure 2-3: The Mohr-Coulomb criterion.

The Mohr-Coulomb criterion can be written in terms of the invariants I1 and J2 and
the Lode angle 03 (Ref. 1, Ref. 9) when the principal stresses are sorted as
123. The yield function then reads

J
F y = --- I 1 sin  + -----2-   1 + sin   cos  –  1 – sin   cos   + ------  – c cos  = 0
1 2
3 3 3

The tensile meridian is defined when =and the compressive meridian when =.

Rearranging terms, the Mohr-Coulomb criterion reads

Fy = J 2 m    + I 1 – k = 0

where

m    = cos   –   6  – 1  3 sin  sin   –   6  ,  = sin   3 , and k = c cos 

In the special case of frictionless material, (  = 0 , =0, k=c), the Mohr-Coulomb

criterion reduces to a Tresca’s maximum shear stress criterion, 13=2k or
equivalently


Fy = J 2 cos   – --- – k = 0
 6

PLASTICITY MODELS FOR SOILS | 21

 Drucker-Prager Criterion
The Mohr-Coulomb criterion causes numerical difficulties when treating the plastic
flow at the corners of the yield surface. The Drucker-Prager model neglects the
influence of the invariant J3 (introduced by the Lode angle) on the cross-sectional
shape of the yield surface. It can be considered as the first attempt to approximate the
Mohr-Coulomb criterion by a smooth function based on the invariants I1 and J2
together with two material constants (which can be related to Mohr-Coulomb’s
coefficients)

Fy = J 2 + I 1 – k = 0

This is sometimes also called the extended von Mises criterion, since it is equivalent to
the von Mises criterion for metals when setting =

Figure 2-4: The Drucker-Prager criterion.

The coefficients in the Drucker-Prager model can be matched to the coefficients in the
Mohr-Coulomb criterion by

2 sin  2 3c cos 
 = -------  -------------------------- and k = ---------------------------
3  3  sin    3  sin  

The symbol ± is related to either matching the tensile meridian (positive sign) or the
compressive meridian (negative sign) of Mohr-Coulomb’s pyramid.

The matching at the tensile meridian (=) comes from setting

m  0  =  3 + sin     2 3 

22 | CHAPTER 2: GEOMECHANICS THEORY

in the Mohr-Coulomb criterion, and the matching at the compressive meridian (=
) from setting

m    3  =  3 – sin     2 3 

Compressive meridian Tensile meridian

Figure 2-5: The Drucker-Prager criterion showing the tensile and compressive meridians
(inner and outer circles), and the Lode angle compared to the cross section of
Mohr-Coulomb criterion in the -plane.

In the special case of frictionless material, (=0, =0, k = 2c  3 ), the Drucker

-Prager criterion reduces to the von Mises criterion

J 2 = 2c  3

When matching Drucker-Prager criterion to Mohr-Coulomb criterion in 2D

plane-strain applications, the parameters are

tan  3c
 = ------------------------------------ and k = ------------------------------------
2 2
9 + 12 tan  9 + 12 tan 

and when matching both criteria in 2D plane-stress applications, the matching

parameters are:

1 2
 = -------  sin  and k = -------  c cos 
3 3

PLASTICITY MODELS FOR SOILS | 23

 DILATATION ANGLE
The Mohr-Coulomb yield criterion is sometimes used with a nonassociated plastic
potential. This plastic potential could be either a Drucker-Prager criterion, or the same
Mohr-Coulomb yield function but with a different slope with respect to the
hydrostatic axis, in which case the angle of internal friction is replaced by the dilatation
angle, which is normally smaller (Ref. 5).

Also, when using a Drucker-Prager criterion matched to a Mohr-Coulomb criterion,

the plastic potential could also be nonassociated, in which case the difference between
the dilatation angle and the angle of internal friction would result in a yield surface and
plastic potential portrayed by two cones with different angles with respect to the
hydrostatic angle.

ELLIPTIC CAP
The Mohr-Coulomb and Drucker-Prager criteria portray a conic yield surface which is
open in the hydrostatic axis direction. Normally, these soil models are not accurate
above a given limit pressure because real-life materials cannot bear infinite loads and
still behave elastically. A simple way to overcome this problem is to add an elliptical
end-cap to these soil models.

The elliptic cap is an elliptic yield surface of a semi-axes as shown in Figure 2-6. The
initial pressure pa (SI units: Pa) denotes the pressure at which the elastic range
circumscribed by either a Mohr-Coulomb pyramid or a Drucker-Prager cone is not
valid any longer, so a cap surface is added. The limit pressure pb gives the curvature of
the ellipse, and denotes the maximum admissible hydrostatic pressure for which the
material starts deforming plastically. Pressures higher than pb are not allowed

pa pb p

Figure 2-6: Elliptic cap model in Haigh–Westergaard coordinate system.

24 | CHAPTER 2: GEOMECHANICS THEORY

Note that the sign convention for the pressure is taken from the Structural Mechanics
Module: positive sign under compression, so pa and pb are positive parameters.
Figure 2-6 shows the cap in terms of the variables

q = 3J 2    and p = – I 1     3

Matsuoka-Nakai Criterion
Matsuoka and Nakai (Ref. 3) discovered that the sliding of soil particles occurs in the
plane in which the ratio of shear stress to normal stress has its maximum value, which
they called the mobilized plane. They defined the yield surface as

F y =  9 + 9 2 I 3 – I 1 I 2 = 0

where the parameter =/nSTP equals the maximum ratio between shear stress
and normal stress in the spatially mobilized plane (STP-plane), and the invariants are
applied over the effective stress tensor (this is the Cauchy stress tensor minus the fluid
pore pressure).

The Matsuoka-Nakai criterion circumscribes the Mohr-Coulomb criterion in dry soils,

when

2 2
 = ----------- tan 
3

and  denotes the angle of internal friction in Mohr-Coulomb criterion.

Figure 2-7: The Matsuoka-Nakai criterion and Mohr-Coulomb criterion in the principal
stress space.

PLASTICITY MODELS FOR SOILS | 25

 Lade-Duncan Criterion
The Lade-Duncan criterion was originally developed to model a large volume of
laboratory sample test data of cohesionless soils. This criterion is defined as

3
F y = kI 3 – I 1 = 0

where I1 and I3 are the first and third stress invariants respectively, and k is a parameter
related to the direction of the plastic strain increment in the triaxial plane. The
parameter k can vary from 27 for hydrostatic stress conditions (1=2=3), up to a
critical value kc at failure. In terms of the invariants I1, J2,and J3, this criterion can be
written as

F y = J 3 – --- I 1 J 2 +  ------ – --- I 1 = 0

1 1 1 3
3  27 k

The Lade-Duncan criterion can be fitted to the compressive meridian of the

Mohr-Coulomb surface by choosing

k =  3 – sin   3 /  cos  2   1 – sin   

with  as the angle of internal friction in Mohr-Coulomb criterion

Lade-Duncan
Matsuoka-Nakai

Mohr-Coulomb

Figure 2-8: Comparing the Mohr-Coulomb, Matsuoka-Nakai, and Lade-Duncan

criteria when matching the tensile meridian.

The Lade-Duncan criterion does not match the Mohr-Coulomb criterion

(nor the Matsuoka-Nakai criterion) at the tensile meridian.

26 | CHAPTER 2: GEOMECHANICS THEORY

Tension Cut-Off
It appears that the Mohr-Coulomb and Drucker-Prager criteria predict tensile
strengths larger than the experimental measurements on soil samples. This discordance
can be mended by the introduction of the Rankine or tension cut-off criterion.

The Rankine criterion states that a material stops deforming elastically when the
biggest principal stress 1 reaches a maximum tensile stress, also called tension cut-off
limit t.

In terms of the principal stress, Rankine criterion reads

F = 1 – t = 0

For soils and clays, the maximum tensile stress can be estimated from the material
parameters, such as the cohesion c and the friction angle  For instance, the tip of the
cone in Mohr-Coulomb criterion is reached when

cos 
 1 = c ------------
sin 

therefore, the tension cut-off should be chosen such as

cos 
 t  c ------------
sin 

The Mohr-Coulomb criterion together with a tension cut-off is sometimes called

modified Mohr-Coulomb criterion (Ref. 19).

Tension cut-off is also available with Concrete material models.

PLASTICITY MODELS FOR SOILS | 27

 Theory for the Cam-Clay Material
In this section:

• About the Cam-clay Material

• Volumetric Elastic Behavior
• Hardening and Softening
• Including Pore Fluid Pressure

About the Cam-clay Material

The Cam-clay material model was developed at the University of Cambridge in the
1970s, and since then it has experienced different modifications. The modified
Cam-clay model is the most commonly used due to the smooth yield surface, and it is
the one implemented in the Geomechanics Module.

The Cam-clay model is a so-called critical state model, where the loading and
unloading of the material follows different trajectories in stress space. The model also
features hardening and softening of clays. Different formulations can be found in
textbooks about these models (see Ref. 13, 14, and 15).

The yield function is written in terms of the variables

q = 3J 2    and p = – I 1     3

Following the Structural Mechanics Module sign convention:

F y = q 2 + M 2  p – p c p = 0

This is an ellipse in p-q plane, with a cross section independent of Lode angle and
smooth for differentiation. Note that p, q and pc are always positive variables.

The material parameter M0 defines the slope of a line in the p-q space called critical
state line, and it can be related to the angle of internal friction  in the
Mohr-Coulomb criterion

6 sin   
M = --------------------------
3 – sin   

28 | CHAPTER 2: GEOMECHANICS THEORY

 q = Mp
q
Critical state line
softening

hardening

Nonlinear elastic region

pc/2 pc p

Figure 2-9: Modified Cam-clay surface in the p-q plane. The ellipse circumscribes a
nonlinear elastic region.

In the Cam-clay model, hardening is controlled by the consolidation pressure pc,

which depends exponentially on the volumetric plastic strain pvol.

– B p  pvol
p c = p c0 e (2-4)

The volumetric plastic strain is available in the variable [Link] and

the consolidation pressure in the variable [Link].

Here, the parameter pc0 is the initial consolidation pressure, and the exponent Bp is a
parameter which depends on the initial void ratio e0, the swelling index , and the
compression index 

1 + e0
B p = ---------------
–

The initial void ratio, the compression index, and the swelling index are all positive
parameters and must fulfill

0     , so B p  0

THEORY FOR THE CAM-CLAY MATERIAL | 29

 The void ratio e is the ratio between pore volume and solid volume. It can
be written in terms of the porosity  as e=  /(1  ).

In this formulation, the compression index  is the slope of the virgin isotropic
consolidation line, and  is the slope of the rebound-reloading line (also called
loading-reloading line) in the e versus ln(p) plane.

e
N



ln(prefN) ln(pC0) ln(p)

Figure 2-10: Slopes of the virgin isotropic consolidation line, and rebound-reloading line
in the e vs. ln(p) plane. The reference void ratio N is measured at the reference pressure
prefN.

If an Initial Stress and Strain feature subnode is added to the Cam-clay

material, the initial consolidation pressure pc0 must be made equal or
bigger than one third of minus the trace of the initial stress tensor,
otherwise the initial stress state is outside the Cam-clay ellipse.

Volumetric Elastic Behavior

The stress-strain relation beyond the elastic range is of great importance in soil
mechanics. For additive decomposition of strains, Cauchy’s stress tensor is written as

 –  0 = C:   –  0 –  th –  p –  c 

Here,  is the Cauchy stress tensor,  is the total strain tensor, 0 and 0 are the initial
stress and strain tensors, th is the thermal strain tensor, p is the plastic strain tensor,
c is the creep strain tensor, and, C is the fourth-order elasticity tensor.

30 | CHAPTER 2: GEOMECHANICS THEORY

For a linear elastic material, the trace of the Cauchy’s stress tensor is linearly related to
the volumetric elastic strain (the trace of the elastic strain tensor) by the elastic bulk
modulus

p = – I 1     3 = p 0 – K elvol

here p0=trace0/3 is the trace of the initial stress tensor 0, and K is the bulk
modulus, a constant parameter independent of the stress or strain.

The elastic volumetric strain is available in the variable [Link].

The modified Cam-clay model introduces a nonlinear relation between stress and
volumetric elastic strain

– B e  elvol 1 + e0
p – p0 = –K0 e with B e = ---------------


and K0 a reference bulk modulus. This formulation gives a tangent bulk modulus
KTBepp0. The reference bulk modulus is calculated from the initial
consolidation pressure pc0, and the void ratio at reference pressure N.

Hardening and Softening

The yield surface for the modified Cam-clay model reads

F y = q 2 + M 2  p – p c p = 0

The associated flow rule (Qp=Fy) and the yield surface written in terms of these two
invariants, Fy(I1, J2), gives a rate equation for the plastic strain tensor calculated from
the derivatives of Fy with respect to the stress tensor 

Q p Q p I 1 Q p J 2
 p =  p ---------- =  p  ---------- -------- + ---------- ---------
·
  I 1  J 2  

Here, p means the plastic multiplier, see Plastic Flow for Small Strains
and Hardening Models.

THEORY FOR THE CAM-CLAY MATERIAL | 31

 The plastic strain rate tensor · p includes both deviatoric and isotropic parts. Note that

I 1   = I and J 2   = dev   

These relations can be used for writing the plastic flow as

1 Q p Q p 3
 p =  p  – --- ---------- I + ---------- ------- dev    =  p  – --- M 2  2p – p c  I + 3dev   
· 1
 3 p q 2q   3 

The trace of the plastic strain rate tensor (the volumetric plastic strain rate · pvol )
then reads

· ·
 pvol = trace   p  =  p  – M 2  2p – p c  

This relation explains the reason why there is isotropic hardening for p>pc/2 and
isotropic softening for p<pc/2. So the volumetric plastic strain can either increase or
decrease.

In the Cam-clay model, the hardening is controlled by the consolidation pressure

variable p c   pvol  as a function of volumetric plastic strain, as written in Equation 2-4.
Hardening introduces changes in the shape of the Cam-clay ellipse, since its major
semiaxis depends on the value of pc.

Including Pore Fluid Pressure

When a pore fluid pressure pf is added to the Cam-clay material, the yield surface is
shifted on the p axis

Fy = q 2 + M 2  p – pf – pc   p – pf  = 0

The quantity ppf is normally regarded as the effective pressure, or effective stress,
which should not be confused with von Mises stress.

32 | CHAPTER 2: GEOMECHANICS THEORY

Failure Criteria for Concrete, Rocks,
a nd O t he r B r i t t l e Mat eri al s
In this section:

• Bresler-Pister Criterion
• Willam-Warnke Criterion
• Ottosen Criterion
• Original Hoek-Brown Criterion
• Generalized Hoek-Brown Criterion

Bresler-Pister Criterion
The Bresler-Pister criterion (Ref. 2) was originally devised to predict the strength of
concrete under multiaxial stresses. This failure criterion is an extension of the
Drucker-Prager criterion to brittle materials and can be expressed in terms of the stress
invariants as

2
Fy = J2 + k1 I1 + k2 I1 + k3 = 0 (2-5)

here, k1, k2, and k3 are material parameters.

This criterion can also be written (Ref. 17) in term of the uniaxial compressive strength
fc and the octahedral normal and shear stresses

 oct = I 1 /3 and  oct = 2J 2 /3

so Equation 2-5 simplifies to

 oct  oct  oct

2
F y = --------- – a + b ---------- – c ---------
- = 0
fc fc f c2

Here, the parameters a, b, and c are obtained from the uniaxial compression, uniaxial
tension, and biaxial compression tests, respectively. The octahedral normal stress oct
is considered positive when tensile, and fc is taken positive.

FAILURE CRITERIA FOR CONCRETE, ROCKS, AND OTHER BRITTLE MATERIALS | 33

 Willam-Warnke Criterion
The Willam-Warnke criterion (Ref. 10) is used to predict failure in concrete and other
cohesive-frictional materials such as rock, soil, and ceramics. Just as the Bresler-Pister
criterion, it depends only on three parameters. It was developed to describe initial
concrete failure under triaxial conditions. The failure surface is convex, continuously
differentiable, and is fitted to test data in the low compression range. The material is
considered perfect elastoplastic (no hardening).

The original “three-parameter” Willam-Warnke failure criterion was defined as

 oct
--- --------- + r      --- – ----  oct – 1 = 0
3 1 1
Fy = (2-6)
5 fc   f t f b 

here, fc is the uniaxial compressive strength, ft is uniaxial tensile strength, and fb is

obtained from the biaxial compression test. All parameters are positive. The octahedral
normal and shear stresses are defined as usual

 oct = I 1 /3 , and  oct = 2J 2 /3

This can be written as

--- r     ------ I 1 – f c = 0
1
Fy = J2 + 5
2  3 

with

 = fb ft /  fb fc – ft fc 

The function r( describes the segment of an ellipse on the octahedral plane when
0      3 . By using the Lode angle , the dimensionless function r( is defined as

2
2r c  r c2 – r t2  cos  + r c  2r t – r c  4  r c2 – r t2  cos  + 5r t2 – 4 r c r t
r    = --------------------------------------------------------------------------------------------------------------------------------------------------------------
2
-
4  r c2 – r t2  cos  +  r c – 2r t  2

Here, the tensile and compressive meridian rt and rc are defined in terms of the positive
parameters fc, fb, and ft:

fb ft
rt = 6
--- ----------------------------
5 2f b f c + f t f c

34 | CHAPTER 2: GEOMECHANICS THEORY

 fb ft
rc = 6
--- -------------------------------------------
5 3f b f t + f b f c – f t f c

The function r( can be interpreted as the friction angle which depends on the Lode
angle (Ref. 10).



rc
r 
rt

Figure 2-11: The deviatoric section of Willam-Warnke failure criterion.

Ottosen Criterion
The Ottosen criterion is a four-parameter failure criterion proposed for short-time
loading of concrete. It corresponds to a smooth convex failure surface with curved
meridians, which is open in the negative (compressive) direction of the hydrostatic
axis. The curve in the pi-plane changes from almost triangular to a more circular shape
with increasing hydrostatic pressure. The criterion is in agreement with experimental
results over a wide range of stress states, including both triaxial tests along the tensile
and the compressive meridian and biaxial tests (Ref. 18).

The Ottosen criterion is commonly written as (Ref. 17, Ref. 18):

a
F y = ---- J 2 +     J 2 + bI 1 – f c = 0
fc

FAILURE CRITERIA FOR CONCRETE, ROCKS, AND OTHER BRITTLE MATERIALS | 35

 In this formulation, the parameters a>0, b>0, k1>0, and k2>0 are dimensionless, and
fc>0 is the uniaxial compressive strength of concrete (with a positive sign). The
function      0 (dimensionless) is defined as


 k cos  1 --- acos  k 2 cos  3   J3  0
 1 3 

 = 

  1 
 k 1 cos  --3- – --3- acos  – k 2 cos  3   J 3  0,


The parameter k1 is called the size factor. The parameter k2 (also called shape factor)
is positive and bounded to 0k21(Ref. 17, Ref. 18).

Typical values for these parameters are obtained by curve-fitting the uniaxial
compressive strength fc, uniaxial tensile strength ft, and from the biaxial and triaxial
data (for instance, a typical biaxial compressive strength of concrete is 16% higher than
the uniaxial compressive strength). The parameters fc, fb, and ft are positive.
TABLE 2-1: TYPICAL PARAMETER VALUES FOR OTTOSEN FAILURE CRITERION ( Ref. 18).
ft/ fc a b k1 k2 lt lc
0.08 1.8076 4.0962 14.4863 0.9914 14.4725 7.7834
0.10 1.2759 3.1962 11.7365 0.9801 11.7109 6.5315
0.12 0.9218 2.5969 9.9110 0.9647 9.8720 5.6979

The compressive and tensile meridians (as defined in the Willam-Warnke criterion) are

1 1
r c = ----- = -------------------
 c   /3 

1 1
r t = ----- = -----------
t   0 

For concrete, the ratio

 c / t = r t /r c

36 | CHAPTER 2: GEOMECHANICS THEORY

normally lies between 0.54~0.58.

The Ottosen criterion is equivalent to Drucker-Prager when a = 0 and

= constant.

Original Hoek-Brown Criterion

The Hoek-Brown criterion is an empirical type of model which is commonly used
when dealing with rock masses of variable quality. The Hoek-Brown criterion is widely
used within civil engineering and is popular because the material parameters can be
estimated based on simple field observations together with knowledge of the uniaxial
compressive strength of the intact rock material. The Hoek–Brown criterion is one of
the few nonlinear criteria widely accepted and used by engineers to estimate the yield
and failure of rock masses. The original Hoek-Brown failure criterion states (Ref. 5)

 1 =  3 + m c  3 + s c2

where 1230 are the principal stresses at failure (as defined in geotechnical
engineering; that is, an absolute value), c is the uniaxial compressive strength of the
intact rock (positive parameter), and m and s are positive material parameters.

If the expression is converted into to the sign convention for principal stresses in the
Structural Mechanics Module, it becomes

 1 =  3 + – m c  1 + s c2

with c, m, and s positive material parameters. (In this case, note that 1 sc/m).

As developed originally, there is no relation between the parameters m and s and the
physical characteristics of a rock mass measured in laboratory tests. However, for intact
rock, s1 and mmi, which is measured in a triaxial test.

For jointed rock masses, 0s1 and mmi. The parameter m usually lies in the
range 5m30 (Ref. 7).
TABLE 2-2: CHARACTERISTIC VALUES FOR DIFFERENT ROCK TYPES

m ROCK TYPE

5 Carbonate rocks, dolomite, limestone

10 Consolidated rocks, mudstone, shale

FAILURE CRITERIA FOR CONCRETE, ROCKS, AND OTHER BRITTLE MATERIALS | 37

 TABLE 2-2: CHARACTERISTIC VALUES FOR DIFFERENT ROCK TYPES

m ROCK TYPE

15 Sandstone
20 Fine grained rocks
25 Coarse grained rocks
The Hoek-Brown criterion can be written in terms of the invariants I1 and J2 and the
Lode angle 03, so

 1
F y = 2 J 2 sin   + --- –  c s – m ------ = 0
 3 c

Generalized Hoek-Brown Criterion

The generalized Hoek-Brown criterion was developed in order to fit the Geological
Strength Index (GSI) classification of isotropic rock masses (Ref. 6). A new
relationship between GSI, m, s and the newly introduced parameter a was developed,
to give a smoother transition between very poor quality rock masses (GSI<25) and
stronger rocks

1 a
 1 –  3 =  c  s – m b ------
c

In terms of the invariants J2 and the Lode angle 0      3 this equals

 1 a
F y = 2 J 2 sin   + --- –  c  s – m b ------ = 0
3 c

where 123 are the principal stresses (using the Structural Mechanics Module
conventions) of the effective stress tensor (this is, the stress tensor minus the fluid pore
pressure).

The positive parameter mb is a reduced value of the material constant mi:

m b = m i exp  --------------------------
GSI – 100
28 – 14D

s and a are positive parameters for the rock mass given by the following relationships:

38 | CHAPTER 2: GEOMECHANICS THEORY

 s = exp  --------------------------
GSI – 100
 9 – 3D 

a = --- + ---  exp  ------------- – exp  --------- 

1 1 – GSI – 20
2 6 15 3

The disturbance factor D was introduced to account for the effects of stress relaxation
and blast damage, and it varies from 0 for undisturbed in-situ rock masses to 1 for very
damaged rock masses.
TABLE 2-3: DISTURBANCE FACTOR IN ROCK MASSES

D DESCRIPTION OF ROCK MASS

0 Undisturbed rock mass

0~0.5 Poor quality rock mass
0.8 Damaged rock mass
1.0 Severely damaged rock mass

FAILURE CRITERIA FOR CONCRETE, ROCKS, AND OTHER BRITTLE MATERIALS | 39

 Elastoplastic Material Theory
In this section:

• Introduction to Small and Large Plastic Strains

• Plastic Flow for Small Strains
• Isotropic Plasticity
• Yield Function
• Hardening Models
• Plastic Flow for Large Strains
• Numerical Solution of the Elastoplastic Conditions
• References for Elastoplastic Materials

Working with the Geomechanics Materials

Introduction to Small and Large Plastic Strains

There are two implementations of plasticity available in COMSOL Multiphysics. One
is based on the additive decomposition of strains, which is the most suitable approach
in the case of small strains, and the other one is based on the multiplicative
decomposition of the deformation gradient, which is more suitable when large plastic
strains occurs.

The additive decomposition of strains is used when small plastic strain is selected as
the plasticity model. The stress-strain relationship is then written as

 –  0 = C:   –  0 –  th –  p –  c 

Here:

•  is the Cauchy stress tensor

•  is the total strain tensor
• 0 and 0 are the initial stress and strain tensors
• th is the thermal strain tensor

40 | CHAPTER 2: GEOMECHANICS THEORY

• p is the plastic strain tensor
• c is the creep strain tensor, and
• C is the fourth-order elasticity tensor.

The elastic strain tensor is computed by removing inelastic and initial strains from the
total strain tensor

 el =  –  0 –  th –  p –  c

In the special case of zero initial strains and no creep or thermal expansion, the elastic
strain tensor simplifies to

 el =  –  p (2-7)

If  or p are large, the additive decomposition might produce incorrect results;

however, the additive decomposition of strains is widely use for metal and soil
plasticity.

When Small plastic strain is selected as the plasticity model for the
Plasticity node, and the Include geometric nonlinearity check box is
selected on the study settings window, the second Piola-Kirchhoff stress
is used instead of the Cauchy stress, and the additive decomposition of
strains is understood as the summation of Green-Lagrange strains.

When large plastic strain is selected as the plasticity model, the total deformation
gradient tensor is multiplicatively decomposed into elastic deformation gradient and
plastic deformation gradient

F = F el F p

The plastic deformation is removed from the total deformation gradient,

–1
F el = FF p

so the elastic Green-Lagrange strain tensor is computed from the elastic

deformation gradient tensor Fel

1 T
 el = ---  F el F el – I 
2

ELASTOPLASTIC MATERIAL THEORY | 41

 and the plastic Green-Lagrange strain tensor is computed from the plastic
deformation gradient tensor

1 T
 p = ---  F p F p – I 
2

As opposed to the small strain formulation described in Equation 2-7, the total, plastic
and elastic Green-Lagrange strain tensors are related as

–T –1
 el = F p   –  p F p

Under multiplicative decomposition, the elastic right Cauchy-Green tensor and the
plastic right Cauchy-Green tensor are defined as

T T
C el = F el F el and C p = F p F p

When Large plastic strain is selected as the plasticity model for the
Plasticity node, the Include geometric nonlinearity check box in the study
settings is automatically selected and becomes unavailable.

Plastic Flow for Small Strains

The flow rule defines the relationship between the increment of the plastic strain
tensor · p and the current state of stress, for a yielded material subject to further
loading. When Small plastic strain is selected as the plasticity model for the Plasticity
node, the direction of the plastic strain increment is defined by

· Q p
 p =  ----------


Here, is a positive multiplier (also called the consistency parameter or plastic
multiplier) which depends on the current state of stress and the load history, and Qp
is the plastic potential.

The “dot” (for · p ) means the rate at which the plastic strain tensor
changes with respect to Qp/. It does not represent a true time
derivative. Some authors call this formulation rate independent
plasticity.

42 | CHAPTER 2: GEOMECHANICS THEORY

The direction of the plastic strain increment · p is perpendicular to the surface (in the
space of principal stresses) defined by the plastic potential Qp.

The plastic multiplier  is determined by the complementarity or Kuhn-Tucker

conditions

  0 , F y  0 and F y = 0 (2-8)

where Fy is the yield function. The yield surface encloses the elastic region defined by
Fy<0. Plastic flow occurs when Fy=0.

If the plastic potential and the yield surface coincide with each other Qp=Fy, the
flow rule is called associated, and the rate in Equation 2-9 is solved together with the
conditions in Equation 2-8.

· F y
 p =  --------- (2-9)


For a non-associated flow rule, the yield function does not coincide with the plastic
potential, and together with the conditions in Equation 2-8, the rate in Equation 2-10
is solved for the plastic potential Qp (often, a smoothed version of Fy).

· Q p
 p =  ---------- (2-10)


The evolution of the plastic strain tensor · p (with either Equation 2-9 or
Equation 2-10, plus the conditions in Equation 2-8) is implemented at Gauss points
in the plastic element elplastic.

Isotropic Plasticity
For isotropic plasticity, the plastic potential Qp is written in terms of at most three
invariants of Cauchy’s stress tensor

Q p    = Q p  I 1    J 2    J 3    

where the invariants of the stress tensor are

I 1    = trace   
1
J 2    = --- dev   :dev   
2
J 3    = det  dev    

ELASTOPLASTIC MATERIAL THEORY | 43

 so that the increment of the plastic strain tensor · p can be decomposed into

Q p Q p I 1 Q p J 2 Q p J 3
 p =  ---------- =   ---------- -------- + ---------- --------- + ---------- ---------
·
 I 1  J 2  J 3 

The increment in the plastic strain tensor · p includes in a general case both deviatoric
and volumetric parts. The tensor · p is symmetric given the following properties

I 1
-------- = I

J (2-11)
---------2 = dev   

J 3
--------- = dev   dev    – 2
--- J 2 I
 3

A common measure of inelastic deformation is the effective plastic strain rate, which
is defined as

· 2
 pe = --- · p :· p (2-12)
3

The trace of the incremental plastic strain tensor, which is called the volumetric plastic
strain rate · pvol , depends only on the reliance of the plastic potential on the first
invariant I1(), sinceJ2/ and J3/ are deviatoric tensors

Q p Q p
 pvol = trace   p  =  trace  ---------- = 3 ----------
· ·
 I 1

For metal plasticity under the von Mises and Tresca criteria, the volumetric plastic
strain rate is always zero because the plastic potential is independent of the invariant
I1. This is known as J2 plasticity. However, this is not the case for most materials
used in geotechnical applications. For instance, a nonzero volumetric plastic strain is
explicitly used in the Cam-Clay material.

Plasticity

44 | CHAPTER 2: GEOMECHANICS THEORY

 The effective plastic strain and the volumetric plastic strain are available in
the variables [Link] and [Link].

Yield Function
When an associated flow rule is applied, the yield function must be smooth, that is,
continuously differentiable with respect to the stress. In COMSOL Multiphysics, the
following form is used:

F y =     –  ys

where ys is the yield stress.

The predefined form of the effective stress is the von Mises stress, which is
commonly used in metal plasticity:

3
 mises = 3J 2    = --- dev   :dev   
2

Other expressions can be defined, such as Tresca stress or another user-defined

expression.

The Tresca effective stress is calculated from the difference between the largest and the
smallest principal stress

 tresca =  1 –  3

A user-defined yield function can by expressed in terms of invariants of the stress tensor
such as the pressure (volumetric stress)

1
p = – --- I 1   
3

the effective (von Mises) stress mises, or other invariants, principal stresses, or stress
tensor components.

Plasticity

ELASTOPLASTIC MATERIAL THEORY | 45

 Hardening Models
The plasticity model implements three different kinds of hardening models for
elastoplastic materials:

• Perfect plasticity (no hardening)

• Isotropic hardening
• Kinematic hardening

PERFECT (OR IDEAL) PLASTICITY

In this case the plasticity algorithm solves either the associated or non-associated flow
rule for the plastic potential Qp

· Q p
 p =  ----------


with the yield function

F y =     –  ys0

In the settings for plasticity you specify the effective stress for the yield function
from a von Mises stress, a Tresca stress, or a user-defined expression.

When Large plastic strain is selected as the plasticity model for the
Plasticity node, either the associate or non-associated flow rule is applied
as written in Equation 2-15.

ISOTROPIC HARDENING
In this case the plasticity algorithm solves either the associated or non-associated flow
rule for the plastic potential Qp

· Q p
 p =  ----------


with the yield function

F y =     –  ys   pe 

where pe is the effective plastic strain. The variable yspe is the yield stress, which
now depends on the effective plastic strain. The yield stress versus the effective plastic
strain can be specified in two different ways — tangent data (linear isotropic

46 | CHAPTER 2: GEOMECHANICS THEORY

hardening), or hardening function data. When using hardening function data, the
hardening curve could also depend on other variables, such as stress or temperature.

Tangent Data (Linear Isotropic Hardening)

In this case, an isotropic tangent modulus ETiso is given. This tangent modulus is
defined as (stress increment / total strain increment), and it relates the hardening to
the effective plastic strain linearly. The yield stress yspe is a linear function of the
effective plastic strain

 ys   pe  =  ys0 +  h   pe 

1 1 1
with  h   pe  = k pe and --- = -------------- – ----
k E Tiso E

here, ys0 is the initial yield stress, and k is the isotropic hardening modulus. A value
for ETiso is entered in the isotropic tangent modulus section for the Plasticity node.
The Young’s modulus E is taken from the linear elastic material. For orthotropic and
anisotropic elastic materials, E represents an average Young’s modulus.

Hardening Function Data

In this case, define the (usually nonlinear) hardening function hpe such that the
yields stress reads

 ys   pe  =  ys0 +  h   pe 

This is the preferred way to define nonlinear hardening models.

The internal variable for the effective plastic strain is named [Link].
The effective plastic strain evaluated at Gauss points is named
[Link], where solid is the name of the interface identifier for the
physics interface.

When Large plastic strain is selected as the plasticity model for the
Plasticity node, either the associate or non-associated flow rule is applied
as written in Equation 2-15.

ELASTOPLASTIC MATERIAL THEORY | 47

 KINEMATIC HARDENING
The algorithm solves either the associated or non-associated flow rule for the plastic
potential Qp

· Q p
 p =  ----------


with the yield function defined as

F y =    –  shift  –  ys0 and  shift = f   p 

Here, ys0 is the initial yield stress, and the effective stress is either the von Mises
stress or a user-defined expression. The stress tensor used in the yield function is shifted
by what is usually called the back stress, shift.

The back stress is generally not only a function of the current plastic strain but also of
its history. In the case of linear kinematic hardening, f p is a linear function of the
plastic strain tensor p, this is also known as Prager’s hardening rule. The
implementation of kinematic hardening assumes a linear evolution of the back stress
tensor with respect to the plastic strain tensor:

2
 shift = --- c p
3

where the work hardening constant is calculated from

--- =  --------------- – ----

1 1 1
c  E Tkin E

The value for ETkin is entered in the kinematic tangent modulus section and the
Young’s modulus E is taken from the linear elastic material model. For orthotropic and
anisotropic elastic materials, E represents an average Young’s modulus.

When kinematic hardening is added, both the plastic potential and the
yield surface are calculated with effective invariants, that is, the invariants
of the tensor defined by the difference between the stress tensor minus
the back-stress,  eff =  –  shift . The invariant of effective deviatoric
tensor is named solid.II2sEff, which is used when a von Mises
plasticity is computed together with kinematic hardening.

48 | CHAPTER 2: GEOMECHANICS THEORY

Plastic Flow for Large Strains
When Large plastic strain is selected as the plasticity model for the Plasticity node, a
multiplicative decomposition of deformation (Ref. 1, Ref. 2, and Ref. 3) is used, and
the associated plastic flow rule can be written as the Lie derivative of the elastic left
Cauchy-Green deformation tensor Bel:

1 
– --- L  B el  =  ------- B el (2-13)
2 

The plastic multiplier  and the yield function  (written in terms of the Kirchhoff
stress tensor ) satisfy the Kuhn-Tucker condition, as done for infinitesimal strain
plasticity

  0 ,   0 and  = 0

The yield function in Ref. 1 and Ref. 2 was written in terms of
Kirchhoff stress and not Cauchy stress because the authors defined the
plastic dissipation with the conjugate energy pair  and d, where d is the
rate of strain tensor.

The Lie derivative of Bel is then written in terms of the plastic right Cauchy-Green rate

·
L  B el  = FC p– 1 F T (2-14)

By using Equation 2-13 and Equation 2-14the either associated or non-associated

plastic flow rule for large strains is written as (Ref. 2)

1 · Q p
– --- FC p– 1 F T =  ---------- B el (2-15)
2 

together with the Kuhn-Tucker conditions for the plastic multiplier  and the yield
function Fy

  0 , F y  0 and F y = 0 (2-16)

For the associated flow rule, the plastic potential and the yield surface coincide with
each other (Qp=Fy), and for the non-associated case, the yield function does not
coincide with the plastic potential.

In COMSOL Multiphysics, the elastic left Cauchy-Green tensor is written in terms of

the deformation gradient and the right Cauchy-Green tensor, so Bel=FCp1F. The

ELASTOPLASTIC MATERIAL THEORY | 49

 plastic flow rule is then solved at Gauss points in the plastic element elplastic for the
inverse of the plastic deformation gradient Fp1, so that the variables in Equation 2-15
are replaced by

· · · T = FF – 1 F – T F T
C p– 1 = F p– 1 F p–T + F p–1 F p– T and B el = F el F el p p

The flow rule then reads

1 · · Q p
– ---  F p– 1 F p– T + F p– 1 F p–T  = F –1 ---------- FF p– 1 F p–T (2-17)
2 

After integrating the flow rule in Equation 2-17, the plastic Green-Lagrange strain
tensor is computed from the plastic deformation tensor

1 T
 p = ---  F p F p – I 
2

and the elastic Green-Lagrange strain tensor is computed from the elastic deformation
gradient tensor Fel=FFp1

1 T
 el = ---  F el F el – I 
2

When Large plastic strain is selected as the plasticity model for the
Plasticity node, the effective plastic strain variable is computed as the true
effective plastic strain (also called Hencky or logarithmic plastic strain).

Numerical Solution of the Elastoplastic Conditions

A backward Euler discretization of the pseudo-time derivative is used in the plastic flow
rule. For small plastic strains, this gives

Q p
 p –  p,old =  ----------


where old denotes the previous time step and  = t , where t is the pseudo-time
step length. For large plastic strains, Equation 2-17 gives

1 T T T Q p T
– --- F  2MM – M old M – MM old  =  ---------- FMM
2 
–1
where M = F p .

50 | CHAPTER 2: GEOMECHANICS THEORY

For each Gauss point, the plastic state variables (p and M, respectively) and the plastic
multiplier, , are computed by solving the above time-discretized flow rule together
with the complementarity conditions

  0 , F y  0 and F y = 0

This is done as follows (Ref. 1):

1 Elastic-predictor: Try the elastic solution  p =  p,old (or M = M old ) and  = 0 .

If this satisfies F y  0 it is done.
2 Plastic-corrector: If the elastic solution does not work (this is F y  0 ), solve the
nonlinear system consisting of the flow rule and the equation F y = 0 using a
damped Newton method.

References for Elastoplastic Materials

1. J. Simo and T. Hughes, Computational Inelasticity, Springer, 1998.

2. J. Simo, “Algorithms for Static and Dynamic Multiplicative Plasticity that Preserve
the Classical Return Mapping Schemes of the Infinitesimal Theory,” Computer
Methods in Applied Mechanics and Engineering, vol. 99, pp. 61–112, 1992.

3. J. Lubliner, Plasticity Theory, Dover, 2008.

4. R. Hill, “A Theory of the Yielding and Plastic Flow of Anisotropic Metals,” Proc.
Roy. Soc. London, vol. 193, pp. 281–297, 1948.

5. N. Ottosen and M. Ristinmaa, The Mechanics of Constitutive Modeling, Elsevier

Science, 2005.

ELASTOPLASTIC MATERIAL THEORY | 51

 Creep and Viscoplasticity
In this section:

• About Creep
• Fundamental Creep Material Models
• Solver Settings for Creep

Working with the Geomechanics Materials

About Creep

In the literature, the terms viscoplasticity and creep are often used
interchangeably to refer to the class of problems related to
rate-dependent plasticity.

Creep is an inelastic time-dependent deformation that occurs when a material is

subjected to stress (typically much less than the yield stress) at sufficiently high
temperatures.

The creep strain rate, in a general case, depends on stress, temperature, and time,
usually in a nonlinear manner:

·
 c = F cr   T t 

It is often possible to separate these effects as shown in this equation:

F cr   t T  = f 1   f 2  T f 3  t 

Experimental data shows three types of behavior for the creep strain rate at constant
stress as function of time. Researchers normally subdivide the creep curve into three
regimes, based on the fact that many different materials show similar responses:

• In the initial primary creep regime (also called transient creep) the creep strain rate
decreases with time to a minimum steady-state value.

52 | CHAPTER 2: GEOMECHANICS THEORY

• In the secondary creep regime the creep strain rate is almost constant. This is also
called steady-state creep.
• In the tertiary creep regime the creep strain increases with time until a failure
occurs.

When this distinction is assumed, the total creep rate can be additively split into
primary, secondary, and tertiary creep rates

·
 c = F cr1 + F cr2 + F cr3

In most cases, Fcr1 and Fcr3 depend on stress, temperature and time, while secondary
creep, Fcr2, depends only on stress level and temperature. Normally, secondary creep
is the dominant process. Tertiary creep is seldom important because it only accounts
for a small fraction of the total lifetime of a given material.

c

1

primary creep
2

secondary creep tertiary creep

1 >2

log time

Figure 2-12: Uniaxial creep as a function of logarithmic time.

Creep

CREEP AND VISCOPLASTICITY | 53

 Fundamental Creep Material Models
Despite the fact that the creep response of a given material is related to its atomic
structure, a macroscopic (continuum mechanics) description is normally appropriate
for modeling scientific and engineering problems.

The fundamental mathematical models available for modeling creep are:

• Creep Potential
• Volumetric Creep
• Deviatoric Creep
• User-Defined Creep

These creep models are contributing subnodes to the Linear Elastic Material and they
can also be combined with Plasticity and Thermal Expansion subnodes for more
advanced models.

• Creep
• Plasticity
• See the Solid Mechanics interface in the Structural Mechanics Module
User’s Guide for details about the Thermal Expansion and Linear
Elastic Material nodes.

CREEP PO TENTIAL
Some authors use a creep potential to describe the secondary creep rate, so that the
creep rate is written in a way similar to the flow rule for plasticity:

· Q cr
 c =  ----------- and   0


Here, Qcr is a user-defined creep potential, which is normally written in terms of

invariants of the stress tensor.

Volumetric creep is obtained when the creep potential depends only on the first
invariant of Cauchy stress tensor, I1, since

Q cr Q cr
----------- = ----------- I
 I 1

This is equivalent to that the creep potential would depend on the pressure pI1.

54 | CHAPTER 2: GEOMECHANICS THEORY

When the creep potential depends only on the second deviatoric invariant of Cauchy
stress tensor, J2, the deviatoric creep model is obtained since

Q cr Q cr
----------- = ----------- dev   
 J 2

This is equivalent to that the creep potential would depend on the effective stress
e3J2.

When (in SI units) the creep potential, Qcr, is given in units of Pa, the rate multiplier
is given in units of 1/s.

VO L U M E T R I C C R E E P
The creep strain rate is calculated by solving the rate equation

· 1
 c = --- F cr I
3

so the creep rate tensor is a diagonal tensor, and the trace of the creep rate tensor, the
volumetric creep strain rate, equals the user input Fcr

·
trace   c  = F cr (2-18)

The creep rate, Fcr, usually depends on the first invariant of Cauchy stress I1 or the
pressure pI1, in addition to the temperature and other material parameters.

Volumetric creep is not generally used to model creep in metals, but it is commonly
used to model creep in soils or other geological materials.

DEVIATORIC CREEP
The creep strain rate is calculated by solving the rate equation

·
 c = F cr n D

Here, nD is a deviatoric tensor coaxial to the stress tensor.

The creep rate, Fcr, normally depends on the second deviatoric invariant of the stress
J2 or the effective or von Mises (effective) stress e, in addition to the temperature
and other material parameters.

The deviatoric tensor nD is defined as

3 dev   
n D = --- -----------------
2 e

CREEP AND VISCOPLASTICITY | 55

 The resulting creep strain rate tensor is also deviatoric, since trace (nD)

·
trace   c  = F cr trace  n D  = 0

Given the property

3
n D :n D = ---
2

the effective creep strain rate equals the absolute value of the user input Fcr

· 2
 ce = --- · c :· c = F cr
3

Deviatoric creep is very popular to model creep in metals and alloys. For example,
Norton’s law is a deviatoric creep model.

The effective creep strain and the effective creep strain rate are available in
the variables [Link] and [Link].

USER-DEFINED CREEP
The creep strain rate is calculated by solving the rate equation

·
 c = F cr

where Fcr is a user-defined symmetric tensor field.

Potential, Volumetric, Deviatoric, or User defined creep do not overwrite

each other. These are all contributing nodes to the Linear Elastic Material
node.

ENERGY DISSIPATION
Since creep is an inelastic process, the dissipated energy density can be calculated by
integrating the creep dissipation rate density (SI unit: W/m3) given by

· ·
W cdr = : c

In case many creep sub-nodes are added to a Linear Elastic Material node, the creep
dissipation rate density is calculated from the total creep strain rate tensor · c .

56 | CHAPTER 2: GEOMECHANICS THEORY

The total energy dissipated by creep in a given volume can be calculated by a volume
integration of the dissipated creep energy density Wc (SI unit: J/m3).

When the Calculate dissipated energy check box is selected, the creep
dissipation rate density is available under the variable [Link] and the
dissipated creep energy density under the variable Wc.

• Creep
• Linear Elastic Material in the Structural Mechanics User’s Guide.

Solver Settings for Creep

Creep is a time-dependent phenomenon.

PHYSICS INTERFACE SETTINGS

In the Model Builder, click the Solid Mechanics node. In the settings window, under
Structural Transient Behavior, select Quasi-static to treat the elastic behavior as
quasi-static (with no mass effects; that is, no second-order time derivatives for the
displacement variables). Selecting this option gives a more efficient solution for
problems where the variation in time is slow when compared to the natural frequencies
of the system.

SOLVER SETTINGS
When Quasi-static is selected on the physics interface settings window, the automatic
solver suggestion changes the method for the Time Stepping from Generalized alpha to
BDF.

For a Fully Coupled node (or Segregated node for multiphysics problems), the default
Nonlinear method under Method and Termination is Automatic (Newton). To get a faster
computation time when the effective strain rate is low or moderate, select Constant
(Newton) as the Nonlinear method instead.

In the COMSOL Multiphysics Reference Manual:

• Studies and Solvers

• Fully Coupled and Segregated

CREEP AND VISCOPLASTICITY | 57

 References for the Geomechanics
Module
1. W.F. Chen and E. Mizuno, Nonlinear Analysis in Soil Mechanics: Theory and
Implementation (Developments in Geotechnical Engineering), 3rd ed., Elsevier
Science, 1990.

2. B. Bresler and K.S. Pister, “Strength of Concrete Under Combined Stresses,” ACI
Journal, vol. 551, no. 9, pp. 321–345, 1958.

3. H. Matsuoka and T. Nakai, “Stress-deformation and Strength Characteristics of Soil

Under Three Different Principal Stresses,” Proc. JSCE, vol. 232, 1974.

4. H. Matsuoka and T. Nakai, “Relationship Among Tresca, Mises, Mohr-Coulomb,

and Matsuoka-Nakai Failure Criteria,” Soils and Foundations, vol. 25, no. 4, pp.123–
128, 1985.

5. H.S. Yu, Plasticity and Geotechnics, Springer, 2006.

6. V. Marinos, P. Marinos, and E. Hoek, “The Geological Strength Index: Applications

and Limitations,” Bull. Eng. Geol. Environ., vol. 64, pp. 55–65, 2005.

7. J. Jaeger, N. G. Cook, and R. Zimmerman, Fundamentals of Rock Mechanics, 4th

ed., Wiley-Blackwell, 2007.

8. G. C. Nayak and O. C. Zienkiewicz, “Convenient Form of Stress Invariants for

Plasticity,” J. Struct. Div. ASCE, vol. 98, pp. 949–954, 1972.

9. A.J. Abbo and S.W. Sloan, “A Smooth Hyperbolic Approximation to the

Mohr-Coulomb Yield Criterion,” Computers and Structures, vol. 54, no. 3, pp. 427–
441, 1995.

10. K.J. Willam and E.P. Warnke, “Constitutive Model for the Triaxial Behavior of
Concrete,” IABSE Reports of the Working Commissions, Colloquium (Bergamo):
Concrete Structures Subjected to Triaxial Stresses, vol. 19, 1974.

11. B.H.G. Brady and E.T. Brown, Rock Mechanics for Underground Mining, 3rd
ed., Springer, 2004.

12. H.A. Taiebat and J.P. Carter, “Flow Rule Effects in the Tresca Model,” Computer
and Geotechnics, vol. 35, pp. 500–503, 2008.

58 | CHAPTER 2: GEOMECHANICS THEORY

13. A. Stankiewicz et al., “Gradient-enhanced Cam-Clay Model in Simulation of Strain
Localization in Soil,” Foundations of Civil and Environmental Engineering, no.7,
2006.

14. D.M. Wood, Soil Behaviour and Critical State Soil Mechanics, Cambridge
University Press, 2007.

15. D.M. Potts and L. Zadravkovic, Finite Element Analysis in Geothechnical

Engineering, Thomas Telford, 1999.

16. W. Tiecheng et al., Stress-strain Relation for Concrete Under Triaxial Loading,
16th ASCE Engineering Mechanics Conference, 2003.

17. W.F. Chen, Plasticity in Reinforced Concrete, McGraw-Hill, 1982.

18. N. Ottosen, “A Failure Criterion for Concrete,” J. Eng. Mech. Division, ASCE,
vol. 103, no. 4, pp. 527–535, 1977.

19. N. Ottosen and M. Ristinmaa, The Mechanics of Constitutive Modelling, Elsevier,

2005.

REFERENCES FOR THE GEOMECHANICS MODULE | 59

60 | CHAPTER 2: GEOMECHANICS THEORY
 3

The Geomechanics Materials

The Geomechanics Module has materials that are used in combination with the
Solid Mechanics interface to account for the plasticity of soils and failure criteria in
rocks, concrete, and other materials used in geotechnical applications.

In this chapter:

• Working with the Geomechanics Materials

61
 Working with the Geomechanics
Materials
In this section:

• Adding a Material to a Solid Mechanics Interface

• Plasticity
• Soil Plasticity
• Concrete
• Rocks
• Cam-Clay Material
• Creep

See Geomechanics Theory for background information about these

materials.

Adding a Material to a Solid Mechanics Interface

Creating a New Model in the COMSOL Multiphysics Reference Manual

In general, to add a node, go to the Physics toolbar, no matter what

operating system you are using. However, to add subnodes as is the case
with those described here, you have to right-click the parent node.

1 Add a Solid Mechanics interface from the Structural Mechanics branch ( ) when
adding a physics interface.

62 | CHAPTER 3: THE GEOMECHANICS MATERIALS

2 In the Model Builder, right-click the Solid Mechanics node ( ) to add a Cam-Clay
Material subnode. Or right-click a Linear Elastic Material node to add Soil Plasticity,
Concrete, and Rocks nodes.

Except for the settings described in this section, see The Solid Mechanics
Interface in the Structural Mechanics Module User’s Guide for details.

Plasticity
Right-click the Linear Elastic Material node to add the Plasticity subnode. It adds the
equations for plasticity in order to make it possible to model elastoplastic materials.

Elastoplastic Material Theory

DOMAIN SELECTION
From the Selection list, choose the domains to add plasticity to a linear elastic material

By default, this node inherits the selection from its parent node, and only
a selection that is a subset of the parent node’s selection can be used.

MODEL INPUTS
Use this section to define model inputs such as the temperature field if the material
model uses a temperature-dependent material property. If no model inputs are
required, this section is empty.

PLASTICITY MODEL
Use this section to define the plastic properties of the material.

Plasticity Model
Select Small plastic strains or Large plastic strains to apply either an additive or
multiplicative decomposition between elastic and plastic strains.

Yield Function F
The Yield function F (SI unit: Pa) defines the limit of the elastic regime Fys0.

WO R K I N G W I T H T H E G E O M E C H A N I C S M A T E R I A L S | 63
 Select a Yield function F criterion—von Mises stress (the default), Tresca stress, or User
defined:

• The default is von Mises stress with associate plastic potential.

• Select Tresca stress to use a Tresca yield criterion. The plastic potential is
nonassociated with the von Mises criterion as flow rule.
• If User defined is selected, enter a different value or expression. Write any expression
in terms of the stress tensor variables or its invariants (SI unit: Pa) in the field.

Also select the Plastic potential Q (SI unit: Pa) related to the flow rule—Associated (the
default), von Mises, or User defined (nonassociated). If User defined (SI unit: Pa) is
selected enter a value in the Q field as required.

Initial Yield Stress

For all yield function criteria, the default Initial yield stress ys0 (SI unit: Pa) uses values
From material and represents the stress level where plastic deformation starts. Select
User defined to enter a different value or expression.

Hardening Model
For all yield function criteria, select the type of plasticity from the Hardening model
list—Isotropic or Perfectly plastic.

• Select Perfectly plastic (ideal plasticity) if it is a material that can undergo plastic
deformation without any increase in yield stress.
• If Isotropic is selected, also choose an Isotropic hardening model—Use tangent data
or Use hardening function data.
- If Use tangent data is selected, the default Isotropic tangent modulus ETiso
(SI unit: Pa) uses values From material (if it exists) or User defined, and so the yield
level ys is modified as hardening occurs. The Young’s modulus E is taken from

64 | CHAPTER 3: THE GEOMECHANICS MATERIALS

 the elastic material properties. The yield level ys is related to the effective plastic
strain pe as
E Tiso
 ys =  ys0 + ------------------------  pe
E Tiso
1 – --------------
E

This is a linear isotropic hardening model, so the yield level increases

proportionally to the effective plastic strain pe.

- If Use hardening function data is selected, the Hardening functionh(pe) (SI unit:
Pa) uses values From material (if it exists) or User defined, the yield level ys is
modified as
 ys =  ys0 +  h   pe 

This definition implies that the hardening function must be zero at zero plastic
strain. In other words, ysys0 whenpe0. With this option it is possible to
fit nonlinear isotropic hardening curves depending on stress level, temperature or
other variables.

Soil Plasticity
Right-click the Linear Elastic Material node to add a Soil Plasticity subnode. The Soil
Plasticity node adds the equations for linear elasticity and plasticity and the interface
for defining the elastic material properties and yield surface. The yield criteria are
described in the theory section:

• Drucker-Prager Criterion
• Mohr-Coulomb Criterion
• Matsuoka-Nakai Criterion
• Lade-Duncan Criterion

DOMAIN SELECTION
From the Selection list, choose the domains to compute the displacements, stresses, and
strains.

MODEL INPUTS
Define model inputs, for example, the temperature field of the material uses a
temperature-dependent material property. If no model inputs are required, this section
is empty.

WO R K I N G W I T H T H E G E O M E C H A N I C S M A T E R I A L S | 65
 SOIL PLASTICITY
Select a Yield criterion—Drucker-Prager, Mohr-Coulomb, Matsuoka-Nakai, or
Lade-Duncan.

Drucker-Prager
If required, select the Match to Mohr-Coulomb criterion check box (see Mohr-Coulomb
Criterion). If this check box is selected, the default values for Cohesion c (SI unit: Pa)
and the Angle of internal friction  (SI unit: rad) are taken From material. If User defined
is selected, then enter other values or expressions.

If required, select the Use dilatation angle in plastic potential check box. If this check
box is selected, then enter a value or expression for the Dilatation angle (SI unit: rad).
Alternatively, select From material. The dilatation angle replaces the angle of internal
friction when defining the plastic potential.

If the Match to Mohr-Coulomb criterion check box is not selected, then the default
Drucker-Prager alpha coefficient (dimensionless) and Drucker-Prager k coefficient k
(SI unit: Pa) are taken From material. Select User defined to enter other values or
expressions.

If required, select the Include elliptic cap check box (see Elliptic Cap). Enter values or
expressions for the semi-axes of the Elliptic cap parameter Pa and Elliptic cap parameter
Pb (SI unit: Pa).

Mohr-Coulomb
The default Angle of internal friction  (SI unit: rad) and Cohesion c (SI unit: Pa) are
taken From material. Select User defined to enter other values or expressions.

If required, select the Use dilatation angle in plastic potential check box. If this check
box is selected, then enter a value or expression for the Dilatation angle (SI unit: rad).
Alternatively, select From material. The dilatation angle replaces the angle of internal
friction when defining the plastic potential.

Under Plastic potential select either Drucker-Prager matched at compressive meridian,

Drucker-Prager matched at tensile meridian, or Associated.

If required, select the Include elliptic cap check box (see Elliptic Cap). Enter values or
expressions for the semi-axes of the Elliptic cap parameter Pa and Elliptic cap parameter
Pb (SI unit: Pa).

66 | CHAPTER 3: THE GEOMECHANICS MATERIALS

Matsuoka-Nakai
If required, select the Match to Mohr-Coulomb criterion check box. If this check box is
selected, the default Angle of internal friction  (SI unit: rad) is taken From material.
Select User defined to enter other values or expressions.

If the Match to Mohr-Coulomb criterion check box is not selected, then the default
Matsuoka-Nakai mu coefficient  (dimensionless) is taken From material. Select User
defined to enter other values or expressions.

Lade-Duncan
If required, select the Match to Mohr-Coulomb criterion check box. If this check box is
selected, then enter a value or expression for the Angle of internal friction 
(SI unit: rad). Alternatively, select From material.

If the Match to Mohr-Coulomb criterion check box is not selected, then the default
Lade-Duncan k coefficient k(dimensionless) is taken From material. Select User defined
to enter other values or expressions.

TE N S I O N C U T - O F F
If required, select the Include tension cut-off check box (see Tension Cut-Off). Enter a
value or expression for the Max tensile stress t (SI unit: N/m2). Use this to constrain
the soil plasticity model with an extra yield surface, which limits the maximum principal
stress.

• Deep Excavation: model library path

Geomechanics_Module/Soil_Models/deep_excavation
• Flexible and Smooth Strip Footing on Stratum of Clay: model library
path Geomechanics_Module/Soil_Models/flexible_footing

Concrete
Right-click the Linear Elastic Material node to add a Concrete subnode. The Concrete
node adds the equations for linear elasticity, and the interface for defining the elastic
material properties and failure surface. The failure criteria are described in the theory
section:

• Bresler-Pister Criterion
• Willam-Warnke Criterion

WO R K I N G W I T H T H E G E O M E C H A N I C S M A T E R I A L S | 67
 • Ottosen Criterion

The settings for Domain Selection, Model Inputs, and Linear Elastic Model
are the same as for the Soil Plasticity.

CONCRETE MODEL
Select a Concrete criterion—Bresler-Pister, Willam-Warnke, or Ottosen.

Bresler-Pister
The defaults for the Uniaxial tensile strength t, Uniaxial compressive strength c, and
Biaxial compressive strength b (SI units: Pa) are taken From material. Select User
defined to enter other values or expressions.

Willam-Warnke
The defaults for the Uniaxial tensile strength t, Uniaxial compressive strengthc, and
Biaxial compressive strength b (SI units: Pa) are taken From material. Select User
defined to enter other values or expressions.

Ottosen
The defaults for the Uniaxial tensile strengthc (SI unit: Pa), Ottosen’s parameters a
and b (dimensionless), Size factor k1 (dimensionless), and Shape factor k2
(dimensionless) are taken From material. Select User defined to enter other values or
expressions.

TE N S I O N C U T - O F F
If required, select the Include tension cut-off check box (see Tension Cut-Off). Enter a
value or expression for the Max tensile stress t (SI unit: N/m2). Use this to constraint
the concrete model with an extra yield surface, which limits the maximum principal
stress.

68 | CHAPTER 3: THE GEOMECHANICS MATERIALS

Rocks
Right-click the Linear Elastic Material node to add a Rocks subnode. The Rocks node
adds the equations for linear elasticity, and the interface for defining the elastic material
properties and failure surface. The failure criterion are described in the theory section:

• Original Hoek-Brown Criterion

• Generalized Hoek-Brown Criterion

The settings for Domain Selection, Model Inputs, and Linear Elastic Model
are the same as for the Soil Plasticity.

ROCK MODEL
Select a Rock criterion—Original Hoek-Brown or Generalized Hoek-Brown.

Original Hoek-Brown
The defaults for the Uniaxial compressive strengthc (SI unit: Pa), Hoek-Brown m
parameter m (dimensionless), and Hoek-Brown s parameter s (dimensionless) are taken
From material. Select User defined to enter other values or expressions.

Generalized Hoek-Brown
The defaults for the Uniaxial compressive strengthc (SI unit: Pa), Geological strength
index GSI (dimensionless), Disturbance factor D (dimensionless), and Intact rock
parameter (dimensionless) mi are taken From material. Select User defined to enter
other values or expressions.

Cam-Clay Material
The Cam-Clay Material node adds the equations and interface for defining the material
properties for the modified Cam-clay material (see Theory for the Cam-Clay Material).
Right-click to add the Thermal Expansion and Initial Stress and Strain subnodes.

DOMAIN SELECTION
From the Selection list, choose the domains to compute the displacements, stresses, and
strains.

WO R K I N G W I T H T H E G E O M E C H A N I C S M A T E R I A L S | 69
 MODEL INPUTS
The Pore fluid pressure pfluid (SI unit: Pa) is user-defined by default. The default value
is 1 atm, but you can change it to another value or expression for the pore fluid
pressure.

CAM-CLAY MATERIAL
To use a mixed formulation by adding the negative mean pressure as an extra
dependent variable to solve for, select the Nearly incompressible material check box.

From the Specify list, define the elastic properties either in terms of Poisson’s ratio or
Shear modulus.

The defaults for the Poisson’s ratio (dimensionless) or the Shear modulus G
(SI unit: Pa), Density (SI unit: kg/m3), Cam-Clay M parameter M (dimensionless),
Swelling index (dimensionless), Compression index  (dimensionless), and Void ratio
at reference pressure N (dimensionless) are taken From material. Select User defined for
any of these parameters to enter other values or expressions.

For the Cam-Clay M parameter you can alternatively select Match to

Mohr-Coulomb criterion which then matches the slope of the virgin
consolidation line to the Angle of internal friction. Then select the Angle
of internal friction  (SI unit: rad) as From material (the default) or User
defined.

Enter a value or expression for the Reference pressure for the parameter N prefN
(SI unit: Pa), the Initial void ratio e0,and the Initial consolidation pressure pc0
(SI unit: Pa).

GEOMETRIC NONLINEARITY
Select the Force linear strains check box as required.

Geometric Nonlinearity Theory for the Solid Mechanics Interface in the

Structural Mechanics Module User’s Guide

Isotropic Compression: model library path

Geomechanics_Module/Verification_Models/isotropic_compression

70 | CHAPTER 3: THE GEOMECHANICS MATERIALS

 Thermal Expansion and Initial Stress and Strain in the Structural
Mechanics Module User’s Guide

Creep
Right-click the Linear Elastic Material node to add a Creep subnode to a model. Use the
Creep node to define the material model.

See Creep and Viscoplasticity for background theory about this feature.

DOMAIN SELECTION
From the Selection list, choose the domains to define.

MODEL INPUTS
From the Temperature T (SI unit: K) list, select an existing temperature variable from
a heat transfer interface (for example, Temperature (ht/sol1)), if any temperature
variables exist, or select User defined to enter a value or expression for the temperature
(the default is 293.15 K).

COORDINATE SYSTEM SELECTION

The Global coordinate system is selected by default. The Coordinate system list contains
any additional coordinate systems that the model includes.

CREEP DATA
Select a Material model—Potential, Volumetric, Deviatoric, or User defined.

• If Potential is selected, enter the Rate multiplier (SI unit: 1/s) (the default is 0
(1/s)) and the Creep potential Qcr (SI unit: N/m2) (the default is 0 N/m2).

WO R K I N G W I T H T H E G E O M E C H A N I C S M A T E R I A L S | 71
 • If Deviatoric or Volumetric is selected, enter the Creep rate Fcr (SI unit: 1/s). The
default is 0 1/s.
• If User defined is selected, enter each element for the symmetric Creep rate tensor Fcr
(SI unit: 1/s). The defaults are 0 1/s.

The Energy Dissipation section is also available on the Linear Elastic Material
settings window as an advanced feature (to display the section, click the
Show button ( ) and select Advanced Physics Options).

See Energy Dissipation in the theory section for more information about
this option.

72 | CHAPTER 3: THE GEOMECHANICS MATERIALS

I n d e x
A absolute values 37 disturbance factor 39
alpha coefficient 66 documentation 7
angle of internal friction 20, 26, 66 dolomite 37
associated flow rule 43 Drucker-Prager criterion 22

B back stress 48 ductile materials 19

biaxial compression 34 E effective creep strain rate 56

biaxial data 36 effective plastic strain rate 44
biaxial tension 68 effective stress tensor 38
Bresler-Pister criterion 33 elastic deformation gradient 41

C calcite 37 elastic Green-Lagrange strain tensor 41

Cam-Clay material (node) 69 elastic right Cauchy-Green tensor 42

Cam-clay model 28 elastic volumetric strain variable 31

carbonate rocks 37 elastoplastic materials 16, 34

Cauchy stress tensor 13 elliptic cap 24

Cayley-Hamilton theorem 13 elplastic 43, 50

ceramics 34 emailing COMSOL 9

circle, Mohr 20 F failure surfaces 16

clays 19 flow rule 42
cohesion 20, 66 fluid pore pressure 25
cohesionless soils 26 friction, angle 20, 26, 66
cohesive-frictional materials 34 frictionless materials 21, 23
complementarity 43
G Gauss points 43, 50
compressive meridians 14, 21, 34–35
generalized Hoek-Brown criterion 38
compressive stresses 12
geological strength index (GSI) 38
concrete 33, 35
concrete (node) 67 H Haigh–Westergaard coordinates 14

consistency parameter 42 hardening function data 47

continuum mechanics 12 hardening models

creep (node) 71 plasticity node settings 64

creep dissipation rate density 56 theory 46

creep strain rate 52 Hencky plastic strain 50

crystal cleavage 37 hexagonal prism 18

Hoek-Brown criterion 37
D deviatoric creep 56
hydrostatic axis 14
deviatoric stress 13, 17
hydrostatic pressure 12, 20
dissipated energy density 56
hydrostatic stress 26

INDEX| 73
 I ideal plasticity 46 mu coefficient 67
initial yield stress 47 multiaxial stress states 33
Internet resources 6
N nearly incompressible material 70
isotropic hardening 46
non-associated flow rule 43
isotropic materials 17
isotropic plasticity 43 O octahedral normal stress 33

isotropic rocks 38 octahedral plane 15

Ottosen criterion 35
K k coefficient 67
over-consolidation pressure 29
kinematic hardening 47
knowledge base, COMSOL 9 P perfectly elastoplastic materials 17

Kuhn-Tucker conditions 43 perfectly plastic hardening 46

plastic deformation gradient 42
L Lade-Duncan criterion 26
plastic element 43, 50
large plastic strain 41
plastic flow rule 50
Lie derivative 49
plastic Green-Lagrange strain 42
limestone 37
plastic multiplier 42
Lode angle 14
plastic potential 42
Hoek-Brown criterion 38
plastic right Cauchy-Green tensor 42
Mohr-Coulomb criterion 21
plasticity (node) 63
Tresca criterion 18
plasticity models 63
Willam-Warnke criterion 34
plasticity, theory 12
logarithmic plastic strain 50
Prager’s hardening rule 48
M marble 37 primary creep 52
Matsuoka-Nakai criterion 25 principal stresses 15
meridians, tensile and compressive 14,
R rate independent plasticity 42
21, 34–35
rate of strain tensor 49
metal plasticity 44
rock mass 39
metals 19, 22
rock types 37
mixed formulation 70
rock yield criteria 20
mobilized planes 25
rocks 37
Model Libraries window 8
rocks (node) 69
model library examples
Cam-Clay material 70 S saturated cohesive soils 19

soil plasticity 67 secondary creep 53

modified Cam-clay model 28 shape factors 36

modified Mohr-Coulomb criterion 27 shear stresses 14, 20, 33

Mohr-Coulomb criterion 20 sign convention 12

MPH-files 8 size factors 36

small plastic strain 40

74 | I N D E X
 soil deformation theory 30 W web sites, COMSOL 9
soil plasticity 14 Willam-Warnke criterion 34
soil plasticity (node) 65 Y yield function 45
soil yield criteria 20 yield functions 17, 20
spatially mobilized planes (STP) 25 yield stress levels 64
steady-state creep 53 yield surface 43
stress invariants 13 yield surfaces 16, 22
T tangent data 46
tangent modulus 47
technical support, COMSOL 9
tensile meridian 14
tensile meridians 21–22, 26, 34–35
tensile normal strains and stresses 12
tension cut-off 27
tertiary creep 53
transient creep 52
Tresca effective stress variable 18
Tresca stress 45
Tresca yield criterion 18
triaxial conditions 34
triaxial data 36

U undrained shear strength 18

uniaxial compression 36–37
uniaxial compressive strength 33
uniaxial tension 17, 68
user community, COMSOL 9

V variables
effective creep strain 56
effective plastic strain 45, 47
elastic volumetric strain 31
Tresca effective stress 18
volumetric plastic strain 29, 45
volumetric creep strain rate 55
volumetric plastic strain rate 44
volumetric plastic strain variable 29
von Mises criterion 17, 44
von Mises stress 45

INDEX| 75
76 | I N D E X