yade-users team mailing list archive

[Xue lihua <question683776@xxxxxxxxxxxxxxxxxxxxx>]

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.