# yade-users team mailing list archive

## [Question #695387]: one dimensional element that interacts with spheres

```New question #695387 on Yade:

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 :

###

import math

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

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

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

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])

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())
,NewtonIntegrator(damping=0.3,label='newton')
]

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
n = 1-vol_s/vol
nbFrictCont=0.
for i in O.interactions:
if i.isReal and i.phys.cohesionBroken:
nbFrictCont+=1
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=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)
O.pause()

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

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

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

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:
o.shape.color[2]=0

O.step()

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)

O.step()

newton.damping=DAMPSHEAR
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)
sys.exit(0)

thanks for the answer and sorry for my english.

Alessandro

--