Data Files for “Probing a Granular Specimen” Example

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