Tunnel Excavation in a Bonded Block Model
Problem Statement
Note
To view this project in 3DEC, use the menu command . 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 vertexvertex 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 fromgeometry
). 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 vertexvertex contacts is suppressed using the command block contact vertexvertex
off. This reduces memory and and run time without having a significant effect on the results (this is discussed further below).
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 
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 crosssectional 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.
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 jointset
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 vertexvertex contacts is suppressed using the command block contact vertexvertex
off. This reduces memory and and run time without having a significant effect on the results.
Table 2 summarizes the contact strength parameters. The contacts are assigned the MohrCoulomb joint constitutive model. The joint strength properties are assigned a Gaussian distribution. A crosssection showing joint cohesion values is shown in Figure 6
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 
The insitu 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.
VertexVertex Contacts
The same model was run with vertexvertex contacts on. This results in approximately 800,000 vertexvertex 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.
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 vertexvertex contacts
block contact vertexvertex off
; generate tet blocks with randomness
block generate fromgeometry set 'box' minedge 0.6 maxedge 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 stiffnorm 1050e9 stiffshear 525e9
; gps on the boundary
block gridpoint group 'top' range posz 10
block gridpoint group 'bottom' range posz 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 velz 0 gradient 0 0 0.005
block gridpoint apply velz 0.05 range group 'bottom'
block gridpoint apply velz 0.05 range group 'top'
model largestrain 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 joinbycontact on
block create brick 30 30 10 10 20 20
; cut out bbm region
block cut jointset dip 90 dipdir 90 origin 10 0 0
block cut jointset dip 90 dipdir 90 origin 10 0 0
block cut jointset origin 0 0 10 range posx 10 10
block cut jointset origin 0 0 2.5 range posx 10 10
block join
block group 'continuum'
; delete bbm block and make bbm
block delete range posx 10 10 posz 2.5 10
geom set 'box'
geometry generate box 10 10 10 10 2.5 10
; suppress formation of vertexvertex contacts
block contact vertexvertex off
; generate tet blocks with randomness
block generate fromgeometry set 'box' minedge 0.6 maxedge 1.4 random
; cut out tunnel
block cut jointset origin 0 0 2.5 jointsetid 99 range posx 5 5
block cut jointset dip 90 dipdir 90 jointsetid 99 origin 5 0 0 range posz 2.5 2.5
block cut jointset dip 90 dipdir 90 jointsetid 99 origin 5 0 0 range posz 2.5 2.5
block group 'tunnel' range posx 5 5 posz 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 stiffnorm 1050e9 stiffshear 525e9 dilation 5
block contact propdist friction 45 deviationgaussian 5
block contact propdist cohesion 2e6 deviationgaussian 0.5e6
block contact propdist tension 9e5 deviationgaussian 0.2e5
block contact mattable default prop stiffnorm 1050e9 stiffshear 525e9 friction 45 dilation 5
; make construction joints elastic
block contact jmodel assign elastic range jointset 99
block contact prop stiffnorm 1050e9 stiffshear 525e9 range jointset 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 ratiox 2 ratioy 2 overburden [overburden]
; boundary conditions
block face apply stresszz [overburden] range posz 20
block gridpoint apply velocityy 0. range positiony 10
block gridpoint apply velocityy 0. range positiony 10
block gridpoint apply velocityx 0. range positionx 30
block gridpoint apply velocityx 0. range positionx 30
block gridpoint apply velocity 0 0 0 range positionz 20
; histories
block his name 'topmiddle' dispz pos 0 0 2.5
model largestrain on
model solve
model save 'tunnelinitial'
Was this helpful? ...  Itasca Software © 2023, Itasca  Updated: Feb 16, 2024 