some question in triaxial test modeling

Asked by Xue on 2019-06-24

Hello everyone:
  the following is my script:

******************************
from yade import pack,utils,plot

from yade import utils

matP = FrictMat()
matP.density=2500#kg/m^3
matP.young = 3e9
matP.poisson = 0.21
matP.frictionAngle = 0.65

matW = FrictMat()
matW.density = 1950#kg/m^3
matW.young = 1.2e9
matW.poisson = 0.33
matW.frictionAngle = 0.23

O.materials.append(matP)
O.materials.append(matW)

#cylinder in triaxial
#cyl=yade.geom.facetCylinder((0,0,0),radius=15,height=60,orientation=Quaternion((0,0,1),0),wallMask=7) # test this code singlely
#O.bodies.append(cyl)
#O.bodies.erase()

#wall1=utils.wall(0,axis=2,sense=1,material=matW)
#wall2=utils.wall(0,axis=2,sense=1,material=matW)
#O.bodies.append(wall1)

#pred=pack.inCylinder((0,0,0),(0,0,60),15)#theunitoflengthiscm
#sp=SpherePack()
#sp=pack.randomDensePack(pred,radius=6,rRelFuzz=3)
#spheres=sp.toSimulation(color=(0,1,1),material=matP)
#former context is concern about how to build the model in cylinder,have some question.

wall1=utils.wall(0,axis=2,sense=1,material=matW)
wall2=utils.wall(0.15,axis=0,sense=-1,material=matW)
wall3=utils.wall(-0.15,axis=0,sense=1,material=matW)
wall4=utils.wall(-0.15,axis=1,sense=1,material=matW)
wall5=utils.wall(0.15,axis=1,sense=-1,material=matW)

O.bodies.append((wall1,wall2,wall3,wall4,wall5))

#pre=pack.inAlignedBox((-0.15,-0.15,0),(0.15,0.15,0.60));
#spheres=pack.randomDensePack(pre,radius=6,rRelFuzz=2.5,material=matP); # in chapter 3, packing algorithms
#O.bodies.append([sphere(c, r) for c, r in spheres]); #id=-1 is the last particle in box.

sp=pack.SpherePack()
sp.makeCloud((-0.15,-0.15,0),(0.15,0.15,0.60),rMean=0.02,rRelFuzz=.006,periodic=True)
sp.toSimulation()

#pre=pack.inAlignedBox((-15,-15,0),(15,15,60));
#spheres=pack.randomDensePack(pre,radius=6,rRelFuzz=2.5,material=matP); # in chapter 3, packing algorithms
#O.bodies.append([sphere(c, r) for c, r in spheres]); #id=-1 is the last particle in box.

sp=pack.SpherePack()
sp.makeCloud((-0.15,-0.15,0),(0.15,0.15,0.60),rMean=0.02,rRelFuzz=.006,periodic=True)
sp.toSimulation()

lpm=O.bodies[-1].id
O.bodies[lpm].state.blockerDOFs='xyzXYZ'
global inipos# global var, attention
inipos=O.bodies[lpm].state.pos[2]
O.bodies.erase(lpm)
wall6= utils.wall(inipos-0.02,axis=2,sense=-1,material=matW)
O.bodies.append(wall6);
lpm=O.bodies[-1].id;
inipos=O.bodies[lpm].state.pos[2];

from yade import timing
yade.timing.reset() #Zero all timing data

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(damping=0.3, label='newton'),
 PyRunner(command='Step()',iterPeriod=1,label='step'),
 PyRunner(command='Save()',iterPeriod=50)
]

def Step():
 from yade import utils, plot
 global stepnum, inipos, deltpos, n
 lpm=O.bodies[-1].id
 platepos=O.bodies[lpm].state.pos[2]
 plateforce=O.forces.f(lpm)[2]
 O.bodies[lpm].state.blockedDOFs='xyzXYZ'
 O.bodies[lpm].state.vel=(0,0,-0.005)
 platepos2=inipos-platepos
 wallx1pos=O.bodies[1].state.pos[0]
 wallx2pos=O.bodies[2].state.pos[0]
 wally1pos=O.bodies[3].state.pos[1]
 wally2pos=O.bodies[4].state.pos[1]
 areaplate=platepos*(wallx1pos-wallx2pos)
 wallx1force=O.forces.f(1)[0]
 wallx2force=O.forces.f(2)[0]
 wally1force=O.forces.f(3)[1]
 wally2force=O.forces.f(4)[1]
 fmw=[abs(wallx1force)+abs(wallx2force)+abs(wally1force)+abs(wally2force)]/(4*areaplate)
 G=(0.5*areaplate)/(matW.young*10*0.0005)
 wav=G*(fmw-500)
 O.bodies[1].state.vel[0]=-wav#可能有错误,正负号
 O.bodies[2].state.vel[0]=wav
 O.bodies[3].state.vel[1]=-wav
 O.bodies[4].state.vel[1]=wav
 delvol=(wallx2pos-wallx1pos)*(wally2pos-wally1pos)*platepos-0.09*inipos#体积变化量
 plot.addData(Pf=plateforce,Pp=platepos2,Pps=platepos2,Pv=delvol,Cof=fmw,**O.energy)
 plot.saveDataTxt('Triaxial-ballast.txt.bz2',vars=('Pp','Pf','Pv','Cof'))
 if stepnum==1:
  O.bodies[lpm].state.vel=(0,0,-0.005)
  if plateforce< 4500: return
  stepnum=stepnum+1
 elif stepnum==2:
  O.bodies[lpm].state.blockedDOFs='xyzXYZ'
  delconf=abs(fwm-500)
  if abs > 100: return
  export.Spheres('/tmp/Triaxial-pre.dat')
  print "wall6pos=%.5f"%(platepos)
  O.pause()

def Save():
 global n
 O.save('Compression'+str(n)+'.yade.gz')
 n=n+1

#plot.plots={'Pp':('Pf',),'Pps':('Pv',)}
plot.plots={'Pp':('Pf',)}
plot.plot()
#plot.saveDataTxt('Compression.txt.bz2',vars('Pp','Pf'))
from yade import qt
O.dt = 0.0005
V = qt.View()
V.screenSize = (550,450)
V.sceneRadius = 20
V.eyePosition = (0.7,0.5,0.1)
V.upVector = (0,0,1)
V.lookAt = (0.15,0.15,0.1)
O.run()
******************************************************

However, when i tested it in yadedaily it showed some problems in below:
   FATAL /data/trunk/pkg/common/InsertionSortCollider.cpp:495 spatialOverlapPeri: Body #0 spans over half of the cell size 0.3 (axis=0, see flag allowBiggerThanPeriod)
FATAL /data/trunk/core/ThreadRunner.cpp:30 run: Exception occured:
/data/trunk/pkg/common/InsertionSortCollider.cpp: Body larger than half of the cell size encountered.

 Besides that, i found that the walls in lateral are different from what i designed in line 35-41 after i added particles in the model.
 So, i would be greatful if you can help me find out the causes of these problems.
 Thanks a lot!

 xue

Question information

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

Hello,
just a comment, do you use both walls and periodicity by intention?
cheers
Jan

Xue (1q12) said : #2

Dear Jan,

First of all, thank you for your reply.

Based on your suggestion, I checked the code and found that there was an error in the setting of Wall. Unfortunately, after my correction, the model still did not achieve the expected results. Specifically, the stress-strain curve I got was wrong. So I suspect that if there's a problem with the 'def step'command?

Best regards

xue

Launchpad Janitor (janitor) said : #3

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