yade-users team mailing list archive

[Akm <question683526@xxxxxxxxxxxxxxxxxxxxx>]

Question #683526 on Yade changed:
https://answers.launchpad.net/yade/+question/683526

    Status: Needs information => Open

Akm gave more information on the question:
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()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.