Tunnel Excavation in a Bonded Block Model

Problem Statement

Note

To view this project in 3DEC, use the menu command Help ► Examples…. Choose “ExampleApplications/BBM_Tunnel” and select “BBM_tunnel.prj” to load. The project’s main data files are shown at the end of this example.

A 10x5 m rectangular tunnel is excavated in fractured crystalline rock at a depth of 1100 m. The area of interest is represented by tetrahedral blocks bonded together, and the rest of the model is elastic continuum. Fragmentation is examined in the area of interest. The example also shows the effects of suppressing the vertex-vertex contacts.

Calibration

A bonded block model (BBM) is composed of many blocks bonded together. The blocks are usually tetrahedra. The properties of the intact rock are usually known from lab studies and these properties are assigned to the zones that make up the blocks. The rock mass properties may be determined from GSI or some other rock classification scheme used to estimate the effect of fractures on the rock strength and stiffness. To create a model of the rock mass, the user needs to specify the joint properties that yield the desired rock mass behavior. This is an inverse problem that usually involves guessing at some properties and then running a test and adjusting the properties until the desired behavior is observed.

Because bonded block models are computationally expensive, usually only a small volume of the entire model is represented this way. The rest of the model is represented by an equivalent continuum. The stiffness of the equivalent continuum should be the same as the macro stiffness of the BBM. The macro stiffness of the BBM is obtained by performing a compression test on a BBM sample. It is important that the block sizes in the calibration test are the same as the blocks that will be used in the field scale model since the joint spacing affects the stiffness.

The elastic properties of the BBM in this model are shown in Table 1. To determine the macro stiffness of the material, a 10x10x20 m sample is created and subjected to an unconfined compression test. The model is created by filling a box with tetrahedral elements (block generate from-geometry). Figure 1 shows the model. The random keyword is used to give some randomness to the block sizes (they will be assigned a uniform size distribution between the minimum and maximum specified edge lengths). This is useful to prevent a very regular packing which may occur when all zones are the same size (Figure 2).

In this example, the formation of vertex-vertex contacts is suppressed using the command block contact vertex-vertex off. This reduces memory and and run time without having a significant effect on the results (this is discussed further below).

Table 1: Model Elastic Properties

Property

Value

Intact Young’s Modulus (GPa)

60

Intact Poisson’s Ratio

0.25

Contact normal stiffness (GPa/m)

1050

Contact shear stiffness (GPa/m)

525

../../../../../_images/bbmtunnel-calibration-blocks.png

Figure 1: 3DEC calibration model showing blocks generated with the random keyword.

../../../../../_images/bbmtunnel-blocks-uniform.png

Figure 2: 3DEC calibration model showing blocks generated without the random keyword.

The UCS test is run by applying velocities to the top and bottom gridpoints and monitoring the vertical stress (sum of reaction forces on the top and bottom gridpoints divided by the cross-sectional area), and the vertical and lateral strains (displacement differences between selected gridpoints on the edges of the model). The data file for this model is shown below (bbm_calibration.dat).

Young’s Modulus is then calculated from vertical stress / vertical strain and Poisson’s Ratio is calculated from horizontal strain / vertical strain. The evolution of Young’s modulus and Poisson’s Ratio in the model are shown in Figure 3 and Figure 4. From these graphs it can be seen that the macro Young’s modulus is approximately 51 GPa and the macro Poisson’s Ratio is approximately 0.23. These values will be used in the equivalent continuum portion of the tunnel model.

../../../../../_images/bbmtunnel-youngs.png

Figure 3: Evolution of Young’s Modulus in the calibration model.

../../../../../_images/bbmtunnel-poisson.png

Figure 4: Evolution of Poisson’s Ratio in the calibration model.

Tunnel Model

For this example, it is assumed that the area of interest is above the tunnel. The model is created by first making a 60x20x40 m block. The region of interest above the proposed tunnel is then cut out with block cut joint-set commands. The cut out block is then deleted and the volume is filled with tetrahedral elements as in the above calibration example. The tunnel is then cut out and the blocks inside of the tunnel are grouped. Finally, the model is zoned with coarser zones in the continuum region. Figure 5 shows the model.

As in the calibration model, the formation of vertex-vertex contacts is suppressed using the command block contact vertex-vertex off. This reduces memory and and run time without having a significant effect on the results.

../../../../../_images/bbmtunnel-model.png

Figure 5: 3DEC tunnel model.

Table 2 summarizes the contact strength parameters. The contacts are assigned the Mohr-Coulomb joint constitutive model. The joint strength properties are assigned a Gaussian distribution. A cross-section showing joint cohesion values is shown in Figure 6

Table 2: Contact Strength Properties

Property

Mean Value

Standard Deviation

Contact dilation (degrees)

5

Contact friction (degrees)

45

5

Contact cohesion (MPa)

2

0.5

Contact tension (MPa)

0.9

0.2

../../../../../_images/bbmtunnel-contacts.png

Figure 6: Cross-section of the model showing joint cohesion values.

The in-situ stress is lithostatic, with the vertical stress equal to the weight of the overburden material and the horizontal stresses equal to twice the vertical.

All boundaries are fixed except the top, on which a vertical stress is applied. Gravity is turned on and the model is solved in large strain.

Tunnel Excavation

Gradual excavation is simulated using the block relax command. After the relaxation, the model is run for another 5000 steps. The model solve command cannot be used here because blocks are falling from the roof and the model never comes to equilibrium.

After 5000 steps, fragments are computed as shown in Figure 7.

../../../../../_images/bbmtunnel-fragments.png

Figure 7: Fragments in the model at the end of the simulation.

Vertex-Vertex Contacts

The same model was run with vertex-vertex contacts on. This results in approximately 800,000 vertex-vertex contacts. The resulting fragments for this model are shown in Figure 8. Some slight differences can be seen, but the overall behavior is similar. This model uses 10% more memory and takes 20% longer to run.

../../../../../_images/bbmtunnel-fragments-vv.png

Figure 8: Fragments in the model at the end of the simulation for a model with vertex-vertex contacts.

Data Files

bbm_calibration.dat

model new
model random 10000

geom set 'box'
geometry generate box -5 5 -5 5 -10 10

; suppress formation of vertex-vertex contacts
block contact vertex-vertex off

; generate tet blocks with randomness
block generate from-geometry set 'box' min-edge 0.6 max-edge 1.4 random

block zone gen edge 1.4

; assign properties
block zone cmodel assign elastic
block zone prop dens 2750 young 60e9 poiss 0.25

block contact jmodel assign elastic
block contact prop stiff-norm 1050e9 stiff-shear 525e9 

; gps on the boundary
block gridpoint group 'top' range pos-z 10
block gridpoint group 'bottom' range pos-z -10

; set up for histories
fish def setup
  global top_gps = list(block.gp.list)(block.gp.isgroup(::block.gp.list,'top'))
  global bottom_gps = list(block.gp.list)(block.gp.isgroup(::block.gp.list,'bottom'))
  
  global area = 10*10  
  global height = 20
  global width = 10 
  global depth = 10
  global gp_left = block.gp.near(-5,0,0)
  global gp_right = block.gp.near(5,0,0)
  global gp_front = block.gp.near(0,-5,0)
  global gp_back = block.gp.near(0,5,0)
  global gp_top = block.gp.near(0,0,10)
  global gp_bottom = block.gp.near(0,0,-10)
end
[setup]

; stress/strain histories
fish def szz
  ; compression positive
  local szz_local = -list.sum(block.gp.force.reaction.z(::top_gps))
  szz_local += list.sum(block.gp.force.reaction.z(::bottom_gps))
  szz_local /= (area*2.0)
  global ezz = (block.gp.disp.z(gp_bottom) - block.gp.disp.z(gp_top))/height
  global young = szz_local / ezz

  local lat_strain1 = (block.gp.disp.x(gp_right) - block.gp.disp.x(gp_left))/width
  local lat_strain2 = (block.gp.disp.y(gp_back) - block.gp.disp.y(gp_front))/depth
  global poiss = 0.5*(lat_strain1+lat_strain2) / ezz
  szz = szz_local
end

fish his szz
fish his ezz
fish his young
fish his poiss

; set up initial velocities to speed analysis
block gridpoint ini vel-z 0 gradient 0 0 -0.005

block gridpoint apply vel-z 0.05 range group 'bottom'
block gridpoint apply vel-z -0.05 range group 'top'

model large-strain off

model cycle 2000

model save 'calibration'

bbm_tunnel.dat

model new
model random 10000

; this prevents unjoining of blocks when new cuts are made
block join-by-contact on

block create brick -30 30 -10 10 -20 20

; cut out bbm region
block cut joint-set dip 90 dip-dir 90 origin -10 0 0 
block cut joint-set dip 90 dip-dir 90 origin 10 0 0 
block cut joint-set origin 0 0 10 range pos-x -10 10
block cut joint-set origin 0 0 -2.5 range pos-x -10 10
block join
block group 'continuum'

; delete bbm block and make bbm
block delete range pos-x -10 10 pos-z -2.5 10

geom set 'box'
geometry generate box -10 10 -10 10 -2.5 10

; suppress formation of vertex-vertex contacts
block contact vertex-vertex off

; generate tet blocks with randomness
block generate from-geometry set 'box' min-edge 0.6 max-edge 1.4 random

; cut out tunnel
block cut joint-set origin 0 0 2.5 jointset-id 99 range pos-x -5 5 
block cut joint-set dip 90 dip-dir 90 jointset-id 99 origin -5 0 0 range pos-z -2.5 2.5
block cut joint-set dip 90 dip-dir 90 jointset-id 99 origin 5 0 0 range pos-z -2.5 2.5
block group 'tunnel' range pos-x -5 5 pos-z -2.5 2.5

; zone
block zone gen edge 2.8 range group 'continuum'
block zone gen edge 1.4

; assign bbm zone properties
block zone cmodel assign elastic
block zone prop dens 2750 young 60e9 poiss 0.25

; assign equivalent continuum zone properties
block zone prop dens 2750 young 51e9 poiss 0.23 range group 'continuum'

; assign contact properties
block contact jmodel assign mohr
block contact prop stiff-norm 1050e9 stiff-shear 525e9 dilation 5
block contact prop-dist friction 45 deviation-gaussian 5
block contact prop-dist cohesion 2e6 deviation-gaussian 0.5e6
block contact prop-dist tension 9e5 deviation-gaussian 0.2e5 

block contact mat-table default prop stiff-norm 1050e9 stiff-shear 525e9 friction 45 dilation 5

; make construction joints elastic
block contact jmodel assign elastic range joint-set 99
block contact prop stiff-norm 1050e9 stiff-shear 525e9 range joint-set 99

; initial conditions
; assume 1100 m below the surface and horizontal = 2*vertical
model gravity 9.81

[overburden = -9.81*2750*1100]
block insitu topo ratio-x 2 ratio-y 2 overburden [overburden]

; boundary conditions
block face apply stress-zz [overburden] range pos-z 20
block gridpoint apply velocity-y 0. range position-y -10
block gridpoint apply velocity-y 0. range position-y 10
block gridpoint apply velocity-x 0. range position-x -30
block gridpoint apply velocity-x 0. range position-x 30
block gridpoint apply velocity 0 0 0 range position-z -20

; histories
block his name 'top-middle' disp-z pos 0 0 2.5

model large-strain on

model solve

model save 'tunnel-initial'