yade-users team mailing list archive

[Yor1 <question295520@xxxxxxxxxxxxxxxxxxxxx>]

New question #295520 on Yade:
https://answers.launchpad.net/yade/+question/295520

Hello everybody !

I try to simulate a simple shear test.
I simulate the simple shear test with 400000 iterations.
The problem is when the simulation is between 1000 and 1100 iterations i obtained this error message:

"FATAL /home/hamdi/YADE/sources/trunk/core/ThreadRunner.cpp:30 run: Exception occured: 
std::bad_alloc"

I don't undersatand the origin of the error.

Best regards!
Jabrane

P.S this is my script


# -*- coding: utf-8 -*-
O=Omega()

from yade import ymport, utils, plot

################## parameters definition
PACKING='X1Y1Z1_1k'
DAMPING=0.5
dtCOEFF=0.5
normalSTRESS=5e5
normalVEL=normalSTRESS/1e8 # 0.001 for 100kPa // optimized for normalVEL=normalSTRESS/1e8?
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
OUT='shearTest_'+PACKING+'_K10_D3000_Y50e9A03F18T45e5C45e6_500kPa_shearVel1'

###################### Import of the sphere assembly
### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
## copy spheres from the packing into the scene
O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),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.3

### loading platens
O.bodies.append(utils.box(center=(xinf+X/2.,yinf-thickness+Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness,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,oversizeFactor*Z/2),material=wallMat,fixed=True))
topPlate=O.bodies[-1]


O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='ss2d3dg'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),

 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=100,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?
    #print yTranslation.velocity, sigmaN
    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 top box
    O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,3*(ysup/4)+2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
    butee=O.bodies[-1]
    O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,3*(ysup/4)+2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))     
    O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True)) 
    O.bodies.append(utils.box(center=(xsup/2,3*(ysup/4)+2*Rmean,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))
 ### add bottom box
    O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,(ysup/4)-2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
    piston1=O.bodies[-1]
    O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,(ysup/4)-2*Rmean,zsup/2),extents=(thickness,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
    piston2=O.bodies[-1]
    O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))
    O.bodies.append(utils.box(center=(xsup/2,(ysup/4)-2*Rmean,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness),material=wallMat,fixed=True,wire=True))
 ### 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.saveDataTxt(OUT)



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


from yade import qt
qt.Controller()
qt.View()
v=qt.Renderer()
v.dispScale=(1,1,1) 

O.step();

ss2d3dg.interactionDetectionFactor=-1.
aabb.aabbEnlargeFactor=-1.

O.run(iterMAX)


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.