Re: [Question #683526]: Interaction problem between two bodies - Direct Shear

```Question #683526 on Yade changed:

Status: Needs information => Open

from yade import pack, plot, qt
O.reset()

''''
#Creating the upper box
O.bodies.append(U_box)
#Creating the lower box
O.bodies.append(L_box)

'''

#Creating the upper box
U_left=utils.box((0,.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_left)

U_right=utils.box((.4,.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_right)

U_front=utils.box((.2,.1,1),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_front)

U_back=utils.box((.2,.1,.6),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_back)

#Creating the lower box
L_left=utils.box((0,-.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_left)

L_right=utils.box((.4,-.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_right)

L_front=utils.box((.2,-.1,1),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_front)

L_back=utils.box((.2,-.1,.6),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_back)

L_bottom =utils.box((.2,-0.05,.8),(.2,0,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_bottom)

#Creating the extension for sliding plane so particles don't fall out
#Give extra .01m in Y dir. so that the two facet walls...
#...at the shear interface do not rub against each other
O.bodies.append(U_box_extension1)
O.bodies.append(U_box_extension2)

#Creating the extension for sliding plane so particles don't fall out
O.bodies.append(L_box_extension1)
O.bodies.append(L_box_extension2)

#Material assignment- normalcohesion = Tensile strength, shearCohesion = Shear strength
sample_material=FrictMat(
O.materials.append(sample_material)

wall_material= FrictMat(
O.materials.append(wall_material)

#Speed of the translation engine
Vel=1

#Sphere Pack Creation
sp=pack.SpherePack()
sp.makeCloud((0,-0.05,0.6),(0.4,0.2,1.0),rMean=0.008,rRelFuzz=0.3,porosity=0.5, num = 6000)
sp.toSimulation(material=sample_material)   #already one cloud of particles here

#Collecting the ids for translation engine run
id_list=[]
#for i in L_box: id_list.append(i.id)    #for facet box
id_list.append(L_left.id)
id_list.append(L_right.id)
id_list.append(L_front.id)
id_list.append(L_back.id)
id_list.append(L_bottom.id)
for i in L_box_extension1: id_list.append(i.id)
for i in L_box_extension2: id_list.append(i.id)

#Simulation
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=False,label="cohesiveIp") ],
[Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
),
NewtonIntegrator(gravity=(0,-1,0),damping=0.3),
GlobalStiffnessTimeStepper(active=True,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8), #using global timestep for faster simulation
# (command='checkUnbalanced()',realPeriod=0.5,label='checker'),
PyRunner(command='addplate()',iterPeriod=10,label='runner') #real period is actually not good
]

O.bodies if isinstance(b.shape,Sphere)])

global plate_b
plate_b=wall(min([b.state.pos[1]-(b.shape.radius) for b in O.bodies if isinstance(b.shape,Sphere)]),axis=1,sense=+1)
#O.bodies.append(plate_b)

#Plate is being added to the simulation to consolidate the sample
#if unbalancedForce()<0.5 and O.iter>15000:
global plate
inst_position1=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies if isinstance(b.shape,Sphere)])
plate=utils.box((.2,inst_position1,.8),(.199,0,.199),dynamic=False,fixed=False,wire=False,color=(0,0,1))#,material=wall_material)
O.bodies.append(plate)
plate.state.vel=(0,-2,0)

#Plate rebounding after consolidation upto a maximum load
print('Max load reached for top plate')
plate.state.vel=(0,0.5,0)

#Plate movement stopped
plate.state.vel=(0,0,0)
runner.command='defineCNLcondition()'

#Removing the consolidation plate and adding a new one instead to apply a permanent force on it.
def defineCNLcondition():
print('Time to add the CNL plate...')
if unbalancedForce()>0.4:   return
global mainplate
O.bodies.erase(plate.id)
print('object erase')
inst_position2=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies if isinstance(b.shape,Sphere)])
#print(inst_position2)
mainplate=utils.box((.2,inst_position2*1.2,.8),(.199,0,.199),fixed=False,wire=False,color=(0,1,0))#,material=wall_material)
# mainplate=utils.box((.2,inst_position2*1.2,.8),(.199,0,.199),fixed=False,wire=False,color=(1,0,0),material='wall_mat')
O.bodies.append(mainplate)
mainplate.state.mass=100
print(mainplate)
print('mainplate appended')
#mainplate.state.blockedDOFs = 'xyzXYZ'
# print(mainplate.id)
print(mainplate.dynamic)
O.forces.setPermF(mainplate.id,(0,-13000,0))
print('Starting to shear')
runner.command='shearprocess()'

#Engines required for plotting and horizontal movement
def shearprocess():
#print(O.forces.permF(mainplate.id))
#plot.saveDataTxt('arun.txt.bz2')
O.engines=O.engines+[TranslationEngine(ids=id_list,translationAxis=[1,0,0], velocity=Vel)]

Norm_force=abs(O.forces.f(mainplate.id)[1])
Shear_disp= abs(L_left.state.pos[0]-L_left.state.refPos[0])
Shear_force = abs(O.forces.f(L_left.id)[0])
Shear_stress = abs((O.forces.f(L_left.id)[0])/0.16)
Norm_stress= abs((O.forces.f(mainplate.id)[1])/0.16)

plot.plots={'SD':('SF', 'SS'), 'NS':('SS'), }
plot.plot()

#O.run()
O.saveTmp()

--