yade-users team mailing list archive

[Philip Ho <question631016@xxxxxxxxxxxxxxxxxxxxx>]

New question #631016 on Yade:
https://answers.launchpad.net/yade/+question/631016

I want to make two-ball clumps from spheres.  Would appreciate some help. Thanks.

#rMean is particles size
utils.readParamsFromTable(rMean=.053,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
from yade.params.table import *
from yade import export,pack, plot

#cylinder
O.materials.append(FrictMat(density=2800,young=3e10,poisson=.4,frictionAngle=radians(22),label='walls'))
O.bodies.append(utils.geom.facetCylinder((.5,.5,.75),0.6,1.5,orientation=Quaternion((1,0,0),0),segmentsNumber=16,wallMask=4,material='walls'))
O.bodies.append(utils.geom.facetCylinder((.5,.5,.75),0.6,1.5,orientation=Quaternion((1,0,0),0),segmentsNumber=16,wallMask=2,material='walls'))

#spheres
O.materials.append(FrictMat(density=2650,young=5e8,poisson=.4,frictionAngle=radians(33),label='spheres'))
sp=pack.SpherePack()
sp.makeCloud((0.15,0.15,0.05),(0.9,0.9,3.6),rMean=rMean,rRelFuzz=.0,seed=1)
sp.toSimulation(material='spheres')

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()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.2),
   # the label creates an automatic variable referring to this engine
   # we use it below to change its attributes from the functions called
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.5*utils.PWaveTimeStep()

def checkUnbalanced():
   # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
   if O.iter<10000: return 
   # the rest will be run only if unbalanced is < .1 (stabilized packing)
   if utils.unbalancedForce()>.1: return 
   # add plate at the position on the top of the packing
   # the maximum finds the z-coordinate of the top of the topmost particle
   #O.bodies.append(utils.wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
   #m=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]
   b = len(O.bodies)
   
   #print "number of bodies",b
   #penetrometer
   O.materials.append(FrictMat(density=2800,young=3e10,poisson=.4,frictionAngle=radians(0),label='shaft'))
   O.bodies.append(utils.geom.facetCylinder((.5,.5,2.5),0.1,2.0,orientation=Quaternion((1,0,0),0),segmentsNumber=6,wallMask=4,material='shaft'))
   O.bodies.append(utils.geom.facetCylinder((.5,.5,2.5),0.1,2.0,orientation=Quaternion((1,0,0),0),segmentsNumber=6,wallMask=1,material='shaft'))
   O.bodies.append(utils.geom.facetCone((.5,.5,1.4),0.1,0.0,0.2,orientation=Quaternion((1,0,0),0),segmentsNumber=6,wallMask=4,material='walls'))
   
   #b = len(O.bodies)
   #print "number of bodies2",b
   #O.pause()
   global plate1,plate2,plate3,plate4,plate5,plate6,plate7,plate8,plate9,plate10,plate11,plate12,plate13,plate14,plate15,plate16,plate17,plate18,plate19,plate20,plate21,plate22,plate23,plate24
   plate1=O.bodies[b]  # the last body is the plate
   plate2=O.bodies[b+1]
   plate3=O.bodies[b+2]
   plate4=O.bodies[b+3]
   plate5=O.bodies[b+4]
   plate6=O.bodies[b+5]
   plate7=O.bodies[b+6]  
   plate8=O.bodies[b+7]
   plate9=O.bodies[b+8]
   plate10=O.bodies[b+9]
   plate11=O.bodies[b+10]
   plate12=O.bodies[b+11]
   plate13=O.bodies[b+12]  
   plate14=O.bodies[b+13]
   plate15=O.bodies[b+14]
   plate16=O.bodies[b+15]
   plate17=O.bodies[b+16]
   plate18=O.bodies[b+17]
   plate19=O.bodies[b+18]  
   plate20=O.bodies[b+19]
   plate21=O.bodies[b+20]
   plate22=O.bodies[b+21]
   plate23=O.bodies[b+22]
   plate24=O.bodies[b+23]
   
   # Wall objects are "fixed" by default, i.e. not subject to forces
   # prescribing a velocity will therefore make it move at constant velocity (downwards)
   plate1.state.vel=(0,0,-.2)
   plate2.state.vel=(0,0,-.2)
   plate3.state.vel=(0,0,-.2)
   plate4.state.vel=(0,0,-.2)
   plate5.state.vel=(0,0,-.2)
   plate6.state.vel=(0,0,-.2)
   plate7.state.vel=(0,0,-.2)
   plate8.state.vel=(0,0,-.2)
   plate9.state.vel=(0,0,-.2)
   plate10.state.vel=(0,0,-.2)
   plate11.state.vel=(0,0,-.2)
   plate12.state.vel=(0,0,-.2)
   plate13.state.vel=(0,0,-.2)
   plate14.state.vel=(0,0,-.2)
   plate15.state.vel=(0,0,-.2)
   plate16.state.vel=(0,0,-.2)
   plate17.state.vel=(0,0,-.2)
   plate18.state.vel=(0,0,-.2)
   plate19.state.vel=(0,0,-.2)
   plate20.state.vel=(0,0,-.2)
   plate21.state.vel=(0,0,-.2)
   plate22.state.vel=(0,0,-.2)
   plate23.state.vel=(0,0,-.2)
   plate24.state.vel=(0,0,-.2)
   # start plotting the data now, it was not interesting before
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   #next time, do not call this function anymore, but the next one (unloadPlate) instead
   #checker.command='unloadPlate()'
   
   checker.command='stopUnloading()'

#def unloadPlate():
   ## if the force on plate exceeds maximum load, start unloading
   #a1 = O.forces.f(plate1.id)[2]
   #a2 = O.forces.f(plate2.id)[2]
   #a3 = O.forces.f(plate3.id)[2]
   #a4 = O.forces.f(plate4.id)[2]
   #a5 = O.forces.f(plate5.id)[2]
   #a6 = O.forces.f(plate6.id)[2]
   #if abs(a1+a2+a3+a4+a5+a6)>maxLoad:
      #plate1.state.vel*=-1
      #plate2.state.vel*=-1
      #plate3.state.vel*=-1
      #plate4.state.vel*=-1
      #plate5.state.vel*=-1
      #plate6.state.vel*=-1
      # next time, do not call this function anymore, but the next one (stopUnloading) instead
      #checker.command='stopUnloading()'

def stopUnloading():
   b1 = O.forces.f(plate1.id)[2]
   b2 = O.forces.f(plate2.id)[2]
   b3 = O.forces.f(plate3.id)[2]
   b4 = O.forces.f(plate4.id)[2]
   b5 = O.forces.f(plate5.id)[2]
   b6 = O.forces.f(plate6.id)[2]
   b7 = O.forces.f(plate7.id)[2]
   b8 = O.forces.f(plate8.id)[2]
   b9 = O.forces.f(plate9.id)[2]
   b10 = O.forces.f(plate10.id)[2]
   b11 = O.forces.f(plate11.id)[2]
   b12 = O.forces.f(plate12.id)[2]
   b13 = O.forces.f(plate13.id)[2]
   b14 = O.forces.f(plate14.id)[2]
   b15 = O.forces.f(plate15.id)[2]
   b16 = O.forces.f(plate16.id)[2]
   b17 = O.forces.f(plate17.id)[2]
   b18 = O.forces.f(plate18.id)[2]
   b19 = O.forces.f(plate19.id)[2]
   b20 = O.forces.f(plate20.id)[2]
   b21 = O.forces.f(plate21.id)[2]
   b22 = O.forces.f(plate22.id)[2]
   b23 = O.forces.f(plate23.id)[2]
   b24 = O.forces.f(plate24.id)[2]
   #if abs(b1+b2+b3+b4+b5+b6)<minLoad:
   #if abs(b1+b2+b3+b4+b5+b6+b7+b8+b9+b10+b11+b12+b13+b14+b15+b16+b17+b18+b19+b20+b21+b22+b23+b24)>maxLoad:
      # O.tags can be used to retrieve unique identifiers of the simulation
      # if running in batch, subsequent simulation would overwrite each other's output files otherwise
      # d (or description) is simulation description (composed of parameter values)
      # while the id is composed of time and process number
      #plot.saveDataTxt(O.tags['d.id']+'.txt')
      #O.pause()
   
   #stopping the cone tip from penetrate the bottom of the cylinder
   #z_check=0.bodies[1214].state.pos[2]-0.2
   #z_check=0.bodies[1215].state.pos[2]-0.2
   #z_check=0.bodies[1216].state.pos[2]-0.2
   #z_check=0.bodies[1217].state.pos[2]-0.2
   #z_check=0.bodies[1218].state.pos[2]-0.2
   #z_check=0.bodies[1219].state.pos[2]-0.2
   #if z_check=0.2
   #O.pause()
   
def addPlotData():
   #if not isinstance(O.bodies[-1].shape,Wall):
      #plot.addData(); return
   Fz1=O.forces.f(plate1.id)[2]
   Fz2=O.forces.f(plate2.id)[2]
   Fz3=O.forces.f(plate3.id)[2]
   Fz4=O.forces.f(plate4.id)[2]
   Fz5=O.forces.f(plate5.id)[2]
   Fz6=O.forces.f(plate6.id)[2]
   Fz7=O.forces.f(plate7.id)[2]
   Fz8=O.forces.f(plate8.id)[2]
   Fz9=O.forces.f(plate9.id)[2]
   Fz10=O.forces.f(plate10.id)[2]
   Fz11=O.forces.f(plate11.id)[2]
   Fz12=O.forces.f(plate12.id)[2]
   Fz13=O.forces.f(plate13.id)[2]
   Fz14=O.forces.f(plate14.id)[2]
   Fz15=O.forces.f(plate15.id)[2]
   Fz16=O.forces.f(plate16.id)[2]
   Fz17=O.forces.f(plate17.id)[2]
   Fz18=O.forces.f(plate18.id)[2]
   Fz19=O.forces.f(plate19.id)[2]
   Fz20=O.forces.f(plate20.id)[2]
   Fz21=O.forces.f(plate21.id)[2]
   Fz22=O.forces.f(plate22.id)[2]
   Fz23=O.forces.f(plate23.id)[2]
   Fz24=O.forces.f(plate24.id)[2]
   Ftot=Fz1+Fz2+Fz3+Fz4+Fz5+Fz6+Fz7+Fz8+Fz9+Fz10+Fz11+Fz12+Fz13+Fz14+Fz15+Fz16+Fz17+Fz18+Fz19+Fz20+Fz21+Fz22+Fz23+Fz24
   #plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=utils.unbalancedForce(),i=O.iter)
   s=plate16.state.pos[2]
   plot.addData(Fz=Ftot,i=O.iter,s=s)
   
   if plate16.state.pos[2]<2.2:
     plot.saveDataTxt('output.txt')
     O.pause()

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

yade.qt.View()
#O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
utils.waitIfBatch()

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