Data files for “Genesis and Testing of a Soft-Bonded Material”Example

doall_A10.dat

model new

; Generate packing
model random 10001
[resolution = 10]   ; approx. number of balls accross specimen width 
[davg  = 2.5e-3]    ; average ball diameter (m)
[dratio = 1.66]     ; ball max/min diameter ratio [-]
[tag = 'A'+ string(resolution)]
program call 'make_brick.dat' 

; Install contact model
model restore [string.build('brick_%1-iso',tag)]
[igap     = 0.0*davg ] 
[emod     = 50.0e9   ]
[kratio   = 2.5      ]
[fric     = 0.5      ]
[ten_mean = 20.0e6  ]
[ten_std  = 0.0e6    ]
[coh_mean = 100.0e6  ]
[coh_std  = 0.0e6    ]
[soft     = 0.0      ]
[cut      = 1.0      ]
[tag += '_G0_E50_KR25_S0_C1']
program call 'bond_sb.dat' 

; Perform Direct Tension Tests
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'tt.dat'
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 0.1]
[tag += '_SR01']
program call 'tt.dat'

; Perform Unconfined Compression Test
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'ucs.dat'
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 0.1]
[tag += '_SR01']
program call 'ucs.dat'

; Perform Direct Tension Test (using domain distortion)
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'tt2.dat'
; Perform Unconfined Compression Test (using domain distortion)
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'ucs2.dat'

; Perform Confined Compression Tests
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[pc = -10.0e3]
[strain_rate = 1.0]
[tag = tag + '_P10_SR1']
program call 'triax.dat'
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[pc = -100.0e3]
[strain_rate = 1.0]
[tag = tag + '_P100_SR1']
program call 'triax.dat'
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[pc = -500.0e3]
[strain_rate = 1.0]
[tag = tag + '_P500_SR1']
program call 'triax.dat'
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[pc = -1000.0e3]
[strain_rate = 1.0]
[tag = tag + '_P1000_SR1']
program call 'triax.dat'
[tag = 'A10_G0_E50_KR25_S0_C1']
model restore [string.build('brick_%1-bonded',tag)]
[pc = -10000.0e3]
[strain_rate = 1.0]
[tag = tag + '_P10000_SR1']
program call 'triax.dat'
program return
Top

doall_A10_S10_C0.dat

model new

[tag = 'A10']
model restore [string.build('brick_%1-iso',tag)]
; Install contact model
[igap     = 0.0*davg ] 
[emod     = 50.0e9   ]
[kratio   = 2.5      ]
[fric     = 0.5      ]
[ten_mean = 20.0e6  ]
[ten_std  = 0.0e6    ]
[coh_mean = 100.0e6  ]
[coh_std  = 0.0e6    ]
[soft     = 10.0     ]
[cut      = 0.0      ]
[tag += '_G0_E50_KR25_S10_C0']
program call 'bond_sb.dat' 

; Perform Direct Tension Test (using domain distortion)
[tag = 'A10_G0_E50_KR25_S10_C0']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'tt2.dat'
; Perform Unconfined Compression Test (using domain distortion)
[tag = 'A10_G0_E50_KR25_S10_C0']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'ucs2.dat'

program return
Top

doall_A10_S20_C0_COH1000.dat

model new

[tag = 'A10']
model restore [string.build('brick_%1-iso',tag)]
; Install contact model
[igap     = 0.0*davg ] 
[emod     = 50.0e9   ]
[kratio   = 2.5      ]
[fric     = 0.5      ]
[ten_mean = 20.0e6  ]
[ten_std  = 0.0e6    ]
[coh_mean = 100.0e7  ]
[coh_std  = 0.0e6    ]
[soft     = 20.0     ]
[cut      = 0.0      ]
[tag += '_G0_E50_KR25_S20_C0_COH1000']
program call 'bond_sb.dat' 

; Perform Direct Tension Test (using domain distortion)
[tag = 'A10_G0_E50_KR25_S20_C0_COH1000']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'tt2.dat'
; Perform Unconfined Compression Test (using domain distortion)
[tag = 'A10_G0_E50_KR25_S20_C0_COH1000']
model restore [string.build('brick_%1-bonded',tag)]
[strain_rate = 1.0]
[tag += '_SR1']
program call 'ucs2.dat'

program return
Top

make_brick.dat

; filename: make_brick.dat
; 
; Generate an isotropic, homogeneous assembly of balls with periodic
; boundary conditions, and settle under low isotropic confinement 
; using a servo-control algorithm to distort the periodic domain 
;
; ===================================================================

; Setup domain extent (aspect ratio 2:1) and boundary conditions
model domain extent [-0.5*resolution*davg] [0.5*resolution*davg] ...
                    [-0.5*resolution*davg] [0.5*resolution*davg] ...
                    [-1.0*resolution*davg] [1.0*resolution*davg] ...
                     condition periodic
                     
; Distribute balls to a target low porosity and set density and 
; local damping attributes
[dmin = 2.0*davg / (1.0 + dratio)]
[dmax = dmin * dratio]
ball distribute porosity 0.36 radius [0.5*dmin] [0.5*dmax] 
ball attribute density 2000.0 damp 0.7

; Choose the contact model
contact cmat default model linear ...
                     method deformability emod 1e8 kratio 2.5

; Monitoring
model history name 'ratio-average' mechanical ratio-average
model history name 'timestep' timestep

; Perform some initial cycles and periodically reset the 
; spin and velocities to zero to relax the inital overlaps,
; then solve to equilibrium 
model cycle 1000 calm 100
model solve ratio-average 5e-3
model save [string.build('brick_%1-pkg',tag)]

; Use a servo-control algorithm to settle the assembly under
; a presecribed low isotropic stress state

; Load utility functions
program call 'servo.fis' suppress

; Parameter definitions (used by the servo-control function)
[servo_do_z = true     ]  
[servo_ts_xx = -1.0e5  ]
[servo_ts_yy = -1.0e5  ]
[servo_ts_zz = -1.0e5  ]
[servo_gain  =  1.0e-4 ]
[servo_tol = 1e-2]
[s_xx = 0.0]
[s_yy = 0.0]
[s_zz = 0.0]

[Lx = domain.max.x() - domain.min.x()]
[Ly = domain.max.y() - domain.min.y()]
[Lz = domain.max.z() - domain.min.z()]
measure create id 1 radius [0.95*0.5*Lx]

; Monitoring
fish history name 's_xx' @s_xx
fish history name 's_yy' @s_yy
fish history name 's_zz' @s_zz
fish history name 'L_x'  @Lx
fish history name 'L_y'  @Ly
fish history name 'L_z'  @Lz

model domain tolerance 1e-3
model solve fish-call 1.0 @servo_domain ...
            fish-halt @servo_halt ...
            ratio-average 1e-3 ...
            cycles 100 and
model save [string.build('brick_%1-iso',tag)]
; ===================================================================
program return
; EOF: make_brick.dat
Top

servo.fis

; filename: servo.fis
; 
; Utility function to define a servo-control algorithm 
; to distort the periodic domain to apply a target confinement
;
; ===================================================================

fish define servo_domain
  global Lx = domain.max.x() - domain.min.x()
  global Ly = domain.max.y() - domain.min.y()
  global Lz = domain.max.z() - domain.min.z()
  local vol = Lx*Ly*Lz
  local m = measure.find(1)
  local mr = math.sqrt((0.5*Lx)^2+(0.5*Ly)^2+(0.5*Lz)^2)
  local mv = (4.0/3.0)*math.pi*mr^3
  measure.radius(m) = mr
  local stress = measure.stress(m) *mv/vol
  global s_xx = stress(1,1)
  global servo_xdiff = servo_ts_xx - s_xx
  global servo_xrate = servo_gain *servo_xdiff 
  domain.strain.rate.xx = servo_xrate
  global s_yy = stress(2,2)
  global servo_ydiff = servo_ts_yy - s_yy
  global servo_yrate = servo_gain *servo_ydiff 
  domain.strain.rate.yy = servo_yrate
  global s_zz = stress(3,3)
  if servo_do_z then 
    global servo_zdiff = servo_ts_zz - s_zz
    global servo_zrate =  servo_gain * servo_zdiff
  endif
  domain.strain.rate.zz = servo_zrate
end

fish define servo_halt
  servo_halt = false
  local delta = math.max(math.abs(servo_xdiff/servo_ts_xx) ...
               ,math.abs(servo_ydiff/servo_ts_yy))
  if servo_do_z then 
    delta = math.max(math.abs(servo_zdiff/servo_ts_zz),delta)
  endif
  if delta < servo_tol then
    servo_halt = true
    domain.strain.rate.xx = 0.0
    domain.strain.rate.yy = 0.0
    domain.strain.rate.zz = 0.0
  endif
end
; ===================================================================
program return
; EOF: servo.fis
Top

bond_sb.dat

ball attribute force-contact  multiply 0.0  ...
               moment-contact multiply 0.0  ...
               velocity multiply 0.0        ...
               spin multiply 0.0            ...
               displacement multiply 0.0
contact cmat default model 'softbond' ...
                   method deformability emod @emod kratio @kratio ...
                   property fric @fric  sb_mode 1                 ...
                            sb_soft @soft sb_cut @cut             
contact cmat proximity @igap
model clean all
contact cmat apply

fish define normal_dist(c,mean,sdev)
; Return a non-negative value from a normal distribution with
; given mean and standard deviation.
;
  loop while 1 # 0
    local _val = mean + (math.random.gauss * sdev)
    if _val >= 0.0 then
      normal_dist = _val
      exit
    end_if
  end_loop
end

contact property sb_coh @normal_dist(@coh_mean,@coh_std)
contact property sb_ten @normal_dist(@ten_mean,@ten_std)

model save [string.build('brick_%1-bonded',tag)]
Top

tt.dat

fish define setup_tt
  ; identify boundary balls
  loop foreach local c contact.list.all
      if contact.offset.z(c) # 0 then
          ball.group(contact.end1(c)) = 'boundary'
          ball.group(contact.end2(c)) = 'boundary'
      endif
  endloop
  
  ; install measurement sphere
  Lx0 = domain.max.x() - domain.min.x()
  Ly0 = domain.max.y() - domain.min.y()
  Lz0 = domain.max.z() - domain.min.z()
  local mr = 0.5*math.min(Lx0,Ly0,Lz0)
  measure.radius(measure.find(1)) = 0.9*mr  
  global btop = ball.near(0.0,0.0,domain.max.z()) 
  global bbot = ball.near(0.0,0.0,domain.min.z()) 
  ball.group(btop,'gage') = 'top'
  ball.group(bbot,'gage') = 'bot'
  global bright = ball.near(domain.max.x(),0.0,0.0) 
  global bleft  = ball.near(domain.min.x(),0.0,0.0) 
  global bfront = ball.near(0.0,domain.max.y(),0.0) 
  global bback  = ball.near(0.0,domain.min.y(),0.0) 
  ball.group(bright,'gage') = 'right'
  ball.group(bleft,'gage')  = 'left'
  ball.group(bfront,'gage') = 'front'
  ball.group(bback,'gage')  = 'back'
  Lx0 = ball.pos.x(bright) - ball.pos.x(bleft)
  Ly0 = ball.pos.y(bfront) - ball.pos.y(bback)
  Lz0 = ball.pos.z(btop)   - ball.pos.z(bbot)
  ; enlarge domain extent
  domain.min.x = 1.5*domain.min.x
  domain.max.x = 1.5*domain.max.x
  domain.min.y = 1.5*domain.min.y
  domain.max.y = 1.5*domain.max.y
  domain.min.z = 1.5*domain.min.z
  domain.max.z = 1.5*domain.max.z  
end
@setup_tt

model domain condition stop
model clean all
contact method bond gap @igap range contact all 

; reset the proximity, delete inactive contacts 
; and clean the model to reduce the number of inactive contacts
contact cmat proximity 0.0
contact delete range contact inactive
model clean all

ball fix velocity spin range group 'boundary'
;contact property sb_ten 1e20 sb_coh 1e20 range group 'boundary' match 1

fish define e_xx
  Lx = ball.pos.x(bright) - ball.pos.x(bleft)
  e_xx = (Lx - Lx0) / Lx0
end
fish define e_yy
  Ly = ball.pos.y(bfront) - ball.pos.y(bback)
  e_yy = (Ly - Ly0) / Ly0
end
fish define e_zz
  Lz = ball.pos.z(btop) - ball.pos.z(bbot)
  e_zz = (Lz - Lz0) / Lz0
end
fish define e_zz
  Lz = ball.pos.z(btop) - ball.pos.z(bbot)
  e_zz = (Lz - Lz0) / Lz0
end
model cycle 1
[s_xx = measure.stress.xx(measure.find(1))]
[s_yy = measure.stress.yy(measure.find(1))]
[s_zz = measure.stress.zz(measure.find(1))]

fish define s_zz
  s_zz = measure.stress.zz(measure.find(1))
end

fish define poiss_ratio
  poiss_ratio = - 0.5*(e_xx+e_yy) / e_zz
end

fish define emod_secant
  emod_secant = s_zz / e_zz
end

history purge
fish history name 'e_xx' @e_xx 
fish history name 'e_yy' @e_yy 
fish history name 'e_zz' @e_zz 
fish history name 'emod_secant' @emod_secant
fish history name 'poiss_ratio' @poiss_ratio 

fish define halt
  halt = false
  s_zz_max = math.max(s_zz,s_zz_max)
  if s_zz < peak_frac*s_zz_max
    halt = 1
  endif
end

ball attribute velocity-z [ 0.5*strain_rate*Lz0] ...
    range position-z 0.0 @Lz0   ...
          group 'boundary'
ball attribute velocity-z [-0.5*strain_rate*Lz0] ...
    range position-z [-Lz0] 0.0 ...
          group 'boundary'

[peak_frac = 0.5 ]
model save [string.build('tt_%1-initial',tag)]
model solve fish-halt @halt ...
            cycles 100 and
model save [string.build('tt_%1-final',tag)]

; ===================================================================
program return
;EOF
Top

ucs.dat

fish define setup_ucs
  ; identify boundary balls
  loop foreach local c contact.list.all
      if contact.offset.z(c) # 0 then
          ball.group(contact.end1(c)) = 'boundary'
          ball.group(contact.end2(c)) = 'boundary'
      endif
  endloop
  
  global zmax = 0.0
  global zmin = 0.0
  local cntmax = 0
  local cntmin = 0
  loop foreach local b ball.groupmap('boundary')
    local z = ball.pos.z(b)
    if z > 0.0 then
        zmax   += z + ball.radius(b)
        cntmax +=1
    else
        zmin   += z - ball.radius(b)
        cntmin +=1    
    endif
  endloop
  zmax /= cntmax
  zmin /= cntmin
    
  ; install measurement sphere
  Lx0 = domain.max.x() - domain.min.x()
  Ly0 = domain.max.y() - domain.min.y()
  Lz0 = domain.max.z() - domain.min.z()
  global mr = 0.5*math.min(Lx0,Ly0,Lz0)
  measure.radius(measure.find(1)) = 0.9*mr
  
  global bright = ball.near(domain.max.x(),0.0,0.0) 
  global bleft  = ball.near(domain.min.x(),0.0,0.0) 
  global bfront = ball.near(0.0,domain.max.y(),0.0) 
  global bback  = ball.near(0.0,domain.min.y(),0.0) 
  ball.group(bright,'gage') = 'right'
  ball.group(bleft,'gage')  = 'left'
  ball.group(bfront,'gage') = 'front'
  ball.group(bback,'gage')  = 'back'
  Lx0 = ball.pos.x(bright) - ball.pos.x(bleft)
  Ly0 = ball.pos.y(bfront) - ball.pos.y(bback)
  
  ; enlarge domain extent
  domain.min.x = 1.5*domain.min.x
  domain.max.x = 1.5*domain.max.x
  domain.min.y = 1.5*domain.min.y
  domain.max.y = 1.5*domain.max.y
  domain.min.z = 1.5*domain.min.z
  domain.max.z = 1.5*domain.max.z
  ;
  
end
@setup_ucs

model domain condition stop
model clean all                           
contact method bond gap @igap range contact all 

; reset the proximity, delete inactive contacts 
; and clean the model to reduce the number of inactive contacts
contact cmat proximity 0.0
contact delete range contact inactive
model clean all

wall generate name 'box' box [-0.5*Lx0]  [0.5*Lx0] ...
                             [-0.5*Ly0]  [0.5*Ly0] ...
                             [-0.5*Lz0]  [0.5*Lz0] ...
                             expand 1.25
wall delete range set name 'boxLeft'                             
wall delete range set name 'boxRight'                             
wall delete range set name 'boxFront'                             
wall delete range set name 'boxBack'                             

contact cmat default type 'ball-facet' ...
                     model softbond ...
                     method deformability emod [emod] kratio @kratio ...
                     property sb_mode 1 fric @fric

model clean all

[wtop = wall.find('boxTop')]
[wbot = wall.find('boxBottom')]
[Lz0 = wall.pos.z(wtop) - wall.pos.z(wbot)]

fish define e_xx
  Lx = ball.pos.x(bright) - ball.pos.x(bleft)
  e_xx = (Lx - Lx0) / Lx0
end
fish define e_yy
  Ly = ball.pos.y(bfront) - ball.pos.y(bback)
  e_yy = (Ly - Ly0) / Ly0
end
fish define e_zz
  Lz = wall.pos.z(wtop) - wall.pos.z(wbot)
  e_zz = (Lz - Lz0) / Lz0
end
model cycle 1
[s_xx = measure.stress.xx(measure.find(1))]
[s_yy = measure.stress.yy(measure.find(1))]
[s_zz = measure.stress.zz(measure.find(1))]

fish define s_zz
  s_zz = measure.stress.zz(measure.find(1))
end

fish define poiss_ratio
  poiss_ratio = - 0.5*(e_xx+e_yy) / e_zz
end

fish define emod_secant
  emod_secant = s_zz / e_zz
end

history purge
fish history name 'e_xx' @e_xx 
fish history name 'e_yy' @e_yy 
fish history name 'e_zz' @e_zz 
fish history name 'emod_secant' @emod_secant
fish history name 'poiss_ratio' @poiss_ratio 

fish define halt
  halt = false
  s_zz_min = math.min(s_zz,s_zz_min)
  if s_zz > peak_frac*s_zz_min
    halt = 1
  endif
end

wall attribute velocity-z [-0.5*strain_rate*Lz0] range set name 'boxTop'
wall attribute velocity-z [ 0.5*strain_rate*Lz0] range set name 'boxBottom'
model save [string.build('ucs_%1-initial',tag)]
[peak_frac = 0.8 ]
model solve fish-halt @halt ...
            cycles 100 and
model save [string.build('ucs_%1-final',tag)]

; ===================================================================
program return
;EOF
Top

tt2.dat

fish define setup_tt
  ; identify boundary balls
  loop foreach local c contact.list.all
      if contact.offset.z(c) # 0 then
          ball.group(contact.end1(c)) = 'boundary'
          ball.group(contact.end2(c)) = 'boundary'
      endif
  endloop
  
  ; install measurement sphere
  Lx0 = domain.max.x() - domain.min.x()
  Ly0 = domain.max.y() - domain.min.y()
  Lz0 = domain.max.z() - domain.min.z()
  local mr = 0.5*math.min(Lx0,Ly0,Lz0)
  measure.radius(measure.find(1)) = 0.9*mr
  global btop = ball.near(0.0,0.0,domain.max.z()) 
  global bbot = ball.near(0.0,0.0,domain.min.z()) 
  global bright = ball.near(domain.max.x(),0.0,0.0) 
  global bleft  = ball.near(domain.min.x(),0.0,0.0) 
  global bfront = ball.near(0.0,domain.max.y(),0.0) 
  global bback  = ball.near(0.0,domain.min.y(),0.0) 
  ball.group(bright,'gage') = 'right'
  ball.group(bleft,'gage')  = 'left'
  ball.group(bfront,'gage') = 'front'
  ball.group(bback,'gage')  = 'back'
  ball.group(btop,'gage') = 'top'
  ball.group(bbot,'gage') = 'bot'
 
  Lx0 = ball.pos.x(bright) - ball.pos.x(bleft)
  Ly0 = ball.pos.y(bfront) - ball.pos.y(bback)
  Lz0 = ball.pos.z(btop) - ball.pos.z(bbot)

  ; enlarge domain extent
  domain.min.x = 1.5*domain.min.x
  domain.max.x = 1.5*domain.max.x
  domain.min.y = 1.5*domain.min.y
  domain.max.y = 1.5*domain.max.y
  domain.min.z = 1.5*domain.min.z
  domain.max.z = 1.5*domain.max.z
  
end
@setup_tt

model domain condition stop
model clean all

contact method bond gap @igap range contact all 

; reset the proximity, delete inactive contacts 
; and clean the model to reduce the number of inactive contacts
contact cmat proximity 0.0
contact delete range contact inactive
model clean all

ball fix velocity-z range group 'boundary'

fish define e_xx
  Lx = ball.pos.x(bright) - ball.pos.x(bleft)
  e_xx = (Lx - Lx0) / Lx0
end

fish define e_yy
  Ly = ball.pos.y(bfront) - ball.pos.y(bback)
  e_yy = (Ly - Ly0) / Ly0
end

fish define e_zz
  Lz = ball.pos.z(btop) - ball.pos.z(bbot)
  e_zz = (Lz - Lz0) / Lz0
end
model cycle 1

[s_xx = measure.stress.xx(measure.find(1))]
[s_yy = measure.stress.yy(measure.find(1))]
[s_zz = measure.stress.zz(measure.find(1))]

fish define s_zz
  s_zz = measure.stress.zz(measure.find(1))
end

fish define poiss_ratio
  poiss_ratio = - 0.5*(e_xx+e_yy) / e_zz
end

fish define emod_secant
  emod_secant = s_zz / e_zz
end

history purge
fish history name 'e_xx' @e_xx 
fish history name 'e_yy' @e_yy 
fish history name 'e_zz' @e_zz 
fish history name 'emod_secant' @emod_secant
fish history name 'poiss_ratio' @poiss_ratio 

fish define halt
  halt = false
  s_zz_max = math.max(s_zz,s_zz_max)
  if s_zz < peak_frac*s_zz_max
    halt = 1
  endif
end

model domain strain-rate-zz @strain_rate

[peak_frac = 0.8 ]
model save [string.build('tt2_%1-initial',tag)]
model solve fish-halt @halt ...
            cycles 100 and
model save [string.build('tt2_%1-final',tag)]

; ===================================================================
program return
;EOF
Top

ucs2.dat

fish define setup_tt
  ; identify boundary balls
  loop foreach local c contact.list.all
      if contact.offset.z(c) # 0 then
          ball.group(contact.end1(c)) = 'boundary'
          ball.group(contact.end2(c)) = 'boundary'
      endif
  endloop
  
  ; install measurement sphere
  Lx0 = domain.max.x() - domain.min.x()
  Ly0 = domain.max.y() - domain.min.y()
  Lz0 = domain.max.z() - domain.min.z()
  local mr = 0.5*math.min(Lx0,Ly0,Lz0)
  measure.radius(measure.find(1)) = 0.9*mr
  
  global bright = ball.near(domain.max.x(),0.0,0.0) 
  global bleft  = ball.near(domain.min.x(),0.0,0.0) 
  global bfront = ball.near(0.0,domain.max.y(),0.0) 
  global bback  = ball.near(0.0,domain.min.y(),0.0) 
  global btop = ball.near(0.0,0.0,domain.max.z()) 
  global bbot = ball.near(0.0,0.0,domain.min.z()) 
  ball.group(bright ,'gage') = 'right'
  ball.group(bleft  ,'gage') = 'left'
  ball.group(bfront ,'gage') = 'front'
  ball.group(bback  ,'gage') = 'back'  
  ball.group(btop   ,'gage') = 'top'
  ball.group(bbot   ,'gage') = 'bot'

; enlarge domain extent
  domain.min.x = 1.5*domain.min.x
  domain.max.x = 1.5*domain.max.x
  domain.min.y = 1.5*domain.min.y
  domain.max.y = 1.5*domain.max.y
  domain.min.z = 1.5*domain.min.z
  domain.max.z = 1.5*domain.max.z
  ;
  Lx0 = ball.pos.x(bright) - ball.pos.x(bleft)
  Ly0 = ball.pos.y(bfront) - ball.pos.y(bback)
  Lz0 = ball.pos.z(btop) - ball.pos.z(bbot)  
end
@setup_tt

model domain condition stop
model clean all

contact method bond gap @igap range contact all 

; reset the proximity, delete inactive contacts 
; and clean the model to reduce the number of inactive contacts
contact cmat proximity 0.0
contact delete range contact inactive
model clean all
model cycle 1
[s_xx = measure.stress.xx(measure.find(1))]
[s_yy = measure.stress.yy(measure.find(1))]
[s_zz = measure.stress.zz(measure.find(1))]

ball fix velocity-z range group 'boundary'

fish define e_xx
  Lx = ball.pos.x(bright) - ball.pos.x(bleft)
  e_xx = (Lx - Lx0) / Lx0
end
fish define e_yy
  Ly = ball.pos.y(bfront) - ball.pos.y(bback)
  e_yy = (Ly - Ly0) / Ly0
end
fish define e_zz
  Lz = ball.pos.z(btop) - ball.pos.z(bbot)
  e_zz = (Lz - Lz0) / Lz0
end

fish define poiss_ratio
  poiss_ratio = - 0.5*(e_xx+e_yy) / e_zz
end

fish define emod_secant
  emod_secant = s_zz / e_zz
end

history purge
fish history name 'e_xx' @e_xx 
fish history name 'e_yy' @e_yy 
fish history name 'e_zz' @e_zz 
fish history name 'emod_secant' @emod_secant
fish history name 'poiss_ratio' @poiss_ratio 

fish define s_zz
  s_zz = measure.stress.zz(measure.find(1))
end

fish define halt
  halt = false
  s_zz_min = math.min(s_zz,s_zz_min)
  if s_zz > peak_frac*s_zz_min
    halt = 1
  endif
end

model domain strain-rate-zz [-strain_rate]

[peak_frac = 0.8 ]
model save [string.build('ucs2_%1-initial',tag)]
model solve fish-halt @halt ...
            cycles 100 and
model save [string.build('ucs2_%1-final',tag)]

; ===================================================================
program return
;EOF
Top

triax.dat

contact method bond gap @igap range contact all 

contact cmat proximity 0.0
contact delete range contact inactive
model clean all

; Bring to istotropic confinement
[Lx0 = domain.max.x() - domain.min.x()]
[Ly0 = domain.max.y() - domain.min.y()]
[Lz0 = domain.max.z() - domain.min.z()]

fish define e_xx
  Lx = domain.max.x() - domain.min.x()
  e_xx = (Lx - Lx0) / Lx0
end
fish define e_yy
  Ly = domain.max.y() - domain.min.y()
  e_yy = (Ly - Ly0) / Ly0
end
fish define e_zz
  Lz = domain.max.z() - domain.min.z()
  e_zz = (Lz - Lz0) / Lz0
end
model cycle 1
history purge
fish history name 'e_xx' @e_xx 
fish history name 'e_yy' @e_yy 
fish history name 'e_zz' @e_zz 

[servo_do_z  = true    ]
[servo_ts_xx = pc      ]
[servo_ts_yy = pc      ]
[servo_ts_zz = pc      ]
[servo_gain  =  1.0e-5 ]
[servo_tol   =  1.0e-5 ]
model solve fish-call 1.0 @servo_domain ...
            fish-halt @servo_halt ...
            ratio-average 1e-5 ...
            cycles 100 and

model save [string.build('triax_%1-iso',tag)]

; perform compression
history purge
[servo_do_z = false    ]
[servo_zrate = -strain_rate]
[servo_gain  =  1.0e-5 ]
[Lx0 = domain.max.x() - domain.min.x()]
[Ly0 = domain.max.y() - domain.min.y()]
[Lz0 = domain.max.z() - domain.min.z()]

fish define emod_secant
  emod_secant = (s_zz - 0.5*(s_xx+s_yy)) / e_zz
end
fish define poiss_ratio
  poiss_ratio = - 0.5*(e_xx+e_yy) / e_zz
end

fish history name 'emod_secant' @emod_secant 
fish history name 'poiss_ratio' @poiss_ratio 

fish define halt
  halt = false
  if e_zz < -emax
    halt = 1
  endif
end

[emax = 2.0e-4]
model solve fish-call 1.0 @servo_domain ...
            fish-halt @halt ...
            cycles 100 and
[ymod  =   (s_zz - 0.5*(s_xx+s_yy)) / e_zz]
[poiss = - 0.5*(e_xx+e_yy) / e_zz]
model save [string.build('triax_%1-elastic',tag)]

fish define halt
  halt = false
  s_zz_min = math.min(s_zz,s_zz_min)
  if s_zz > peak_frac*s_zz_min
    halt = 1
    exit
  endif
  if e_zz < -emax
    halt = 1
    exit
  endif  
end

[emax = 15.0e-3]
[peak_frac = 0.8 ]
model solve fish-call 1.0 @servo_domain ...
            fish-halt @halt ...
            cycles 100 and
model save [string.build('triax_%1-final',tag)]
; ===================================================================
program return
;EOF
Top