Genesis and Testing of a Soft-Bonded Material

Problem Statement

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.

The following example illustrates genesis and testing of a bonded-particle model using the Soft-Bond contact model. The system composed of balls is first generated in a periodic domain, then settled under low confining stress. This is achieved by distorting the periodic domain dimensions in order to achieve the prescribed stress state, using a servo-algorithm implemented as a FISH function. The soft-bond contact model is then installed at existing contacts, while the contact forces are reset. From this initial, stress-free system, the following tests are performed: a direct tension test, an unconfined compression test and a confined compression test. The methodologies used at each stage of the modeling process are discussed, and the response of the material is investigated.

Numerical Model

Specimen Genesis

The initial packing is generated by calling the file make_brick.dat, which uses the FISH functions defined in servo.fis. The specimen is generated in the entire model domain, with periodic boundary conditions in all directions, to avoid boundary effects. The function servo_domain implements the servo-control algorithm procedure, which aims at updating the domain extent in each direction in order to achieve and maintain a prescribed stress. The stress in the specimen is monitored using a measurement sphere. Note that this sphere is extended to encompass the entire domain, and the value of the stress that it computes is scaled to use the actual volume of the domain.

At this stage, the objective is to generate a dense, well connected packing. Since the contact model and properties will be overridden in a subsequent stage, a simple linear contact model is used. The friction coefficient is set to zero, and the deformability of the contacts is set to a value reasonably large with respect to the prescribed target stress, to optimize the computation time while preventing large overlaps. Note that increasing the target stress, or using frictional contacts would modify the geometry of the packing, and therefore its subsequent behavior.

The file make_brick.dat is partially parameterized for convenience. The file doall_A10.dat demonstrates how it can be executed, and the parameters chosen in this example:

The average ball diameter is set to 2.5mm, with a maximum to minimum diameter ratio of 1.66, to prevent local arrangement of the balls into crystalline structures. The parameter resolution used in this example denotes the target number of balls across the lateral dimension of the specimen, which is set to have a 2:1 aspect ratio.

Contact Model and Properties

Once the initial packing has been created, the file bond_sb.dat is used to install the soft-bond contact model and prepare the specimen for further testing. The contact forces are reset to zero, as well as the translational and rotational velocities of the balls and their accumulated forces/moments and displacements, and the contact model is replaced with the soft-bond model. Please refer to the documentation of this contact model for a detailed description of its behavior and properties.

The file bond_sb.dat is partially parameterized for convenience. The file doall_A10.dat demonstrates how it can be executed, and the parameters chosen in this example follow:

The contact effective modulus is set to 50GPa, with a normal-to-shear stiffness ratio of 2.5 and a friction coefficient of 0.5. The contact tensile strength is set to 20MPa and the cohesion to 100MPa. Note that although the file bond_sb.dat is set to allow for a normal distribution of the strength parameters to be set, this functionality is not used in this example and all contacts are assigned the same strength parameters. Finally, the softening coefficient is set to 0.0 and the softening tensile strength factor is set to unity. With these values, softening is prohibited, and the model behaves essentially similar to the linear parallel-bond contact model with respect to failure.

Direct Tension Test

Two methods to simulate a direct tension test are described and compared below. The first method identifies balls on the two opposite ends of the specimen in the vertical direction, and fixes their velocity. The second method distorts the domain outward in the vertical direction.

Method 1: Using Ball Fixity

The file tt.dat is used to perform the direct tension test by fixing the velocities of a top and a bottom layer of balls. It is partially parameterized for convenience, and the file doall_A10.dat demonstrates how it can be executed, and the parameters chosen in this example follow:

The test is performed twice, with different values of the prescribed strain-rate (1.0 and 0.1). The stress-strain responses are shown in the figures below.

../../../../../_images/p3d-example-softbonded-tt1.png

Figure 1: Stress vs Strain for a direct tension test with a strain rate of 1.0.

../../../../../_images/p3d-example-softbonded-tt2.png

Figure 2: Stress vs Strain for a direct tension test with a strain rate of 0.1.

Both curves show the typical response of a brittle material, with a Young’s modulus of about 29GPa and a maximal tensile strength of about 8MPa in both cases. The strain rate has a small impact on the results: with the largest strain rate, the peak strength is slightly larger and reached at a slightly larger strain, with more pronounced oscillations at the beginning of the test.

Method 2: Using Domain Distortion

A similar test is performed by expanding the domain extent in the vertical direction. The file used to apply this test is tt2.dat, which is called by doall_A10.dat:

The stress-strain curve obtained in this case (strain-rate of 1.0) is shown in the figure below. The response agrees quantitatively well with that obtained with the first method (Figure Figure #exa-sbnd-tt1)

../../../../../_images/p3d-example-softbonded-tt3.png

Figure 3: Stress vs Strain for a direct tension test with a strain rate of 1.0 (method 2).

Unconfined Compression Test

Two methods to simulate an unconfined compression test are described and compared below. The first method uses walls on the two opposite ends of the specimen in the vertical direction, and fixes their velocity. The second method distorts the domain inward in the vertical direction at a constant rate.

Method 1: Using Wall Platens

The file ucs.dat is used to perform the unconfined compression test using wall platens. It is partially parameterized for convenience, and the file doall_A10.dat demonstrates how it can be executed, and the parameters chosen in this example follow:

The test is performed twice, with different values of the prescribed strain-rate (1.0 and 0.1). The stress-strain responses are shown in the figures below.

../../../../../_images/p3d-example-softbonded-ucs1.png

Figure 4: Stress vs Strain for an unconfined compression test with a strain rate of 1.0.

../../../../../_images/p3d-example-softbonded-ucs2.png

Figure 5: Stress vs Strain for an unconfined compression test with a strain rate of 0.1.

The peak compressive strength is about 32MPa, i.e. about 4 times the peak strength in tension, as typically observed for e.g. linear parallel-bonded models.

Method 2: Using Domain Distortion

A similar test is performed by contracting the domain extent in the vertical direction. The file used to apply this test is ucs2.dat, which is called by doall_A10.dat:

The stress-strain curve obtained in this case (strain-rate of 1.0) is shown in the figure below. The response agrees quantitatively well with that obtained with the first method (Figure 4)

../../../../../_images/p3d-example-softbonded-ucs3.png

Figure 6: Stress vs Strain for an unconfined compression test with a strain rate of 1.0 (method 2).

Confined Compression Test

The file triax.dat is used to perform confined compression tests at various confining stress. It operates with periodic boundary conditions, and uses the servo-algorithm introduced above to apply the target confinement, then to maintain the lateral stress while contracting the specimen in the vertical direction at a constant strain-rate. It is partially parameterized for convenience, and the file doall_A10.dat demonstrates how it can be executed, and the parameters chosen in this example (strain-rate of 1.0 and a confining stress of 10kPa, 100kPa, 500kPa, 1MPa and 10MPa):

The stress-strain curves obtained are shown in the figure below, which shows a transition to a more ductile response at large confinement.

../../../../../_images/p3d-example-softbonded-triax.png

Figure 7: Stress vs Strain for an confined compression tests with a strain rate of 1.0 and increasing confinement.

Results and Discussion

The files discussed above demonstrate an effective set of procedures that can be used to generate a bonded-particle model, and perform direct-tension, unconfined compression and confined compression tests. In this example, the soft-bond contact model is used. The results shown above all use a set of parameters that forbid softening, and the resulting behavior has been shown to be that of a brittle material, with a compressive-to-tensile strength ratio of about 4.0, which is the value typically obtained with linear parallel-bonded materials, but is typically much low than the ratio most rocks exhibit.

A second set of tests is performed by executing the file doall_A10_S10_C0.dat. In this case, softening is activated by setting a value of 10.0 for the softening factor, and a value of 0.0 for the tensile strength softening factor. This means that once the tensile strength is reached at a specific contact, this contact enters a softening regime with a softening stiffness that is 1/10th of the loading stiffness, and breaks only if the tensile force reaches 0.0. The stress-strain curves obtained for a direct tension and unconfined compression tests are shown below.

../../../../../_images/p3d-example-softbonded-S10-tt.png

Figure 8: Stress vs Strain for a direct tension with a softening parameter of 10.

../../../../../_images/p3d-example-softbonded-S10-ucs.png

Figure 9: Stress vs Strain for an unconfined compression with a softening parameter of 10.

Under direct tension, the stress-strain curve exhibits a clear softening before reaching the peak strength, which is of about 14MPa (compared to 8MPa with no softening). Under compression, the peak strength is reached at about 100MPa, which leads to a compressive-to-tensile strength ratio of about 7.

A third set of tests is performed by executing the file doall_A10_S20_C0_COH1000.dat. In this case, softening is activated by setting a value of 20.0 for the softening factor, a value of 0.0 for the tensile strength softening factor, and the bond cohesion is increased to forbid shear failure.

../../../../../_images/p3d-example-softbonded-S20-tt.png

Figure 10: Stress vs Strain for a direct tension with a softening parameter of 20 and no shear failure allowed.

../../../../../_images/p3d-example-softbonded-S20-ucs.png

Figure 11: Stress vs Strain for an unconfined compression with a softening parameter of 20 and no shear failure allowed.

With these parameters, the peak strength increases to a value of about 200MPa in compression, and a value of about 17MPa in tension, i.e. a compressive-to-tensile strength ratio of about 12.0.

Data Files

doall_A10.dat

model new
model large-strain on
; 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

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

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

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

; 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]
; Load utility functions
program call 'servo.fis' suppress

; 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

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

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)
; Compute 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)]

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' matches 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

[peak_frac = 0.5 ]
fish define halt
  halt = false
  global 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'

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

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 

[peak_frac = 0.8 ]
fish define halt
  halt = false
  global 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

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 

[peak_frac = 0.8 ]
fish define halt
  halt = false
  global 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

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

[peak_frac = 0.8]
fish define halt
  halt = false
  global 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

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

[emax = 2.0e-4]
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)]

[peak_frac = 0.8 ]
fish define halt
  halt = false
  global 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

Endnote