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.
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.
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.
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.
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.
Figure 5: Results of the UCS test on the flat-jointed model.
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.
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
Potyondy, D. O., and P. A. Cundall. “A Bonded-Particle Model for Rock,” Int. J. Rock Mech. & Min. Sci., 41(8), 1329-1364 (2004).
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
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Aug 13, 2024 |
