New question #275284 on Yade:
https://answers.launchpad.net/yade/+question/275284
Hi all,
There seems to be a bug in yade when simulating problems with a periodic cell
length of the order of the diameter (typically ~2d): some particles are
completely overlapping each other, and behave as a unique particle.
>From what I tested, the bug is independent from:
- the particles material properties (stiffness, viscous damping, friction
angle, density),
- the contact law (can be also pbserved with
Law2_ScGeom_FrictPhys_CundallStrack with a time step ~10^{-6})
- the time step
- the diameter of the particles
- the number of particle layers
I tested on yadedaily 1.10.0-72-d9ab58c~precise, yade from source git-3e1e44a,
on two different computer running with Ubuntu 14.04.2 LTS.
Below the message, you will find a simplified script that reproduce the
problem. It is just a gravity deposition with or without lateral walls
(lateralWalls = 1 or 0), prescribing the width and length of the periodic cell,
and the number of particle layers that you want to obtain.
The bug arises when the length or width of the periodic cell is inferior to
about 2.2d. Adding lateral walls (lateralWalls = 1 in the script) solves the
problem. For periodic cell length/width of 5d or upper, I never had any
problems.
Are you able to reproduce the bug ? If yes, do you have any idea of why that is
happening ?
Thanks in advance !
Raphael
from yade import pack, plot
import math
##
## Main parameters of the simulation
##
diameterPart = 6e-3 #Diameter of the particles, in meter
densPart = 2500 #Density of the particles, kg/m3
lengthCell = 2 #Streamwise length of the periodic cell, in diameter
widthCell = 10. #Spanwise length of the periodic cell, in diameter
Nlayer = 5 #nb of layers of particle, in diameter
#Option to put lateral walls and break the biperiodicity
lateralWalls = 0
###########################
## PARAMETERS DECLARATION##
###########################
#Geometrical configuration, define useful quantities
height = 10*Nlayer*diameterPart #heigth of the periodic cell
length = lengthCell*diameterPart #length of the stream
width = widthCell*diameterPart #width of the stream
gravityVector = Vector3(0.,0.0,-9.81)
groundPosition = height/4.0
#Particles contact law parameters
maxPressure = densPart*0.6*Nlayer*diameterPart*abs(gravityVector[2]) #Estimated
max particle pressure ("hydrostatic")
normalStiffness = maxPressure*diameterPart*1e4 #Evaluate the minimal normal
stiffness to be in the rigid particle limit (cf Roux and Combe 2002)
youngMod = normalStiffness/diameterPart #Young modulus of the particles from
the stiffness wanted.
#Material
#O.materials.append(ViscElMat(en=0.5, et=0., young=youngMod, poisson=0.5,
density=densPart, frictionAngle=0.4, label='Mat'))
O.materials.append(FrictMat(young=youngMod, poisson=0.5, density=densPart,
frictionAngle=0.4, label='Mat'))
#Definition of the semi-periodic cell
O.periodic = True
O.cell.setBox(length,width,height)
#To be compatible with lateral walls
leftLimitY = 0.0
rightLimitY = width
centerLimitY = width/2.0
#If lateralWalls, add the walls and increase the cell size to avoid particle
touching the walls to appear on the other side
if lateralWalls:
#Warn the user
print '\nlateralWalls option activated: mono-periodic boundary
condition only !\n'
#Increase the cell size
O.cell.setBox(length,width+2*diameterPart,height)
#Modify accordingly the position of the center of the cell and the wall
right and left position
leftLimitY = diameterPart
rightLimitY = width+diameterPart
centerLimitY = diameterPart + width/2.0
#Define the wall and add to the simulation
sidePlaneL = box(center=
(length/2.0,leftLimitY,height/2.0),extents=(2000,0,height*10),fixed=True,wire =
True,color = (1,0,0), material = 'Mat')
sidePlaneR = box(center=
(length/2.0,rightLimitY,height/2.0),extents=(2000,0,height*10.0),fixed=True,wire=True,
material = 'Mat',color = (0,0,1))
O.bodies.append([sidePlaneR,sidePlaneL])
# Ground reference and Wall on the side
lowPlane = box(center=
(length/2.0,centerLimitY,groundPosition),extents=(200,200,0),fixed=True,wire=False,color
= (0.,1.,0.),material = 'Mat') #Build a plane to have a reference for the eyes
O.bodies.append(lowPlane) #add to simulation
#Create a loose cloud of particle inside the cell
partCloud = pack.SpherePack()
partVolume = pi/6.*pow(diameterPart,3.)
partNumber = int(Nlayer*0.6*diameterPart*length*width/partVolume) #Volume of
beads
#Define the deposition height considering that the packing realised by make
cloud is 0.2
depositionHeight = Nlayer*0.6/0.1*diameterPart #Consider that the packing
realised by make cloud is 0.2
#Create a cloud of partNumber particles PHF change to 0.1
partCloud.makeCloud(minCorner=(0,leftLimitY,groundPosition+diameterPart+1e-4),maxCorner=(length,rightLimitY,groundPosition+depositionHeight),rRelFuzz=0.,
rMean=diameterPart/2.0, num = partNumber)
#Send this packing to simulation with material Mat
partCloud.toSimulation(material='Mat')
#Evaluate the deposition time considering the free-fall time of the highest
particle to the ground (overestimation)
depoTime = sqrt(depositionHeight*2/abs(gravityVector[2]))
#########################
#### SIMULATION LOOP#####
#########################
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod
= True),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
# [Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
# [Ip2_FrictMat_FrictMat_FrictPhys()],
# [Law2_ScGeom_FrictPhys_CundallStrack()]
,label = 'interactionLoop'),
PyRunner(command='gravityDeposition(depoTime)',virtPeriod = 0.05,label
= 'packing'),
GlobalStiffnessTimeStepper(defaultDt = 1e-4, viscEl =
True,timestepSafetyCoefficient = 0.7, label = 'GSTS'),
NewtonIntegrator(damping=0.2, gravity=gravityVector, label='newton')
]
#save the initial configuration to be able to recharge the simulation starting
configuration easily
O.saveTmp()
#run
#O.run()
def gravityDeposition(lim):
if O.time<lim : return
else:
packing.virtPeriod = 0.2
O.pause()
return
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
