yade-users team mailing list archive

[littleFish <question674229@xxxxxxxxxxxxxxxxxxxxx>]

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

Dear all,
I have run a very simple script to prepare a packing with 10,000 particles. But when this process is finished with unbalancedForce lower than 0.01, I find there are more than 2,000 particles out of the total 10,000 particles without any interactions with other particles or walls ('float particles' for which O.bodies[i.id].intrs()==[]). I wonder if this is normal or there is something wrong with my simulation. The codes are attached below and thanks for your consideration.

from yade import pack, qt, plot
import time
import numpy as np

############################################
###   DEFINING VARIABLES AND MATERIALS   ###
############################################

ts = time.time()
pressure = -1e5
r = 0.005 
n = 10000 
size0 = 0.3
mn,mx=Vector3(0,0,0),Vector3(size0,size0,size0)


## create materials for spheres and plates 
O.materials.append(CohFrictMat(density=2650,young=6e8,poisson=.8,frictionAngle=0.5,isCohesive=True,normalCohesion=5e10,shearCohesion=5e10,momentRotationLaw=False,label='spheres'))
O.materials.append(CohFrictMat(density=0,young=8.8e10,poisson=.8,frictionAngle=0.,isCohesive=False,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=r,rRelFuzz=0.4,periodic=False,num=n,seed=1)
sp.toSimulation(material='spheres')

print(len(O.bodies))
############################
###   DEFINING ENGINES   ###
############################

triax=TriaxialStressController(
	## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
	## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
	#maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
	#finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
	thickness = 0,
	## switch stress/strain control using a bitmask. What is a bitmask, huh?!
	## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
	## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
	## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
	## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
	stressMask = 7,
	internalCompaction=False, # If true the confining pressure is generated by growing particles
        goal1 = pressure,
        goal2 = pressure,
        goal3 = pressure,
        max_vel = 1,
)

newton=NewtonIntegrator(gravity=(0,0,0), damping=.1)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),],
		[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionOnNewContacts=False,label='cohesiveIp'),],
		[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),]
	),
	## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	#TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
	newton,
        PyRunner(command='addPlotData()',iterPeriod=1000),
        PyRunner(command='check()',iterPeriod=1000),
        PyRunner(command='savePack()',iterPeriod=1000),
]

def addPlotData():
   
   pp = utils.porosity()  #this is the porosity of the cell. 
   ee = pp / (1-pp)  #this is the void ratio of the 3D cell.  
   a = [i for i in O.bodies if i.intrs()==[]]
   plot.addData(unbalanced=unbalancedForce(),i=O.iter,   
                s11 = -getStress()[0,0],
                s22 = -getStress()[1,1],
                s33 = -getStress()[2,2],
                void = ee,
                Num = len(a)
   )

plot.plots={'i':('unbalanced',),'i ':('s11','s22','s33'),
   ' i ':('void',),' i':('Num')
}
# show the plot
plot.plot()

yade.qt.Controller(), yade.qt.View()

def savePack():
   O.save('packing_100kPa_Tmp.yade.gz')

def check():
  unb=unbalancedForce()
  if unb<0.01 and abs(pressure-triax.meanStress)/-pressure<0.001:
    print('Packing generated!')
    print('time: '+str((time.time()-ts)/60.)+'min')
    O.save('packing_100kPa.yade.gz')
    O.pause()

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