# Some Problems about Direct Shear Test

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

################## 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.density=DENS

wallMat = FrictMat()
wallMat.young=YOUNG
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
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]
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))
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))
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

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
plot.saveGnuplot(Output)

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

################# graphical intervace
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:
Expired
For:
Assignee:
No assignee Edit question
Last query:
2019-09-10
2019-09-25
 Jan Stránský (honzik) said on 2019-09-09: #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).
###
O.bodies.append(ymport.textExt('/tmp/compressed.txt',format='x_y_z_r'))
###
?

cheers
Jan

 Xue (1q12) said on 2019-09-10: #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 on 2019-09-10: #3

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

 Launchpad Janitor (janitor) said on 2019-09-25: #4

This question was expired because it remained in the 'Needs information' state without activity for the last 15 days.