Some Problems about Direct Shear Test

Asked by Xue on 2019-09-09

Hello, everyone. I'm trying to simulate the direct shear process with the following program[1]. Before running the script, I generated particles in another script. But when I run the following script, the system can't find the particles. But when I input the commands of the script sentence by sentence in the Yade program (before importing the particles), I can import the particles normally. Is there someone can tell me where the problem is?
[1]the script:
from yade import ymport, utils, plot
from yade import ymport, utils, plot
from yade import pack,export
from yade import polyhedra_utils
from yade import geom,utils,plot
from yade import qt,ymport,timing
from yade import *

################## parameters definition
DAMPING=0.5
dtCOEFF=0.5
normalSTRESS=5e5
normalVEL=normalSTRESS/1e8 # 0.001 for 100kPa , i dont know the reason
shearVEL=1*normalVEL # try different values to ensure quasi-static conditions
intR=1.263 #1.263 for X1Y1Z1_5k
DENS=3000
YOUNG=50e9
FRICT=18
ALPHA=0.3
TENS=45e5
COH=45e6
iterMAX=400000
###################### Import of the sphere assembly
### material definition
sphereMat = FrictMat()
sphereMat.young=YOUNG
sphereMat.frictionAngle=radians(FRICT)
sphereMat.density=DENS

wallMat = FrictMat()
wallMat.young=YOUNG
wallMat.frictionAngle=radians(0)
wallMat.density=DENS

O.materials.append(sphereMat)
O.materials.append(wallMat)

## copy spheres from the packing into the scene
O.bodies.append(ymport.textExt('/tmp/compressed.txt',format='x_y_z_r',material=sphereMat))
## preprocessing to get dimensions
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
## initial surface
S0=X*Z
## spheres factory
R=0
Rmax=0
numSpheres=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
  o.shape.color=(0.7,0.5,0.3)
  numSpheres+=1
  R+=o.shape.radius
  if o.shape.radius>Rmax:
   Rmax=o.shape.radius
  Rmean=R/numSpheres
###################### creation of shear box
thickness=Y/100
oversizeFactor=1.05
O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,3*(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
butee=O.bodies[-1]
### add top box
O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,3*(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
O.bodies.append(utils.box(center=(xinf+X/2,3*(Y/4)+yinf,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
O.bodies.append(utils.box(center=(xinf+X/2,3*(Y/4)+yinf,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
### add bottom box
O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
piston1=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
piston2=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2,(Y/4)+yinf,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
O.bodies.append(utils.box(center=(xinf+X/2,(Y/4)+yinf,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
### loading platens,it is a xz plate?
O.bodies.append(utils.box(center=(xinf+X/2.,yinf-thickness+Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness/2,oversizeFactor*Z/2),material=wallMat,fixed=True))
bottomPlate=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2.,ysup+thickness-Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness/2,oversizeFactor*Z/2),material=wallMat,fixed=True))
topPlate=O.bodies[-1]

###################### Engines
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
  ),
 GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),timestepSafetyCoefficient=dtCOEFF),
 TranslationEngine(ids=[topPlate.id],translationAxis=(0.,-1.,0.),velocity=0.,label='yTranslation'),
 PyRunner(iterPeriod=1,initRun=True,command='servoController()',label='servo'),
 NewtonIntegrator(damping=DAMPING,gravity=(0,0,0),label='damper'),
 PyRunner(iterPeriod=iterMAX/1000,initRun=True,command='dataCollector()'),
]
###################### Engines definition ( servoController() and dataCollector() )
shearing=False
sigmaN=0
tau=0
Fs1=0
Fs2=0
Xdispl=0
px0=0
Ydispl=0
py0=topPlate.state.pos[1]
prevTranslation=0
n=0

def servoController():
 global px0, py0, sigmaN, n, Fn1, Fn2, shearing, butee, piston1, piston2
 Fn1=abs(O.forces.f(topPlate.id)[1])
 Fn2=abs(O.forces.f(bottomPlate.id)[1])
 sigmaN=Fn1/(S0-2*Xdispl*Z) # necessary?
 if shearing==False:
  if yTranslation.velocity<normalVEL:
   yTranslation.velocity+=normalVEL/1000
  if sigmaN>(0.9*normalSTRESS):
   yTranslation.velocity=normalVEL*((normalSTRESS-sigmaN)/normalSTRESS)
 if shearing==False and abs((normalSTRESS-sigmaN)/normalSTRESS)<0.001 and unbalancedForce()<0.01:
  yTranslation.velocity=0
  n+=1
 if n>1000 and abs((sigmaN-normalSTRESS)/normalSTRESS)<0.001:
  print 'stress on joint plane = ', utils.forcesOnPlane((X/2,Y/2,Z/2),(0,1,0))/S0

### add engine and initialisation
O.engines=O.engines[:6]+ [TranslationEngine(ids=[piston1.id,piston2.id],translationAxis=(-1.,0.,0.),velocity=0.)]+O.engines[6:]
px0=piston1.state.pos[0]
py0=topPlate.state.pos[1]
shearing=True
print 'shearing! || iteration=', O.iter
if shearing==True:
 #yTranslation.velocity=0. # if you want a constant normal diaplacement control
 yTranslation.velocity=50*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
 if O.engines[6].velocity<shearVEL:
  O.engines[6].velocity+=(shearVEL/1000)

def dataCollector():
 global Xdispl, Ydispl, tau, Fs1, Fs2
 if shearing:
  Fs1=abs(O.forces.f(butee.id)[0])
  Fs2=abs(O.forces.f(piston1.id)[0])
  Xdispl=abs(piston1.state.pos[0]-px0)
 tau=Fs1/(S0-2*Xdispl*Z)
 Ydispl=abs(topPlate.state.pos[1]-py0)
 numberOfCohesiveBonds=0
 for i in O.interactions:
  if i.phys.isCohesive:
   numberOfCohesiveBonds+=1
 yade.plot.addData({'t':O.time,'i':O.iter,'Xdispl':Xdispl,'sigmaN':sigmaN,'Fn1':Fn1,'Fn2':Fn2,'tau':tau,'Fs1':Fs1,'Fs2':Fs2,'Ydispl':Ydispl,'unbF':utils.unbalancedForce(),'Ek':utils.kineticEnergy(),'bonds':numberOfCohesiveBonds})
 plot.saveGnuplot(Output)

# defines window plot
plot.plots={'i':('sigmaN','tau')}

################# graphical intervace
from yade import qt
qt.Controller()
qt.View()
v=qt.Renderer()
v.dispScale=(1,1,1) # displacements are scaled for visualization purpose

################# to manage interaction detection factor during the first timestep
O.step();
################# initializes the interaction detection factor to its default value (new contacts will be created between strictly contacting particles only)
ss2d3dg.interactionDetectionFactor=-1.
aabb.aabbEnlargeFactor=-1.

################# simulation starts here
O.run(iterMAX)

Question information

Language:
English Edit question
Status:
Needs information
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-09-10
Last reply:
2019-09-10
Jan Stránský (honzik) said : #1

Hello,

> the system can't find the particles.

what are the symptoms (errors etc.)?

Try to create a MWE [1]
M = minimal (focusing as much as possible on the problematic piece of code, leaving many irrelevant lines out, not discouraging us to try your script).
maybe for your case:
###
from yade import ymport
O.bodies.append(ymport.textExt('/tmp/compressed.txt',format='x_y_z_r'))
###
?

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Xue (1q12) said : #2

Hello, Jan!
Thank you very much for your help!
The problem I encountered was that if I opened Yade in the terminal and then entered the command line by line, after typing O.bodies.append (Symort.textExt ('/tmp/compressed.txt', format='x_y_z_r'), I would be able to generate pre-prepared particles. Besides that, i have confirmed that the positions of boxes are correct. But if I type: yade-j 12 xx.py (xx.py is my script file) directly into the terminal, It will report an error which means that no particles can be found. (For some reason, please forgive me for not presenting the error information to you right now.)
I suspect that the operation of importing particles is correct, but there was an error in the simulation of the direct shear process, so I presented all my code.

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

please read [1] and provide all necessary information.
If you get an error, please paste the complete error message

> the system can't find the particles.

what are the symptoms (errors etc.)?
what is "the system"?
what is "particles"? the saved file? data in the file? ... (the error message could have already answer it)

> it will report an error

see above.

thanks
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Can you help with this problem?

Provide an answer of your own, or ask Xue for more information if necessary.

To post a message you must log in.