Data Files for “TunnelBBM” Example

geometry.dat

; Create the rigid block geometry. Initially this will just 
; be a tetrahedral packing that will later be cut. 
model new

; Set the random seed so we get the same results each time
model random 10000

; Specify the domain extent
model domain extent -10 10 

; Create a geometry set to construct rigid blocks to fill the 
; interior of the set.  
geometry set 'blocks' 
[lb = vector(-6.,-1.5,-6.)]
[ub = vector(6.,1.5,6.)]
geometry generate box [comp.x(lb)] [comp.x(ub)] ... 
                      [comp.y(lb)] [comp.y(ub)] ...
                      [comp.z(lb)] [comp.z(ub)]

; Make geometry nodes throughout the box and give them 
; weights (as extra variables) to construct the rblock tetrahedra. 
; This will help to eliminate the crystaline mesh using the INTERNAL keyword. 
[maxEdge = 0.5]
geometry node extra 1 [maxEdge]
fish define addExtra
    spacing = maxEdge
    std = spacing/2.5
    num = (ub - lb) / spacing + vector(1,1,1)
    gs = geom.set.find("blocks")
    loop i(1,int(comp.x(num))+1)
        xloc = comp.x(lb) + (i-1)*spacing
        loop j(1,int(comp.y(num))+1)
            yloc = comp.y(lb) + (j-1)*spacing
            loop k(1,int(comp.z(num))+1)
                zloc = comp.z(lb) + (k-1)*spacing
                gn = geom.node.create(gs,vector(xloc,yloc,zloc))
                val = (math.random.uniform() - 0.5)*2.0*std + spacing
                ;loop while val < spacing / 3. 
                ;    val = (math.random.uniform() - 0.5)*2.0*std + spacing
                ;endloop
                geom.node.extra(gn,1) = val
            endloop
        endloop
    endloop
end
[addExtra]
rblock construct from-geometry 'blocks' maximum-edge [maxEdge] internal

; Create the tunnel as a geometry set. The tunnel is composed 
; of a cylinder and a box
geometry set 'tunnel-box'
geometry generate box -2 2 -1.5 1.5 -2 2
geometry set 'tunnel-arc'
geometry generate cylinder base 1 -1.5 1  radius 1.0 axis 0 1.0 0 ...
                                          height 3.0 resolution 0.4 cap false
geometry generate cylinder base -1 -1.5 1 radius 1.0 axis 0 1.0 0 ... 
                                          height 3.0 resolution 0.4 cap false

; Import the geometry polygons as fractures for cutting
fracture import from-geometry 'tunnel-box' 
fracture import from-geometry 'tunnel-arc'

; Cut the rigid blocks, making sure not to cut those not intersected 
; by the tunnel boundary
rblock cut fractures 'tunnel-arc' throughgoing group 'supportBlock' ... 
                                  intersects-box 1 2 -1.5 1.5 1 2 ...
                                  sliver-remove false sliver-group 'sliver' slot 'sliver'
rblock cut fractures 'tunnel-arc' throughgoing group 'supportBlock' ... 
                                  intersects-box -2 -1 -1.5 1.5 1 2 ...
                                  sliver-remove false sliver-group 'sliver' slot 'sliver'
rblock cut fractures 'tunnel-box' throughgoing group 'supportBlock' ... 
                                  intersects-box -1 1 -1.5 1.5 -2 2 ...
                                  sliver-remove false sliver-group 'sliver' slot 'sliver'
rblock cut fractures 'tunnel-box' throughgoing group 'supportBlock' ... 
                                  intersects-box -2 2 -1.5 1.5 -2 1 ...
                                  sliver-remove false sliver-group 'sliver' slot 'sliver'

; Put the tunnel blocks in a group - shrink the range to exactly the 
; tunnel dimensions
rblock group 'tunnel' range geometry-space 'tunnel-box' ...
                            count odd position-x -1 1
rblock group 'tunnel' range geometry-space 'tunnel-box' ...
                            count odd position-z -2 1
rblock group 'tunnel' range geometry-space 'tunnel-arc' ...
                            count odd

; Scale the rigid blocks so that they overlap initially. This
; is important as the contact position for the bonds is based 
; on the original overlap volume. This contact position is 
; incrementally updated when bonded. 
rblock scale relative 1.01 keep-inertial

; Find the rblocks at the boundary and assign a group. This will 
; be used later for fixing the out-of-plane boundaries
fish define getBoundary
    loop foreach rb rblock.list
        pmin = rblock.facet.vertex.pos.y(rb,1,1)
        pmax = rblock.facet.vertex.pos.y(rb,1,1)
        loop i(1,rblock.facet.num(rb))
            pmin = math.min(rblock.facet.vertex.pos.y(rb,i,1),pmin)
            pmax = math.max(rblock.facet.vertex.pos.y(rb,i,1),pmax)
            pmin = math.min(rblock.facet.vertex.pos.y(rb,i,2),pmin)
            pmax = math.max(rblock.facet.vertex.pos.y(rb,i,2),pmax)
            pmin = math.min(rblock.facet.vertex.pos.y(rb,i,3),pmin)
            pmax = math.max(rblock.facet.vertex.pos.y(rb,i,3),pmax)
        endloop
        if pmin < comp.y(lb)
            rblock.group(rb,'boundary') = 'boundary'
        endif
        if pmax > comp.y(ub)
            rblock.group(rb,'boundary') = 'boundary'
        endif
    endloop
end
@getBoundary

; Remove the geometry set and fractures
geometry delete
fracture delete

; Save the model 
model save 'BBM_Geometry'
Top

stress.dat

; Install the initial stress state in the specimen. 

; Restore the initial geometry created in the previous
; step.
model new
;model restore 'BBM_Geometry_Half'
model restore 'BBM_Geometry'

; Specify the material properties. The SimpleBBM example
; can be used to calibrate a specimen.  
[density = 2650.0]
[young = 50.0e9]
[poisson = 0.25]
[tensileStrength = 2.e7]
[cohesionFactor = 2.4]
[emod = 7.25e10] 

; Elevation used as the top of the model
[topOfModel = -1200.0]

; Set the rigid block damping and density
rblock attribute density [density] damp 0.7

; Use the contact area
rblock contact-resolution update-area true

; Specify the softbond model properties. Initially
; all bonded contacts have infinite strength. This is 
; useful when equilibrating the model. 
contact cmat default model softbond method deformability emod [emod] ...
                           kratio 2.0 ...
                           property sb_ten 1e100 sb_coh 1e100 ...
                           sb_fa 30.0 fric 0.5 sb_mode 1
                                    
; Create zones 
zone create brick size 10 3 10  point 0 -15 -1.5 -15  point 1 15 -1.5 -15  ...
                                point 2 -15  1.5 -15  point 3 -15 -1.5 15
                                
; Densify the zones near to the rigid blocks
zone densify local maximum-length (1,1,1) range position-z -9 9 position-x -9 9 

; Assign the zone material properties 
zone cmodel assign elastic
zone property density [density] young [young] poisson [poisson]

; Turn on gravity
model gravity 0 0 -9.8

; Initialize zone stresses  
zone initialize-stresses overburden [(topOfModel+15.0)*density*9.81] ...
                         ratio 3.0 1.6 direction-x 1.0 0 0 

; Group the zones that will be nulled where the rigid blocks reside
zone group 'nullZones' range position-z -6 6 position-x -6 6 
zone cmodel assign null range group 'nullZones'

; Attach - necessary as the zones have been densified
zone attach by-face 

; Wrap the null zone edges with wall facets 
wall-zone create name 'wrapper' range position-x -6.5 6.5 ...
                                      position-z -6.5 6.5 position-y -1.4 1.4 

; Have the stiffness be computed 
wall-zone compute-stiffness true

; Identify gridpoints for roller boundary conditions 
zone gridpoint group 'z' slot 'boundary-z' range position-z -100 -15 ...
                                                 position-z 15 100 union
zone gridpoint group 'x' slot 'boundary-x' range position-x -100 -15 ...
                                                 position-x 15 100 union
zone gridpoint group 'y' slot 'boundary-y' range position-y -100 -1.5 ...
                                                 position-y 1.5 100 union

; Fix the rigid blocks and zones
rblock fix velocity spin
zone gridpoint fix velocity

; Create the contacts and scale the timestep
model mechanical timestep scale
model cyc 1

; Ensure that the contact areas are unchanged
contact method area

; Inhibit contacts, checking if they should be reactivated every
; 1000 cycles
rblock contact-resolution inhibit-contacts inhibit-interval 1000

; Bond the contacts where rigid blocks overlap
contact method bond range contact active

; Now erode the rigid blocks, removing edge-edge and vertex-vertex.
; This is important as it will greatly increase computational efficiency.
; These contacts do not contribute to the response of the bonded material. 
rblock erode relative skip-errors 0.1

; Specify the stress
rblock tractions overburden [(topOfModel+12.0)*density*9.81] ...
                         ratio 3.0 1.6 direction-x 1.0 0 0 contact-active false

; Free the rigid blocks and zones
rblock free velocity spin
zone gridpoint free velocity

; Fix the boundary rblocks and zones
rblock fix velocity-y range group 'boundary'
zone gridpoint fix velocity-z range group 'z' slot 'boundary-z'
zone gridpoint fix velocity-x range group 'x' slot 'boundary-x'
zone gridpoint fix velocity-y range group 'y' slot 'boundary-y'

; Save the initial state
model save 'BBM_Initial'

; Solve to an average ratio
model solve ratio-average 1e-5

; Get closer to the target stress state
fish define reachTarget
    loop local i(1,2)
        command
            zone cmodel assign elastic range group 'nullZones'
            zone initialize-stresses overburden ...
                                     [(topOfModel+15.0)*density*9.81] ...
                                     ratio 3.0 1.6 direction-x 1.0 0 0 
            zone cmodel assign null range group 'nullZones'
model cycle 1
model solve ratio-average 1e-5
        endcommand
    endloop
end
@reachTarget

; Solve to a lower average ratio
model solve ratio-average 1e-5

; Save the equilibrated model
model save 'BBM_Equil'
Top

excavate.dat

; Cut out the tunnel 
model new

; Restore the equilibrated model
model restore 'BBM_Equil'

; Delete the tunnel rigid blocks
rblock delete range group 'tunnel'

; Delete the slivers
rblock delete range group 'sliver' slot 'sliver'

; Reset the velocities of the blocks and zones
model calm 

[fixMap = rblock.groupmap('supportBlock')]

; Fix all dof of the the support blocks
rblock fix velocity spin range group 'supportBlock' 

; Get the model back into equilibrium
model cycle 1
model solve ratio-average 1e-5

; For each rigid block in the apply group we need to set the 
; applied force and also, if it isn't bonded to anything, need
; to activate all of the contacts
fish define computeSupportForce
    loop foreach local rb fixMap
        if rblock.isbonded(rb) == false
            loop foreach local con rblock.contactmap.all(rb)
                contact.inhibit(con) = false
            endloop
        endif
        rblock.force.app(rb) = vector(0,0,0)
        rblock.force.app(rb) = -1.0 * rblock.force.unbal(rb)
        rblock.extra(rb,1) = rblock.force.app(rb)
    endloop
end
@computeSupportForce

; Free the support blocks and equilibrate
rblock free velocity spin range group 'supportBlock' 
model cycle 1
model solve ratio-average 1e-5

; Reset the displacement of the blocks and zones
rblock attribute displacement 0 0 0
zone gridpoint initialize displacement 0 0 0

; FISH function used to define the tensile and shear strength. It
; returns a non-negative value from a normal distribution with
; given mean and standard deviation. It must take a contact point
; as the first argument.
fish define normalDist(c,mean,sdev)
    loop while true
        local val = mean + (math.random.gauss * sdev)
            if val >= 0.0 then
                normalDist = val
                exit
            endif
    endloop
end

; Assign the normal distribution of tensile and shear strengths
contact property sb_ten @normalDist [tensileStrength] [tensileStrength/5.0] ...
                                    range contact type 'rblock-rblock' ...
                                    contact active
contact property sb_coh @normalDist [tensileStrength * cohesionFactor] ...
                                    [(tensileStrength*cohesionFactor)/5.0] ...
                                    range contact type 'rblock-rblock' ...
                                    contact active

; Create a DFN for crack tracking
[shearNetwork = dfn.create("shear")]
[tensileNetwork = dfn.create("tensile")]
[diskArray = array.create(5)]
[diskArray(1) = "disk"]
[trackCracks = false]

; Define a fish callback to track bond breakages. 
fish define onBondBreak(arr)
    con = arr(1)
    if trackCracks == true
        ; Do not track cracks that happen in the periodic space
        if math.mag(contact.offset(con)) == 0.0
            local mode = arr(2)
            diskArray(2) = contact.pos(con)
            diskArray(3) = math.sqrt(contact.prop(con,'sb_area') / math.pi)
            local norm = contact.normal(con)
            diskArray(4) = math.dip.from.normal(norm)/math.degrad
            diskArray(5) = math.ddir.from.normal(norm)/math.degrad
            if diskArray(5) < 0.0
                diskArray(5) = 360.0+diskArray(5)
            endif
            if mode == 1
                fracture.create(tensileNetwork,diskArray)
            else
                fracture.create(shearNetwork,diskArray)
            endif
        endif
    endif
end

; Add this callback to the bond_break event
fish callback add @onBondBreak event bond_break   

; Reset the velocities of the blocks and zones
model calm

; Solve to equilibrium with the correct tensile and shear 
; strengths installed
model cycle 1
model solve ratio-average 1e-5

; Reset the block and zone displacements for the last time
rblock attribute displacement 0 0 0
zone gridpoint initialize displacement 0 0 0

; Start tracking cracks
[trackCracks = true]

; Save the model
model save 'BBM_Tunnel' 


; FISH halt used to relax the support force. Basically you specify 
; the current percentage, final percentage, increment to reduce the
; percentage by, and the target average ratio. As cycling continues 
; the halt checks for the average ratio smaller than the target value. 
; If not smaller then continue to cycle otherwise reduce the current 
; percentage by the increment. This reduction continues until the final
; percentage is achieved. 
[curPercent = 100]
[endPercent = 30]
[increment = 1]
[targetAveRatio = 5.e-4]
fish define halt
    halt = 0
    currentAveRatio = math.max(zone.mech.ratio.avg,rblock.mech.ratio.avg)
    if curPercent <= endPercent
        if currentAveRatio <= targetAveRatio 
            halt = 1
        endif
        exit
    endif
    if currentAveRatio <= targetAveRatio
        curPercent = math.max(0.0,float(curPercent) - float(increment))
        loop foreach local rb fixMap
            rblock.force.app(rb) = rblock.extra(rb,1) * curPercent / 100.0
        endloop
    endif
end

; Record the current percentage as a history and solve
fish history name 1 @curPercent
history interval 1
model solve fish-halt @halt

; Save the model where the support stress is ramped to 30 percent of the
; original value
model save 'BBM_30_Percent'

; End at 5 percent with an increment of 0.5 percent, solve, and save the model 
; where the support stress is ramped to 5 percent of the original value
[endPercent = 5]
[increment = 0.5]
model solve fish-halt @halt
model save 'BBM_5_Percent'

; End at 0 percent with an increment of 0.01 percent, solve, and save the 
; model where the support stress is ramped to 0 percent of the original value
[endPercent = 0]
[increment = 0.01]
model solve fish-halt @halt
model save 'BBM_0_Percent'

; Register the rigid block contacts and compute the fragments
fragment register rblock-rblock
fragment compute
Top