Data Files for “Probing a Granular Specimen” Example
- StrainUtilities.p2fis
- StressUtilities.p2fis
- MakeSpecimen.p2dat
- CompactSpecimen.p2dat
- LoadUnload.p2dat
- BiaxialTest.p2dat
- MeasureStress.p2dat
- PureShear.p2dat
- SimpleShear.p2dat
- Doall.p2dat
StrainUtilities.p2fis
fish define ini_mstrain(sid)
command
ball attribute displacement multiply 0.0
endcommand
global mstrains = matrix(2,2)
global mp = measure.find(sid)
end
fish define accumulate_mstrain
global msrate = measure.strain.rate(mp)
global mstrains = mstrains + msrate * global.timestep
global xxmstrain = mstrains(1,1)
global xymstrain = mstrains(1,2)
global yxmstrain = mstrains(2,1)
global yymstrain = mstrains(2,2)
end
fish define ini_gstrain(ly)
local pos = vector(0.0,-0.5*ly)*0.75
global gb1 = ball.near(pos)
ball.group(gb1) = 'gauge_balls'
pos = vector( 0.0,-0.5*ly)*0.25
global gb2 = ball.near(pos)
ball.group(gb2) = 'gauge_balls'
pos = vector( 0.0, 0.5*ly)*0.25
global gb3 = ball.near(pos)
ball.group(gb3) = 'gauge_balls'
pos = vector(0.0, 0.5*ly)*0.75
global gb4 = ball.near(pos)
ball.group(gb4) = 'gauge_balls'
global gl12_0 = ball.pos(gb2) - ball.pos(gb1)
global gl13_0 = ball.pos(gb3) - ball.pos(gb1)
global gl14_0 = ball.pos(gb4) - ball.pos(gb1)
end
fish define xygstrain
local gl12 = ball.pos(gb2) - ball.pos(gb1)
global xygstrain12 = (comp.x(gl12) - comp.x(gl12_0))/ comp.y(gl12_0)
local gl13 = ball.pos(gb3) - ball.pos(gb1)
global xygstrain13 = (comp.x(gl13) - comp.x(gl13_0))/ comp.y(gl13_0)
local gl14 = ball.pos(gb4) - ball.pos(gb1)
global xygstrain14 = (comp.x(gl14) - comp.x(gl14_0))/ comp.y(gl14_0)
xygstrain = (xygstrain12 + xygstrain13 + xygstrain14) /3.0
end
StressUtilities.p2fis
fish define wsxx
wsxx = 0.5*(wall.force.contact.x(wp_left) ...
- wall.force.contact.x(wp_right))/ wly
end
fish define wsyy
wsyy = 0.5*(wall.force.contact.y(wp_bot) ...
- wall.force.contact.y(wp_top) )/ wlx
end
fish define compute_averagestress
global asxx = 0.0
global asxy = 0.0
global asyx = 0.0
global asyy = 0.0
loop foreach local contact contact.list("ball-ball")
local cforce = contact.force.global(contact)
local cl = ball.pos(contact.end2(contact)) ...
- ball.pos(contact.end1(contact))
asxx = asxx + comp.x(cforce)*comp.x(cl)
asxy = asxy + comp.x(cforce)*comp.y(cl)
asyx = asyx + comp.y(cforce)*comp.x(cl)
asyy = asyy + comp.y(cforce)*comp.y(cl)
endloop
asxx = - asxx / (wlx*wly)
asxy = - asxy / (wlx*wly)
asyx = - asyx / (wlx*wly)
asyy = - asyy / (wlx*wly)
end
fish define compute_spherestress(rad)
command
contact group 'insphere' remove
contact group 'insphere' range circle radius @rad
endcommand
global ssxx = 0.0
global ssxy = 0.0
global ssyx = 0.0
global ssyy = 0.0
loop foreach contact contact.groupmap("insphere","ball-ball")
local cf = contact.force.global(contact)
local cl = ball.pos(contact.end2(contact)) ...
- ball.pos(contact.end1(contact))
ssxx = ssxx + comp.x(cf) *comp.x(cl)
ssxy = ssxy + comp.x(cf) *comp.y(cl)
ssyx = ssyx + comp.y(cf) *comp.x(cl)
ssyy = ssyy + comp.y(cf) *comp.y(cl)
endloop
local vol = (math.pi*rad^2)
ssxx = - ssxx / vol
ssxy = - ssxy / vol
ssyx = - ssyx / vol
ssyy = - ssyy / vol
end
MakeSpecimen.p2dat
; fname: make_specimen.p2dat
;
; Generate a dense granular assembly within a box
;
; ============================================================================
; Load utility FISH functions for later use
program call 'StrainUtilities.p2fis' suppress
program call 'StressUtilities.p2fis' suppress
; Define domain extent and default contact model and properties in the CMAT
; note that since friction is not set (and defaults to zero), shear forces
; won't develop
model domain extent -0.1 0.1
contact cmat default model linear method deformability emod 1.0e8 kratio 2.5
; Generate a box-wall comprised of 4 independant walls. Each wall is expanded
; in preparation for upcoming motion, and global FISH variables that hold
; pointers to the walls are created for convenience.
wall generate name 'vessel' box -0.035 0.035 expand 1.5
[wp_left = wall.find('vesselLeft')]
[wp_right = wall.find('vesselRight')]
[wp_bot = wall.find('vesselBottom')]
[wp_top = wall.find('vesselTop')]
; FISH functions to get vessel dimensions
fish define wlx
wlx = wall.pos.x(wp_right) - wall.pos.x(wp_left)
end
fish define wly
wly = wall.pos.y(wp_top) - wall.pos.y(wp_bot)
end
; Generate a cloud of overlapping balls with a target porosity
; and assign density and local damping attributes
model random 10001
ball distribute porosity 0.2 radius 0.0006 0.0009 box -0.035 0.035
ball attribute density 2500.0 damp 0.7
;define ball and wall friction property
ball property 'fric' @ballFriction
wall property 'fric' @wallFriction
; solve to equilibr'ium
model cycle 1000 calm 10
model solve ratio-average 1e-5
model calm
; identify floaters using FISH function defined in make_utilities.p2fis
fish define identify_floaters
loop foreach local ball ball.list
ball.group.remove(ball,'floaters')
local contactmap = ball.contactmap(ball)
local size = map.size(contactmap)
if size <= 1 then
ball.group(ball) = 'floaters'
endif
endloop
end
@identify_floaters
@ini_gstrain(@wly)
program return
; ============================================================================
; eof: make_specimen.p2dat
CompactSpecimen.p2dat
[ly0 = wly]
[lx0 = wlx]
[v0 = wlx*wly]
[txx = -5.0e5]
[tyy = -5.0e5]
wall servo activate on force-x [ txx*wly] ...
velocity-max 0.1 range set name 'vesselRight'
wall servo activate on force-x [-txx*wly] ...
velocity-max 0.1 range set name 'vesselLeft'
wall servo activate on force-y [ tyy*wlx] ...
velocity-max 0.1 range set name 'vesselTop'
wall servo activate on force-y [-tyy*wlx] ...
velocity-max 0.1 range set name 'vesselBottom'
fish define servo_walls
wall.servo.force.x(wp_right) = txx*wly
wall.servo.force.x(wp_left) = -txx*wly
wall.servo.force.y(wp_top) = tyy*wlx
wall.servo.force.y(wp_bot) = -tyy*wlx
end
fish callback add @servo_walls 9.0
fish history name 41 @wsxx
fish history name 42 @wsyy
model orientation-tracking on
model calm
[tol = 5e-3]
fish define stop_me
if math.abs((wsyy - tyy)/tyy) > tol
exit
endif
if math.abs((wsxx - txx)/txx) > tol
exit
endif
if mech.solve("ratio-average") > 1e-6
exit
endif
stop_me = 1
end
ball attribute displacement multiply 0.0
model solve fish-halt @stop_me
measure create id 1 rad [0.4*(math.min(lx0,ly0))]
[porosity = measure.porosity(measure.find(1))]
@compute_spherestress([0.4*(math.min(lx0,ly0))])
@compute_averagestress
@identify_floaters
return
LoadUnload.p2dat
;define ball and wall friction property
ball property 'fric' @ballFriction
wall property 'fric' @wallFriction
[rate = 0.01]
wall attribute velocity-y [-rate*wly] range set name 'vesselTop'
wall attribute velocity-y [ rate*wly] range set name 'vesselBottom'
wall servo activate off range set name 'vesselTop' set name 'vesselBottom' union
model calm
[ly0 = wly]
[lx0 = wlx]
[wsyy0 = wsyy]
[wexx = 0.0]
[weyy = 0.0]
[wevol = 0.0]
fish define wexx
wexx = 2.0*(wlx - lx0) / (wlx+lx0)
end
fish define weyy
weyy = 2.0*(wly - ly0) / (wly+ly0)
end
fish define wevol
wevol = wexx + weyy
end
fish history name 51 @weyy
fish history name 52 @wexx
fish history name 53 @wevol
history purge
[stop_load = 0]
[target = 5e-4]
fish define stop_load
if -weyy >= target then
stop_load = 1
endif
end
model solve fish-halt @stop_load
[young_modulus = (wsyy-wsyy0)/weyy]
[poisson_ratio = -wexx/weyy]
[shear_modulus = 0.5*young_modulus/(1+poisson_ratio)]
model save 'biaxial-load'
wall attribute velocity-y [ rate*wly] range set name 'vesselTop'
wall attribute velocity-y [-rate*wly] range set name 'vesselBottom'
[stop_unload = 0]
fish define stop_unload
if weyy >= 0.0 then
stop_unload = 1
endif
end
model solve fish-halt @stop_unload
BiaxialTest.p2dat
;define ball and wall friction property
ball property 'fric' @ballFriction
wall property 'fric' @wallFriction
[ly0 = wly]
[lx0 = wlx]
[wexx = 0.0]
[weyy = 0.0]
[wevol = 0.0]
fish define wexx
wexx = (wlx - lx0) / lx0
end
fish define weyy
weyy = (wly - ly0) / ly0
end
fish define wevol
wevol = wexx + weyy
end
fish history name 51 @wexx
fish history name 52 @weyy
fish history name 53 @wevol
history purge
[rate = 0.2]
wall servo activate off range set name 'vesselTop'
wall servo activate off range set name 'vesselBottom'
wall attribute velocity-y [-rate*wly] range set name 'vesselTop'
wall attribute velocity-y [ rate*wly] range set name 'vesselBottom'
[stop_me = 0]
[target = 0.075]
fish define stop_me
if weyy <= -target then
stop_me = 1
endif
end
ball attribute displacement multiply 0.0
model calm
model solve fish-halt @stop_me
program return
MeasureStress.p2dat
;restore specimen
;measure create id 1 radius [0.4*math.min(wlx,wly)]
;list measure
;set echo off
; call 004_stress_utilities.p2fis
;set echo on
@compute_averagestress
;@compute_wallstress
@compute_spherestress([0.3*math.min(wlx,wly)])
return
PureShear.p2dat
wall servo activate off
fish callback remove @servo_walls 1.0
;define ball and wall friction property
ball property 'fric' @ballFriction
wall property 'fric' @wallFriction
[ly0 = wly]
[lx0 = wlx]
[defxy = 0]
[tauxy = 0]
fish define shear_strain_stress
x_strain = (wlx-lx0)
y_strain = (wly-ly0)
defxy = x_strain/ly0 + y_strain/lx0
tauxy = (wsxx-wsyy)/2
end
[load = 1]
[rate = 0.001]
[xfactor = 0.0]
[yfactor = 0.0]
fish define cyclic_shear
wall.vel.x(wp_left) = xfactor*rate ; left wall
wall.vel.x(wp_right) = -1.0 * xfactor*rate ; right wall
wall.vel.y(wp_top) = -1.0 * yfactor*rate ; top wall
wall.vel.y(wp_bot) = yfactor*rate ; bottom wall
if defxy < 0.0005
if load = 1
xfactor = 1.0
yfactor = -2.0
exit
else
if defxy < -0.0005
load = 1
exit
endif
exit
endif
else
xfactor = -1.0
yfactor = 2.0
load = 0
endif
end
fish history name 51 @tauxy
fish history name 61 @defxy
fish callback add @shear_strain_stress 2.0
fish callback add @cyclic_shear 3.0
model orientation-tracking on
model solve time 0.2
SimpleShear.p2dat
wall servo activate off
fish callback remove @servo_walls 1.0
;define ball and wall friction property
ball property 'fric' @ballFriction
wall property 'fric' @wallFriction
[rate = 1.0]
wall attribute velocity-x [ 1.0*rate*wly] range set name 'vesselTop'
wall attribute rotation-center [-0.5*wlx] [-0.5*wly] ...
spin [-1.0*rate] range set name 'vesselLeft'
wall attribute rotation-center [0.5*wlx] [-0.5*wly] ...
spin [-1.0*rate] range set name 'vesselRight'
[defxy = 0]
fish define shear_strain
defxy = defxy + rate*global.timestep
end
@ini_mstrain(1)
fish history name 11 @xxmstrain
fish history name 12 @xymstrain
fish history name 13 @yxmstrain
fish history name 14 @yymstrain
measure history name 31 stress-xx id 1
measure history name 32 stress-xy id 1
measure history name 33 stress-yx id 1
measure history name 34 stress-yy id 1
fish history name 61 @asxy
fish history name 51 @defxy
@ini_gstrain(@wly)
fish history name 21 @xygstrain
fish history name 22 @xygstrain12
fish history name 23 @xygstrain13
fish history name 24 @xygstrain14
fish callback add @shear_strain 2.0
fish callback add @accumulate_mstrain 11.0
fish callback add @compute_averagestress 43.0
model orientation-tracking on
model solve time 0.01
Doall.p2dat
;---------BIAXIAL TEST ON A LOOSE SAMPLE
model new
;define ball and wall friction values to generate a loose sample
[ballFriction = 0.3 ]
[wallFriction = 0.0 ]
[filename = 'loose' ]
program call 'MakeSpecimen'
model save ['specimen'+filename]
program call 'CompactSpecimen'
model save ['biaxial-iso'+filename]
program call 'BiaxialTest'
model save ['biaxial-final'+filename]
;---------BIAXIAL TEST ON A DENSE SAMPLE
model new
;define ball and wall friction values to generate a dense sample
[ballFriction = 0.0 ]
[wallFriction = 0.0 ]
[filename = 'dense' ]
program call 'MakeSpecimen'
model save ['specimen'+filename]
program call 'CompactSpecimen'
model save ['biaxial-iso'+filename]
[ballFriction = 0.3 ]
program call 'BiaxialTest'
model save ['biaxial-final'+filename]
;---------ESTIMATION OF ELASTIC PARAMETERS - LOAD UNLOAD TEST
[filename = 'dense']
model restore ['biaxial-iso'+filename]
[ballFriction = 0.3 ]
program call 'LoadUnload'
model save 'biaxial-unload'
;---------SIMPLE SHEAR TEST ON A DENSE SAMPLE
[filename = 'dense']
model restore ['biaxial-iso'+filename]
[ballFriction = 0.3 ]
[wallFriction = 0.3 ]
program call 'SimpleShear'
model save 'simpleshear-final'
;---------PURE SHEAR TEST ON A DENSE SAMPLE
[filename = 'dense']
model restore ['biaxial-iso'+filename]
[ballFriction = 0.3 ]
[wallFriction = 0.3 ]
program call 'PureShear'
model save 'pureshear-final'
program return
Was this helpful? ... | PFC 6.0 © 2019, Itasca | Updated: Nov 19, 2021 |