one dimensional element that interacts with spheres

Asked by alessandro solari on 2021-02-05

Goodmorning,

i need help with yade DEM simulation.
i need to know how to insert a one dimensional element (like a cylinder) that interacts with the spheres (that simulated the soil).
i have the script of direct shear test and i would like to know how and where to add this one dimensional element in the script.
 i write the script below :

###

from yade import pack,plot,export
import math

sp=pack.SpherePack()
O.periodic=True

# dimensions of sample (fixed by particle size such as L/D~15)
RADIUS=0.001
length=0.06
height=length/3
width=length
thickness=RADIUS

# friction angles
compFRIC=40. # during compaction (controls porosity)
FRIC=40. # during shear

# boundary conditions
PI=25e2 # sample preparation: pressure applied for the isotropic compaction
SN=25e3 # normal stress
RATE=0.1 # shear velocity (top plate)

# simulation control
DAMPSHEAR=0.3
ITER=4e4
VTK=20
OUT='PTx'

####

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=1620,young=230e6,poisson=0.5,frictionAngle=radians(compFRIC),normalCohesion=0,shearCohesion=0,label='sphereMat'))

upBox = utils.box(center=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])

sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.001,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

#print 'volRatio=',volRatio

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,PeriTriaxController(dynCell=True,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(damping=0.3,label='newton')
 ,PyRunner(command='dataRecorder()',iterPeriod=100,label='recData',dead=False)
 ,VTKRecorder(fileName=OUT+'.',iterPeriod=1,skipNondynamic=1,recorders=['spheres','colors','velocity','bstresses','stress','boxes','intr','boxes'],label='saveSolid',dead=True)
]

def dataRecorder():
 h=vol=vol_s=nb_s=0.
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 for o in O.bodies:
   if isinstance(o.shape,Sphere) and o.shape.color[0]!=1:
    nb_s+=1
    vol_s += 4.*pi/3.*(o.shape.radius)**3
 n = 1-vol_s/vol
 nbFrictCont=0.
 for i in O.interactions:
  if i.isReal and i.phys.cohesionBroken:
   nbFrictCont+=1
 plot.addData(
  iter=O.iter
  ,contactStress01=((contactStress[0,1])/1000),contactStress11=(contactStress[1,1])
  ,xW=(O.bodies[0].state.pos[0]*1000)
  ,height=h
  ,volume=vol
  ,porosity=n
 )

def triaxDone():
 global phase
 volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 vol_s=Rmean=Rmax=nbSph=0
 Rmin=1e6
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
   if o.shape.radius>Rmax: Rmax=o.shape.radius
   if o.shape.radius<Rmin: Rmin=o.shape.radius
   vol_s += 4.*pi/3.*(o.shape.radius)**3
 Rmean=Rmean/nbSph
 n = 1-vol_s/vol
 print ('DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n)
 print ('Changing contact properties now')
 tt=TriaxialCompressionEngine()
 tt.setContactProperties(FRIC)
 triax.dead=True
 O.pause()

#### Initialization
print ('SAMPLE PREPARATION!')

recData.dead=True # uncomment to record what is happening during stress initialization
O.run(600000,1)
saveSolid.dead=False
#saveSolid.fileName=OUT+'_isoConfined.'
O.step()
saveSolid.dead=True
O.save(OUT+'_initialState.yade')

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1])

#### Applying normal stress
print ('NORMAL LOADING! iter=',O.iter)

recData.dead=True
stage=0
stiff=fnPlaten=currentSN=0.
def servo():
 global stage,stiff,fnPlaten,currentSN
 if stage==0:
  currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  unbF=unbalancedForce()
  #print 'SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF
  boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   print ('Normal stress =',currentSN,' | unbF=',unbF)
   ## the following computes the stiffness of the plate (used for stress control of the top plate)
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   print ('DONE! iter=',O.iter)
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  #boundaryVel=copysign(abs(0.5*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt),O.forces.f(0)[1]-fnDesired)
  boundaryVel=copysign(min(RATE,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
  O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]

O.run(35000,1)
print ('STABILIZING! iter=',O.iter)
O.run(10000,1)
# coloring particles to see vertical stripes in material
dxi=4*(2*RADIUS) # can be a function of cell length dxi=O.cell.hSize[0,0]/5.
n=int(length/dxi)
xmin=1e6
for i in range(0,n):
 for o in O.bodies:
  if o.id>1:
   if o.state.pos[0]<xmin: xmin=o.state.pos[0]
   if (o.state.pos[0]>=i*dxi) and (o.state.pos[0]<((i+0.5)*dxi)):
    o.shape.color[1]=1
for o in O.bodies:
 if (o.state.pos[0]>(xmin+0.5*O.cell.hSize[0,0]-RADIUS)) and (o.state.pos[0]<(xmin+0.5*O.cell.hSize[0,0]+3*RADIUS)):
  o.shape.color[2]=0

saveSolid.dead=True
#saveSolid.fileName=OUT+'_normalLoaded.'
O.step()
saveSolid.dead=True
O.save(OUT+'_normallyLoaded.yade')
if recData.dead==False: plot.saveDataTxt(OUT)

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1])

#### preparing for shearing
Gl1_Sphere.stripes=1
print ('Gluing spheres to boundary walls')
gluedSpheres=[]

## gluing particles in contact with the walls
for i in O.interactions:
  if i.isReal:
   if isinstance(O.bodies[i.id1].shape,Box):
    O.bodies[i.id2].shape.color[0]=1
    gluedSpheres += [O.bodies[i.id2]]
   if isinstance(O.bodies[i.id2].shape,Box):
    O.bodies[i.id1].shape.color[0]=1
    gluedSpheres += [O.bodies[i.id1]]

for i in O.interactions:
  if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and O.bodies[i.id2].shape.color[0]==1 ):
   O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   i.phys.initCohesion=True

print ('nb of glued spheres=',len(gluedSpheres))

#### shearing
print ('SHEARING! iter=',O.iter)

for i in O.interactions:
  if i.isReal:
   if isinstance(O.bodies[i.id1].shape,Sphere):
    O.bodies[i.id1].state.blockedDOFs='XYZ'
    O.bodies[i.id1].state.angVel=(0,0,0)

saveSolid.dead=True
O.step()
saveSolid.dead=True
O.save(OUT+'_shearInit.yade')

recData.dead=False
newton.damping=DAMPSHEAR
saveSolid.dead=True
saveSolid.iterPeriod=int(ITER/VTK)
shearVel=0
iterShear=O.iter
while 1:
 O.run(int(10),1)
 if shearVel<RATE:
  shearVel+=(RATE/100.)
 # only top wall moves
 O.bodies[0].state.vel[0]=shearVel
 ## top and bottom walls move
 #O.bodies[0].state.vel[0]=0.5*shearVel
 #O.bodies[1].state.vel[0]=-0.5*shearVel
 if ( O.iter >= (int(iterShear+ITER)) ):
  print ('iter=',O.iter,' -> FINISHED!')
  plot.saveDataTxt(OUT)
  O.save(OUT+'_sheared.yade')
  sys.exit(0)

thanks for the answer and sorry for my english.

Alessandro

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2021-02-05
Last reply:
2021-02-08
Karol Brzezinski (kbrzezinski) said : #1

Hi Alessandro,

I am not sure what you mean by one dimensional element, but I guess that you are interested in "flat" elements such as walls and facets. If so, please look at the examples in the tutorial [1]. There are some examples of spheres interacting with box made of facets.
If you are specifically interested in building cylinder from facets, you can use this function [2].
If you want someone to look into your script, please try to simplify it as much as possible (see description of MWE - minimum working example [3]). Sometimes, preparing MWE allows you to find solution before even asking the question ;)

Cheers,
Karol

[1] https://yade-dem.org/doc/tutorial-examples.html

[2] https://yade-dem.org/doc/yade.geom.html?highlight=facetcylinder#yade.geom.facetCylinder

[3] https://www.yade-dem.org/wiki/Howtoask

Can you help with this problem?

Provide an answer of your own, or ask alessandro solari for more information if necessary.

To post a message you must log in.