incrememntal compression with engine

Hello !

I try to simulate triaxial compression with confinement of 1 MPa test on 3 steps:
1) isotropic loading sigma1=sigma2=sigma3=1MPa
2) loading the specimen by incrementing sigma1=1 MPa every 200 iterations until peak
3) in the post-peak stage, i load the specimen with strain Rate = 0.01 m/s

The problem is the code take in consideration the first and the third step. In fact, when the isotropic loading finish, i have a loading with the strain rate and not with the incrementation of sigma1.

Jabrane

This is the code.

from yade import ymport, utils , plot
import math

PACKING='X1Y2Z1_20k'
OUT=PACKING+'_compressionTest'

DAMP=0.4
saveData=100
iterMax=2000000
saveVTK=10000

confinement=-1e6
delta_stress=-1e6
stress_max=-200e6
strainRate=-0.01

intR=1.4450546
DENS=4000
YOUNG=65e9
FRICT=10
ALPHA=0.4
TENS=8e6
COH=160e6

def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))

O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/numSpheres

O.reset()

mn,mx=Vector3(xinf+0.1*Rmean,yinf+0.1*Rmean,zinf+0.1*Rmean),Vector3(xsup-0.1*Rmean,ysup-0.1*Rmean,zsup-0.1*Rmean)
walls=utils.aabbWalls(oversizeFactor=1.5,extrema=(mn,mx),thickness=min(X,Y,Z)/100.,material=wallMat)
wallIds=O.bodies.append(walls)

O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),

        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
        NewtonIntegrator(damping=DAMP,label="newton"),
        PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk'),
 TriaxialStressController(internalCompaction=False,label='triax'),
]

O.step();
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

O.engines = O.engines[:6]+[TriaxialStressController(stressMask=7,goal1=confinement,goal2=confinement,goal3=confinement,max_vel=0.01,internalCompaction=False,label='triax')]

while 1:
  if confinement==0:
    O.run(1000,True) # to stabilize the system
    break
  O.run(100,True)
  unb=unbalancedForce()
  meanS=abs(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'unbalanced force:',unb,' mean stress: ',meanS
  if unb<0.005 and abs(meanS-abs(confinement))/abs(confinement)<0.001:
    O.run(1000,True) # to stabilize the system
    e10=triax.strain[0]
    e20=triax.strain[1]
    e30=triax.strain[2]
    break

O.bodies[wallIds[2]].mat.frictionAngle=radians(30)
O.bodies[wallIds[3]].mat.frictionAngle=radians(30)

O.engines = O.engines[:6]+[TriaxialStressController(stressMask=7,goal1=confinement,goal2=confinement,goal3=confinement,max_vel=0.1,internalCompaction=False,label='triax')]
for i in range(0,int(-1e-6*stress_max)):
 if ( abs(triax.goal2) < abs(stress_max) ):
  O.run(200,True)
  triax.goal2+=delta_stress

O.engines = O.engines[:6]+[TriaxialStressController(stressMask=5,goal1=confinement,goal2=strainRate,goal3=confinement,max_vel=0.01,internalCompaction=False,label='triax')]
O.run(200000,True)

O.run(iterMax)