Interaction problem between two bodies

In an effort to increase the sustainability of this knowledge base, please copy and paste the code here instead of using a link. Thank you :-)

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

''''
#Creating the upper box
U_box=geom.facetBox((.2,.1,.8),(.2,.1,.2),wallMask=51)
O.bodies.append(U_box)
#Creating the lower box
L_box=geom.facetBox((.2,-.1,.8),(.2,.1,.2),wallMask=55)
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
U_box_extension1=geom.facetBox((-.2,.5,.8),(.2,.5,.2),wallMask=4)
O.bodies.append(U_box_extension1)
U_box_extension2=geom.facetBox((.6,.5,.8),(.2,.5,.2),wallMask=4)
O.bodies.append(U_box_extension2)

#Creating the extension for sliding plane so particles don't fall out
L_box_extension1=geom.facetBox((.6,-.5001,.8),(.2,.5,.2),wallMask=8)
O.bodies.append(L_box_extension1)
L_box_extension2=geom.facetBox((-.2,-.5001,.8),(.2,.5,.2),wallMask=8)
O.bodies.append(L_box_extension2)

#Material assignment- normalcohesion = Tensile strength, shearCohesion = Shear strength
sample_material=FrictMat(
    young=50e6 ,poisson=0.25 ,density=1800 ,frictionAngle=radians(45),label='sample_mat')
O.materials.append(sample_material)

wall_material= FrictMat(
    young=4e9 ,poisson=0.25 ,density=5000 ,frictionAngle=radians(50) ,label='wall_mat')
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
]
maxLoad=1e5
minLoad=1e1

inst_position_lower=min([b.state.pos[1]-(b.shape.radius) for b in 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)
print('Plate in -Y added')

#Plate is being added to the simulation to consolidate the sample
def addplate():
    #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)
        print('Plate added')
        plate.state.vel=(0,-2,0)
        runner.command='unloadPlate()'

#Plate rebounding after consolidation upto a maximum load
def unloadPlate():
    if abs(O.forces.f(plate.id)[1])>maxLoad:
        print('Max load reached for top plate')
        plate.state.vel=(0,0.5,0)
        runner.command='stopUnloadingPlate()'

#Plate movement stopped
def stopUnloadingPlate():
    if abs(O.forces.f(plate.id)[1])<minLoad:
        print('Unloading 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)]
    O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=10)]

def addPlotData():
    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.addData(NF=Norm_force, SD=Shear_disp, SF=Shear_force, SS=Shear_stress, NS=Norm_stress)

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

#O.run()
O.saveTmp()

Can anyone please tell me where I have gone wrong with the code? The 'dynamic' property is deprecated as listed in the Yade documentation. So I am using the 'Fixed' parameter to specify the body's movement. I have tried several ways but when I apply a permanent force on the main plate, the plate just passes through the spheres without any interaction. I am unsure about why this is happening. Help would be very much appreciated.

Hello,

I have tried your code with Yade 2018.02b and it seems to work as expected (one wall is compressing the sample of spheres). Consider upgrading Yade (1.20.0 is almost 5 years old).

cheers
Jan

Thanks, Jan. I guess that the problem is with the version of Yade. I will upgrade it appropriately.