Rock Testing

Introduction

Note

The project file for this example may be viewed/run in PFC.[1] The data files used are shown at the end of this example.

This example demonstrates how to simulate laboratory testing of rock core samples; specifically an unconfined compression test and a direct tension test. In the laboratory, a direct tension test is difficult to perform, but in the idealized world of numerical models, this test is possible and provides a more reliable estimate of rock tensile strength than the indirect Brazilian test.

Two different models are compared: a sample bonded with linear parallel bonds and a sample that uses the flat joint contact model. It is well known that parallel bonded models in PFC produce unrealistically low ratios of compressive to tensile strength [Potyondy2004]. This problem can be overcome by using the flat joint contact model in which the contacts have a finite length and therefore particle rotations are resisted [Potyondy2012].

PFC2D Model

Unbonded Sample

First an unbonded assembly with the Linear model assigned to all contacts is created following a similar procedure as in the tutorial “Generating a Bonded Assembly.” The sample is 10 cm long and 5 cm in diameter, conforming to typical laboratory core samples. The walls are created such that they extend past the edge of the sample to allow for changes in the sample size during loading (see {”make_sample.p2dat” in 2D; “make_sample.p3dat” in 3D}).

The contact stiffnesses are set using contact methods instead of properties. The contact methods perform small calculations to automatically set the properties. So, for example, in 3D, the normal stiffness at a contact is proportional to the effective Young’s modulus and the particle radius as described here. We can therefore set the contact deformability property ‘emod’ (in stress units) and the linear stiffnesses will automatically be calculated based on the particle radii. In this example, it is assumed that stress quantities are in Pa. The unbonded sample is shown below.

../../../../../../../_images/p2d-tuto-rocktest-unbonded.png

Figure 1: Unbonded material model.

The parallel bonded sample is created using the commands in {”parallel_bonded.p2dat” in 2D; “parallel_bonded.p3dat” in 3D}. The parallel bond stiffnesses are again set using methods. Be careful with parallel bonds in that both the linear stiffness (contact method deform) and the bond material stiffness (contact method pb_deform) are set.

When the contact type is changed from linear to parallel bond, the linear stiffnesses are not inherited and must be set again.

One other thing to note is the use of lin-mode 1

when resetting the linear forces to 0. This ensures that contact normal forces will be generated based on relative displacements (incrementally), rather than absolute overlaps between entities.

Unconfined Compressive Strength

A UCS test is performed by simply moving the top wall downward and the bottom wall upward as shown in the file {”ucs.p2dat” in 2D; “ucs.p3dat” in 3D}. The only trick is recording the axial stress and strain. A small FISH function is used to sum all of the vertical forces acting on the walls and then dividing by the sample width to get the axial stress. The strain is simply calculated from the wall displacement divided by the initial sample length. The initial sample width and length are determined at the start of the test. Note that this will introduce a small error in the stress calculation because the sample width will actually change during the test and this is not accounted for in the function.

A set of FISH functions that are very useful when simulating rock tests are given in {”fracture.p2fis” in 2D; “fracture.p3fis” in 3D}. These functions track bond breakages, and for each bond breakage, a fracture is added to a DFN. The cracks can then be easily plotted during and/or after the test. In addition, the functions delineate fragments in the sample. A fragment is a piece of the model that is not connected to any other pieces with intact bonds.

The sample is loaded until the axial stress falls below 70% of the peak stress. The result of the UCS test on the parallel bonded sample is shown in the next figure and the axial stress versus the axial strain plotted on top of the sample. The stress and strain are calculated assuming compressive stress is negative, but in the next figure, the plot is reversed to show compression positive in keeping with rock mechanics conventions. The UCS for this sample is 33.9 MPa. Note that running this example on different computers may produce slightly different results due to the randomness inherent in the particle packing.

../../../../../../../_images/p2d-tuto-rocktest-parallel-ucs.png

Figure 2: Parallel bonded material after the UCS test. Cracks are shown in black and red lines indicating tensile and shear failure, respectively. Axial stress vs. axial strain (positive compression) is plotted as a green line.

Direct Tension Test

A direct tension test is simulated on the same numerical sample as was used in the UCS test (see {”tension.p2dat” in 2D; “tension.p3dat” in 3D}). To perform the test, the walls are deleted and a row of particles at the top and bottom of the sample are identified as “grips” (figure below). The top grip is moved upward and the bottom grip downward to apply tensile loading.

../../../../../../../_images/p2d-tuto-rocktest-parallel-grips.png

Figure 3: Particle grips used to apply tensile loading are shown in green and red.

Since the walls have been deleted, the stress and strain must be measured in a different way. Stress is measured by creating a single measurement circle in the center of the sample. The stress is computed by averaging the contact forces within the circle as described here.

The strain is computed by observing the displacement of two gage particles. A particle at the top of the model and another at the bottom. The relative displacement of the two gage particles divided by the initial distance between them gives the axial strain.

As with the UCS test, the model is loaded until the axial stress is 70% of the peak. The stress-strain response is shown in Figure 4 with the final sample configuration. Remember we are assuming compression positive, so the tensile strength is the most negative value in the plot. The peak tensile stress is 6.8 MPa. This yields a ratio of compressive to tensile strength of ≈5, considerably lower than the expected value of 10–20.

../../../../../../../_images/p2d-tuto-rocktest-parallel-tt.png

Figure 4: Parallel bonded material after the direct tension test plotted with axial stress versus axial strain (compression positive).

Flat-Joint Contact Model

It is well known that the parallel bonded model yields unrealistically low ratios of compressive to tensile strength and unrealistically low friction angles for rock [Potyondy2004]. A likely reason for this is that round particles are used to simulate angular grains. There is not enough resistance to rotation, especially after bonds have broken. The flat-joint contact model overcomes this deficiency by simulating a contact as a {flat line in 2D; disk in 3D} composed of many subcontacts such that rotation is resisted (see the “Flat-Joint Model” section for a detailed description).

The UCS and direct tension tests were repeated with a flat-jointed material (shown in {”flat_joint.p2dat” in 2D; “flat_joint.p3dat” in 3D}). The same properties were assigned as for the parallel bonded model. Stress-strain curves and failed samples for the two tests are shown in the next two figures. These plots show a peak compressive stress of 70.5 MPa and tensile strength of 7.1 MPa, yielding a ratio of compressive to tensile strength of ≈10.

Interestingly, the tensile strength is almost the same as for the parallel bonded material, but the compressive strength is much higher. The mechanism of tensile failure is approximately the same in both samples; the breaking of a few bonds in tension and rapid propagation of the fracture across the sample. In compression, the failure is much more complicated and involves some combination of axial splitting and shear failure along en-echelon arrays of tensile cracks. The flat-joint model shows many more cracks, but the strength is higher, possibly because the cracks cannot so easily coalesce to cause failure because particle rotations are somewhat suppressed.

../../../../../../../_images/p2d-tuto-rocktest-fj-ucs.png

Figure 5: Results of the UCS test on the flat-jointed model.

../../../../../../../_images/p2d-tuto-rocktest-fj-tt.png

Figure 6: Results of the direct tension test on the flat-jointed model.

PFC3D Model

The 3D model is essentially the same with a slightly larger ball size. The sample shape in 3D is a cylinder. The final state for the flat-jointed 3D model in the UCS test is shown below.

../../../../../../../_images/p3d-tuto-rocktest-fj-ucs.png

Figure 7: 3D flat-joint model after UCS test.

Discussion

This example shows how to simulate a UCS test and direct tension test and shows how the flat-joint contact model gives a more realistic compressive to tensile strength ratio than the parallel bonded model. This example is really a simplified version of what is provided in the material-modeling support package. This package creates and defines particular instances of parallel-bonded and flat-jointed materials, and also provides crack monitoring as well as facilities to perform standard rock-mechanics tests (including direct-tension and compression) upon these materials.

References

[Potyondy2004] (1,2)

Potyondy, D. O., and P. A. Cundall. “A Bonded-Particle Model for Rock,” Int. J. Rock Mech. & Min. Sci., 41(8), 1329-1364 (2004).

[Potyondy2012]

Potyondy, D.O. (2012) “A Flat-Jointed Bonded-Particle Material for Hard Rock,” paper ARMA 12-501 in Proceedings of 46th U.S. Rock Mechanics/Geomechanics Symposium, Chicago, USA, 24–27 June 2012.

Data Files

make_sample.dat (2D)

; fname: make_sample.p2dat
;
;  Create unbonded sample
;=========================================================================
model new
model large-strain on
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 deformability 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 automatic
model calm

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

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

parallel_bonded.dat (2D)

; 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 deformability emod 1.0e9 kratio 1.0

; set stiffness of bond material
contact method pb_deformability emod 1.0e9 kratio 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 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

ucs.dat (2D)

model large-strain on
; fname: ucs.pdat
;
;  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
[global peak_fraction = 0.7]
model solve fish-halt loadhalt_wall
fish list [peak_stress]

;=========================================================================
; eof: ucs.pdat

tension.dat (2D)

model large-strain on
; fname: tension.dat
;
;  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
[global peak_fraction = 0.7]
model cycle 1000 ; required, else initial loading ends prematurely - DJB 5/19
model solve fish-halt loadhalt_meas
fish list [peak_stress]

;=========================================================================
; eof: tension.dat

flat_joint.dat (2D)

; 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 
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 deformability emod 1.0e9 kratio 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_nr 4 fj_track on ...
    method deformability emod 1.0e9 kratio 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

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 fragment ids
;
;=============================================================================

fish define add_crack(entries)
    global crack_num
    global frag_time
    global xmin,xmax,ymin,ymax
    local contact    = entries(1)
    local mode       = entries(3)
    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.p2fis

make_sample.dat (3D)

; fname: make_sample.dat
;
;  Create unbonded sample
;=============================================================================
model new
model large-strain on
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 automatic
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 radius 0.02 not

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

parallel_bonded.dat (3D)

; fname: parallel_bonded.dat
;
;  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 deformability emod 1.0e9 kratio 1.0

; set stiffness of bonds using methods
contact method pb_deformability emod 1.0e9 kratio 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.dat

ucs.dat (3D)

model large-strain on
; fname: ucs.dat
;
;  Perform an unconfined compressive strength test on a bonded sample
;=============================================================================
program call 'ss_wall.fis' suppress
program call 'fracture.fis' suppress

; 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 vibrations
model cycle 1000

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

;=============================================================================
; eof: ucs.dat

tension.dat (3D)

model large-strain on
; fname: tension.dat
;
;  Perform a direct tension test on a bonded sample
;=============================================================================

program call 'ss_gage.fis' suppress
program call 'fracture.fis' suppress

; 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
fish list [peak_stress]

;=============================================================================
; eof: tension.dat

flat_joint.dat (3D)

; 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 
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 deformability emod 1.0e9 kratio 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 deformability emod 1.0e9 kratio 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.dat

doall_parallel.dat (3D)

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

; restore unbonded sample and bond it
program call 'parallel_bonded'

program log-file 'parallel3d.log'
program log on

; run ucs test
program call 'ucs'
model save 'parallel_ucs'

; run tension test
model restore 'parallel_bonded'
program call 'tension'
model save 'parallel_tension'

program log off

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

doall_fj.dat (3D)

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

; restore unbonded specimen and install flatjoint model
program call 'flat_joint'

program log-file 'fj3d.log'
program log on

; run ucs test
program call 'ucs'
model save 'fj_ucs'

; run tension test
model restore 'flatjoint'
program call 'tension'
model save 'fj_tension'

program log off

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

fracture.fis (3D)

; fname: fracture.fis
;
; 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(3)
    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 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
    global crack_num
    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
    global frag_time
    global xmin,xmax,ymin,ymax,zmin,zmax
    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 zmin = domain.min.z()

end

; ============================================================================
; eof: fracture.fis

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)
;
    global peak_fraction
    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

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 & 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)
;
    global peak_fraction
    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

Endnote