yade-users team mailing list archive

[Weimin Song <question376480@xxxxxxxxxxxxxxxxxxxxx>]

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

Weimin Song gave more information on the question:
Code above has bugs. Below is the modified code:

readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=2e8,minLoad=10)


from yade.params.table import *

from yade import pack, plot

aggregate=O.materials.append(FrictMat(density=3e3, young=30e9,poisson=.3,frictionAngle=.5,label="aggregate"))
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))


sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation(material=aggregate)
#sp.toSimulation()

wallMat=O.materials.append(FrictMat(density=4e3, young=500e9, poisson=0.1, frictionAngle=0.5, label="wallMat"))
wall=O.bodies.append(wall((.5,.5,1),axis=2,sense=-1, color=(0,0,1), material=wallMat))
     
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.5*PWaveTimeStep()

def checkUnbalanced():
   if O.iter<10000: return 
   if unbalancedForce()>0.2: return 
#   O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1, material=wallMat))
   global plate        # without this line, the plate variable would only exist inside this function
   plate=O.bodies[-1]  # the last particles is the plate
   plate.state.vel=(0,0,-.4)
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=50)]
   checker.command='unloadPlate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>maxLoad:
      plate.state.vel*=-1
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[2])<minLoad:
      plot.saveDataTxt('1.txt')
      O.pause()
      ######remove the facet and append the box
      O.bodies.erase(wall)
      boxMat=O.materials.append(FrictMat(density=4e3, young=500e9, poisson=0.1, frictionAngle=0.5, label="boxMat"))
      O.bodies.append(box((.5,.5,1.05),(0.1,0.2,0.05),orientation=Quaternion((1, 0, 0), 0), dynamic=True, fixed=False, wire=False, color=(0.5,0.5,0.5), highlight=False, material=boxMat))
      O.engines=[
         ForceResetter(),
         # sphere, facet, box
         InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
         InteractionLoop(
         # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
            [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack()]
         ),
         NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
         PyRunner(command='checkUnbalanced1()',realPeriod=2,label='checker1'),
      ]
   
      O.dt=.05*PWaveTimeStep() 

      def checkUnbalanced1():
         global box        # without this line, the box variable would only exist inside this function
         box=O.bodies[-1]  # the last particles is the box
         box.state.vel=(0,0,-1)
         O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=50)]
         checker.command='unloadBox()'

      def unloadBox():
         if abs(O.forces.f(box.id)[2])>maxLoad:
            box.state.vel*=-1
            checker.command='stopUnloading1()'
            
            
      def stopUnloading1():
         if abs(O.forces.f(plate.id)[2])<minLoad:
            plot.saveDataTxt('2.txt')
            O.pause()
   
def addPlotData():
   if not isinstance(O.bodies[-1].shape,Wall):
      plot.addData(); return
   Fz=O.forces.f(plate.id)[2]
   plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots={'i':('unbalanced',),'w':('Fz',)}
plot.plot()


O.run()
waitIfBatch()

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