Data Files for “Rock Testing” Example

make_sample.p2.dat

; fname: make_sample.p2dat
;
;  Create unbonded sample
;==============================================================================
model new
model title 'Testing Bonded Particle Model'

; Set the domain extent
model domain extent -0.05 0.05 -0.1 0.1 condition destroy

contact cmat default model linear method deform emod 1.0e9 kratio 0.0 
contact cmat default property dp_nratio 0.5 

; create walls that extend past the edges of the sample
wall create vertices -0.03,0.05 0.03,0.05 id 1
wall create vertices -0.03,-0.05 0.03,-0.05 id 2
wall create vertices -0.025,-0.06 -0.025,0.06 id 3
wall create vertices 0.025,-0.06 0.025,0.06 id 4

model random 10001
ball distribute porosity 0.1 radius 0.5e-3 0.75e-3 box -0.025 0.025 -0.05 0.05
ball attribute density 2500 damp 0.7

; Calm the system
model cycle 1000 calm 10
; Solve the system to a target limit (here the average force ratio)
; Use density scaling to quickly reach equilibrium
model mechanical timestep scale
model solve ratio-average 1e-4
model mechanical timestep auto
model calm

; delete side walls
wall delete range id 3
wall delete range id 4

model save 'unbonded'
;==============================================================================
; eof: make_sample.p2dat
Top

parallel_bonded.p2dat

; fname: parallel_bonded.p2dat
;
;  Create parallel bonded sample
;==============================================================================
model restore 'unbonded' 

contact model linearpbond range contact type 'ball-ball'
contact method bond gap 0.5e-4

; set linear stiffness
contact method deform emod 1.0e9 krat 1.0

; set stiffness of bond material
contact method pb_deform emod 1.0e9 krat 1.0

; set bond strengths
contact property pb_ten 10.0e6 pb_coh 50.0e6 pb_fa 0.0

; set some damping at the contacts
contact property dp_nratio 0.5

; set ball-ball friction to non-zero value
contact property fric 0.577 range contact type 'ball-ball' 

; Reset ball displacement 
ball attribute displacement multiply 0.0

; Set the linear force to 0.0 and force a reset of the linear contact forces. 
contact property lin_force 0.0 0.0 lin_mode 1
ball attribute force-contact multiply 0.0 moment-contact multiply 0.0 
model cycle 1
model solve ratio-average 1e-5
model save 'parallel_bonded' 

;==============================================================================
; eof: parallel_bonded.p2dat
Top

ucs.p2dat

; fname: ucs.p2dat
;
;  Perform an unconfined compressive strength test on a bonded sample
;==============================================================================
program echo off
program call 'ss_wall.fis' 
program call 'fracture.p2fis' 
program echo on

; set up global parameters for stress and strain measurement
@setup_wall

; apply loading by moving top and bottom walls
wall attribute velocity-y -0.1 range id 1
wall attribute velocity-y  0.1  range id 2

; apply a small amount of damping
ball attribute damp 0.1

; record histories
fish history name 1 @axial_stress_wall
fish history name 2 @axial_strain_wall

; record cracks
@track_init
fish history name 3 @crack_num

; run some steps to get past initial vibrations
model cycle 1000

; run the test until stress falls below 70% of the peak
fish set peak_fraction = 0.7
model solve fish-halt @loadhalt_wall
fish list @peak_stress

;==============================================================================
; eof: ucs.p2dat
Top

tension.p2dat

; fname: tension.p2dat
;
;  Perform a direct tension test on a bonded sample
;==============================================================================

program echo off
program call 'ss_gage.fis'
program call 'fracture.p2fis'
program echo on

; set up global parameters for strain measurement
@setup_gage

; delete the walls
wall delete

; identify balls to form top and bottom grips
ball group 'top_grip' range position-y 0.049 0.05
ball group 'bottom_grip' range position-y -0.049 -0.05

; apply loading by moving top and bottom platens
ball fix velocity-y range group 'top_grip'
ball attribute velocity-y 0.01 range group 'top_grip'
ball fix velocity-y range group 'bottom_grip'
ball attribute velocity-y -0.01 range group 'bottom_grip'

; apply a small amount of damping
ball attribute damp 0.1

; add a measurement circle to record stress
; location defaults to 0,0
measure create radius 0.02

; record histories
measure history name 1 stress-yy id 1
fish history name 2 @axial_strain_gage

; record cracks
@track_init
fish history name 3 @crack_num

; run the test until stress falls below 70% of the peak
fish set peak_fraction = 0.7
model cycle 1000
model solve fish-halt @loadhalt_meas
fish list @peak_stress

;==============================================================================
; eof: tension.p2dat
Top

flat_joint.p2dat

; fname: flat_joint.p2dat
;
;  Create sample with flatjoint contacts
;==============================================================================
model restore 'unbonded'

; here we want all existing contacts to assume flat_joint properties
contact model flatjoint range contact type 'ball-ball'
contact property fj_track on

; bond ball-ball contacts
contact method bond gap 0.5e-4 range contact type 'ball-ball'

; set linear stiffness
contact method deform emod 1.0e9 krat 1.0

; set bond strengths
contact property fj_ten 10.0e6 fj_coh 50.0e6 fj_fa 0.0

; set ball-ball friction to non-zero value
contact property fj_fric 0.577 range contact type 'ball-ball'

; Set new contacts to be flatjoint contacts (unbonded)
contact cmat default model flatjoint ...
    property fj_n 4 fj_track on ...
    method deformability emod 1.0e9 krat 1.0

; Reset ball displacement 
ball attribute displacement multiply 0.0

; Set the linear force to 0.0 and force a reset of the linear contact forces. 
contact property lin_force 0.0 0.0 lin_mode 1
ball attribute force-contact multiply 0.0 moment-contact multiply 0.0 
model cycle 1
model solve ratio-average 1e-5

model save 'flatjoint'

;==============================================================================
; eof: flat_joint.p2dat
Top

fracture.p2fis

; fname: fracture.p2fis
;
; Simple environment to track fragmentation in a BPM.
; Track LinearPBond model "bond_change" events and turn them into fractures.
; Use Fragment logic and Ball Result logic to record fragemnt ids 
;
;=============================================================================

fish define add_crack(entries)
    local contact    = entries(1)
    local mode       = entries(2)
    local frac_pos   = contact.pos(contact)
    local norm       = contact.normal(contact)
    local dfn_label  = 'crack'
    local frac_size
    local bp1 = contact.end1(contact)
    local bp2 = contact.end2(contact)
    local ret = math.min(ball.radius(bp1),ball.radius(bp2))
          ;contact.method(contact,'pb_radius')
    frac_size = ret
    local inDir = vector(-comp.y(norm),comp.x(norm))
    local vert1 = frac_pos + inDir * frac_size
    local vert2 = frac_pos - inDir * frac_size
    local arg = array.create(4)
    arg(1) = 'vertices'
    arg(2) = 2
    arg(3) = vert1
    arg(4) = vert2
    crack_num = crack_num + 1
    
    if mode = 1 then 
        ; failed in tension
        dfn_label = dfn_label + '_tension'
    else if mode = 2 then
        ; failed in shear
        dfn_label = dfn_label + '_shear'
    endif
        
    global dfn = dfn.find(dfn_label)
    if dfn = null then
        dfn = dfn.create(dfn_label)
    endif
    local fnew = fracture.create(dfn,arg)
    fracture.prop(fnew,'age')  = mech.time.total
    fracture.extra(fnew,1) = bp1
    fracture.extra(fnew,2) = bp2
    crack_accum += 1
    if crack_accum > 50
        if frag_time < mech.time.total
            frag_time = mech.time.total
            crack_accum = 0
            command 
                fragment compute
            endcommand
            ; go through and update the fracture positions
            loop for (local i = 0, i < 2, i = i + 1)
                local name = 'crack_tension'
                if i = 1
                    name = 'crack_shear'
                endif
                dfn = dfn.find(name)
                if dfn # null
                    loop foreach local frac dfn.fracturelist(dfn)
                        local ball1 = fracture.extra(frac,1)
                        local ball2 = fracture.extra(frac,2)
                        if ball1 # null
                            if ball2 # null
                                local len=fracture.len(frac)/2.0
                                local pos=(ball.pos(ball1)+ball.pos(ball2))/2.0
                                if comp.x(pos)-len > xmin
                                    if comp.x(pos)+len < xmax
                                        if comp.y(pos)-len > ymin
                                            if comp.y(pos)+len < ymax
                                                fracture.pos(frac) = pos
                                            endif
                                        endif
                                    endif
                                endif
                            endif
                        endif
                    endloop
                endif
            endloop
        endif
    endif
end

fish define track_init
    command
        fracture delete
        fragment clear
        fragment register ball-ball
    endcommand
  ; activate fishcalls
    command
        fish callback remove @add_crack
        fish callback add @add_crack event bond_break
    endcommand
    ; reset global variables
    global crack_accum = 0
    global crack_num = 0
    global track_time0 = mech.time.total
    global frag_time = mech.time.total
    global xmin = domain.min.x()
    global ymin = domain.min.y()
    global xmax = domain.max.x()
    global ymax = domain.max.y()

end

;=============================================================================
; eof: fracture.p3fis
Top

make_sample.p3dat

; fname: make_sample.p3dat
;
;  Create unbonded sample
;=============================================================================
model new
model title 'Testing Bonded Particle Model'

; Set the domain extent
model domain extent -0.05 0.05 -0.05 0.05 -0.1 0.1 condition destroy

contact cmat default model linear method deformability emod 1.0e9 kratio 0.0
contact cmat default property dp_nratio 0.5

; create walls extending past the edges of the sample
wall generate id 1 plane dip 0 dip-direction 0   position 0 0 0.04
wall generate id 2 plane dip 0 dip-direction 0   position 0 0 -0.04
wall generate id 3 plane dip 90 dip-direction 90 position -0.025 0 0
wall generate id 4 plane dip 90 dip-direction 90 position 0.025 0 0
wall generate id 5 plane dip 90 dip-direction 0  position 0 -0.025 0
wall generate id 6 plane dip 90 dip-direction 0  position 0 0.025 0

model random 10002
ball distribute porosity 0.2 radius 1.0e-3 1.5e-3 ...
                box -0.025 0.025 -0.025 0.025 -0.04 0.04
ball attribute density 2500 damp 0.7

; Calm the system
model cycle 1000 calm 10
; Solve the system to a target limit (here the average force ratio)
; Use density scaling to quickly reach equilibrium
model mechanical timestep scale
model solve ratio-average 1e-4
model mechanical timestep auto
model calm

; Delete side walls
; Be careful to include the keyword 'walls' or else facets will be deleted
wall delete walls range id 3 6

; make cylinder
ball delete range cylinder end-1 0 0 -0.04 end-2 0 0 0.04 rad 0.02 not

model save 'unbonded'
;=============================================================================
; eof: make_sample.p3dat
Top

parallel_bonded.p3dat

; fname: parallel_bonded.p3dat
;
;  Create parallel bonded sample
;==============================================================================
model restore 'unbonded'

contact model linearpbond range contact type 'ball-ball'
contact method bond gap 2.0e-4

; set linear stiffness using methods
contact method deform emod 1.0e9 krat 1.0

; set stiffness of bonds using methods
contact method pb_deform emod 1.0e9 krat 1.0

; set bond strengths
contact property pb_ten 10.0e6 pb_coh 50.0e6 pb_fa 0.0

; set some damping at the contacts
contact property dp_nratio 0.5

; set ball-ball friction to non-zero value
contact property fric 0.577 range contact type 'ball-ball'

; Reset ball displacement 
ball attribute displacement multiply 0.0

; Set the linear force to 0.0 and force a reset of the linear contact forces. 
contact property lin_force 0.0 0.0 0.0 lin_mode 1
ball attribute force-contact multiply 0.0 moment-contact multiply 0.0 
model cycle 1
model solve ratio-average 1e-5

model save 'parallel_bonded'

;==============================================================================
; eof: parallel_bonded.p3dat
Top

ucs.p3dat

; fname: ucs.p3dat
;
;  Perform an unconfined compressive strength test on a bonded sample
;==============================================================================
echo off
call 'ss_wall.fis'
call 'fracture.p3fis'
echo on

; set up global parameters for stress and strain measurement
@setup_wall

; apply loading by moving top and bottom walls
wall attribute velocity-z -0.15 range id 1
wall attribute velocity-z  0.15 range id 2

; apply a small amount of damping
ball attribute damp 0.1

; record histories
fish history name 1 @axial_stress_wall
fish history name 2 @axial_strain_wall

; record cracks
@track_init
fish history name 3 @crack_num

; cycle a few steps to get past initial vibations
model cyc 1000

; run the test until stress falls below 70% of the peak
[peak_fraction = 0.7]
model solve fish-halt @loadhalt_wall
list @peak_stress

;==============================================================================
; eof: ucs.p3dat
Top

tension.p3dat

; fname: tension.p3dat
;
;  Perform a direct tension test on a bonded sample
;==============================================================================

echo off
call 'ss_gage.fis'
call 'fracture.p3fis'
echo on

; set up global parameters for strain measurement
@setup_gage

; delete the walls
wall delete

; identify balls to form top and bottom grips
ball group 'top_grip' range position-z 0.035 0.04
ball group 'bottom_grip' range position-z -0.035 -0.04

; apply loading by moving top and bottom platens
ball fix velocity-z range group 'top_grip'
ball attribute velocity-z 0.03 range group 'top_grip'
ball fix velocity-z range group 'bottom_grip'
ball attribute velocity-z -0.03 range group 'bottom_grip'

; apply a small amount of damping
ball attribute damp 0.1

; add a measurement circle to record stress
; location defaults to 0,0
measure create radius 0.02

; record histories
measure history name 1 stress-zz id 1
fish history name 2 @axial_strain_gage

; record cracks
@track_init
fish history name 3 @crack_num

; run the test until stress falls below 70% of the peak
[peak_fraction = 0.7]
model solve fish-halt @loadhalt_meas
list @peak_stress

;==============================================================================
; eof: tension.p3dat
Top

flat_joint.p3dat

; fname: flat_joint.p3dat
;
;  Create sample with flatjoint contacts
;==============================================================================
model restore 'unbonded'

; here we want all existing contacts to assume flat_joint properties
contact model flatjoint range contact type 'ball-ball'
contact property fj_track on

; bond ball-ball contacts
contact method bond gap 2.0e-4 range contact type 'ball-ball'

; set stiffness of bonds using methods
contact method deform emod 1.0e9 krat 1.0

; set bond strengths
contact property fj_ten 10.0e6 fj_coh 50.0e6 fj_fa 0.0

; set ball-ball friction to non-zero value
contact property fj_fric 0.577 range contact type 'ball-ball'

; Set new ball-ball contacts to be flatjoint contacts (unbonded)
contact cmat default model flatjoint ...
    property fj_track on ...
    method deform emod 1.0e9 krat 1.0 ...
    type 'ball-ball'

; Reset ball displacement 
ball attribute displacement multiply 0.0
; Set the linear force to 0.0 and force a reset of the linear contact forces. 
contact property lin_force 0.0 0.0 0.0 lin_mode 1
ball attribute force-contact multiply 0.0 moment-contact multiply 0.0 
model cycle 1
model solve ratio-average 1e-5

model save 'flatjoint'

;==============================================================================
; eof: flat_joint.p3dat
Top

doall_parallel.p3dat

; fname: doall_parallel.p3dat
;
;  Create and test parallel bonded material
;==============================================================================

; restore unbonded sample and bond it
call 'parallel_bonded.p3dat'

log-file 'parallel3d.log'
log on

; run ucs test
call 'ucs.p3dat'
model save 'parallel_ucs'

; run tension test
model restore 'parallel_bonded'
call 'tension.p3dat'
model save 'parallel_tension'

log off

;==============================================================================
; eof: doall_parallel.p3dat
Top

doall_fj.p3dat

; fname: doall_fj.p3dat
;
;  Create and test flat joint material
;==============================================================================

; restore unbonded specimen and install flatjoint model
call 'flat_joint.p3dat'

log-file 'fj3d.log'
log on

; run ucs test
call 'ucs.p3dat'
model save 'fj_ucs'

; run tension test
model restore 'flatjoint'
call 'tension.p3dat'
model save 'fj_tension'

log off

;==============================================================================
; eof: doall_fj.p3dat
Top

fracture.p3fis

; fname: fracture.p3fis
;
; Simple environment to track fragmentation in a BPM.
; Track LinearPBond model "bond_change" events and turn them into fractures.
; Use Fragment logic and Ball Result logic to record fragemnt ids 
;
;=============================================================================

fish define add_crack(entries)
    local contact    = entries(1)
    local mode       = entries(2)
    local  frac_pos   = contact.pos(contact)
    local norm       = contact.normal(contact)
    local dfn_label  = 'crack'
    local  frac_size
    local  bp1 = contact.end1(contact)
    local  bp2 = contact.end2(contact)
    local  ret = math.min(ball.radius(bp1),ball.radius(bp2))
    frac_size = ret
    local arg = array.create(5)
    arg(1) = 'disk'
    arg(2) = frac_pos
    arg(3) = frac_size
    arg(4) = math.dip.from.normal(norm)/math.degrad
    arg(5) = math.ddir.from.normal(norm)/math.degrad
    if arg(5) < 0.0
      arg(5) = 360.0+arg(5)
    end_if
    crack_num = crack_num + 1
    
    if mode = 1 then 
        ; failed in tension
        dfn_label = dfn_label + '_tension'
    else if mode = 2 then
        ; failed in shear
        dfn_label = dfn_label + '_shear'
    endif
        
    global dfn = dfn.find(dfn_label)
    if dfn = null then
        dfn = dfn.create(dfn_label)
    endif
    local fnew = fracture.create(dfn,arg)
    fracture.prop(fnew,'age')  = mech.time.total
    fracture.extra(fnew,1) = bp1
    fracture.extra(fnew,2) = bp2
    crack_accum += 1
    if crack_accum > 50
        if frag_time < mech.time.total
            frag_time = mech.time.total
            crack_accum = 0
            command 
                fragment compute
            endcommand
            ; go through and update the fracture positions
            loop for (local i = 0, i < 2, i = i + 1)
                local name = 'crack_tension'
                if i = 1
                    name = 'crack_shear'
                endif
                dfn = dfn.find(name)
                if dfn # null
                    loop foreach local frac dfn.fracturelist(dfn)
                        local ball1 = fracture.extra(frac,1)
                        local ball2 = fracture.extra(frac,2)
                        if ball1 # null
                            if ball2 # null
                                local len=fracture.diameter(frac)/2.0
                                local pos=(ball.pos(ball1)+ball.pos(ball2))/2.0
                                if comp.x(pos)-len > xmin
                                    if comp.x(pos)+len < xmax
                                        if comp.y(pos)-len > ymin
                                            if comp.y(pos)+len < ymax
                                                if comp.z(pos)-len > zmin
                                                    if comp.z(pos)+len < zmax
                                                        fracture.pos(frac)= pos
                                                    end_if
                                                end_if
                                            endif
                                        endif
                                    endif
                                endif
                            endif
                        endif
                    endloop
                endif
            endloop
        endif
    endif
end

fish define track_init
    command
        fracture delete
        model results clear-map
        fragment clear
        fragment register ball-ball
    endcommand
  ; activate fishcalls
    command
        fish callback remove @add_crack event bond_break
        fish callback add @add_crack event bond_break
    endcommand
    ; reset global variables
    global crack_accum = 0
    global crack_num = 0
    global track_time0 = mech.time.total
    global frag_time = mech.time.total
    global xmin = domain.min.x()
    global ymin = domain.min.y()
    global xmax = domain.max.x()
    global ymax = domain.max.y()
    global zmin = domain.min.z()
    global zmax = domain.max.z()

end

;=============================================================================
; eof: fracture.p3fis
Top

ss_wall.fis

; fname: ss_wall.fis
;
;  Fish functions to record stress and strain in wall-loaded core samples
;========================================================================
fish define setup_wall
;
; Find the walls and get initial sample dimensions
;
  global wp_top = wall.find(1) ; assume wall 1 is the top wall
  global wp_bottom = wall.find(2) ; assume wall 2 is the bottom wall
  global vertical_direction = global.dim
  global sample_height = wall.pos(wp_top,vertical_direction) ...
                         - wall.pos(wp_bottom,vertical_direction)
  
  ; assume x and y are approximately the same in 3D
  local xmin = 1.0e12
  local xmax = -1.0e12
  loop foreach bp ball.list
    local ball_xmin = ball.pos.x(bp) - ball.radius(bp)
    xmin = math.min(xmin,ball_xmin)
    local ball_xmax = ball.pos.x(bp) + ball.radius(bp)
    xmax = math.max(xmax,ball_xmax)
  end_loop
  
  local diameter_ = xmax - xmin
  if global.dim = 2
    global cross_sectional_area = diameter_
  else
    ; assume cylindrical sample in 3D
    cross_sectional_area = math.pi*0.25*diameter_*diameter_
  end_if
  
end
========================================================
fish define axial_stress_wall
;
; Compute axial stress (positive tension) using walls
;
; Assumes global variables sample_height and sample_width have been set
;
  
  local force1 = -wall.force.contact(wp_top,vertical_direction)
  local force2 = wall.force.contact(wp_bottom,vertical_direction)
  axial_stress_wall = 0.5*(force1+force2)/cross_sectional_area
end
========================================================
fish define axial_strain_wall
;
; Compute axial strain (positive tension) using walls
;
; Assumes global variable sample_width has been set
;
  axial_strain_wall = 2.0*wall.disp(wp_top,vertical_direction)/sample_height
end
==========================================================
fish define loadhalt_wall
;
; Function used to stop a test when axial stress decreases
; to some fraction of the peak
; Assumes axial_stress_wall is a function that returns axial stress
; 
; INPUT: peak_fraction - fraction of peak stress that dictates 
;                        the stopping of the test (global float)
;
  loadhalt_wall = 0
  local abs_stress = math.abs(axial_stress_wall)
  global peak_stress = math.max(abs_stress,peak_stress)
  if abs_stress < peak_stress*peak_fraction
    loadhalt_wall = 1
  end_if
end
;=============================================================================
; eof: ss_wall.fis
Top

ss_gage.fis

; fname: ss_gage.fis
;
;  Fish functions to record strain using gage balls
========================================================
fish define setup_gage
; 
; Find gage particles to record axial strain when there are no walls
; and set initial sample height
;

  global vertical_direction = global.dim
  
  ; assume x and y are approximately the same in 3D
  local bottom_ = 1.0e12
  local top_ = -1.0e12
  loop foreach bp ball.list
    top_ = math.max(ball.pos(bp,vertical_direction),top_)
    bottom_ = math.min(ball.pos(bp,vertical_direction),bottom_)
  end_loop
  if global.dim = 2
    local top_loc = vector(0.0,top_)
    local bottom_loc = vector(0.0,bottom_)
  else
    top_loc = vector(0.0,0.0,top_)
    bottom_loc = vector(0.0,0.0,bottom_)
  end_if
  
  global gage_top = ball.near(top_loc)
  global gage_bottom = ball.near(bottom_loc)
  global sample_height = ball.pos(gage_top,vertical_direction) ...
                         - ball.pos(gage_bottom,vertical_direction)
end
========================================================
fish define axial_strain_gage
;
; Compute axial strain using gage balls (positive tension)
;
; Assumes global variables gage_top, 
; gage_bottom and sample_height have been set
;
  axial_strain_gage = (ball.disp(gage_top,vertical_direction) ...
                     - ball.disp(gage_bottom,vertical_direction))/sample_height
end
==========================================================
fish define loadhalt_meas
;
; Function used to stop a test when axial stress decreases
; to some fraction of the peak
; Assumes measurement circle 1 is installed
; 
; INPUT: peak_fraction - fraction of peak stress that dictates 
;                        the stopping of the test (global float)
;
  loadhalt_meas = 0
  local mp = measure.find(1)
  local abs_stress = math.abs(meas.stress(mp,global.dim,global.dim))
  global peak_stress = math.max(abs_stress,peak_stress)
  if abs_stress < peak_stress*peak_fraction
    loadhalt_meas = 1
  end_if
end
;=============================================================================
; eof: ss_gage.fis
Top