yade-users team mailing list archive

[Clinton Schramm <question294531@xxxxxxxxxxxxxxxxxxxxx>]

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

Clinton Schramm gave more information on the question:
****Ammended Code


readParamsFromTable(rMean=.02,rRelFuzz=.003,maxLoad=1e5,minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *
numspherespack1=1000
numspherespack2 = 173
Zextentbox = 0.200 #was .200
width = 0.154
length = 0.441
cylinderheight = 0.130
cylinderradius = 0.05
nosides = 10

#Set materials for spheres, and box
O.materials.append(FrictMat(young=78.71e6,poisson=0.27,frictionAngle=.5934,density=2600,label='spheres')) #Need to find youngs modulus
O.materials.append(FrictMat(young=30e9,poisson=0.2,frictionAngle=0,label='walls'))

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((0,0,0),(width/2,length/2,Zextentbox),wallMask=31,material='walls'))   #(centre),(x,y,z) Old Box with correct dimensions   
TBM=O.bodies.append(geom.facetCylinder(center=(-0.077,.1555,-.100),radius=cylinderradius,height=cylinderheight,orientation=Quaternion((1,0,0),3.1415926/2),segmentsNumber=nosides,wallMask=7,angleRange=None,closeGap=False,wire=False,highlight=False))

sp=pack.SpherePack()
sp.makeCloud((-0.077,0.180,-0.200),(0.077,-0.200,-0.155),rMean=rMean,rRelFuzz=rRelFuzz,num=numspherespack2) #Cloud for under tunnel
#May need to make cloud larger tos full up the model
sp.makeCloud((-0.077,0.180,-0.05),(0.077,-0.180,.55),rMean=rMean,rRelFuzz=rRelFuzz,num=numspherespack1) #Move cloud in +z 400
sp.toSimulation(material='spheres',color=(0,1,1)) 

O.engines=[
   ForceResetter(),
   # sphere, facet, wall
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),	
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
   SnapshotEngine(iterPeriod=300,fileBase='/home/clinton/Desktop/Video/1/Centri-',label='cent'),
   # the label creates an automatic variable referring to this engines
   PyRunner(command='checkUnbalanced1()',iterPeriod = 200,label='checker'),
   #PyRunner(command='Changecolor()',iterPeriod = 200),
]
O.dt=.5*PWaveTimeStep()

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time

# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced1():	
   # 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<5000: return 
   # the rest will be run only if unbalanced is < .1 (stabilized packing)
   if 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(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1,material='walls'))
   global plate        # without this line, the plate variable would only exist inside this function
   plate=O.bodies[-1]  # the last particles is the plate
   # Wall objects are "fixed" by default, i.e. not subject to forces
   # prescribing a velocity will therefore make it move at constant velocity (downwards)
   plate.state.vel=(0,0,-.7)
   # 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='unloadPlate1()'

def unloadPlate1():
   # if the force on plate exceeds maximum load, start unloading
   if abs(O.forces.f(plate.id)[2])>maxLoad :
      plate.state.vel*=-3
      # next time, do not call this function anymore, but the next one (stopUnloading) instead
      checker.command='stopUnloading()'
      
def stopUnloading():

   if abs(O.forces.f(plate.id)[2])<minLoad:
	plate.state.vel*=0
	areacylinder = (.5*nosides*cylinderradius**2*(sin(degrees(360)/10)))/2
	Volcylinder = areacylinder*cylinderheight
       	height = max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])+Zextentbox
	print 'Height of heighest sphere is: ',height
	v1=Volcylinder+sum([4./3.*pi*(b.shape.radius)**3. for b in O.bodies if isinstance(b.shape,Sphere)])
	v2=width*length*height
	Voidratio=(v2-v1)/v2*100.
	print'The Void Ratio:',Voidratio,'%','This is the second one'
	
    	checker.command='checkUnbalancedtunnelmovement()'
	
def checkUnbalancedtunnelmovement():
	O.engines=O.engines+[TranslationEngine(translationAxis=[0,1,0],velocity=0.3,ids=TBM)]
	# start plotting the data now, it was not interesting before
   	#O.engines=O.engines+[PyRunner(command='addPlotDataForceTunnel()',iterPeriod=200)]
	# we use it below to changefor s in O.bodies if isinstance(s.shape,Sphere):
	global CylinderforceX,  CylinderforceY, CylinderforceZ
	O.engines=O.engines+[PyRunner(command='Changecolor()',iterPeriod=200)]
	
   	CylinderforceX = utils.sumForces([10,14,16,22,26,30,34,36,42,46],(1,0,0))
	CylinderforceY = utils.sumForces([10,14,16,22,26,30,34,36,42,46],(0,1,0))
	CylinderforceZ = utils.sumForces([10,14,16,22,26,30,34,36,42,46],(0,0,1))
 	
	Cylindernormalforce = sumFacetNormalForces([10,14,16,22,26,30,34,36,42,46],1)
	#print 'The sum of normal forces are: ',Cylindernormalforce
   
	plot.addData(CylinderforceX=CylinderforceX,tx=(O.bodies[20].state.pos[1]-O.bodies[20].state.refPos[1]),CylinderforceY=CylinderforceY,ty=(O.bodies[20].state.pos[1]-O.bodies[20].state.refPos[1]),CylinderforceZ=CylinderforceZ,tz=(O.bodies[20].state.pos[1]-O.bodies[20].state.refPos[1]),Cylindernormalforce=Cylindernormalforce,t=(O.bodies[20].state.pos[1]-O.bodies[20].state.refPos[1]),i=O.iter)
   	# start plotting the data now, it was not interesting before
   	#O.engines=O.engines+[PyRunner(command='addPlotDataForceTunnel()',iterPeriod=200)]
   	# next time, do not call this function anymore, but the next one (unloadPlate) instead
   	#O.pause()

#plot.saveDataTxt(O.tags['d.id']+'.txt')

def Changecolor():   
   for s in O.bodies:
	s.shape.color=scalarOnColorScale(s.state.displ().norm(),0,2)
	#print 0.1*s.state.displ()
		# for changing color based on the movement of sheres. its attributes from the functions called
 	#print 'Change'  


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)
   
#Plotting the force on the tunnel face
def addPlotDataForceTunnel():
  print 'test => ',CylinderforceX
  plot.addData(CylinderforceX,tx=(O.bodies[20].state.pos[1]-O.bodies[20].state.refPos[1]),i=O.iter)


# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots={'tx':('CylinderforceX',),'ty':('CylinderforceY',),'tz':('CylinderforceZ',),'t':('Cylindernormalforce',)}
#plot.plots={'i':('unbalanced',),'w':('Fz',),'tx':('Facex',)}

plot.plot()


#O.saveTmp()
from yade import qt
v=qt.View()
v.eyePosition=(-1,0.3,0.25);v.upVector=(0,0,1); v.lookAt=(0,0,0); v.axes=False; v.sceneRadius=1.9

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