period triaxial test and apply shear stress

Asked by jsonscript on 2020-02-24

Hello, everyone, my program runs out of some problems, I want to consult.

My code is divided into three parts, generate particles in the packing,consolidate the packing to reach the required confining pressure,shear the packing to have the required initial shear stress, respectively.

When I applied the shear stress in the third part of the code , I applied the shear speed O.cell.velGrad[0,2]=0.01 to reach the initial shear stress, but the shear strain(gamma) was always 0, and the speed(.cell.velGrad[0,2]) printed on the terminal was always 0,
I don't know why.

Here is my code
############################################
### generate particles in the packing ###
############################################
from yade import plot
from yade import pack,qt

O.periodic=True
O.cell.setBox(.016,.016,.016)
radius=0.4e-3
O.materials.append(FrictMat(young=70e9,poisson=0.3,frictionAngle=0.13,density=2650*1e4,label='spheres'))
sp=pack.SpherePack()
sp.makeCloud(minCorner=(0,0,0),maxCorner=(.016,.016,.016),psdSizes=[0.00045,0.000525,0.0006,0.000675,0.00075,0.000825,0.0009],psdCumm=[0,0.167,.333,.5,.667,0.833,1],num=10000,distributeMass=True,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
V_sphere=utils.getSpheresVolume()#Compute the total volume of spheres in the simulation, mask parameter is considered
V_cell=O.cell.hSize[0,0] * O.cell.hSize[1,1] * O.cell.hSize[2,2]
print 'void ratio:',(V_cell-V_sphere)/V_sphere,'dt:',utils.PWaveTimeStep()
O.save('gen-ball.xml')
qt.View()

############################################
### consolidate the packing to reach the required confining pressure ###
############################################
from yade import plot,qt
from yade.pack import *
O=Omega()
O.load('gen-ball.xml')
radius=0.4e-3
for b in O.bodies:
  b.mat.frictionAngle=0.119

V_sphere=utils.getSpheresVolume() #the volume of the sphere
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=False)]
 ),
 PeriTriaxController(dynCell=True,maxUnbalanced=2e-3,relStressTol=2e-5,goal=[-1e5,-1e5,-1e5],stressMask=7,globUpdate=5,maxStrainRate=[0.1,0.1,0.1],doneHook='triaxDone()',label='triax'),
 NewtonIntegrator(damping=.2),
]

def triaxDone():
 print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
 O.pause()

O.dt=0.4*utils.PWaveTimeStep()
O.run()
O.wait()

for b in O.bodies:
  b.mat.frictionAngle=0.46364

for i in O.interactions:
  i.phys.tangensOfFrictionAngle=tan(0.46364)

O.engines[3].maxStrainRate=[0,0,0]
O.engines[3].maxUnbalanced=1e-6
O.engines[3].relStressTol=1e-3 #Relative stress tolerance

for i in O.interactions:
  un=i.geom.penetrationDepth
  i.phys.unMax=un
O.save('Phase1.xml.bz2')

#######################################
### shear the packing to have the required initial shear stress ###"
#######################################
from yade.pack import *

O=Omega()

O.load('Phase1.xml.bz2')

O.engines=O.engines[0:5]
O.engines[3].maxStrainRate=[0.01,0.01,0.01]
O.engines[3].maxUnbalanced=2e-4
O.run(1,True)

V_sphere=utils.getSpheresVolume() #the volume of the sphere

Pn=10000
cycle=0
g0=O.cell.hSize[0,2]*100.0/O.cell.hSize[2,2]#shear strain is used in %

def triaxDone():
 print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
 O.pause()

O.cell.velGrad[0,2]=0.01
O.dt=utils.PWaveTimeStep()

time=0
count=0
count1=0
zero=1
stop=0
while 1:
  #O.run(time*10,True)
  O.run(100,True)
  count=count+1
  count1=count1+1
  gamma=O.cell.hSize[0,2]*100/O.cell.hSize[2,2]-g0
  stress=utils.getStress()
  sigmazx=stress[0,2]/-1000
  sigmazz=stress[2,2]/1000.0
  if abs(sigmazx)>26:
    count=0
    count1=0
    O.cell.velGrad[0,2]=0
    O.save('Circular-x25-y0-m.xml.bz2')
    break

  print 'Speed:',O.cell.velGrad[0,2],'time:',O.realtime,'cycle:',cycle,'szx:',sigmazx,'stop:',stop

Thanks for your help

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-02-24
Last query:
2020-02-24
Last reply:
2020-02-24
Best Jan Stránský (honzik) said : #1

Hello,

because of implementation reasons, you have to set the velGrad as a Matrix3 at once (I can provide links if wanted). So
O.cell.velGrad[0,2]=Matrix3(0,0,0.01, 0,0,0, 0,0,0)
or
O.cell.velGrad = O.cell.velGrad + Matrix3(0,0,0.01, 0,0,0, 0,0,0)
should fix the problem

cheers
Jan

Jan Stránský (honzik) said : #3

> I think it should be O.cell.velGrad=Matrix3(0,0,0.01, 0,0,0, 0,0,0).

yes, of course, sorry for the typo
Jan

Hopefully the solution will appear more clearly in the documentation of velGrad after [1].
Nevertheless it was already clearly indicated there that assigning a single component would not work (just type "O.cell.velGrad?" after turning O.periodic True ); never a bad idea to check documentation of a parameter when you have problems with it.
B

[1] https://gitlab.com/yade-dev/trunk/-/commit/d7c2611a3cd5e78d1ecdaf231089b20d8936136d