Geomechanics Module

Application Library Manual

Geomechanics Module Application Library Manual
© 1998–2017 COMSOL
Protected by U.S. Patents listed on [Link]/patents, and U.S. Patents 7,519,518; 7,596,474;
7,623,991; 8,457,932; 8,954,302; 9,098,106; 9,146,652; 9,323,503; 9,372,673; and 9,454,625. Patents
pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement ([Link]/comsol-license-agreement) and may be used or copied only under the terms
of the license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, Capture the Concept, COMSOL Desktop,
LiveLink, and COMSOL Server 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]/trademarks.
Version: COMSOL 5.3

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]/product-download
• Product Updates: [Link]/support/updates
• COMSOL Blog: [Link]/blogs
• Discussion Forum: [Link]/community
• Events: [Link]/events
• COMSOL Video Gallery: [Link]/video
• Support Knowledge Base: [Link]/support/knowledgebase

Part number: CM021802

 Created in COMSOL Multiphysics 5.3

Block Verification

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Introduction
This example shows how to set up a compression test on a prestressed soil sample. Due to
a simple stress state, it is possible to determine the vertical yield stress analytically. The soil
sample is modeled with soil plasticity and the Mohr-Coulomb criterion.

Model Definition
In this example, we consider a block of soil of 1 m length on each side. The soil is pressed
from the sides by boundary loads in the x and y directions, and from the top by a
prescribed displacement.
Prescribed displacement

Boundary load

Roller

1m

Figure 1: Dimensions, boundary conditions, and loads for the test.

ELASTIC PROPERTIES
The soil properties are taken from a standard clay.

• Young’s modulus, E = 207 MPa, and Poisson’s ratio ν = 0.3.

SOIL PLASTICITY
• Cohesion c = 70 kPa, and angle of internal friction φ = 30° .
• Use the Mohr-Coulomb criterion, with non-associated flow rule, and Drucker-Prager
matched at compressive meridian as plastic potential.

2 | B L O C K VE R I F I C A T I O N
CONSTRAINTS AND LOADS
• The test is based on uniaxial compression. Fix the normal displacement at the lower, left,
and right boundaries with a roller boundary condition (these are the boundaries at
x = 0, y = 0, and z = 0).
• On the other two vertical walls (x = 1 m and y = 1 m), apply boundary loads of 300 kPa
and 200 kPa.
• Prescribe the in-situ stress via the External Stress node.
• The soil sample is subjected to a loading at the top boundary. By means of a parametric
sweep, gradually increase the displacement up to 8 mm.

Results and Discussion

The soil cube experiences a homogeneous stress state. This is seen in the in the default plot
group shows the von Mises stress.

Figure 2: Effective stress and deformation in the soil sample after applying 8 mm displacement
from the top.

From the Mohr circle, the Mohr-Coulomb criterion can be written in terms of the biggest
and smallest principal stress:

3 | B L O C K VE R I F I C A T I O N
 1 1
--- ( σ 1 – σ 3 ) + --- ( σ 1 + σ 3 ) sin φ – c cos φ = 0
2 2

Since stress in the y direction is the largest principal stress and the stress in z direction is
the smallest principle stress at the onset of yielding, its analytical value can be obtained.
Manipulation of the above formula gives

2c cos φ – σ YY ( 1 + sin φ )
σ ZZ = --------------------------------------------------------------- (1)
( sin φ – 1 )

The stresses history together with the analytical value of the stress in the z direction at the
onset of yielding is shown in Figure 3. Plasticity is reached after deflection of about
3.5 mm.

Figure 3: This plot shows how the soil sample behaves elastically until it reaches the yield surface
at the compressive meridian.

Application Library path: Geomechanics_Module/Verification_Examples/

block_verification

4 | B L O C K VE R I F I C A T I O N
Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Stationary.
6 Click Done.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

Disp 0 0
X_stress -3e5[Pa] -3E5 Pa In-situ stress, xx component
Y_stress -2e5[Pa] -2E5 Pa In-situ stress, yy component
Z_stress -1e5[Pa] -1E5 Pa In-situ stress, zz component

In-situ stresses are set with negative sign to fit the structural mechanics convention
which assumes negative stresses in compression, and positive in tension.

GEOMETRY 1

Block 1 (blk1)
1 On the Geometry toolbar, click Block.
2 In the Settings window for Block, click Build All Objects.
The geometry consists in a simple unit block.

5 | B L O C K VE R I F I C A T I O N
SOLID MECHANICS (SOLID)

Linear Elastic Material 1

In the Model Builder window, expand the Component 1 (comp1)>Solid Mechanics (solid)
node, then click Linear Elastic Material 1.

Soil Plasticity 1
1 On the Physics toolbar, click Attributes and choose Soil Plasticity.
2 In the Settings window for Soil Plasticity, locate the Soil Plasticity section.
3 From the Yield criterion list, choose Mohr-Coulomb.
The Mohr-Coulomb criterion is used as yield surface. Use the non-associated Drucker-
Prager plastic potential matched to the Mohr-Coulomb model at the compressive
meridian.

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click
Linear Elastic Material 1.

External Stress 1
1 On the Physics toolbar, click Attributes and choose External Stress.
2 In the Settings window for External Stress, locate the External Stress section.
3 From the list, choose Diagonal.
4 In the Sext table, enter the following settings:

X_stress 0 0
0 Y_stress 0
0 0 Z_stress

Roller 1
1 On the Physics toolbar, click Boundaries and choose Roller.
2 Select Boundaries 1–3 only.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 6 only.
3 In the Settings window for Boundary Load, locate the Force section.

6 | B L O C K VE R I F I C A T I O N
4 Specify the FA vector as

X_stress x
0 y
0 z

Boundary Load 2
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 5 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

0 x
Y_stress y
0 z

Prescribed Displacement 1
1 On the Physics toolbar, click Boundaries and choose Prescribed Displacement.
2 Select Boundary 4 only.
3 In the Settings window for Prescribed Displacement, locate the Prescribed Displacement
section.
4 Select the Prescribed in z direction check box.
5 In the u 0z text field, type -Disp.
Roller conditions are assumed on Boundaries 1, 2, and 3. Boundaries 4 and 5 are
loaded, in-situ stresses are defined, and finally a prescribed displacement is applied at the
top boundary.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.

7 | B L O C K VE R I F I C A T I O N
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 207e6 Pa Basic
Poisson’s ratio nu 0.3 1 Basic
Density rho 2000 kg/m³ Basic
Cohesion cohesion 70e3 Pa Mohr-Coulomb
Angle of internalphi 30[deg] rad Mohr-Coulomb
internal friction

MESH 1

Mapped 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
More Operations>Mapped.
2 Select Boundary 4 only.

Size
1 In the Model Builder window, right-click Mesh 1 and choose Swept.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Extra coarse.
4 Click Build All.

STUDY 1

Step 1: Stationary
Set up an auxiliary continuation sweep for the Disp parameter.

1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

Disp 8[mm]*range(0,0.05,1)

Use continuation to model the displacement at the top surface. The prescribed
displacement goes from 0 to 8 mm in 20 steps.

8 | B L O C K VE R I F I C A T I O N
6 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
The default plot shows the von Mises stress for the final step.

1 In the Model Builder window, under Results click Stress (solid).

2 On the Stress (solid) toolbar, click Plot.
3 Click the Zoom Extents button on the Graphics toolbar.

Add a 1D plot to show the evolution of stress tensor components versus the displacement
at the top surface.

1D Plot Group 2
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Stresses vs Displacement in the
Label text field.
3 Click to expand the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Stress vs vertical displacement.
5 Locate the Plot Settings section. Select the x-axis label check box.
6 In the associated text field, type Vertical displacement (mm).
7 Select the y-axis label check box.
8 In the associated text field, type Stress (kPa).

Point Graph 1
1 Right-click Stresses vs Displacement and choose Point Graph.
2 Select Point 1 only.
3 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Model>Component 1>
Solid Mechanics>Stress>Stress tensor (spatial frame)>[Link] - Stress tensor,
x component.
4 Locate the y-Axis Data section. From the Unit list, choose kPa.
5 Locate the x-Axis Data section. From the Parameter list, choose Expression.
6 In the Expression text field, type Disp*1000.
7 Click to expand the Coloring and style section. Click to expand the Legends section.
Select the Show legends check box.

9 | B L O C K VE R I F I C A T I O N
8 From the Legends list, choose Manual.
9 In the table, enter the following settings:

Legends
sxx stress

10 On the Stresses vs Displacement toolbar, click Plot.

Point Graph 2
1 Right-click Results>Stresses vs Displacement>Point Graph 1 and choose Duplicate.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Model>Component 1>
Solid Mechanics>Stress>Stress tensor (spatial frame)>[Link] - Stress tensor,
y component.
3 On the Stresses vs Displacement toolbar, click Plot.
4 Locate the Legends section. In the table, enter the following settings:

Legends
syy stress

Point Graph 3
1 Right-click Point Graph 1 and choose Duplicate.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Model>Component 1>
Solid Mechanics>Stress>Stress tensor (spatial frame)>[Link] - Stress tensor,
z component.
3 On the Stresses vs Displacement toolbar, click Plot.
4 Locate the Legends section. In the table, enter the following settings:

Legends
szz stress

Point Graph 4
1 Right-click Point Graph 1 and choose Duplicate.
The analytical yield level is given by Equation 1.
2 In the Settings window for Point Graph, locate the y-Axis Data section.
3 In the Expression text field, type (2*[Link]*cos([Link])-
Y_stress*(1+sin([Link])))/(sin([Link])-1).

10 | B L O C K VE R I F I C A T I O N
4 Locate the Coloring and Style section. Find the Line style subsection. From the Line list,
choose Dashed.
5 From the Color list, choose Black.
6 Locate the Legends section. In the table, enter the following settings:

Legends
theoretical yield stress

7 On the Stresses vs Displacement toolbar, click Plot.

11 | B L O C K VE R I F I C A T I O N
12 | B L O C K VE R I F I C A T I O N
 Created in COMSOL Multiphysics 5.3

Concrete Beam With Reinforcement Bars

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Introduction
Concrete structures almost always contain reinforcements in the shape of steel bars
(rebars). In COMSOL Multiphysics, individual rebars can be modeled by adding a Truss
interface to the Solid Mechanics interface used for the concrete beam. The solid mesh for
the concrete and the rebar mesh can be independent of each other, since the displacements
are mapped from within the solids onto the rebars.

Model Definition
This example shows how to include steel reinforcement that is much smaller than the
geometrical dimensions of the concrete structure. The truss interface is used to model the
steel reinforcements instead of a 3D solid. This removes the necessity to model the bars
with small mesh which saves computational time. The geometry of the concrete beam is
given in Figure 1.

Boundary load

Rigid connector

Rigid connector
20 cm
4m

30 cm

Figure 1: The concrete beam is 30 cm in width, 20 cm in height, and 4 meters in length. Due
to symmetry, only half of its width is modeled.

In the example, most dimensions such as height, width, and length of the concrete
structure are parameterized. The number of rebar layers is also given by a parameter, and
the number of rebars per layer is calculated from the spacing in width dimensions and the
minimal distance from the lateral faces. In this example, six steel bars 10 mm in diameter

2 | CONCRETE BEAM WITH REINFORCEMENT BARS

are placed in four parallel layers along the concrete beam. See Figure 2.

Figure 2: A mapped mesh of 6 by 6 elements is swept through the length of the concrete beam.
One hundred elements are used for each reinforcement bar.

In this example the effect of the gravity load is simulated. Moreover a deflection due to a
vertical boundary load is ramped up to 20 kN/m2 by means of a parametric sweep.

Results and Discussion

Three different studies are done. The first study models the concrete beam as an isotropic
elastic material, the second study adds the reinforcements bars, and the third study
includes the effect of plastic deformation in the concrete, modeled using the Ottosen

3 | CONCRETE BEAM WITH REINFORCEMENT BARS

criterion. Figure 3 shows the comparison for the vertical displacement of the three studies.

Figure 3: Deflection along the top surface of the concrete beam due to gravity and external
load.

The simulation shows how force is transferred from the concrete beam to its steel
reinforcement bars. Figure 4 shows effective stresses in the linear elastic model, Figure 5
shows the stress distribution in the reinforced linear elastic concrete. Compare both figures
with each other and notice the change in stress level of the concrete once the bars are
added. Figure 7 shows axial stresses in the bars. It is clear that the compressive stresses are
experienced in the upper bars and tensile in the lower bars. Figure 8 shows the von Mises
stresses in the concrete beam with Ottosen criterion.

4 | CONCRETE BEAM WITH REINFORCEMENT BARS

Figure 4: von Mises stress in a linear elastic beam.

Figure 5: von Mises stress in a linear elastic beam after adding the reinforcement bars.

5 | CONCRETE BEAM WITH REINFORCEMENT BARS

Figure 6: Axial force in the reinforcements bars.

Figure 7: Axial stress in the steel rebars.

6 | CONCRETE BEAM WITH REINFORCEMENT BARS

Figure 8: von Mises stress in the reinforced beam after adding the Ottosen criterion for the
concrete.

Notes About the COMSOL Implementation

Since steel reinforcements are relatively thin compared to the concrete structures, it is
assumed that they are only capable of transmitting axial forces. The bending stiffness of
each bar does not contribute much to the overall total bending stiffness of the section,
therefore the reinforcement bars are modeled with truss elements instead of beam
elements.

In civil engineering it is also common practice that the rebars are pretensioned, but this
effect is not included in the example. However it can easily be incorporated by adding
initial strain in the trusses.

In this example, the concrete is “glued” to the steel rebars, so bonding effects are not
included.

The connection technique used in this example works well as long as the total stiffness of
the rebars is smaller than the stiffness contribution from the concrete. Also, the size of the
solid elements should be significantly larger than the physical volume occupied by the
rebars passing through them. A very refined mesh would actually show stresses in the
solids that increase without bounds where the rebars are attached.

7 | CONCRETE BEAM WITH REINFORCEMENT BARS

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

Application Library path: Geomechanics_Module/Tutorials/concrete_beam

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 In the Select Physics tree, select Structural Mechanics>Truss (truss).
5 Click Add.
6 Click Study.
7 In the Select Study tree, select Preset Studies for Selected Physics Interfaces>Stationary.
8 Click Done.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 Click Load from File.
4 Browse to the model’s Application Libraries folder and double-click the file
concrete_beam_parameters.txt.

GEOMETRY 1

Block 1 (blk1)
1 On the Geometry toolbar, click Block.

8 | CONCRETE BEAM WITH REINFORCEMENT BARS

2 In the Settings window for Block, locate the Size and Shape section.
3 In the Width text field, type length.
4 In the Depth text field, type width/2.
5 In the Height text field, type height.

Bézier Polygon 1 (b1)

1 On the Geometry toolbar, click More Primitives and choose Bézier Polygon.
2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 2, set x to length.
5 In row 1, set y to (bars_across_width-1)/2*width_spacing.
6 In row 2, set y to (bars_across_width-1)/2*width_spacing.
7 In row 1, set z to layer_spacing_first.
8 In row 2, set z to layer_spacing_first.
9 Locate the Selections of Resulting Entities section. Click New.
10 In the New Cumulative Selection dialog box, type bars_inhalf in the Name text field.
11 Click OK.

Array 1 (arr1)
1 On the Geometry toolbar, click Transforms and choose Array.
2 In the Settings window for Array, locate the Input section.
3 From the Input objects list, choose bars_inhalf.
4 Locate the Size section. In the y size text field, type floor(bars_across_width/2).
5 Locate the Displacement section. In the y text field, type -width_spacing.

Bézier Polygon 2 (b2)

1 On the Geometry toolbar, click More Primitives and choose Bézier Polygon.
2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 2, set x to length.
5 In row 1, set z to layer_spacing_first.
6 In row 2, set z to layer_spacing_first.
7 Locate the Selections of Resulting Entities section. Click New.
8 In the New Cumulative Selection dialog box, type bars_midplane in the Name text field.

9 | CONCRETE BEAM WITH REINFORCEMENT BARS

9 Click OK.

Array 2 (arr2)
1 On the Geometry toolbar, click Transforms and choose Array.
2 Select the objects arr1(1,3,1), arr1(1,1,1), b2, and arr1(1,2,1) only.
3 In the Settings window for Array, locate the Size section.
4 In the z size text field, type bar_layers.
5 Locate the Displacement section. In the z text field, type layer_spacing.

Mirror 1 (mir1)
1 On the Geometry toolbar, click Transforms and choose Mirror.
2 Select all the six bars.

3 In the Settings window for Mirror, locate the Point on Plane of Reflection section.
4 In the z text field, type height/2.
5 Click Build All Objects.
6 Locate the Input section. Select the Keep input objects check box.
7 Click Build All Objects.

Form Union (fin)

1 In the Model Builder window, under Component 1 (comp1)>Geometry 1 click
Form Union (fin).

10 | CONCRETE BEAM WITH REINFORCEMENT BARS

2 In the Settings window for Form Union/Assembly, locate the Form Union/Assembly section.
3 From the Action list, choose Form an assembly.
4 Clear the Create pairs check box.

Bézier Polygon 2 (b2)

Bars in plane symmetry must be created only in case of odd number of bars. Add a if
condition to the creation of the first bar.

If 1 (if1)
1 On the Geometry toolbar, click Programming and choose Add Before Selected>If.
2 In the Settings window for If, locate the If section.
3 In the Condition text field, type mod(bars_across_width,2)==1.

End If 1 (endif1)
1 On the Geometry toolbar, click Programming and choose Add After Selected>End If.
2 In the Settings window for End If, click Build All Objects.

DEFINITIONS

Union 1
1 On the Definitions toolbar, click Union.
2 In the Settings window for Union, type bars in the Label text field.
3 Locate the Geometric Entity Level section. From the Level list, choose Edge.
4 Locate the Input Entities section. Under Selections to add, click Add.
5 In the Add dialog box, In the Selections to add list, choose bars_inhalf and bars_midplane.
6 Click OK.

To make the displacements in the beam available for the bars, use a general extrusion
operator.

General Extrusion 1 (genext1)

1 On the Definitions toolbar, click Component Couplings and choose General Extrusion.
2 Click in the Graphics window and then press Ctrl+A to select all domains.
3 In the Settings window for General Extrusion, locate the Destination Map section.
4 In the x-expression text field, type X.
5 In the y-expression text field, type Y.
6 In the z-expression text field, type Z.
7 Locate the Source section. From the Source frame list, choose Material (X, Y, Z).

11 | CONCRETE BEAM WITH REINFORCEMENT BARS

Explicit 1
On the Definitions toolbar, click Explicit.

SOLID MECHANICS (SOLID)

Linear Elastic Material 1

To model the failure of the material, add a material model to the Solid Mechanics interface.

Concrete 1
1 On the Physics toolbar, click Attributes and choose Concrete.
2 In the Settings window for Concrete, locate the Concrete Model section.
3 From the Concrete criterion list, choose Ottosen.

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundary 2 only.

Rigid Connector 1
1 On the Physics toolbar, click Boundaries and choose Rigid Connector.
2 In the Settings window for Rigid Connector, type Rigid Connector Left in the Label
text field.
3 Select Boundary 1 only.
4 Locate the Prescribed Displacement at Center of Rotation section. Select the
Prescribed in x direction check box.
5 Select the Prescribed in y direction check box.
6 Select the Prescribed in z direction check box.
7 Locate the Prescribed Rotation section. From the By list, choose Constrained rotation.
8 Select the Constrain rotation around x-axis check box.
9 Select the Constrain rotation around z-axis check box.

Rigid Connector Left 1

1 Right-click Rigid Connector Left and choose Duplicate.
2 In the Settings window for Rigid Connector, type Rigid Connector Right in the Label
text field.
3 Select Boundary 6 only.

Gravity 1
1 On the Physics toolbar, click Domains and choose Gravity.

12 | CONCRETE BEAM WITH REINFORCEMENT BARS

2 In the Settings window for Gravity, locate the Domain Selection section.
3 From the Selection list, choose All domains.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 4 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

0 x
0 y
-force_area*para z

TR U S S ( T R U S S )
1 In the Model Builder window, under Component 1 (comp1) click Truss (truss).
2 In the Settings window for Truss, locate the Edge Selection section.
3 From the Selection list, choose bars.
Set the truss discretization to quadratic to fit with the solid. To do so, you first have to
enable the discretization section.

4 In the Model Builder window’s toolbar, click the Show button and select Discretization in
the menu.
5 In the Model Builder window, click Truss (truss).
6 In the Settings window for Truss, click to expand the Discretization section.
7 From the Displacement field list, choose Quadratic.

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Truss (truss) click
Linear Elastic Material 1.

Cross Section Data 1

1 On the Physics toolbar, click Attributes and choose Plasticity.
Set the cross section area to half for the bars on midplane.
2 In the Model Builder window, under Component 1 (comp1)>Truss (truss) click
Cross Section Data 1.
3 In the Settings window for Cross Section Data, locate the Cross Section Data section.
4 In the A text field, type pi*(diam_bar/2)^2*(0.5+0.5*Y>0.1[mm]).

13 | CONCRETE BEAM WITH REINFORCEMENT BARS

Prescribed Displacement 1
1 On the Physics toolbar, click Edges and choose Prescribed Displacement.
Use the general extrusion operator to prescribe the displacements of the bars.

2 In the Settings window for Prescribed Displacement, locate the Edge Selection section.
3 From the Selection list, choose bars_inhalf.
4 Locate the Prescribed Displacement section. Select the Prescribed in x direction check
box.
5 In the u 0x text field, type genext1(u).
6 Select the Prescribed in y direction check box.
7 In the u 0y text field, type genext1(v).
8 Select the Prescribed in z direction check box.
9 In the u 0z text field, type genext1(w).

Because the bar displacements are prescribed, the feature Straight Edge Constraint 1 should
be disabled.

Straight Edge Constraint 1

In the Model Builder window, under Component 1 (comp1)>Truss (truss) right-click
Straight Edge Constraint 1 and choose Disable.

Gravity 1
1 On the Physics toolbar, click Edges and choose Gravity.
2 In the Settings window for Gravity, locate the Edge Selection section.
3 From the Selection list, choose bars_inhalf.

ADD MATERIAL
1 On the Home toolbar, click Add Material to open the Add Material window.
2 Go to the Add Material window.
3 In the tree, select Built-In>Concrete.
4 Click Add to Component in the window toolbar.

MATERIALS

Concrete (mat1)
1 In the Model Builder window, under Component 1 (comp1)>Materials click
Concrete (mat1).
2 In the Settings window for Material, locate the Material Contents section.

14 | CONCRETE BEAM WITH REINFORCEMENT BARS

3 In the table, enter the following settings:

Property Name Value Unit Property

group
Uniaxial compressive strength sigmauc 20e6 Pa Yield stress
parameters
Ottosen a parameter aOttosen 1.3 1 Ottosen
Ottosen b parameter bOttosen 3.2 1 Ottosen
Size factor k1Ottosen 11.8 1 Ottosen
Shape factor k2Ottosen 0.98 1 Ottosen

ADD MATERIAL
1 Go to the Add Material window.
2 In the tree, select Built-In>Structural steel.
3 Click Add to Component in the window toolbar.
4 On the Home toolbar, click Add Material to close the Add Material window.

MATERIALS

Structural steel (mat2)

1 In the Model Builder window, under Component 1 (comp1)>Materials click
Structural steel (mat2).
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Edge.
4 From the Selection list, choose bars.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Initial yield stress sigmags 100[MPa] Pa Elastoplastic material model
Isotropic tangent Et 20[GPa] Pa Elastoplastic material model
modulus

TR U S S ( T R U S S )

Gravity 1
Because gravity is already applied on the concrete volume occupied by steel, substract
concrete density contribution from steel gravity.

15 | CONCRETE BEAM WITH REINFORCEMENT BARS

1 In the Model Builder window, under Component 1 (comp1)>Truss (truss) click Gravity 1.
2 In the Settings window for Gravity, locate the Gravity section.
3 Specify the g vector as

0 x
0 y
-g_const*([Link]/[Link]) z

MESH 1

Edge 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
More Operations>Edge.
2 In the Settings window for Edge, locate the Edge Selection section.
3 From the Selection list, choose bars.

Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Edge 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Distribution section.
3 In the Number of elements text field, type 100.

Mapped 1
1 In the Model Builder window, right-click Mesh 1 and choose More Operations>Mapped.
2 Select Boundary 1 only.

Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Mapped 1 and choose Distribution.
2 Select Edges 1 and 4 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 6.

Distribution 1
1 In the Model Builder window, right-click Mesh 1 and choose Swept.
2 Right-click Swept 1 and choose Distribution.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 40.
5 In the Model Builder window, click Mesh 1.

16 | CONCRETE BEAM WITH REINFORCEMENT BARS

6 In the Settings window for Mesh, click Build All.
7 Click the Transparency button on the Graphics toolbar.
8 Click Go to Default View.
The mesh should look like the one in Figure 2.
9 Click the Transparency button on the Graphics toolbar again to restore the default
transparency state.

STUDY 1
The first study solves only the linear elastic problem in the concrete beam without the
reinforcement bars.

1 In the Model Builder window, click Study 1.

2 In the Settings window for Study, type Without Bars in the Label text field.

WITHOUT BARS

Step 1: Stationary
1 In the Model Builder window, under Without Bars click Step 1: Stationary.
2 In the Settings window for Stationary, locate the Physics and Variables Selection section.
3 Select the Modify physics tree and variables for study step check box.
4 In the Physics and variables selection tree, select Component 1 (comp1)>
Solid Mechanics (solid)>Linear Elastic Material 1>Concrete 1.
5 Click Disable.
6 In the Physics and variables selection tree, select Component 1 (comp1)>Truss (truss).
7 Click Disable.
8 Click to expand the Study extensions section. Locate the Study Extensions section. Select
the Auxiliary sweep check box.
9 Click Add.
10 Click to select row number 1 in the table.
11 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

para range(0,0.1,1)

12 On the Home toolbar, click Compute.

17 | CONCRETE BEAM WITH REINFORCEMENT BARS

RESULTS

Mirror 3D 1
Add a mirror data set to plot the entire beam.

1 On the Results toolbar, click More Data Sets and choose Mirror 3D.
2 In the Settings window for Mirror 3D, locate the Plane Data section.
3 From the Plane list, choose ZX-planes.

Stress (solid)
1 In the Model Builder window, under Results click Stress (solid).
2 In the Settings window for 3D Plot Group, type Stress Without Bars in the Label text
field.
3 Locate the Data section. From the Data set list, choose Mirror 3D 1.
4 Locate the Plot Settings section. Clear the Plot data set edges check box.

Surface 1
1 In the Model Builder window, expand the Results>Stress Without Bars node, then click
Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 From the Unit list, choose MPa.
4 On the Stress Without Bars toolbar, click Plot.
5 Click Go to Default View.

ADD STUDY
Add a second study to solve the model with the reinforcement bars.

1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select Study tree, select Preset Studies.
4 In the Select Study tree, select Preset Studies>Stationary.
5 Click Add Study in the window toolbar.

STUDY 2

Step 1: Stationary
1 On the Home toolbar, click Add Study to close the Add Study window.
2 In the Model Builder window, under Study 2 click Step 1: Stationary.

18 | CONCRETE BEAM WITH REINFORCEMENT BARS

3 In the Settings window for Stationary, locate the Physics and Variables Selection section.
4 Select the Modify physics tree and variables for study step check box.
5 In the Physics and variables selection tree, select Component 1 (comp1)>
Solid Mechanics (solid)>Linear Elastic Material 1>Concrete 1.
6 Click Disable.
7 Locate the Study Extensions section. Select the Auxiliary sweep check box.
8 Click Add.
9 Click to select row number 1 in the table.
10 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

para range(0,0.1,1)

Solution 2 (sol2)
1 On the Study toolbar, click Show Default Solver.
This problem is better solved fully coupled.
2 In the Model Builder window, expand the Solution 2 (sol2) node.
3 Right-click Stationary Solver 1 and choose Fully Coupled.
4 In the Model Builder window, click Study 2.
5 In the Settings window for Study, type With Bars in the Label text field.
6 On the Study toolbar, click Compute.

RESULTS

Mirror 3D 2
1 On the Results toolbar, click More Data Sets and choose Mirror 3D.
2 In the Settings window for Mirror 3D, locate the Data section.
3 From the Data set list, choose With Bars/Solution 2 (sol2).
4 Locate the Plane Data section. From the Plane list, choose ZX-planes.
5 In the Y-coordinate text field, type -1e-10.

Stress (solid)
The first default plot shows the von Mises stress, Figure 5. This result can be compared to
the result without reinforcement bars, Figure 4.

1 In the Model Builder window, under Results click Stress (solid).

19 | CONCRETE BEAM WITH REINFORCEMENT BARS

2 In the Settings window for 3D Plot Group, type Stress With Bars in the Label text field.
3 Locate the Data section. From the Data set list, choose Mirror 3D 2.
4 Locate the Plot Settings section. Clear the Plot data set edges check box.

Surface 1
1 In the Model Builder window, expand the Results>Stress With Bars node, then click
Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 From the Unit list, choose MPa.
4 On the Stress With Bars toolbar, click Plot.
5 Click the Zoom Extents button on the Graphics toolbar.

Force (truss)
The second default plot shows the force in bars, Figure 6.

1 In the Model Builder window, under Results click Force (truss).

2 In the Settings window for 3D Plot Group, type Force With Bars (truss) in the Label
text field.
3 Locate the Data section. From the Data set list, choose Mirror 3D 2.
4 Locate the Plot Settings section. Clear the Plot data set edges check box.

Line
The force of rebars in the midplane must be multiplied by 2.

1 In the Model Builder window, expand the Results>Force With Bars (truss) node, then click
Line.
2 In the Settings window for Line, locate the Expression section.
3 In the Expression text field, type (1+(y==0))*[Link].
4 Locate the Coloring and Style section. In the Tube radius expression text field, type
diam_bar/2.

5 Select the Radius scale factor check box.

6 In the associated text field, type 1.
7 From the Color table list, choose WaveLight.
8 Select the Symmetrize color range check box.
9 On the Force With Bars (truss) toolbar, click Plot.
10 Click Go to Default View.

20 | CONCRETE BEAM WITH REINFORCEMENT BARS

Stress (truss)
The third default plot shows the axial stress in bars, Figure 7.

1 In the Model Builder window, under Results click Stress (truss).

2 In the Settings window for 3D Plot Group, type Stress With Bars (truss) in the Label
text field.
3 Locate the Data section. From the Data set list, choose Mirror 3D 2.
4 Locate the Plot Settings section. Clear the Plot data set edges check box.

Line
1 In the Model Builder window, expand the Results>Stress With Bars (truss) node, then
click Line.
2 In the Settings window for Line, locate the Expression section.
3 From the Unit list, choose MPa.
4 Locate the Coloring and Style section. In the Tube radius expression text field, type
diam_bar/2.

5 Select the Radius scale factor check box.

6 In the associated text field, type 1.
7 From the Color table list, choose WaveLight.
8 Select the Symmetrize color range check box.
9 On the Stress With Bars (truss) toolbar, click Plot.
10 Click Go to Default View.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select Study tree, select Preset Studies.
4 In the Select Study tree, select Preset Studies>Stationary.
5 Click Add Study in the window toolbar.

STUDY 3

Step 1: Stationary
1 On the Home toolbar, click Add Study to close the Add Study window.
Set up an auxiliary continuation sweep for the para parameter.
2 In the Model Builder window, under Study 3 click Step 1: Stationary.

21 | CONCRETE BEAM WITH REINFORCEMENT BARS

3 In the Settings window for Stationary, locate the Study Extensions section.
4 Select the Auxiliary sweep check box.
5 Click Add.
6 Click to select row number 1 in the table.
7 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

para range(0,0.1,1)

8 In the Model Builder window, click Study 3.

9 In the Settings window for Study, locate the Study Settings section.
10 Clear the Generate default plots check box.

Solution 3 (sol3)
1 On the Study toolbar, click Show Default Solver.
2 In the Model Builder window, expand the Solution 3 (sol3) node.
3 Right-click Stationary Solver 1 and choose Fully Coupled.
4 In the Model Builder window, under Study 3>Solver Configurations>Solution 3 (sol3)>
Stationary Solver 1 click Parametric 1.
5 In the Settings window for Parametric, click to expand the Continuation section.
6 From the Predictor list, choose Constant to improve the convergence for the elastoplastic
case.
7 In the Model Builder window, click Study 3.
8 In the Settings window for Study, type With Bars and Ottosen in the Label text field.
9 On the Study toolbar, click Compute.

RESULTS

Mirror 3D 3
1 On the Results toolbar, click More Data Sets and choose Mirror 3D.
2 In the Settings window for Mirror 3D, locate the Data section.
3 From the Data set list, choose With Bars and Ottosen/Solution 3 (sol3).
4 Locate the Plane Data section. From the Plane list, choose ZX-planes.
5 In the Y-coordinate text field, type -1e-10.

22 | CONCRETE BEAM WITH REINFORCEMENT BARS

Duplicate the first von Mises stress plot group to compare results with or without the
failure behavior.

Stress With Bars 1

1 In the Model Builder window, under Results right-click Stress With Bars and choose
Duplicate.
2 In the Settings window for 3D Plot Group, type Stress With Bars and Ottosen in
the Label text field.
3 Locate the Data section. From the Data set list, choose Mirror 3D 3.
4 On the Stress With Bars and Ottosen toolbar, click Plot.
5 Click the Zoom Extents button on the Graphics toolbar.

Do the same with force and stress in reinforcement bars.

Force With Bars (truss) 1

1 In the Model Builder window, under Results right-click Force With Bars (truss) and
choose Duplicate.
2 In the Settings window for 3D Plot Group, locate the Data section.
3 From the Data set list, choose Mirror 3D 3.
4 In the Label text field, type Force With Bars and Ottosen (truss).
5 On the Force With Bars and Ottosen (truss) toolbar, click Plot.

Stress With Bars (truss) 1

1 In the Model Builder window, under Results right-click Stress With Bars (truss) and
choose Duplicate.
2 In the Settings window for 3D Plot Group, type Stress With Bars and Ottosen
(truss) in the Label text field.

3 Locate the Data section. From the Data set list, choose Mirror 3D 3.
4 On the Stress With Bars and Ottosen (truss) toolbar, click Plot.

Add a plot group to visualize the plastic zone.

3D Plot Group 8
1 On the Home toolbar, click Add Plot Group and choose 3D Plot Group.
2 In the Settings window for 3D Plot Group, locate the Data section.
3 From the Data set list, choose Mirror 3D 3.
4 In the Label text field, type Plastic Region.

23 | CONCRETE BEAM WITH REINFORCEMENT BARS

Surface 1
1 Right-click Plastic Region and choose Surface.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Model>Component 1>Solid Mechanics>
Strain>[Link] - Effective plastic strain.
3 Locate the Expression section. In the Expression text field, type [Link]>0.
4 Locate the Coloring and Style section. Clear the Color legend check box.
5 Click to expand the Quality section. From the Resolution list, choose Finer.
6 From the Smoothing list, choose None.

Plastic Region
1 Right-click Results>Plastic Region>Surface 1 and choose Deformation.
2 In the Model Builder window, under Results click Plastic Region.
3 In the Settings window for 3D Plot Group, click to expand the Title section.
4 Locate the Plot Settings section. Clear the Plot data set edges check box.
5 On the Plastic Region toolbar, click Plot.
6 Click Go to Default View.

To compare the deflection of the beam for the three models, proceed as follows.

1D Plot Group 9
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Deflection in the Label text field.
3 Locate the Data section. From the Parameter selection (para) list, choose Last.
4 Click to expand the Title section. From the Title type list, choose Manual.
5 In the Title text area, type Deflection of the beam.
6 Locate the Plot Settings section. Select the x-axis label check box.
7 In the associated text field, type Position on X axis.
8 Select the y-axis label check box.
9 In the associated text field, type Deflection (mm).

Line Graph 1
1 Right-click Deflection and choose Line Graph.
2 In the Settings window for Line Graph, locate the Data section.
3 From the Data set list, choose Without Bars/Solution 1 (sol1).

24 | CONCRETE BEAM WITH REINFORCEMENT BARS

4 From the Parameter selection (para) list, choose Last.
5 Select Edge 5 only.
6 Locate the y-Axis Data section. From the Unit list, choose mm.
7 Click Replace Expression in the upper-right corner of the y-axis data section. From the
menu, choose Model>Component 1>Solid Mechanics>Displacement>
Displacement field (material and geometry frames)>w - Displacement field, Z component.
8 Locate the x-Axis Data section. From the Parameter list, choose Expression.
9 In the Expression text field, type X.
10 Click to expand the Coloring and style section. Locate the Coloring and Style section. In
the Width text field, type 2.
11 Click to expand the Legends section. Select the Show legends check box.
12 From the Legends list, choose Manual.
13 In the table, enter the following settings:

Legends
Linear elastic model

Line Graph 2
1 Right-click Results>Deflection>Line Graph 1 and choose Duplicate.
2 In the Settings window for Line Graph, locate the Data section.
3 From the Data set list, choose With Bars/Solution 2 (sol2).
4 Locate the Legends section. In the table, enter the following settings:

Legends
Linear elastic model with bars

Line Graph 3
1 Right-click Line Graph 1 and choose Duplicate.
2 In the Settings window for Line Graph, locate the Data section.
3 From the Data set list, choose With Bars and Ottosen/Solution 3 (sol3).
4 Locate the Legends section. In the table, enter the following settings:

Legends
Ottosen model with bars

5 On the Deflection toolbar, click Plot.

25 | CONCRETE BEAM WITH REINFORCEMENT BARS

6 Click the Zoom Extents button on the Graphics toolbar.

26 | CONCRETE BEAM WITH REINFORCEMENT BARS

 Created in COMSOL Multiphysics 5.3

D eep E xcav at i on

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Introduction
There are two ways to model an excavation in COMSOL Multiphysics, both of which
include a parametric sweep. One option involves removing the soil one step at a time by
means of a a sweep of the geometry. As the soil is removed, the support it supplies is
removed as well, subjecting the retaining wall to soil stresses from the non excavated side.
The other option is to start with the already excavated geometry, and simulate the soil
removal by modifying a boundary load. The boundary load applies a force on the
excavation side of the retaining wall, equal to (and therefore balancing) the in-situ stresses
on the non excavated side, for the part of the wall that is below the virtual excavation
depth.

This deep excavation example is inspired by a benchmark exercise specified by a working

group of the German Society for Geotechnics (Ref. 1 and Ref. 2). In the example, a
26 meter excavation is modeled by means of a parametric sweep. As the excavation
deepens, three struts are activated using a ramp function and boolean expressions. As the
excavation reaches their depths, the struts are activated as long as the horizontal wall
deflection is greater than what is allowed to [Link] interaction between the soil and the
retaining wall is modeled with identity pairs.
Retaining wall
Struts

20 m
Excavation Soil, upper layer
Symmetry

Soil, lower layer

30 m

90 m

Figure 1: Dimensions and boundary conditions for the deep excavation example.

2 | DEEP EXCAVATION
Model Definition
In this example, a Drucker-Prager criterion is used for studying the soil plasticity, and the
retaining wall is made of a linear elastic material. The following material parameters are
used:

SOIL, UPPER LAYER

• Young’s modulus E = 20 MPa, Poisson’s ratio ν = 0.3, and density ρ = 1900 kg/m3.
• Cohesion c = 0 Pa, and angle of internal friction φ = 35°.

SOIL, LOWER LAYER

• Young’s modulus E = 60 MPa, Poisson’s ratio ν = 0.3, and density ρ = 1900 kg/m3.
• Cohesion c = 0 Pa, and angle of internal friction φ = 35°.

STRUTS
• Young’s modulus E = 200 GPa, length l = 30 m and cross sectional area A = 15 cm2.

RETA INING WA LL
• Young’s modulus E = 30 GPa, Poisson’s ratio ν = 0.15, and density ρ = 2400 kg/m3.

CONSTRAINTS AND LOADS

• The bottom soil layer is supported by a rigid and perfectly rough base. Apply therefore
fixed constraint on the lower horizontal boundary.
• Due to symmetry, model only the right half of the domain. Use the symmetry boundary
condition at the left vertical boundary.
• Add in-situ stresses with the External Stress feature. Note that the stress caused by
gravity is compressive, which in the convention used in the Structural Mechanics
Module means negative sign.
• The in-situ stresses account for the local stress prior excavating. You can define a
different value in the vertical and horizontal directions. The horizontal in-situ stress,
X_stress, is also used on the boundary load applied to the retaining wall.
• The boundary load on the retaining wall is gradually decreased in order to simulate the
excavation steps. At the bottom of the excavation it ensures zero initial displacement.
This strategy avoids re-meshing the excavated volume and thus saves memory.
• Add struts that can be activated and deactivated according to the excavation depth by
means of boolean expressions.

3 | DEEP EXCAVATION
• The struts are active after a maximum horizontal deflection is reached. Use a ramp
function to restrict the allowed wall deflection U_max to 25 mm, and to supply the axial
force that the struts apply to the wall.
• The boundaries between the retaining wall and the soil are modeled with extrusion
operators. These operators constrain the normal displacement between the retaining
wall and the soil to stay in contact, while the tangential displacement is unconstrained.
• The axial stiffness of the struts is estimated from their cross sectional area A, the length
l, and the material Young’s modulus E as

A
S = E ----
l

Results and Discussion

In order to improve convergence on the soil-wall boundary, apply a mapped mesh on the
retaining wall domain, as shown in Figure 2.

Figure 2: A mapped mesh on the retaining wall domain gives a better convergence.

4 | DEEP EXCAVATION
Figure 3 shows the soil deformation and the wall deflection after excavating 26 meters.
Figure 4 shows that the different properties of the two soil layers have an impact on the
plastic deformation.

5 | DEEP EXCAVATION
Figure 3: Deformation of the soil and the retaining wall after excavating 26 meters of soil.

Figure 4: Plastic deformation after excavating 26 meters of soil.

6 | DEEP EXCAVATION
Figure 5: Wall deflection as a function of depth for different excavation steps (color lines).

7 | DEEP EXCAVATION
The struts—placed 4.8 m, 9.3 m, and 14.35 m below the initial surface (Ref. 1)—help to
increase the overall stiffness of the retaining wall. In Figure 5, a maximum allowed
displacement of 2.5 cm constrains the horizontal displacement of the retaining wall.

Figure 6: Surface settlement as a function of the distance from the wall. The different color
lines show the settlement for different excavation stages (0 to 26 m).

As the excavation progresses, it is possible to observe the subsidence on the unexcavated

region. As expected, the subsidence increases as the excavation progresses, and decreases
with the distance to the retaining wall, see Figure 6.

References
1. H.F. Schweiger, Benchmarking in Geotechnics 1. Technical Report CGG IR006 2002,
Institute for Soil Mechanics and Foundation Engineering, Graz University of Technology,
Austria.

2. D. Potts and L. Zdravkovic, Finite Element Analysis in Geotechnical Engineering,

Thomas Telford Publishing, 2001.

8 | DEEP EXCAVATION
Application Library path: Geomechanics_Module/Soil/deep_excavation

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Stationary.
6 Click Done.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.

9 | DEEP EXCAVATION
3 In the table, enter the following settings:

Name Expression Value Description

X_stress -24e3[Pa] -2.4E4 Pa in-situ stress, xx
component
Y_stress -35e3[Pa] -3.5E4 Pa in-situ stress, yy
component
Z_stress -24e3[Pa] -2.4E4 Pa in-situ stress, zz
component
E_struts 2e5[MPa] 2E11 Pa Young’s modulus of
struts
A_struts 15[cm^2] 0.0015 m² Cross section area of
struts
l_struts 30[m] 30 m Length of struts
S_struts E_struts*A_struts/ 1E7 N/m Stiffness of struts
l_struts
U_max -25[mm] -0.025 m Allowed wall
deflection
Stage_1 -4.8[m] -4.8 m First excavation
step, first strut
Stage_2 -9.3[m] -9.3 m Second excavation
step, second strut
Stage_3 -14.35[m] -14.35 m Third excavation
step, third strut
Depth 0 0 Excavation depth
(parameter)

In-situ stresses are set with negative sign to fit the structural mechanics convention.
That convention assumes negative stresses in compression and positive in tension.

Ramp 1 (rm1)
1 On the Home toolbar, click Functions and choose Global>Ramp.
2 In the Settings window for Ramp, locate the Parameters section.
3 In the Location text field, type -U_max.
4 In the Slope text field, type S_struts/U_max.
5 Click to expand the Smoothing section. In the Size of transition zone text field, type
0.0001.

6 Select the Smooth at start check box.

10 | DEEP EXCAVATION
GEOMETRY 1
Create the geometry.

To simplify this step, you can insert a prepared geometry sequence. On the Geometry
toolbar, click Insert Sequence. Browse to the application’s Application Library folder and
double-click the file deep_excavation.mph. Click Build All on the Geometry toolbar.
Then, continue with the instruction after the geometry plot below.

Otherwise, proceed with the following instructions to create the geometry from scratch:

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 90.
4 In the Height text field, type 60.
5 Locate the Position section. In the y text field, type -60.
6 Right-click Rectangle 1 (r1) and choose Build Selected.

Rectangle 2 (r2)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 90.
4 In the Height text field, type 20.
5 Locate the Position section. In the y text field, type -20.
6 Right-click Rectangle 2 (r2) and choose Build Selected.

Union 1 (uni1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Union.
2 Click in the Graphics window and then press Ctrl+A to select both objects.
3 Right-click Union 1 (uni1) and choose Build Selected.

Rectangle 3 (r3)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 30.
4 In the Height text field, type 30.
5 Locate the Position section. In the y text field, type -30.

11 | DEEP EXCAVATION
6 Right-click Rectangle 3 (r3) and choose Build Selected.

Rectangle 4 (r4)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 0.8.
4 In the Height text field, type 30.
5 Locate the Position section. In the x text field, type 30.
6 In the y text field, type -30.
7 Right-click Rectangle 4 (r4) and choose Build Selected.

Difference 1 (dif1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Difference.
2 Select the object uni1 only.
3 In the Settings window for Difference, locate the Difference section.
4 Find the Objects to subtract subsection. Select the Active toggle button.
5 Select the objects r3 and r4 only.
6 Select the Keep input objects check box.
7 Right-click Difference 1 (dif1) and choose Build Selected.

Delete Entities 1 (del1)

1 Right-click Geometry 1 and choose Delete Entities.
2 In the Settings window for Delete Entities, locate the Entities or Objects to Delete section.
3 From the Geometric entity level list, choose Object.
4 Under Selection, click Clear Selection.
5 Click the Select Box button on the Graphics toolbar.
6 Select the objects r3 and uni1 only.
7 Right-click Component 1 (comp1)>Geometry 1>Delete Entities 1 (del1) and choose
Build Selected.

Bézier Polygon 1 (b1)

1 On the Geometry toolbar, click Primitives and choose Bézier Polygon.
2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 1, set x to 30 and y to Stage_1+1.

12 | DEEP EXCAVATION
5 In row 2, set x to 30 and y to Stage_1.
6 Right-click Bézier Polygon 1 (b1) and choose Build Selected.

Bézier Polygon 2 (b2)

1 On the Geometry toolbar, click Primitives and choose Bézier Polygon.
2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 1, set x to 30 and y to Stage_2+1.
5 In row 2, set x to 30 and y to Stage_2.
6 Right-click Bézier Polygon 2 (b2) and choose Build Selected.

Bézier Polygon 3 (b3)

1 On the Geometry toolbar, click Primitives and choose Bézier Polygon.
2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 1, set x to 30 and y to Stage_3+1.
5 In row 2, set x to 30 and y to Stage_3.
6 Right-click Bézier Polygon 3 (b3) and choose Build Selected.

Union 2 (uni2)
1 On the Geometry toolbar, click Booleans and Partitions and choose Union.
2 Select the objects b1, b2, b3, and r4 only.
3 Right-click Union 2 (uni2) and choose Build Selected.

Form Union (fin)

1 In the Model Builder window, under Component 1 (comp1)>Geometry 1 click
Form Union (fin).
2 In the Settings window for Form Union/Assembly, locate the Form Union/Assembly section.
3 From the Action list, choose Form an assembly.

13 | DEEP EXCAVATION
4 Right-click Component 1 (comp1)>Geometry 1>Form Union (fin) and choose
Build Selected.

DEFINITIONS

Explicit 1
1 On the Definitions toolbar, click Explicit.
2 In the Settings window for Explicit, type Wall diaphragm in the Label text field.
3 Locate the Input Entities section. From the Geometric entity level list, choose Boundary.
4 Select Boundary 20 only.

Explicit 2
1 On the Definitions toolbar, click Explicit.
2 In the Settings window for Explicit, type Wall soil in the Label text field.
3 Locate the Input Entities section. From the Geometric entity level list, choose Boundary.
4 Select Boundaries 5 and 6 only.

General Extrusion 1 (genext1)

1 On the Definitions toolbar, click Component Couplings and choose General Extrusion.
2 In the Settings window for General Extrusion, locate the Source Selection section.

14 | DEEP EXCAVATION
3 From the Geometric entity level list, choose Boundary.
4 Click Paste Selection.
5 In the Paste Selection dialog box, type 12 in the Selection text field.
6 Click OK.
7 In the Settings window for General Extrusion, locate the Source section.
8 From the Source frame list, choose Material (X, Y, Z).
9 Locate the Destination Map section. In the X-expression text field, type X.
10 In the Y-expression text field, type Y.

General Extrusion 2 (genext2)

1 On the Definitions toolbar, click Component Couplings and choose General Extrusion.
2 In the Settings window for General Extrusion, locate the Source Selection section.
3 From the Geometric entity level list, choose Boundary.
4 From the Selection list, choose Wall diaphragm.
5 Locate the Source section. From the Source frame list, choose Material (X, Y, Z).
6 Locate the Destination Map section. In the X-expression text field, type X.
7 In the Y-expression text field, type Y.

SOLID MECHANICS (SOLID)

Linear Elastic Material 1

In this model, expected displacements are small compared to the geometry. Click Force
linear strains to use the small strain formulation.

1 In the Model Builder window, expand the Solid Mechanics (solid) node, then click
Linear Elastic Material 1.
2 In the Settings window for Linear Elastic Material, locate the Geometric Nonlinearity
section.
3 Select the Force linear strains check box.

Soil Plasticity 1
1 Right-click Component 1 (comp1)>Solid Mechanics (solid)>Linear Elastic Material 1 and
choose Soil Plasticity.
2 In the Settings window for Soil Plasticity, locate the Soil Plasticity section.
3 Select the Match to Mohr-Coulomb criterion check box.
4 Select Domains 1 and 2 only.

15 | DEEP EXCAVATION
External Stress 1
1 Right-click Linear Elastic Material 1 and choose External Stress.
2 In the Settings window for External Stress, locate the External Stress section.
3 From the list, choose Diagonal.
4 In the Sext table, enter the following settings:

X_stress 0 0
0 Y_stress 0
0 0 Z_stress

Symmetry 1
1 In the Model Builder window, right-click Solid Mechanics (solid) and choose
More Constraints>Symmetry.
2 Select Boundary 1 only.

Fixed Constraint 1
1 Right-click Solid Mechanics (solid) and choose Fixed Constraint.
2 Select Boundary 2 only.

Roller 1
1 Right-click Solid Mechanics (solid) and choose Roller.
2 Select Boundaries 9 and 10 only.

Prescribed Displacement 1
1 Right-click Solid Mechanics (solid) and choose Prescribed Displacement.
2 Select Boundary 4 only.
3 In the Settings window for Prescribed Displacement, locate the Prescribed Displacement
section.
4 Select the Prescribed in y direction check box.
5 In the u 0y text field, type genext1(v).

Prescribed Displacement 2
1 Right-click Solid Mechanics (solid) and choose Prescribed Displacement.
2 In the Settings window for Prescribed Displacement, locate the Boundary Selection
section.
3 From the Selection list, choose Wall soil.

16 | DEEP EXCAVATION
4 Locate the Prescribed Displacement section. Select the Prescribed in x direction check
box.
5 In the u 0x text field, type genext2(u).

Boundary Load 1
1 Right-click Solid Mechanics (solid) and choose Boundary Load.
2 Select Boundaries 3, 8, and 19 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

0 x
Y_stress y

Boundary Load 2
1 Right-click Solid Mechanics (solid) and choose Boundary Load.
2 Click the Select Box button on the Graphics toolbar.
3 Select Boundaries 11 and 13–18 only.
4 In the Settings window for Boundary Load, locate the Force section.
5 Specify the FA vector as

-X_stress*(y<Depth) x
0 y

Boundary Load 3
1 Right-click Solid Mechanics (solid) and choose Boundary Load.
2 In the Settings window for Boundary Load, type Strut_1 in the Label text field.
3 Select Boundary 17 only.
4 Locate the Force section. From the Load type list, choose Total force.
5 Specify the Ftot vector as

-rm1(-u[1/m])*(Depth<Stage_1) x
0 y

Boundary Load 4
1 Right-click Solid Mechanics (solid) and choose Boundary Load.
2 In the Settings window for Boundary Load, type Strut_2 in the Label text field.
3 Select Boundary 15 only.

17 | DEEP EXCAVATION
4 Locate the Force section. From the Load type list, choose Total force.
5 Specify the Ftot vector as

-rm1(-u[1/m])*(Depth<Stage_2) x
0 y

Boundary Load 5
1 Right-click Solid Mechanics (solid) and choose Boundary Load.
2 In the Settings window for Boundary Load, type Strut_3 in the Label text field.
3 Select Boundary 13 only.
4 Locate the Force section. From the Load type list, choose Total force.
5 Specify the Ftot vector as

-rm1(-u[1/m])*(Depth<Stage_3) x
0 y

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, type Soil, upper layer in the Label text field.
3 Locate the Geometric Entity Selection section. From the Selection list, choose Manual.
4 Click Clear Selection.
5 Select Domain 2 only.
6 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 20e6 Pa Basic
Poisson’s ratio nu 0.3 1 Basic
Density rho 1900 kg/m³ Basic
Cohesion cohesion 0 Pa Mohr-Coulomb
Angle of internalphi 35[deg] rad Mohr-Coulomb
internal friction

Material 2 (mat2)
1 Right-click Materials and choose Blank Material.

18 | DEEP EXCAVATION
2 In the Settings window for Material, type Soil, lower layer in the Label text field.
3 Select Domain 1 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 60e6 Pa Basic
Poisson’s ratio nu 0.3 1 Basic
Density rho 1900 kg/m³ Basic
Cohesion cohesion 0 Pa Mohr-Coulomb
Angle of internalphi 35[deg] rad Mohr-Coulomb
internal friction

Material 3 (mat3)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Retaining wall in the Label text field.
3 Select Domain 3 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 30e9 Pa Basic
Poisson’s ratio nu 0.15 1 Basic
Density rho 2400 kg/m³ Basic

MESH 1

Mapped 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Mapped.
Use a mapped mesh inside the wall.
2 In the Settings window for Mapped, locate the Domain Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domain 3 only.

Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Mapped 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Boundary Selection section.

19 | DEEP EXCAVATION
3 From the Selection list, choose Wall diaphragm.
4 Locate the Distribution section. In the Number of elements text field, type 60.

Distribution 2
1 Right-click Mapped 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Boundary Selection section.
3 Click Paste Selection.
4 In the Paste Selection dialog box, type 12 in the Selection text field.
5 Click OK.
6 In the Settings window for Distribution, locate the Distribution section.
7 In the Number of elements text field, type 2.

Size 1
1 In the Model Builder window, right-click Mesh 1 and choose Free Triangular.
2 Right-click Free Triangular 1 and choose Size.
3 In the Settings window for Size, locate the Element Size section.
4 From the Predefined list, choose Finer.

Distribution 1
1 Right-click Free Triangular 1 and choose Distribution.
2 Select Boundary 6 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 40.

Distribution 2
1 Right-click Free Triangular 1 and choose Distribution.
2 Select Boundary 5 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 20.

Distribution 3
1 Right-click Free Triangular 1 and choose Distribution.
2 Select Boundary 4 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 3.

20 | DEEP EXCAVATION
5 Click Build All.
Click the Zoom Box button on the Graphics toolbar and then use the mouse to zoom in
on the contact zone where the mesh is the densest.
The mesh should look like the one in Figure 2.

STUDY 1

Step 1: Stationary
Set up an auxiliary continuation sweep for the ’Depth’ parameter.

1 In the Model Builder window, under Study 1 click Step 1: Stationary.

2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

Depth range(0,-2,-26)

Solution 1 (sol1)
1 On the Study toolbar, click Show Default Solver.
2 In the Model Builder window, expand the Solution 1 (sol1) node.
3 In the Model Builder window, expand the Study 1>Solver Configurations>
Solution 1 (sol1)>Stationary Solver 1 node, then click Parametric 1.
4 In the Settings window for Parametric, click to expand the Continuation section.
5 From the Predictor list, choose Constant.
Restrict the step size to 0.5 to better capture the development of plasticity.
6 Select the Tuning of step size check box.
7 In the Initial step size text field, type 0.5.
8 In the Maximum step size text field, type 0.5.
9 On the Study toolbar, click Compute.

RESULTS

Stress (solid)
The default plot shows the von Mises stress for the last value of the parameter. Modify and
rename this plot group to display the displacement.

21 | DEEP EXCAVATION
1 In the Model Builder window, under Results click Stress (solid).
2 In the Settings window for 2D Plot Group, type Displacement in the Label text field.

Surface 1
1 In the Model Builder window, expand the Results>Displacement node, then click
Surface 1.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Component 1>Solid Mechanics>
Displacement>[Link] - Total displacement.
3 On the Displacement toolbar, click Plot.
4 Click the Zoom Extents button on the Graphics toolbar.

2D Plot Group 2
1 On the Home toolbar, click Add Plot Group and choose 2D Plot Group.
2 In the Settings window for 2D Plot Group, type Plastic Region in the Label text field.

Surface 1
1 Right-click Plastic Region and choose Surface.
2 In the Settings window for Surface, locate the Expression section.
3 In the Expression text field, type [Link]>0.
4 On the Plastic Region toolbar, click Plot.

1D Plot Group 3
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Wall Deflection in the Label text field.
3 Click to expand the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Wall deflection (mm).
5 Locate the Plot Settings section. Select the x-axis label check box.
6 In the associated text field, type Horizontal displacement (mm).
7 Select the y-axis label check box.
8 In the associated text field, type Depth(m).
9 Click to expand the Legend section. From the Position list, choose Upper left.

Line Graph 1
1 Right-click Wall Deflection and choose Line Graph.
2 In the Settings window for Line Graph, locate the Selection section.

22 | DEEP EXCAVATION
3 From the Selection list, choose Wall diaphragm.
4 Locate the y-Axis Data section. In the Expression text field, type y.
5 Locate the x-Axis Data section. From the Parameter list, choose Expression.
6 In the Expression text field, type u.
7 From the Unit list, choose mm.
8 Click to expand the Legends section. Select the Show legends check box.
9 On the Wall Deflection toolbar, click Plot.

1D Plot Group 4
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Surface Settlement in the Label text
field.
3 Locate the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Surface settlement.
5 Locate the Plot Settings section. Select the x-axis label check box.
6 In the associated text field, type Distance from the diaphragm (m).
7 Select the y-axis label check box.
8 In the associated text field, type Vertical displacement (mm).
9 Locate the Legend section. From the Position list, choose Lower right.

Line Graph 1
1 Right-click Surface Settlement and choose Line Graph.
2 Select Boundary 8 only.
3 In the Settings window for Line Graph, locate the y-Axis Data section.
4 In the Expression text field, type v.
5 From the Unit list, choose mm.
6 Locate the Legends section. Select the Show legends check box.
7 On the Surface Settlement toolbar, click Plot.

23 | DEEP EXCAVATION
24 | DEEP EXCAVATION
 Created in COMSOL Multiphysics 5.3

Flexible and Smooth Strip Footing on Stratum of

Clay

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Model Definition
A typical verification example for geotechnical problems is a shallow stratum layer of clay,
see Figure 1. In the example, a vertical load is applied on the clay stratum, and the static
response as well as the collapse load are of interest.
Boundary load
Infinite element domain Symmetry axis

1.57 m
3.66 m

7.32 m

Figure 1: Dimensions, boundary conditions, and pressure load for the stratum of clay.

A N A L Y S I S TY P E

Yield Surface
Assume plane-strain conditions, and model the clay with soil plasticity and Drucker-Prager
criterion.

The yield surface, F, for the Drucker-Prager criterion is given by

F = J2 + α I1 – k = 0

where I1 is the first stress invariant and J2 is the second deviatoric stress invariant.

The first stress invariant is defined using the trace of Cauchy stress tensor:

I 1 = trace ( σ )

The second stress invariant is defined by

1 2 2
I 2 = --- ( I 1 – trace ( σ ) )
2

2 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

The second deviatoric stress invariant can be expressed using the first and the second stress
invariants:

1 2
J 2 = --- I 1 – I 2
3

If two-dimensional plane-strain conditions prevail, the Drucker-Prager criterion matches

the Mohr-Coulomb criterion. For this case the material parameters α and k are given by
the Cohesion c and the Angle of internal friction φ (Ref. 1)

tan φ
α = -----------------------------------------
2
( 9 + 12 tan φ )

3c
k = -----------------------------------------
2
( 9 + 12 tan φ )

Drucker-Prager criterion is the default choice for the Soil Plasticity feature, and the check
box Match to Mohr-Coulomb criterion applies the aforementioned matching of the material
parameters.

Under Soil Plasticity it is also possible to use Mohr-Coulomb criterion

1 1
F = --- ( σ max – σ min ) + --- ( σ max + σ min ) sin φ – c cos φ = 0 ,
2 2

where σmax and σmin are the biggest and smallest principal stresses. The Mohr-Coulomb
criterion defines an irregular hexagon pyramid in the principal stress space. Since this yield
function gives rise to singularities in the derivatives of the yield function, the use of a non-
associated flow rule with a Drucker-Prager plastic potential is chosen. This is done in the
Plastic potential list, with the option Drucker-Prager matched at compressive meridian.

Flow Rule
The flow rule defines the relation between the plastic strain increment in a given direction
and the current level of stress in the same direction. The relation reads

∂Q
dε ij = dλ ----------
∂σ ij

where dλ is the plastic multiplier and Q is the plastic potential. If the yield surface, F, and
the plastic potential, Q, are identical, that is, if F = Q, then it is called an associated flow
rule, otherwise it is called non-associated flow rule.

3 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

SOIL PROPERTIES
• Young’s modulus, E = 207 MPa, and Poisson’s ratio ν = 0.3.
• Cohesion c = 69 MPa, and angle of internal friction φ = 20 degrees.

CONSTRAINTS AND LOADS

• The clay layer is supported by a rigid and perfectly rough base. Apply therefore a fixed
constraint on the lower horizontal boundary.
• Model only the left half of the domain due to symmetry reasons. Use symmetry
boundary condition at the right vertical boundary.
• The stratum is subjected to a footing that you can consider to be flexible and smooth.
The width of the strip footing is 3.14 m. Gradually increase the footing pressure until
the clay layer reaches the collapse load.

INFINITE ELEMENT DOMAIN

• In order to mimic an infinite layer of soil, add an Infinite Element Domain. The scaling
1e3*[Link] means that the spatial variables in this domain are scaled
thousand times the typical geometry length.
• The left vertical boundary is perfectly smooth and a can be assumed to be of the roller
type.
This example is adapted from Ref. 2.

Results and Discussion

You can study in Figure 2 the load-displacement curves for both the Mohr-Coulomb and
the Drucker-Prager criteria. The figure shows the applied footing pressure versus the
centerline displacement (directly beneath the footing’s center) in the y direction. The lines
show the load-displacement curves for Mohr-Coulomb and Drucker-Prager criteria.

4 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

The example uses the SI unit system. The curves are identical up to 300 kPa because the
whole domain is still within the elastic region. From that point when the pressure increases
the behavior diverges. Both curves reach the collapse load at approximately 1.1 MPa.

Figure 2: Footing pressure versus vertical displacement for Mohr-Coulomb and Drucker-
Prager material models.

The development of the plasticity in the soil is shown in Figure 3.

Notes About the COMSOL Implementation

Both Mohr-Coulomb or Drucker-Prager criterion are predefined in the Soil Plasticity
subfeature (the default Mohr-Coulomb criterion uses a nonassociated flow rule, with a
plastic potential implemented as Drucker-Prager matched at compressive meridian). The
yield function contains the first stress invariant, I1, and the second deviatoric stress
invariant, J2, as well as the material property constants.

A suitable modeling technique in a case where the relation between the applied load and
the displacement is highly nonlinear, is to use an algebraic equation that controls the
applied pressure so that the model reaches the desired displacement increments. This is

5 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

implemented using a Global Equation, and the parametric solver steps up the desired
vertical displacement.

Figure 3: Evolution of the effective plastic strain on the clay layer during the parametric
loading.

References
1. W.F. Chen and E. Mizuno, Nonlinear Analysis in Soil Mechanics, Elsevier, 1990.

2. A. Mar. How To Undertake Finite Element Based Geotechnical Analysis, NAFEMS,

2002.

6 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

Application Library path: Geomechanics_Module/Soil/flexible_footing

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 7.32.
4 In the Height text field, type 3.66.
5 Right-click Rectangle 1 (r1) and choose Build Selected.
Add a rectangle to model the infinite element domain.

Rectangle 2 (r2)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 7.32*0.1.
4 In the Height text field, type 3.66.
5 Locate the Position section. In the x text field, type -7.32*0.1.
6 Right-click Rectangle 2 (r2) and choose Build Selected.

7 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

Point 1 (pt1)
1 On the Geometry toolbar, click Primitives and choose Point.
2 In the Settings window for Point, locate the Point section.
3 In the x text field, type 7.32-1.57.
4 In the y text field, type 3.66.
5 Right-click Point 1 (pt1) and choose Build Selected.

Form Union (fin)

In the Model Builder window, under Component 1 (comp1)>Geometry 1 right-click
Form Union (fin) and choose Build Selected.

DEFINITIONS

Infinite Element Domain 1 (ie1)

1 On the Definitions toolbar, click Infinite Element Domain.
The infinite element is scaled by a factor of 1000.
2 Select Domain 1 only.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

para 0 0 Prescribed displacement

DEFINITIONS
Use an integration coupling operator to evaluate the displacement in the center of the
applied pressure (Point 7).

Integration 1 (intop1)
1 On the Definitions toolbar, click Component Couplings and choose Integration.
2 In the Settings window for Integration, locate the Source Selection section.
3 From the Geometric entity level list, choose Point.
4 Select Point 7 only.

8 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

Variables 1
1 On the Definitions toolbar, click Local Variables.
2 In the Settings window for Variables, locate the Variables section.
3 In the table, enter the following settings:

Name Expression Unit Description

footing_pressure pressure[Pa] Footing pressure
vc intop1(v) m Displacement

SOLID MECHANICS (SOLID)

Linear Elastic Material 1

In the Model Builder window, expand the Component 1 (comp1)>Solid Mechanics (solid)
node.

Soil Plasticity 1
1 Right-click Linear Elastic Material 1 and choose Soil Plasticity.
2 In the Settings window for Soil Plasticity, locate the Soil Plasticity section.
3 Select the Match to Mohr-Coulomb criterion check box.

Soil Plasticity 2
1 Right-click Linear Elastic Material 1 and choose Soil Plasticity.
2 In the Settings window for Soil Plasticity, locate the Soil Plasticity section.
3 From the Yield criterion list, choose Mohr-Coulomb.

Fixed Constraint 1
1 In the Model Builder window, right-click Solid Mechanics (solid) and choose
Fixed Constraint.
2 Select Boundaries 2 and 5 only.

Symmetry 1
1 Right-click Solid Mechanics (solid) and choose More Constraints>Symmetry.
2 Select Boundary 8 only.

Roller 1
1 Right-click Solid Mechanics (solid) and choose Roller.
2 Select Boundary 1 only.

9 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

Boundary Load 1
1 Right-click Solid Mechanics (solid) and choose Boundary Load.
2 Select Boundary 7 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

0 x
pressure y

5 In the Model Builder window’s toolbar, click the Show button and select
Advanced Physics Options in the menu.

Global Equations 1
1 Right-click Solid Mechanics (solid) and choose Global>Global Equations.
2 In the Settings window for Global Equations, locate the Global Equations section.
3 In the table, enter the following settings:

Name f(u,ut,utt,t) (1) Initial value Initial value Description

(u_0) (1) (u_t0) (1/s)
pressure vc-para 0 0

4 Locate the Units section. Find the Dependent variable quantity subsection. From the list,
choose Pressure (Pa).
5 Find the Source term quantity subsection. From the list, choose Displacement field (m).

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 207e6 Pa Basic
Poisson’s ratio nu 0.3 1 Basic

10 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

Property Name Value Unit Property group
Cohesion cohesion 69e3 Pa Mohr-Coulomb
Angle of internal internalphi 20[deg] rad Mohr-Coulomb
friction

MESH 1
1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2 In the Settings window for Mesh, locate the Mesh Settings section.
3 From the Sequence type list, choose User-controlled mesh.

Size
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Finer.

Free Triangular 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click
Free Triangular 1.
2 In the Settings window for Free Triangular, locate the Domain Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domain 2 only.

Mapped 1
1 In the Model Builder window, right-click Mesh 1 and choose Mapped.
In the Infinite Element Domain, use a mapped mesh to improve convergence.
2 In the Settings window for Mapped, locate the Domain Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domain 1 only.
5 Click Build All.

The first study is parametric and solves the model assuming a Drucker-Prager criterion.
The parameter represents the vertical displacement in the center of the applied pressure
(Point 7). It runs from 0 to -32 mm with a step size of 0.5 mm.

STUDY 1
1 In the Model Builder window, click Study 1.

11 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

2 In the Settings window for Study, type Drucker-Prager criterion in the Label text
field.

DRUCKER-PRAGER CRITERION

Step 1: Stationary
1 In the Model Builder window, under Drucker-Prager criterion click Step 1: Stationary.
2 In the Settings window for Stationary, locate the Physics and Variables Selection section.
3 Select the Modify physics tree and variables for study step check box.
4 In the Physics and variables selection tree, select Component 1 (comp1)>
Solid Mechanics (solid)>Linear Elastic Material 1>Soil Plasticity 2.
5 Click Disable.
Set up an auxiliary continuation sweep for the ’para’ parameter.
6 Click to expand the Study extensions section. Locate the Study Extensions section. Select
the Auxiliary sweep check box.
7 Click Add.
8 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

para range(0,-5e-4,-32e-3)

Solution 1 (sol1)
1 On the Study toolbar, click Show Default Solver.
2 In the Model Builder window, expand the Drucker-Prager criterion>Solver Configurations
node.
3 In the Model Builder window, expand the Solution 1 (sol1) node.
4 In the Model Builder window, expand the Drucker-Prager criterion>Solver Configurations>
Solution 1 (sol1)>Stationary Solver 1 node, then click Parametric 1.
5 In the Settings window for Parametric, click to expand the Continuation section.
For convergence reasons, set the predictor of the parametric solver as constant.
6 From the Predictor list, choose Constant.
7 On the Study toolbar, click Compute.

12 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

ROOT
The second study is also parametric and solves the model assuming a Mohr-Coulomb
criterion. Again, the parameter represents the vertical displacement in the center of the
applied pressure (Point 7), this time running from 0 to -24.5 mm with a step size of
0.5 mm.

ADD STUDY
1 On the Study toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select Study tree, select Preset Studies>Stationary.
4 Click Add Study in the window toolbar.
5 On the Study toolbar, click Add Study to close the Add Study window.

STUDY 2
1 In the Model Builder window, click Study 2.
2 In the Settings window for Study, type Mohr-Coulomb criterion in the Label text field.

MOHR-COULOMB CRITERION

Step 1: Stationary
1 In the Model Builder window, under Mohr-Coulomb criterion click Step 1: Stationary.
2 In the Settings window for Stationary, locate the Physics and Variables Selection section.
3 Select the Modify physics tree and variables for study step check box.
4 In the Physics and variables selection tree, select Component 1 (comp1)>
Solid Mechanics (solid)>Linear Elastic Material 1>Soil Plasticity 1.
5 Click Disable.
Set up an auxiliary continuation sweep for the ’para’ parameter.
6 Locate the Study Extensions section. Select the Auxiliary sweep check box.
7 Click Add.
8 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

para range(0,-5e-4,-24.5e-3)

13 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

Solution 2 (sol2)
1 On the Study toolbar, click Show Default Solver.
For convergence reasons, change the relative tolerance from 1e-3 to 1e-4 and set the
predictor of the continuation solver as constant.
2 In the Model Builder window, expand the Mohr-Coulomb criterion>Solver Configurations
node.
3 In the Model Builder window, expand the Solution 2 (sol2) node, then click
Stationary Solver 1.
4 In the Settings window for Stationary Solver, locate the General section.
5 In the Relative tolerance text field, type 1e-4.
6 In the Model Builder window, expand the Mohr-Coulomb criterion>Solver Configurations>
Solution 2 (sol2)>Stationary Solver 1 node, then click Parametric 1.
7 In the Settings window for Parametric, locate the Continuation section.
8 From the Predictor list, choose Constant.
9 On the Study toolbar, click Compute.

RESULTS
Remove the infinite domain from the data set for plotting.

Drucker-Prager criterion/Solution 1 (sol1)

In the Model Builder window, expand the Results>Data Sets node, then click Drucker-
Prager criterion/Solution 1 (sol1).

Selection
1 On the Results toolbar, click Selection.
2 In the Settings window for Selection, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domain 2 only.

Mohr-Coulomb criterion/Solution 2 (sol2)

In the Model Builder window, under Results>Data Sets click Mohr-Coulomb criterion/
Solution 2 (sol2).

Selection
1 On the Results toolbar, click Selection.
2 In the Settings window for Selection, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Domain.

14 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

4 Select Domain 2 only.

Stress (solid)
The default plots show the von Mises stress at the final step for each study.

1 In the Model Builder window, under Results click Stress (solid).

2 In the Settings window for 2D Plot Group, type Stress, Drucker-Prager criterion
in the Label text field.

Stress (solid) 1
1 In the Model Builder window, under Results click Stress (solid) 1.
2 In the Settings window for 2D Plot Group, type Stress, Mohr-Coulomb criterion in
the Label text field.

Create a plot that displays the part of the model that undergoes plasticity.

2D Plot Group 3
1 On the Results toolbar, click 2D Plot Group.
2 In the Settings window for 2D Plot Group, type Plastic Region in the Label text field.

Surface 1
1 Right-click Plastic Region and choose Surface.
2 In the Settings window for Surface, locate the Expression section.
3 In the Expression text field, type [Link]>0.
4 Locate the Coloring and Style section. From the Color table list, choose Traffic.

Deformation 1
1 Right-click Results>Plastic Region>Surface 1 and choose Deformation.
2 In the Settings window for Deformation, locate the Scale section.
3 Select the Scale factor check box.
4 In the associated text field, type 10.
5 On the Plastic Region toolbar, click Plot.
6 Click the Zoom Extents button on the Graphics toolbar.

Finally, follow the steps below to reproduce the graph in Figure 2.

1D Plot Group 4
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Footing Pressure vs Displacement
in the Label text field.

15 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

3 Click to expand the Legend section. From the Position list, choose Lower right.

Point Graph 1
1 Right-click Footing Pressure vs Displacement and choose Point Graph.
2 In the Settings window for Point Graph, locate the Data section.
3 From the Data set list, choose Drucker-Prager criterion/Solution 1 (sol1).
4 Select Point 7 only.
5 Click Replace Expression in the upper-right corner of the y-axis data section. From the
menu, choose Component 1>Definitions>Variables>footing_pressure - Footing pressure.
6 Locate the y-Axis Data section. Select the Description check box.
7 In the Expression text field, type abs(footing_pressure).
8 Click Replace Expression in the upper-right corner of the x-axis data section. From the
menu, choose Component 1>Definitions>Variables>vc - Displacement.
9 Locate the x-Axis Data section. Select the Description check box.
10 In the Expression text field, type abs(vc).
11 Click to expand the Legends section. Select the Show legends check box.
12 From the Legends list, choose Manual.
13 In the table, enter the following settings:

Legends
Drucker-Prager

Point Graph 2
1 In the Model Builder window, under Results right-click Footing Pressure vs Displacement
and choose Point Graph.
2 In the Settings window for Point Graph, locate the Data section.
3 From the Data set list, choose Mohr-Coulomb criterion/Solution 2 (sol2).
4 Select Point 7 only.
5 Click Replace Expression in the upper-right corner of the y-axis data section. From the
menu, choose footing_pressure - Footing pressure.
6 Locate the y-Axis Data section. Select the Description check box.
7 In the Expression text field, type abs(footing_pressure).
8 Click Replace Expression in the upper-right corner of the x-axis data section. From the
menu, choose vc - Displacement.
9 Locate the x-Axis Data section. Select the Description check box.

16 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

10 In the Expression text field, type abs(vc).
11 Locate the Legends section. Select the Show legends check box.
12 From the Legends list, choose Manual.
13 In the table, enter the following settings:

Legends
Mohr-Coulomb

Footing Pressure vs Displacement

1 In the Model Builder window, under Results click Footing Pressure vs Displacement.
2 On the Footing Pressure vs Displacement toolbar, click Plot.
This graph shows the footing pressure as a function of the vertical displacement. As you
can see, both curves reach the same failure level.

17 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY

18 | FLEXIBLE AND SMOOTH STRIP FOOTING ON STRATUM OF CLAY
 Created in COMSOL Multiphysics 5.3

Isotropic Compression with Cam-Clay Material

Model

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Introduction
Isotropic compression is a common exercise in soil testing. In this example a modified
Cam-Clay model is examined and in more particular the relation between the void ratio
and the logarithm of the pressure.

Model Definition
In this example, a soil sample is placed inside a cylinder with 10 cm diameter and 10 cm
height, see Figure 1. Due to the symmetry, the model is solved in 2D axial symmetry.

A boundary load produces isotropic compression conditions.

Boundary load
Axial symmetry

10 cm Boundary load

5 cm

Figure 1: Dimensions, boundary conditions, and boundary load for the isotropic compression
test.

CAM-CLAY MATERIAL PROPERTIES

• Density ρ = 2400 kg/m3, Poisson’s ratio ν = 0.2, angle of internal friction φ = 30°,
swelling index κ = 0.013, compression index λ = 0.032, void ratio N = 0.7 at a reference
pressure prefN = 100 kPa, initial void ratio e0 = 0.6646, and initial consolidation
pressure pC0 = 400 kPa.

CONSTRAINTS AND LOADS

• The left boundary is the axis of symmetry, a roller condition is applied at the lower
boundary, and a boundary load is applied on the right and upper boundaries.

2 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

• An hydrostatic initial pressure equal to p0 = 200 kPa is applied through the Initial Stress
and Strain feature.
• The boundary load is applied in three steps. First the pressure increase from p0 to 2.6p0,
then the sample is unloaded until 1.8p0, and finally the pressure increases again until
3.4p0.

In order to reproduce the analytical results of Ref. 1, the load is controlled in a parametric
analysis.

Results and Discussion

The produced void ratio versus pressure is shown in Figure 2. Note that the log operator
is implemented in base “e” and not in base “10”.

Figure 2: Void ratio as a function of the logarithm of the pressure in an isotropic compression
test.

For the parametric sweep between 0.05 and 0.25 (the boundary load ranges from 240 kPa
to 400 kPa) the curve follows the slope defined by the swelling index κ.

Once the user-defined consolidation pressure is reached (pC0 = 400 kPa), the sample of
soil behaves plastically, and the curve follows the slope defined by the compression index λ.

3 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

During the unloading and reloading of the soil (between the parameters 0.4 and 0.8), the
curve in Figure 2 follows the elastic slope defined by the swelling index κ.

Finally, the soil is compressed between the parameters 0.8 and 1, and it undergoes plastic
deformation until reaching its final stage at a void ratio e0 = 0.6387 and pressure
p = 680 kPa.

As expected, once the Cam-Clay ellipse is reached (p = pC0, parameter 0.25), the soil
sample starts deforming plastically. Isotropic hardening expands the major semi-axis of the
ellipse, with the expansion given by the increase in consolidation pressure, see Figure 3.
During unloading-reloading steps (between the parameter values 0.4 and 0.8), the
consolidation pressure is kept constant.

Figure 3: Increase in consolidation pressure due to isotropic hardening.

Reference
1. W.F. Chen and E. Mizuno, Nonlinear Analysis in Soil Mechanics, Elsevier, 1990.

Application Library path: Geomechanics_Module/Verification_Examples/

isotropic_compression

4 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D Axisymmetric.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Stationary.
6 Click Done.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

param 0 0 Parameter
p0 200[kPa] 2E5 Pa Initial pressure

DEFINITIONS

Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Local>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 In the table, enter the following settings:

t f(t)
0 1
0.2 1.8
0.4 2.6
0.6 1.8

5 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

t f(t)
0.8 2.6
1 3.4
4 Click Plot.
The interpolation function is used to define the boundary load. The boundary load first
compresses the soil sample, then it relaxes, and finally it compresses the sample again
3.4 times the initial load.

GEOMETRY 1

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 50e-3.
4 In the Height text field, type 100e-3.
5 Click Build All Objects.

SOLID MECHANICS (SOLID)

Cam-Clay Material 1
1 On the Physics toolbar, click Domains and choose Cam-Clay Material.
2 Click in the Graphics window and then press Ctrl+A to select all domains.
A drained case is assumed.
3 In the Settings window for Cam-Clay Material, locate the Model Input section.
4 In the pfluid text field, type 0.
5 Locate the Cam-Clay Material section. From the M list, choose Match to Mohr-
Coulomb criterion.
6 In the e0 text field, type 0.6646.
7 In the pC0 text field, type 400[kPa].
The consolidation pressure is twice the initial pressure. It means that first steps will
occur in the elastic regime.
The void ratio at consolidation pressure ec0 is deduced from the void ratio at reference
pressure N, the reference pressure prefN, and the initial consolidation pressure pc0 by
the following expression:

6 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

 ec0 = N + λ ln(prefN /pc0)

The initial void ratio e0 is then calculated from the void ratio at consolidation pressure
ec0, the initial consolidation pressure pc0 and the initial pressure p0.
e0 = ec0 + κ ln(pc0 /p0)

Initial Stress and Strain 1

1 On the Physics toolbar, click Attributes and choose Initial Stress and Strain.
2 In the Settings window for Initial Stress and Strain, locate the Initial Stress and Strain
section.
3 In the S0 table, enter the following settings:

-p0 0 0
0 -p0 0
0 0 -p0

In the initial state the soil is subjected to an hydrostatic pressure of p0 = 200 kPa.

Roller 1
1 On the Physics toolbar, click Boundaries and choose Roller.
2 Select Boundary 2 only.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundaries 3 and 4 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 From the Load type list, choose Pressure.
5 In the p text field, type int1(param)*p0.
The load function goes from p0 to 3.4*p0 with an unloading/reloading loop.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.

7 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

3 In the table, enter the following settings:

Property Name Value Unit Property group

Poisson’s ratio nu 0.2 1 Basic
Density rho 2400 kg/m³ Basic
Swelling index kappaSwelling 0.013 1 Cam-Clay
Compression index lambdaComp 0.032 1 Cam-Clay
Void ratio at reference Nvoid 0.7 1 Cam-Clay
pressure
Angle of internal friction internalphi 30[deg] rad Mohr-Coulomb

MESH 1
1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2 In the Settings window for Mesh, click Build All.

STUDY 1

Step 1: Stationary
Set up an auxiliary continuation sweep for the param parameter.

1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

param range(0.05,0.05,1)

The parametric solver is set to go from 0.05 to 1 with 20 steps.

6 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
The first default plot shows the von Mises stress at the final state in 2D.

Modify this plot group as follows to show the consolidation pressure.

1 In the Model Builder window, under Results click Stress (solid).

8 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

2 In the Settings window for 2D Plot Group, type Consolidation Pressure in the Label
text field.

Surface 1
1 In the Model Builder window, expand the Results>Consolidation Pressure node, then click
Surface 1.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Component 1>Solid Mechanics>
Material properties>[Link] - Consolidation pressure.
3 Locate the Expression section. From the Unit list, choose kPa.
4 On the Consolidation Pressure toolbar, click Plot.

Stress, 3D (solid)
The second default is a revolution in 3D of the first default plot.

Modify this plot group as follows to show the consolidation pressure in 3D.

1 In the Model Builder window, under Results click Stress, 3D (solid).

2 In the Settings window for 3D Plot Group, type Consolidation Pressure, 3D in the
Label text field.

Surface 1
1 In the Model Builder window, expand the Results>Consolidation Pressure, 3D node, then
click Surface 1.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Component 1>Solid Mechanics>
Material properties>[Link] - Consolidation pressure.
3 Locate the Expression section. From the Unit list, choose kPa.
4 On the Consolidation Pressure, 3D toolbar, click Plot.

1D Plot Group 3
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Void Ratio in the Label text field.
3 Click to expand the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Void ratio vs. log of the pressure during an
isotropic compression test..

5 Locate the Plot Settings section. Select the x-axis label check box.
6 In the associated text field, type Pressure with log scale.

9 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

7 Select the y-axis label check box.
8 In the associated text field, type Void ratio.

Point Graph 1
1 Right-click Void Ratio and choose Point Graph.
2 Select Point 1 only.
3 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid Mechanics>
Material properties>[Link] - Void ratio.
4 Locate the x-Axis Data section. From the Parameter list, choose Expression.
5 Click Replace Expression in the upper-right corner of the x-axis data section. From the
menu, choose Component 1>Solid Mechanics>Stress>[Link] - Pressure.
6 Locate the x-Axis Data section. In the Expression text field, type log([Link]).
7 On the Void Ratio toolbar, click Plot.

Point Graph 2
1 Right-click Results>Void Ratio>Point Graph 1 and choose Duplicate.
2 In the Settings window for Point Graph, locate the Data section.
3 From the Data set list, choose Study 1/Solution 1 (sol1).
4 From the Parameter selection (param) list, choose From list.
5 In the Parameter values (param) list, select 0.05.
6 Click to expand the Coloring and style section. Locate the Coloring and Style section. In
the Width text field, type 5.
7 Click to expand the Legends section. Select the Show legends check box.
8 From the Legends list, choose Manual.
9 In the table, enter the following settings:

Legends
param=0.05 -> p = 2.4e5 Pa

10 Duplicate Point Graph 2 five times and fill the settings according to the following table:

Name Param value Legend

Point Graph 2 0.05 param=0.05 -> p = 2.4e5 Pa
Point Graph 3 0.25 param=0.25 -> p = 4e5 Pa
Point Graph 4 0.4 param=0.4 -> p = 5.2e5 Pa

10 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

Name Param value Legend
Point Graph 5 0.6 param=0.6 -> p = 3.6e5 Pa
Point Graph 6 0.85 param=0.85 -> p = 5.6e5 Pa
Point Graph 7 1 param=1 -> p = 6.8e5 Pa

1D Plot Group 4
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Consolidation Pressure vs param
in the Label text field.
3 Locate the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Consolidation pressure vs Parameter value.

Point Graph 1
1 Right-click Consolidation Pressure vs param and choose Point Graph.
2 Select Point 1 only.
3 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid Mechanics>
Material properties>[Link] - Consolidation pressure.
4 Locate the y-Axis Data section. From the Unit list, choose kPa.
5 On the Consolidation Pressure vs param toolbar, click Plot.

11 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL

12 | ISOTROPIC COMPRESSION WITH CAM-CLAY MATERIAL MODEL
 Created in COMSOL Multiphysics 5.3

Tri axi al Test

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Introduction
The triaxial test is one of the most common tests used in laboratory for soil testing. The
soil sample is normally placed inside a rubber membrane and then compressed at constant
radial pressure.

In this example, a vertical displacement and a confinement pressure are applied on the
sample. The static response and the collapse load for various confinement pressures are
studied. This example is adapted from Ref. 1.

The material is modeled with the soil plasticity feature and the Drucker-Prager criterion.
The analysis can be simplified by considering the axial symmetry of the example.

Model Definition
A cylindrical apparatus with a 10 cm diameter and a 20 cm height, presses the soil sample
from the top by a prescribed a displacement. A flexible membrane contains the soil radially,
allowing changes in radial forces by controlling the surrounding pressure.
Prescribed displacement
Axial symmetry

20 cm Boundary load,
confinement pressure

5 cm

Figure 1: Dimensions, boundary conditions, and boundary load for the triaxial apparatus.

SOIL PROPERTIES
The soil proprieties are taken from a standard clay.

• Young’s modulus, E = 2.5 MPa, and Poisson’s ratio ν = 0.3.

2 | T R I A X I A L TE S T
• Cohesion c = 12 kPa, and angle of internal friction φ = 30° .
• Add a soil plasticity feature with Drucker-Prager criterion and select Match to Mohr-
Coulomb criterion. Add an Initial Stress and Strain node for including the confinement
pressure Pc as pore fluid pressure.

CONSTRAINTS AND LOADS

• Model only the right half of the domain due to the axial symmetry. Use the Axial
Symmetry boundary condition at the left vertical boundary.
• The soil sample is supported by a rigid and perfectly rough base. Use a Fixed Constraint
at the lower horizontal boundary.
• The right vertical boundary is modeled by a confinement pressure, add there a
Boundary Load. Simulate the influence of the confinement pressure through a
parametric sweep with values 5, 20, and 35 kPa.
• The soil sample is subjected to a loading at the top. Use a Prescribed Displacement
boundary, and gradually increase the vertical displacement up to 15 mm, by means of a
second parametric sweep.

Results and Discussion

After loading the soil sample it is possible to observe the distribution of plastic strains.
Most of the sample suffers from plastic deformation, only a minor part of it remains in the
elastic region, see Figure 2. The red zone shows the volume of soil that undergoes plastic
deformation. The 3D plot of the effective plastic strain is obtained from a revolution of the
2D axisymmetric data set.

3 | TR I A X I A L TE S T
Figure 2: Effective plastic strain in the soil sample.

In order to plot the additional loading stress on the soil sample caused by the prescribed
top displacement, integrate the z component of the reaction force over the top surface
with

F res =  σz 2πr dr
where the factor 2πr comes from the integration in cylindrical coordinates. This integral
is computed by the integration operator intop1, which applies the method Summation
over nodes for calculating the reaction force in the vertical direction.

F res = intop1 ( [Link] )

The extra loading stress (SI unit: Pa) is calculated by the difference between the resulting
stress and the confinement (radial) pressure

F res
σ a = – ----------2- – Pc
πR

here, R is radius of the triaxial apparatus.

4 | T R I A X I A L TE S T
The extra loading stress is the resulting stress on the porous matrix. It is shown in Figure 3
for the three different confinement pressures.

Figure 3: The extra loading stress for different confinement pressures.

Reference
1. D. Potts and L. Zdravkovic, Finite Element Analysis in Geotechnical Engineering,
Thomas Telford Publishing, 2001.

Application Library path: Geomechanics_Module/Verification_Examples/

triaxial_test

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

5 | TR I A X I A L TE S T
MODEL WIZARD
1 In the Model Wizard window, click 2D Axisymmetric.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select Preset Studies>Stationary.
6 Click Done.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

Pc 5[kPa] 5000 Pa Confinement pressure
Disp 1[mm] 0.001 m Prescribed displacement
D 100[mm] 0.1 m Diameter of the sample
H 200[mm] 0.2 m Height of the sample

GEOMETRY 1
1 In the Model Builder window, expand the Component 1 (comp1)>Geometry 1 node, then
click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type D/2.
4 In the Height text field, type H.
5 Click Build All Objects.

6 | T R I A X I A L TE S T
DEFINITIONS

Integration 1 (intop1)
1 On the Definitions toolbar, click Component Couplings and choose Integration.
2 In the Settings window for Integration, locate the Source Selection section.
3 From the Geometric entity level list, choose Boundary.
4 Select Boundary 3 only.
5 Locate the Advanced section. From the Method list, choose Summation over nodes.
The integral operator intop1 is used for calculating the resulting strength over the
whole soil sample. Since the reaction force is defined only at nodes, choose summation
over nodes as integration method.

Variables 1
1 On the Definitions toolbar, click Local Variables.
2 In the Settings window for Variables, locate the Variables section.
3 In the table, enter the following settings:

Name Expression Unit Description

Fres intop1([Link]) Resulting strength

SOLID MECHANICS (SOLID)

Linear Elastic Material 1

In the Model Builder window, expand the Component 1 (comp1)>Solid Mechanics (solid)
node, then click Linear Elastic Material 1.

Soil Plasticity 1
1 On the Physics toolbar, click Attributes and choose Soil Plasticity.
2 In the Settings window for Soil Plasticity, locate the Soil Plasticity section.
3 Select the Match to Mohr-Coulomb criterion check box.

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click
Linear Elastic Material 1.

External Stress 1
1 On the Physics toolbar, click Attributes and choose External Stress.
2 In the Settings window for External Stress, locate the External Stress section.
3 In the Sext text field, type -Pc.

7 | TR I A X I A L TE S T
Fixed Constraint 1
1 On the Physics toolbar, click Boundaries and choose Fixed Constraint.
2 Select Boundary 2 only.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 4 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

-Pc r
0 z

Prescribed Displacement 1
1 On the Physics toolbar, click Boundaries and choose Prescribed Displacement.
2 Select Boundary 3 only.
3 In the Settings window for Prescribed Displacement, locate the Prescribed Displacement
section.
4 Select the Prescribed in r direction check box.
5 Select the Prescribed in z direction check box.
6 In the u 0z text field, type -Disp.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 2.5[MPa] Pa Basic
Poisson’s ratio nu 0.3 1 Basic
Cohesion cohesion 12[kPa] Pa Mohr-Coulomb
Angle of internalphi 30[deg] rad Mohr-Coulomb
internal friction

8 | T R I A X I A L TE S T
MESH 1
Create mapped mesh with element size of 2 mm.

Size 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Mapped.
2 Right-click Mapped 1 and choose Size.
3 In the Settings window for Size, locate the Element Size section.
4 Click the Custom button.
5 Locate the Element Size Parameters section. Select the Maximum element size check box.
6 In the associated text field, type 2[mm].
7 Click Build All.

STUDY 1
Use a parametric sweep to test 3 different confinement pressures.

Parametric Sweep
1 On the Study toolbar, click Parametric Sweep.
2 In the Settings window for Parametric Sweep, locate the Study Settings section.
3 Click Add.
4 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

Pc 5[kPa] 20[kPa] 35[kPa]

Step 1: Stationary
1 In the Model Builder window, under Study 1 click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
Use an auxiliary continuation sweep for the Disp parameter to gradually increase the
displacement.

3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

Disp 15[mm]*range(0,0.05,1)

9 | TR I A X I A L TE S T
6 On the Study toolbar, click Compute.

RESULTS

Stress (solid)
The first default plot shows the von Mises stress in 2D.

Stress, 3D (solid)
The second default plot shows the von Mises stress in 3D obtained by a revolution of the
2D axisymmetric data set.

3D Plot Group 3
1 On the Home toolbar, click Add Plot Group and choose 3D Plot Group.
2 In the Settings window for 3D Plot Group, type Plastic Region in the Label text field.

Surface 1
1 Right-click Plastic Region and choose Surface.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Model>Component 1>Solid Mechanics>
Strain>[Link] - Effective plastic strain.
3 Locate the Expression section. Select the Description check box.
4 In the Expression text field, type [Link]>0.
5 On the Plastic Region toolbar, click Plot.
6 Click the Zoom Extents button on the Graphics toolbar.
The 3D image of the effective plastic strain is obtained by first generating a revolution
of the 2D axisymmetric data set.
Create a 1D plot to show the extra loading stress for different confinement pressures.

1D Plot Group 4
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Extra Loading Stress in the Label text
field.
3 Locate the Data section. From the Data set list, choose Study 1/
Parametric Solutions 1 (sol2).

Global 1
1 Right-click Extra Loading Stress and choose Global.
2 In the Settings window for Global, locate the y-Axis Data section.

10 | TR I A X I A L TE S T
3 In the table, enter the following settings:

Expression Unit Description

-Fres/(pi*(D/2)^2)-Pc kPa Extra loading stress

4 Locate the x-Axis Data section. From the Parameter list, choose Expression.
5 In the Expression text field, type Disp/H.
6 Select the Description check box.
7 In the associated text field, type Global axial strain.
8 On the Extra Loading Stress toolbar, click Plot.

Extra Loading Stress

1 In the Model Builder window, under Results click Extra Loading Stress.
2 In the Settings window for 1D Plot Group, click to expand the Legend section.
3 From the Position list, choose Lower right.
4 On the Extra Loading Stress toolbar, click Plot.

11 | TR I A X I A L TE S T
12 | TR I A X I A L TE S T
 Created in COMSOL Multiphysics 5.3

Tu n n el E xc av at i on

This model is licensed under the COMSOL Software License Agreement 5.3.
All trademarks are the property of their respective owners. See [Link]/trademarks.
Introduction
This example simulates the behavior of the soil during a tunnel excavation. The surface
settlement and the width of the plastic region around the tunnel are important parameters
required to predict the necessary reinforcements during the excavation. This verification
example was adapted from Ref. 1 and Ref. 2.

In order to calculate in situ stresses, use two study steps. In the first study compute the
stress state of the soil before the excavation of the tunnel. In the second study compute
the elastoplastic behavior once the soil is removed. This requires incorporation of the stress
response calculated in the first step.

In order to speed up the calculation consider the soil in the first step as elastic and in the
second step add a soil plasticity material model Drucker-Prager. The example is solved in
2D plane strain.

Model Definition
The geometry consists of a soil layer that is 45 m deep and 90 m wide. A tunnel of 10 m
in diameter is placed at the symmetry axis, 20 m below the surface. A bed rock, 45 m
below the surface, constrains the displacement in the vertical direction, and a roller
boundary is used for simulating the infinite extension of the soil in the lateral direction.

20 m
Symmetry

45 m

5m

90 m

Figure 1: Dimensions and boundary conditions for the tunnel excavation example.

2 | TU N N E L E X C A V A T I O N
SOIL PROPERTIES
The soil properties are adapted from Ref. 2.

• Young’s modulus, E = 12 MPa and Poisson’s ratio ν = 0.495.

• Cohesion c = 130 kPa and angle of internal friction φ = 30° .
• Use the Drucker-Prager criterion and match the material parameters to the Mohr-
Coulomb criterion.

CONSTRAINTS AND LOADS

• At the lower boundary fix the displacement with a Fixed Constraint boundary.
• Use Symmetry at the left boundary and Roller at the right boundary.
• Keep the default Free boundary at the top. Keep a Free boundary on the tunnel’s wall
after removing the soil from this domain.
• Add a Gravity node to account for gravity effects.

3 | TU N N E L E X C A V A T I O N
Results and Discussion
The first plot (Figure 2) shows the stress distribution due to gravity. The roller and
symmetry boundaries create an homogeneous vertical variation.

Figure 2: The von Mises stress in the soil layer before the excavation of the tunnel.

The second plot (Figure 3) shows the stress distribution after excavating the tunnel. In situ
stresses are taken from the first step. Note the increase in the effective stress around the
tunnel.

4 | TU N N E L E X C A V A T I O N
Figure 3: The von Mises stress in the soil layer after excavation of the tunnel.

In the second step, beside removing the tunnel domain, a soil plasticity feature is added.
In Figure 4 the region that experience plasticity is shown.

5 | TU N N E L E X C A V A T I O N
Figure 4: The plastic deformation in the zone near the tunnel after the excavation.

The horizontal displacement and the subsidence of the top surface due to the excavation
is shown in Figure 5 and Figure 6.

6 | TU N N E L E X C A V A T I O N
Figure 5: The horizontal displacement at the top surface

7 | TU N N E L E X C A V A T I O N
Figure 6: The surface subsidence.

References
1. D. Potts and L. Zdravkovic, Finite Element Analysis in Geotechnical Engineering,
Thomas Telford Publishing, 2001.

2. H. Schweiger, “Results from Numerical Benchmark Exercises in Geotechnics,” Proc.

5th European Conference in Numerical Methods in Geotechnical Engineering,
pp. 305–314, 2002.

Application Library path: Geomechanics_Module/Soil/tunnel_excavation

Modeling Instructions
From the File menu, choose New.

8 | TU N N E L E X C A V A T I O N
NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
5 Click Add.
6 Click Study.
7 In the Select Study tree, select Preset Studies for Selected Physics Interfaces>Stationary.
8 Click Done.

GEOMETRY 1

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 90.
4 In the Height text field, type 45.
5 Locate the Position section. In the y text field, type -45.
6 Right-click Rectangle 1 (r1) and choose Build Selected.

Circle 1 (c1)
1 On the Geometry toolbar, click Primitives and choose Circle.
2 In the Settings window for Circle, locate the Size and Shape section.
3 In the Radius text field, type 5.
4 In the Sector angle text field, type 180.
5 Locate the Position section. In the y text field, type -20.
6 Locate the Rotation Angle section. In the Rotation text field, type 270.
7 Right-click Circle 1 (c1) and choose Build Selected.

Form Union (fin)

In the Model Builder window, under Component 1 (comp1)>Geometry 1 right-click
Form Union (fin) and choose Build Selected.

9 | TU N N E L E X C A V A T I O N
Set the first step with full geometry and linear elastic material.

SOLID MECHANICS (SOLID)

Linear Elastic Material 1

Since in this example, the Poisson’s ratio is 0.495, select the Nearly incompressible check-
box to avoid locking effects.

1 In the Model Builder window, expand the Component 1 (comp1)>Solid Mechanics (solid)
node, then click Linear Elastic Material 1.
2 In the Settings window for Linear Elastic Material, locate the Linear Elastic Material
section.
3 Select the Nearly incompressible material check box.

Symmetry 1
1 In the Model Builder window, right-click Solid Mechanics (solid) and choose
More Constraints>Symmetry.
2 Click the Select Box button on the Graphics toolbar.
3 Select Boundaries 1 and 3–5 only.

Fixed Constraint 1
1 Right-click Solid Mechanics (solid) and choose Fixed Constraint.
2 Select Boundary 2 only.

Roller 1
1 Right-click Solid Mechanics (solid) and choose Roller.
2 Select Boundary 7 only.

Gravity 1
1 Right-click Solid Mechanics (solid) and choose Volume Forces>Gravity.
2 In the Settings window for Gravity, locate the Domain Selection section.
3 From the Selection list, choose All domains.

Set the second step with modified geometry and elasto-plastic model with yield criterion
according to Drucker-Prager.

SOLID MECHANICS 2 (SOLID2)

1 In the Model Builder window, under Component 1 (comp1) click
Solid Mechanics 2 (solid2).
2 Select Domain 1 only.

10 | TU N N E L E X C A V A T I O N
Linear Elastic Material 1
1 In the Model Builder window, expand the Solid Mechanics 2 (solid2) node, then click
Linear Elastic Material 1.
2 In the Settings window for Linear Elastic Material, locate the Linear Elastic Material
section.
3 Select the Nearly incompressible material check box.

Soil Plasticity 1
1 Right-click Component 1 (comp1)>Solid Mechanics 2 (solid2)>Linear Elastic Material 1
and choose Soil Plasticity.
2 In the Settings window for Soil Plasticity, locate the Soil Plasticity section.
3 Select the Match to Mohr-Coulomb criterion check box.

External Stress 1
1 Right-click Linear Elastic Material 1 and choose External Stress.
2 In the Settings window for External Stress, locate the External Stress section.
3 From the Sext list, choose Second Piola-Kirchhoff stress (solid/lemm1).

Symmetry 1
1 In the Model Builder window, right-click Solid Mechanics 2 (solid2) and choose
More Constraints>Symmetry.
2 Select Boundaries 1 and 5 only.

Fixed Constraint 1
1 Right-click Solid Mechanics 2 (solid2) and choose Fixed Constraint.
2 Select Boundary 2 only.

Roller 1
1 Right-click Solid Mechanics 2 (solid2) and choose Roller.
2 Select Boundary 7 only.

Gravity 1
1 Right-click Solid Mechanics 2 (solid2) and choose Volume Forces>Gravity.
2 In the Settings window for Gravity, locate the Domain Selection section.
3 From the Selection list, choose All domains.

11 | TU N N E L E X C A V A T I O N
MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young’s modulus E 12e6 Pa Basic
Poisson’s ratio nu 0.495 1 Basic
Density rho 2000 kg/m³ Basic
Cohesion cohesion 130e3 Pa Mohr-Coulomb
Angle of internalphi 30[deg] rad Mohr-Coulomb
internal friction

MESH 1

Size 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Free Triangular.
2 Right-click Free Triangular 1 and choose Size.
3 In the Settings window for Size, locate the Element Size section.
4 From the Predefined list, choose Finer.

Distribution 1
1 Right-click Free Triangular 1 and choose Distribution.
2 Select Boundaries 8 and 9 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 12.
5 Click Build All.

Use two stationary study steps. The first one is used to compute the initial stress state of
the soil. The second step is used to compute the elasto-plastic deformation due to the
excavation of the tunnel.

12 | TU N N E L E X C A V A T I O N
STUDY 1

Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, locate the Physics and Variables Selection section.
3 In the table, clear the Solve for check box for the Solid Mechanics 2 interface.

Step 2: Stationary 2
1 On the Study toolbar, click Study Steps and choose Stationary>Stationary.
2 In the Settings window for Stationary, locate the Physics and Variables Selection section.
3 In the table, clear the Solve for check box for the Solid Mechanics interface.
4 On the Study toolbar, click Compute.

RESULTS

Stress (solid)
The first default plot shows the von Mises stress in the soil at the initial state.

1 In the Model Builder window, under Results click Stress (solid).

2 In the Settings window for 2D Plot Group, type Stress, before excavation in the
Label text field.

Surface 1
1 In the Model Builder window, expand the Results>Stress, before excavation node, then
click Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 Select the Description check box.

Deformation
1 In the Model Builder window, expand the Surface 1 node, then click Deformation.
2 In the Settings window for Deformation, locate the Scale section.
3 Select the Scale factor check box.
4 In the associated text field, type 1.
5 On the Stress, before excavation toolbar, click Plot.
6 Click the Zoom Extents button on the Graphics toolbar.

Stress (solid2)
The second default plot shows the von Mises stress in the soil after the tunnel excavation.

13 | TU N N E L E X C A V A T I O N
1 In the Model Builder window, under Results click Stress (solid2).
2 In the Settings window for 2D Plot Group, type Stress, after excavation in the
Label text field.

Surface 1
1 In the Model Builder window, expand the Results>Stress, after excavation node, then
click Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 Select the Description check box.

Deformation
1 In the Model Builder window, expand the Surface 1 node, then click Deformation.
2 In the Settings window for Deformation, locate the Scale section.
3 Select the Scale factor check box.
4 In the associated text field, type 1.
5 On the Stress, after excavation toolbar, click Plot.

2D Plot Group 3
1 On the Home toolbar, click Add Plot Group and choose 2D Plot Group.
Use this plot group to show the plastic zone after the excavation of the tunnel.
2 In the Settings window for 2D Plot Group, type Plastic Region in the Label text field.

Surface 1
1 Right-click Plastic Region and choose Surface.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner of
the Expression section. From the menu, choose Component 1>Solid Mechanics 2>Strain>
[Link] - Effective plastic strain.
3 Locate the Expression section. Select the Description check box.
4 In the Expression text field, type [Link]>0.
This is a boolean expression which is 1 in the plastic region and 0 elsewhere.

Deformation 1
1 Right-click Results>Plastic Region>Surface 1 and choose Deformation.
2 In the Settings window for Deformation, click Replace Expression in the upper-right
corner of the Expression section. From the menu, choose Component 1>
Solid Mechanics 2>Displacement>u2,v2 -
Displacement field (material and geometry frames).

14 | TU N N E L E X C A V A T I O N
3 Locate the Scale section. Select the Scale factor check box.
4 In the associated text field, type 1.
5 On the Plastic Region toolbar, click Plot.

1D Plot Group 4
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Horizontal Displacement in the Label
text field.
3 Click to expand the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Horizontal displacement at surface.
5 Locate the Plot Settings section. Select the x-axis label check box.
6 In the associated text field, type Distance from tunnel axis (m).
7 Select the y-axis label check box.
8 In the associated text field, type Horizontal displacement (mm).

Line Graph 1
1 Right-click Horizontal Displacement and choose Line Graph.
2 Select Boundary 6 only.
3 In the Settings window for Line Graph, locate the y-Axis Data section.
4 In the Expression text field, type u2.
5 From the Unit list, choose mm.
6 On the Horizontal Displacement toolbar, click Plot.

1D Plot Group 5
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Vertical Displacement in the Label
text field.
3 Locate the Title section. From the Title type list, choose Manual.
4 In the Title text area, type Surface settlement.
5 Locate the Plot Settings section. Select the x-axis label check box.
6 In the associated text field, type Distance from tunnel axis (m).
7 Select the y-axis label check box.
8 In the associated text field, type Vertical displacement (mm).

15 | TU N N E L E X C A V A T I O N
Line Graph 1
1 Right-click Vertical Displacement and choose Line Graph.
2 Select Boundary 6 only.
3 In the Settings window for Line Graph, locate the y-Axis Data section.
4 In the Expression text field, type v2.
5 From the Unit list, choose mm.
6 On the Vertical Displacement toolbar, click Plot.

16 | TU N N E L E X C A V A T I O N