Sequential Tunnel excavation

Updated 02/2024

Sequential Tunnel Excavation

Created By: Roozbeh Geraili Mikola, PhD, PE
Email: adonis4geo@[Link]
Web Page: [Link]
2m

11 m 23 m

15 m

120 m

This tutorial illustrates the use of ADONIS for the analysis of the construction of a SEM/NATM
tunnel. The SEM/NATM is a technique in which ground exposed by excavation is stabilized with
shotcrete to form a temporary lining.

Select Unit
Since we will be using Metric (Pa) units for this tutorial, let’s select “SI: m-Pa-N/m3” from
Calculation Setting under Setting button.

Select:  

Create Geometry
To create the model geometry select the Create Geometry button, then select Draw Rectangle
and enter the first- and second-points’ coordinates (-60,-45 and 60,0) and press Apply.

Select:  

Select Draw Line to draw material lines as listed below:

1
 Sequential Tunnel excavation
Updated 02/2024

Type Start Point End Point

-60,-2 60,-2
Line
-60,-25 60,-25

Select Draw Arc to draw material lines as listed below:

Type Start Point Middle Point End Point Num. of Segments

3.88,-20.23 4.69,-19.62 5.25,-18.77 6

4.95,-13.87 0.00,-11.00 -4.95,-13.87 25

-5.25,-18.77 -4.69,-19.62 -3.88,-20.23 6

-3.88,-20.23 0.00,-21.00 3.88,-20.23 20

Arc 5.25,-18.77 5.57,-17.55 5.63,-16.29 6

5.63,-16.29 5.42,-15.05 4.95,-13.87 6

-5.25,-18.77 -5.57,-17.55 -5.63,-16.29 6

-5.63,-16.29 -5.42,-15.05 -4.95,-13.87 6

-5.63,-16.29 0.00,-17.94 5.63,-16.29 25

Figure 1- Model geometry.

Mesh
Now let’s generate the finite element (FE) mesh. In ADONIS, meshing is a simple two-step
process. First you must discretize the boundaries, and then the Mesh can be generated. You can
also configure various Mesh Setup parameters before generating the mesh.
Select the Discretize/Mesh button in command panel or toolbar or the Mesh menu, then select
Mesh Property button then select the element type. Currently only two types of FE mesh is
available 1) 3-Node triangular and 2) 4-Node quadrilateral elements. For this tutorial 3-Node
triangular element is selected. Then select the Discretize/Mesh button in command panel or
toolbar or the Mesh menu, then select Discretize/Mesh button. Enter 1.0 as Max. Edge Size in
order to create finer meh size. Now press Discretize then press Generate Mesh buttons to
discretize the boundaries and generate FE mesh.

2
 Sequential Tunnel excavation
Updated 02/2024

Select:  

The finite element mesh is generated. ADONIS adopts “Nodal Mixed Discretization (NMD)
methodology to overcome this problem. NMD is introduced by Detourney & Dzik (2006) as an
improvement of Mixed Discretization (MD) technique.
Notes:
1. NMD is only available for T3 element (not Q4).
2. NMD is activated by default.

Assign Soil Properties

To create a new material, select the Assign Material/Excavate button from command panel,
toolbar or menu, then select Soil/Rock tab. From drop down select Mohr-Coulomb criteria then
click on Add Soil/Rock Material button.

Select: 

Now enter the following soil and rock properties for each layer as shown in Figures 2.
Table 1- Material properties.

Young Friction Tensile

Material Material Density Poisson’s Cohesion
Modulus Angle Strength K0
ID Name (kg/m3) Ratio (Pa)
(Pa) (deg.) (Pa)

1 Fill 1950 1.0x107 0.30 1.0x103 29 0.0 0.5

2 Soil 2200 1.0x108 0.25 3.9x104 35 0.0 0.6

3 Rock 2400 5.0x108 0.25 1.0x105 40 1.0x104 0.65

3
 Sequential Tunnel excavation
Updated 02/2024

Figure 2- Create Mohr-Coulomb properties for each material.

To assign the created material to the proper region, select the Fill material (i.e. created material)
from the material list and then use mouse and click on top layer in the main window and click on
Apply Button, or press Enter (Figure 2). Repeat the sequence to create material for the rest of
layers. Figure 3 shows the model configuration after assigning the materials.

Figure 3- Model configuration after assigning the materials.

Boundary Condition
In order to fix (i.e. zero displacement) the external boundary, select Initial Condition then select
Apply Boundary Condition then from the ”Keyword List” click on xyfix under Fix keyword. Then
Use the mouse to select the nodes located on bottom of the model. When finished, click on Apply
Button, or press Enter.

Select:  

Then select xfix and use the mouse to select the nodes located on left and right sides of the
model. The boundary condition has been applied as shown in Figure 4.

4
 Sequential Tunnel excavation
Updated 02/2024

Figure 4- Model configuration after applying boundary condition.

Initialize the Stress Condition

In order to initialize the vertical and horizontal stress in the model select Initial Condition then
select Apply Initial Condition then from the ”Keyword List” click on syy under Element Stress.

Select: 

In addition to constant stress, the Initial command will also take the keywords xvar and yvar
which allow the initial stress to vary linearly over the specified range. The stresses vary linearly
with distance from the global coordinate origin of x=0, y=0:
value_modified = value + vx × x + vy × y
where value is the stress component at origin and (x,y) is position of the gauss point in the global
coordination. For example, the vertical stress at the top and bottom of the Fill layer should be 0
and γFill×h (e.g. 1950×9.81×2.0=38259 Pa). To specify such stress gradient, the initial stress
and yvar values should be 0 and 17003 (i.e. value_modified = 0 + 0 × x + γFill × y) as shown in
Figure 5. The same procedure will be applied for the initialization of the syy, sxx and szz for both
Clay and Sand layers. Table 2 summarizes the values required for each initialization command.
The yvar value for sxx and szz should be applied based on each layer’s K0 value as shown in
Table 2. Please make sure that proper range is used for each layer before click on Apply.

Note: Compression stresses have a negative sign, in accordance with the general sign
convention for internal stresses in ADONIS.

5
 Sequential Tunnel excavation
Updated 02/2024

Figure 5- Stress initialization dialog box, initializing vertical stress for each layer.
Table 2- List of parameters for stress initialization.

Stress Type Initial Value (Pa) Y-Variation Value (Pa) X-Limit (m) Y-Limit (m)

syy 0 19130 -60,60 -2,0

syy 4905 21582 -60,60 -25,-2

syy 53955 23544 -60,60 -40,-25

sxx 0 9564.8 -60,60 -2,0

sxx 2943 12949 -60,60 -25,-2

sxx 35071 15304 -60,60 -40,-25

szz 0 9564.8 -60,60 -2,0

szz 2943 12949 -60,60 -25,-2

szz 35071 15304 -60,60 -40,-25

Before solving the mode to equilibrium let’s specify the gravitational gradient value. For that,
select Setting then Gravity Setting and set the magnitude of gravity in y-direction to 9.81 m/s2.
Note: Downward gravity is positive.

Select:  

Now click on Solve the Model button from toolbar, menu or command panel then select Solve
radio button from the Run Type and click on Run button. The ADONIS Compute engine will
proceed in running the analysis. Model should reach equilibrium after several steps.

Select: 

6
 Sequential Tunnel excavation
Updated 02/2024

After reaching to equilibrium let’s zero out the displacement in the model to get ready for the next
stages. For this, select Initial Condition then select Apply Initial Condition then from the
”Keyword List” click on xydisp (means both x- and y- displacements) under Nodal Displacement
keyword and enter 0.0 in XY-Displacement value box then press Apply. This way you can reset
both x- and y- displacements to zero.

Select:  

Perform the Excavation Sequence

A staged construction calculation is needed in which the tunnel lining is activated and the soil
clusters inside the tunnel are deactivated. The calculation phases are Plastic analyses, Staged
construction. The three-dimensional arching effect is emulated by using the so-called β-method.
The idea is that the initial stresses pk acting around the location where the tunnel is to be
constructed are divided into a part (1 - β) pk that is applied to the unsupported tunnel and a part
Deconfinement method that is applied to the supported tunnel. To define the calculation process,
follow these steps:
Stage 1: Excavate top heading
To excavate the top heading, select the Assign Material/Excavate button from command
panel, toolbar or menu, then select Excavate from the material name list. Then use mouse
and select the top heading region and press Apply.
Stage 2: Perform stress relaxation
To perform stress relaxation in ADONIS, one can use the “solve relax” command, which
is available under Run Type in the Calculation panel. In order to perform stress relaxation
let’s select the Solve Relax and modify the parameters. For this tutorial the assumption
is that 60% of initial stress in the excavated soil should disappear before temporary
support installation (so the remaining 40% is to be considered later). To achieve this the
relax factor should be set to 0.6. The number of steps for stress relaxation to gradually be
applied to the tunnel boundary is set to 250 by default. Then using mouse select the
boundary that stress relaxation would be applied at (i.e. excavation boundary) or directly
specify the range (i.e. xlim=[-8.5,8.5] and ylim=[-22.0,-8.0]) and click on Run button to
start the cycle.
Notes: Relax step can be increased to have more smooth gradual stress relaxation but
250 step seem to be working for most of the cases.

7
 Sequential Tunnel excavation
Updated 02/2024

Select: 

Figure 6- Use “solve relax” command to relax stress on the boundary of tunnel.
Stage 3: Install shotcrete support for top heading
In this stage 100 mm thick shotcrete (Table 3) will be installed at the crown of tunnel. To
do that, select the Structures button from command panel, toolbar or menu, then select
Beam element. Select Edge from Type group-box and after that Range from selection
mode. Change the interface type to Both Sides from the dropdown. Now use mouse to
select the range (or type xlim=[-8.0,8.0] and ylim=[-19.0,-10.0]) that shotcrete support
should be installed and press Apply (Figure 7). Then click on Beam Element Properties
and assign the shotcrete properties as shown in Figure 8. Then select Solve the Model
button from toolbar, menu or command panel and solve the model until reaching to
equilibrium (i.e. select Solve and press Run button).

Note: In ADONIS, the beam element formulation is a plane-stress formulation. If the beam
is representing a structure that is continuous in the direction perpendicular to the analysis
plane (e.g., a concrete tunnel lining), the value specified for E should be divided by (1 −
ν2) to account for plane-strain conditions.
Table 3- Shotcrete properties.

Thickness (m) Young Modulus (Pa) Poisson’s Ratio

0.1 15.0x109 0.2

8
 Sequential Tunnel excavation
Updated 02/2024

Figure 7- Install beam elements on tunnel’s crown.

Figure 8- Assign shotcrete properties.

Stage 4: Excavate Bench
To excavate the bench, select the Assign Material/Excavate button from command
panel, toolbar or menu, then select Excavate from the material name list. Then use mouse
and select the bench region and press Apply.
Stage 5: Install shotcrete support for bench
Repeat the sequence explained at stage 3 and install the beam element at the excavated
part of the tunnel (i.e. bench region) as shown in Figure 9. There is no need for specifying
the beam properties since similar shotcrete parameters will be used for the bench. Then
select Solve the Model button from toolbar, menu or command panel and solve the model
until reaching to equilibrium (i.e. select Solve and press Run button).

Figure 9- Install beam elements on tunnel’s bench.

Post Processing

9
 Sequential Tunnel excavation
Updated 02/2024

After reaching to equilibrium in the final stage, the tools under the Plot Geometry/Results in the
command panel can be used to plot the results. An extensive set of options is available for plotting
the results of the finite element analysis generated by ADONIS. These include contouring of the
nodal and gauss point variables, strain and stress values, forces in the structural elements and
highlighting of plastic zones. Lines or color fills can be selected for all contour plots. To plot
horizontal displacement contour, click on Plot Geometry/Results button then click on Plot
Contour, under Displacement select ydisp and press Apply.

Select:  

Now in order to plot the axial force in the beam elements click on Plot Structural Data and select
axialforce under Beam and press Apply.

Select:  

Figure 10- Vertical displacement contour and beam axial force plots.
The same process can be followed to plot the vertical displacement contour and beam elements
moment force as shown in Figures 10 and 11.

10
 Sequential Tunnel excavation
Updated 02/2024

Figure 11- Horizontal displacement contour and beam moment force plots.

Model’s Script
The commands for this tutorial are listed below:
// create new/fresh model
newmodel()

// set unit
set("unit","stress-pa")

// create geometry
rect("startPoint",-60,-40,"endPoint",60,0)
line("startPoint",-60,-2,"endPoint",60,-2)
line("startPoint",-60,-25,"endPoint",60,-25)

arc("startPoint",3.88,-20.23,"midPoint",4.69,-19.62,"endPoint",5.25,-18.77,"numSeg",6)
arc("startPoint",4.95,-13.87,"midPoint",0.00,-11.00,"endPoint",-4.95,-13.87,"numSeg",25)
arc("startPoint",-5.25,-18.77,"midPoint",-4.69,-19.62,"endPoint",-3.88,-20.23,"numSeg",6)
arc("startPoint",-3.88,-20.23,"midPoint",0.00,-21.00,"endPoint",3.88,-20.23,"numSeg",20)
arc("startPoint",5.25,-18.77,"midPoint",5.57,-17.55,"endPoint",5.63,-16.29,"numSeg",6)
arc("startPoint",5.63,-16.29,"midPoint",5.42,-15.05,"endPoint",4.95,-13.87,"numSeg",6)
arc("startPoint",-5.25,-18.77,"midPoint",-5.57,-17.55,"endPoint",-5.63,-16.29,"numSeg",6)
arc("startPoint",-5.63,-16.29,"midPoint",-5.42,-15.05,"endPoint",-4.95,-13.87,"numSeg",6)
arc("startPoint",-5.63,-16.29,"midPoint",0.00,-17.94,"endPoint",5.63,-16.29,"numSeg",25)

// discretize/mesh
discretize("maxedge",1.0)
gmsh("maxedge",1.0,"elemtype","T3","useNMD","on")

// create new materials

material("create","Mohr-
Coulomb","matid",1,"matname","Fill","density",1950,"shear",3.84615e+06,"bulk",8.33333e+06,"coh",1
e3,"fric",29,"dil",0,"tens",0)
material("create","Mohr-
Coulomb","matid",2,"matname","Soil","density",2200,"shear",4e+07,"bulk",6.66667e+07,"coh",3.9e4,"
fric",35,"dil",0,"tens",0)
material("create","Mohr-
Coulomb","matid",3,"matname","Rock","density",2400,"shear",2e+08,"bulk",3.33333e+08,"coh",1e5,"fr
ic",40,"dil",0,"tens",1.0e4)

// assign materials
material("assign","matid",1,"region",0.0,-1.0)
material("assign","matid",2,"region",0.0,-9.0)
material("assign","matid",2,"region",0.0,-13.0)
material("assign","matid",2,"region",0.0,-18.5)
material("assign","matid",3,"region",0.0,-32.0)

// assign boundary condition

applybc("xfix","xlim",-60.1,-59.9,"ylim",-40.1,0.1)
applybc("xfix","xlim",59.9,60.1,"ylim",-40.1,0.1)
applybc("xyfix","xlim",-60.1,60.1,"ylim",-40.1,-39.9)

// initialize stress condition

11
 Sequential Tunnel excavation
Updated 02/2024

var grav = 9.81

var dsyy1 = 1950*grav
var isyy1 = 0.0
var dsyy2 = 2200*grav
var isyy2 = dsyy2*2.0 - dsyy1*2.0
var dsyy3 = 2400*grav
var isyy3 = dsyy3*(2.0+23.0) - dsyy1*2.0 - dsyy2*23.0
var k01 = 0.5
var k02 = 0.6
var k03 = 0.65
initial("syy",isyy1,"yvar",dsyy1,"xlim",-60,60,"ylim",-2,0)
initial("syy",isyy2,"yvar",dsyy2,"xlim",-60,60,"ylim",-25,-2)
initial("syy",isyy3,"yvar",dsyy3,"xlim",-60,60,"ylim",-40,-25)
initial("sxx",k01*isyy1,"yvar",k01*dsyy1,"xlim",-60,60,"ylim",-2,0)
initial("sxx",k02*isyy2,"yvar",k02*dsyy2,"xlim",-60,60,"ylim",-25,-2)
initial("sxx",k03*isyy3,"yvar",k03*dsyy3,"xlim",-60,60,"ylim",-40,-25)
initial("szz",k01*isyy1,"yvar",k01*dsyy1,"xlim",-60,60,"ylim",-2,0)
initial("szz",k02*isyy2,"yvar",k02*dsyy2,"xlim",-60,60,"ylim",-25,-2)
initial("szz",k03*isyy3,"yvar",k03*dsyy3,"xlim",-60,60,"ylim",-40,-25)

// apply gravity force

set("gravity",0,grav)

// solve to initialize
solve()

// reset the displacement to zero

initial("xydisp",0)

// excavate the top region (sand)

excavate("region",0.0,-13.0)

// solve with 60% stress relaxation (40% of load remains)

solve("relax","relaxFactor",0.6,"relaxStep",250,"xlim",-8.5,8.5,"ylim",-22.0,-8.0)

// install shotcrete
var beam_nu = 0.2; // shotcrete poisson's ratio
var beam_th = 0.1; // shotcrete thickness
var beam_b = 1.0; // out of plane length
var beam_area = beam_th * beam_b; // shotcrete area
var beam_I = (beam_b*[Link](beam_th,3))/12.0; // shotcrete moment of inertia
var beam_ymod = 15e9 / ([Link](beam_nu,2)); // E should be divided by (1 − nu^2) to account
for plane-strain conditions.
structure("drawbeam","beamid",1,"xlim",-6.5,6.5,"ylim",-16.3,-10.5)
structure("material","beamid",1,"area",beam_area,"I",beam_I,"ymod",beam_ymod)

// solve to equilibrium
solve()

// excavate bench
excavate("region",0.0,-19.5)

// install shotcrete
structure("drawbeam","beamid",1,"xlim",-6.5,6.5,"ylim",-21.0,-16.3)
// solve to equilibrium
solve()

// plot vertical displacement contour and beam axial forces

plot("contour","ydisp")
plot("struc","beam","axialforce")

// plot horizontal displacement contour and beam moment forces

tab("plot")
plot("contour","xdisp")
plot("struc","beam","moment")

12