New question #683776 on Yade:
https://answers.launchpad.net/yade/+question/683776
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)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
