Tunnel excavation under in-situ stress#

In-situ stress is an important issue that significantly influence the damage distribution and surrounding rock stability in underground excavation. In this example, we will show how to create the excavation model, how to get stress equlibrium before excavation and how to use step-by-step excavation in OpenFDEM.

Figure 1: Geometry of the Rock and Excavation#

1 Main Steps#

Initialize the model
Import mesh
Create groups
Assign material properties
Assign boundaries
Setup output
Run model

2 Codes#

Warning

Please note that, all the setting in the example is based on the developer’s computer. You may need to change these details based on your own environment.

To start with programming, create a new empty input file or copy it from the existing examples (later add the .of extension). Begin writing the following commands:

2.1 Initialize the model#

1 Create a new run and clean up the old memories.

of.new

Set up the folder name to save your results. Here the folder name is set to be “result”.

of.set.result "result"

Set the number of cores to be used for running this model in parallelization.

of.set.omp 15

2.2 Create the geometry and the mesh#

See Figure 1, we will create a square rock model with a 0.05 m radius hole at the center of the rock.

2.3 Assign Material Properties#

The material properties of this model is shown as the table below:

Parameter	Value
Continuum Triangular Elements
model	elastic
density (\(kg/m^3\))	2700
E (Pa)	30e9
\(\nu\)	0.3
damp	1.0

Contact Material Properties
model	MC
friction	0.3

Note

Check element materal, and contact materal to see more materials.

Set the material properties as shown below.

of.mat.element 'default' elastic den 2700 E 30e9 v 0.3 damp 1.0
of.mat.contact 'default' MC fric 0.3

2.4 Group the Edges#

Group the four boundaries of the model for applying the in-situ stress and boundary conditions.

of.group.edge 'bottom' range plane.on [0.0 0.0 1.0 0.0]
of.group.edge 'up' range plane.on [0.0 1.0 1.0 1.0]
of.group.edge 'left' range plane.on [0.0 0.0 0.0 1.0]
of.group.edge 'right' range plane.on [1.0 0.0 1.0 1.0]

2.5 Assign Boundary Conditions#

For this example, 30 MPa compressive stress in the xx direction and 5 MPa compressive stresses in the yy direction are assigned to the rock.

of.boundary.element.stress [rock] xx -30e6 xy 0.0 yy -5e6

To quickly balance the model, absorbing boundary conditions are assigned to four boundaries in the model.

of.boundary.edge.viscous 'up'  normal shear
of.boundary.edge.viscous 'bottom'  normal shear
of.boundary.edge.viscous 'right'  normal shear
of.boundary.edge.viscous 'left'  normal shear

2.6 Set the Outputs#

Set the output interval to be every 2000 steps and output all fields variables and fracture variables.

of.history.pv.interval 2000
of.history.pv.field default
of.history.pv.fracture default

In this step, model will first achieve the equilibrium without excavation under in-situ stresses state. The result of velocity or kinematic ratio should be less than 0.00001.

of.step 50000

Insert CZM after achieve the equlibrium of in-situ stresses

Parameter	Value
Cohesive Material Properties
model	EM
tension (Pa)	1e6
cohesion (Pa)	20e6
friction	0.8
GI (\(J/m^2\))	2
GII (\(J/m^2\))	40

of.mesh.insert 'default'
of.mat.cohesive 'default' EM ten 1e6 coh 20e6 fric 0.8 GI 2 GII 40 

One step excavation method. Remove the elements in the “excavation” group.

of.boundary.excavation 'excavation'

Clear the absorbing boundaries

of.boundary.edge.clear 'up'
of.boundary.edge.clear 'bottom'
of.boundary.edge.clear 'right'
of.boundary.edge.clear 'left'

Add pressures at the external boundaries

of.boundary.edge.pressure 'right'  normal 30e6  
of.boundary.edge.pressure 'up'  normal 5e6  
of.boundary.edge.pressure 'bottom'  normal 5e6  
of.boundary.edge.pressure 'left' normal 30e6 

The program will run 500000 steps in total. In other words, it will output 25 files for reference.

of.step 500000

Finalize the model and clear all the temporary memories.

of.finalize

Save the notepad and double click the .of file to run the program.

3 Run the Program#

When you run the program, you can first check the mesh that was created by Gmsh as shown in Figure 2. If the mesh has a good quality, you can close the window to continue run the program.

Mesh of the sample — Figure 2: Mesh created by Gmsh#

Check the boundary conditions here. Absorbing boundary conditions are on for both normal and shear directions. In-situ stresses are applied to the boundaries.

Triangle elements — Figure 3: Check the mesh setting from command window#

4 Results#

To hide hole of excavation, threshold filter can be used.

Stress XX results.

Stress YY results.

In-situ Stress YY — Figure 5: In-situ Stress XX#

5 Full Script#

Apply_Insitu_Stress.of (click to download from Gitlab)

# initialization, this command is to clear the memory in your last run, it is not 
# mandatory but strongly recommend.
of.new

# Set the path folder of your output results, //result// in the same path of your input 
# file will be created by default
of.set.result "result"

# Set the number of cores you want to use, the parallization will be automatically 
# turned on when you use the following command, otherwise the serilization will be used 
# by default
of.set.omp 15


#################################### create mesh #######################################
of.geometry.square 'rock' [0 1 0 1]
of.geometry.cut.circle 'hole' 'rock' [0.5 0.5 0.05 100]

of.geometry.mesh.size 'default' 0.03
of.geometry.mesh auto

of.group.element 'excavation' range circle.in [0.5 0.5 0.05]

################################ assign material parameters ############################
# Assign material for matrix.
of.mat.element 'default' elastic den 2700 E 30e9 v 0.3 damp 1.0
# assign materials for contacts
of.mat.contact 'default' MC fric 0.3

##################################### create groups ####################################
# OpenFDEM can manually group the nodes, elements, cohesive elements and edges by using 
# the region of box, circle and plane.
of.group.edge 'bottom' range plane.on [0.0 0.0 1.0 0.0]
of.group.edge 'up' range plane.on [0.0 1.0 1.0 1.0]
of.group.edge 'left' range plane.on [0.0 0.0 0.0 1.0]
of.group.edge 'right' range plane.on [1.0 0.0 1.0 1.0]

################################### assign boundaries ##################################
# Set the in-situ stresses at the boundaries
of.boundary.element.stress [rock] xx -30e6 xy 0.0 yy -5e6

# add absorbing boundary to quickly balance the model
of.boundary.edge.viscous 'up'  normal shear
of.boundary.edge.viscous 'bottom'  normal shear
of.boundary.edge.viscous 'right'  normal shear
of.boundary.edge.viscous 'left'  normal shear

####################################### set output #####################################
# Set the interval of writing ParaView results.
of.history.pv.interval 2000
# Output all results by default.
of.history.pv.field default
of.history.pv.fracture default

##################################### execute model ####################################

# Run the in-situ stress equilibrium step. Energy equilibrium is achieved when the 
# maximum velocity is less than 1e-5 or the kinematic ratio is less than 1e-5 or the 
# mechanical ratio is less than 1e-5
of.step 50000

# insert CZM after achieve the equlibrium of in-situ stresses
of.mesh.insert 'default'
of.mat.cohesive 'default' EM ten 1e6 coh 20e6 fric 0.8 GI 2 GII 40 
of.boundary.excavation 'excavation'

# Clear the absorbing boundaries
of.boundary.edge.clear 'up'
of.boundary.edge.clear 'bottom'
of.boundary.edge.clear 'right'
of.boundary.edge.clear 'left'

# Add pressures at the external boundaries
of.boundary.edge.pressure 'right'  normal 30e6  
of.boundary.edge.pressure 'up'  normal 5e6  
of.boundary.edge.pressure 'bottom'  normal 5e6  
of.boundary.edge.pressure 'left' normal 30e6 

# Run the excavation step
of.step 500000

# Finalize the model and clear the memory
of.finalize