Example: Tunnel Through a Fractured Rock (Voronoi)

Problem Statement

Note

The project file for this example is available to be viewed/run in in FLAC2D.[1]. The project’s main data files are shown at the end of this example.

This example demonstrates the methodology for generating a voronoi fracture network in a FLAC2D model and its effect on the stability of a tunnel with an arched roof. The tunnel is assumed to be 232 m below the surface. Liner support is then added to try to stabilize the tunnel in the fractured rock mass. The rock blocks (between the joints) are assumed to be elastic. Rock properties and joint properties are shown in Table 1 and the properties of the shotcrete liner are shown in Table 2. Note that the liner is modeled as elastic beams. The strengths in Table 2 are used to generate moment-thrust diagrams at the end of the example. The liner is assumed to be perfectly bonded to the rock.

Table 1: Rock and Joint Properties

Property	Value
Rock Density	2700 kg/m3
Rock Young’s Modulus	30 GPa
Rock Poisson’s Ratio	0.22
Joint Normal Stiffness	22.5 GPa/m
Joint Shear Stiffness	1.5 GPa/m
Joint Friction	30 degrees
Joint Residual Friction	15 degrees
Joint Cohesion	0.8 MPa
Joint Tension	0.1 MPa

Table 2: Shotcrete Liner Properties

Property	Value
Density	2400 kg/m3
Young’s Modulus	20 GPa
Poisson’s Ratio	0.15
Thickness	0.1 m
Tensile Yield Strength	50 MPa
Compressive Yield Strength	50 MPa

Model Building

The model is built using i Sketch as described below.

• Create a new Sketch Set. Select the Tunnel Wizard and create a tunnel as shown in Figure 1. Note that the edges of the tunnel are automatically assigned a group name (Tunnel Surface). This will be useful later when we add a liner.

• Draw a box around the tunnel going from -6,-6 to 6,6. Thiswill define the fractured region.

• Draw a bigger box from -25,-25 to 25,25. This will define the outer boundary. The model should look like this:

• Select the Voronoi Block tool. Under Domain, select Closed Polygon, slick the Select button and click inside of the inner box. Fill out the rest of the dialog as shown:

• You can click Preview to see the fractures before generating. To actually create the nodes and edges in Sketch, click OK. The model should appear as shown below. Note that the Voronoi block edges are automatically assigned a group name “Voronoi-1” in slot “voronoi”. This will be useful later when we create joints.

• In the zoning menu, select Auto Size and set Zone Length to 0.25. Set the minimum number of zones per edge to 2. This is necessary to properly separate the faces when we turn the Voronoi edges into joints.

• We don’t need this much detail away from the tunnel, so select the four outer edges and set the zone length for these to 1.

• Mesh all Polygons

• Click inside the tunnel and assign this blocks the group name “Tunnel”.

• Create the zones and save the state.

Figure 6: Zoned model.

Initial State

Create a data file and add a command to restore the model created in Sketch. Next you will want to skin the model using the command zone face skin. This automatically assigns group names to the faces on the outside of the model. This could also be done by clicking a button in the Model Pane. It is important that this step be performed before separating all the faces to make joints, otherwise all of the new faces will also be automatically be assigned group names.

To create the joints, use the commands below:

zone joint configure
zone joint create by-slot 'voronoi' separate distinct-geometry

This automatically separates all of the faces with a group name assigned in slot “voronoi”. A special keyword distinct-geometry is required because the voronois have a non-manifold topology (think of a T-junction). This is not required for separation when each joint is distinct (i.e. not intersecting).

Note

The macro stiffness of the fractured part of the model will be less than the stiffness of the rock, due to the presence of the joints. In a real project you should consider reducing the stiffness of the unfractured material so that it matches the macro stiffness of the fractured part of the model.

After specifying joint and rock properties, initial conditions and boundary conditions, the model is solved. Displacements are small as expected (< 1 mm).

Excavation

Excavation is performed with the zone relax excavate command. The resulting displacements are shown in Figure 7. The maximum displacement is just over 3.5 mm in the bottom of the tunnel. Joint shear displacements are shown in Figure 8.

Figure 7: Displacements in the model after tunnel excacavation.

Figure 8: Joint shear displacements in the model after tunnel excacavation.

Liner

The model is now rerun and a liner is installed after relaxing the excavation by 80%. This follows the usual method for simulating a tunnel excavation with support in 2D, where the supporting effect of the face cannot be considered directly. See Multistage Tunnel Excavation and Support (FLAC2D) for details on how to construct a ground reaction curve and obtain the desired percentage of relaxation before liner installation.

In this example beam elements are used instead of liners since we are assuming that the shotcrete is perfectly bonded to the rock. Note that beams are assumed to have a short length in the out-of-plane direction (unlike liners, which are assumed to be continuous), therefore to make beams behave in a plane-strain manner, you need to divide the Young’s modulus by \(1-\nu^2\). For the moment of inertia, we assume a rectangular cross section, so the MOI is \(ab^3/12\) where a is the out-of-plane dimension (1) and b is the thickness (0.1).

In this example, a liner isn’t really necessary because the rock is stable without one, but it is possible that for a different random configuration of voronois, there could conceivably be wedges forming and blocks falling. For this reason, the liner is only added to the sides and the top of the tunnel.

The resulting displacements are shown in Figure 9 and the liner force magnitudes are shown in Figure 10. A moment-thrust diagram for the liner elements is shown in Figure 11. As expected, the liner elements are well within the a factor of safety of 1.2.

Figure 9: Displacements in the model after tunnel excacavation with an installed liner.

Figure 10: Force magnitudes in the liner.

Figure 11: Force magnitudes in the liner.

Data Files

initial.dat

model new
program call 'sketch'

; skin model to get face groups BEFORE separating faces
; this could be done in the model pane
zone face skin

; create zone joints and set properties
zone joint configure
zone joint create by-slot 'voronoi' separate distinct-geometry

contact property stiffness-normal 22.5e9 stiffness-shear 1.5e9 ...
  friction 30 friction-residual 15 ...
  cohesion 0.8e6 tension 0.1e6
  
; rock model and properties
zone cmodel assign elastic
zone property density 2700 young 30e9 poisson 0.22 

; boundary conditions
zone face apply velocity 0 0 range group 'East' or 'West'
zone face apply velocity 0 0 range group 'Top' or 'Bottom'

; Initial conditions.  Top of tunnel is at a depth of 232.65 m
;    Top of model is 23 m above the top of the tunnel
; Assume horizontal stresses are equal to vertical
model gravity 9.81
zone initialize-stresses ratio 1 overburden [-(232.65-23)*9.81*2700]

; solve it
model large-strain off
model solve convergence 1

model save 'initial'

excavate.dat

model restore 'initial'

zone gridpoint initialize displacement 0,0

zone history displacement position 0,2

zone relax excavate range group 'Tunnel'

model solve convergence 1

model save 'excavate'

liner.dat

model restore 'initial'

zone gridpoint initialize displacement 0,0

zone history displacement position 0,2

; relax to 0.2 of initial stress
zone relax excavate minimum 0.2 range group 'Tunnel'

model solve convergence 1

; delete tunnel zones and add liner
zone delete range group 'Tunnel'

; only add to the sides and top
structure beam create by-zone-face range group 'Tunnel Surface' position-y -2.49 2

; use beams so we don't have to specify coupling properties
; just be careful to divide modulus by 1-poiss
structure beam property density 2.4E3 young [2E10/(1-0.15^2)] poisson 0.15 ...
  cross-sectional-area 0.1 moi [1*0.1^3/12] 

model solve convergence 1

model save 'liner'

Endnote

[1]

This example is listed under the present document title in the i Examples dialog (use the menu command Help -> Examples… to access).