# 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,s,color=(0,0,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies.state.pos-O.bodies.state.pos)*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.state.pos-O.bodies.state.pos
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!=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.state.pos*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.state.pos-O.bodies.state.pos)*O.cell.hSize[0,0]*O.cell.hSize[2,2])
h=O.bodies.state.pos-O.bodies.state.pos
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)/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Normal stress (contacts) = ',getStress((O.bodies.state.pos-O.bodies.state.pos)*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)/(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.state.vel=boundaryVel
if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
stage+=1
fnPlaten=O.forces.f(0)
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.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)-fnDesired)/stiff/O.dt),O.forces.f(0)-fnDesired)
boundaryVel=copysign(min(RATE,abs(0.333*(O.forces.f(0)-fnDesired)/stiff/O.dt)),O.forces.f(0)-fnDesired)
O.bodies.state.vel=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<xmin: xmin=o.state.pos
if (o.state.pos>=i*dxi) and (o.state.pos<((i+0.5)*dxi)):
o.shape.color=1
for o in O.bodies:
if (o.state.pos>(xmin+0.5*O.cell.hSize[0,0]-RADIUS)) and (o.state.pos<(xmin+0.5*O.cell.hSize[0,0]+3*RADIUS)):
o.shape.color=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)/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Normal stress (contacts) = ',getStress((O.bodies.state.pos-O.bodies.state.pos)*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=1
gluedSpheres += [O.bodies[i.id2]]
if isinstance(O.bodies[i.id2].shape,Box):
O.bodies[i.id1].shape.color=1
gluedSpheres += [O.bodies[i.id1]]

for i in O.interactions:
if i.isReal and ( O.bodies[i.id1].shape.color==1 and O.bodies[i.id2].shape.color==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.state.vel=shearVel
## top and bottom walls move
#O.bodies.state.vel=0.5*shearVel
#O.bodies.state.vel=-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 on 2021-02-08: #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 . 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 .
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 ). Sometimes, preparing MWE allows you to find solution before even asking the question ;)

Cheers,
Karol

## 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.