bug with periodic cell of the order of the particles diameter

Asked by Raphaël Maurin

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

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jérôme Duriez (jduriez) said :
#1

Hi Raphael,

I gave a try launching your script (after correcting the material definition ;-) ) but I still do not get exactly your issue, unfortunately. In case others are in the same case than me, would it be possible for you to isolate a simpler configuration when this issue arises ?

Revision history for this message
Raphaël Maurin (raphael-maurin) said :
#2

Hi Jerome,

I will try to be more explicit on the way to reproduce the bug, and try to put a script that works directly ! =)
Taking the script below, which is a simple gravity deposition with bi-periodic boundary condition + a possibility to break the periodicity with lateral walls when putting the flag lateralWalls equal to 1.
- You run the script exactly as it is. Open the 3D viewer. Put run. After the gravity deposition, it pauses. If you press run again, it will make a simulation increment of 0.2s. When you do it, you will see (at least I see!) that the particles takes less and less space at the bottom of the sample. They actually completely overlap each other for no reason.
- You can redo exactly the same with the flag lateralWalls to 1, and you will see that the problem is gone.
- You can play with the script, and vary the width or the length of the system. It happens only if the width or the length of the periodic cell is of the order of 2 or lower.

I hope it is more clear ! Can you first tell me if you can reproduce the bug ?

Thanks !

Raphael

new script:
from yade import pack, plot
import math

diameterPart = 6e-3 #Diameter of the particles, in meter

length = 10*diameterPart #length of the stream
width = 2*diameterPart #width of the stream
Nlayer = 5 #nb of layers of particle, in diameter

#Option to put lateral walls and break the biperiodicity
lateralWalls = 0

height = 10*Nlayer*diameterPart #heigth of the periodic cell
groundPosition = height/4.0

O.materials.append(ViscElMat(en=0.5, et=0., young=5e6, poisson=0.5, density=2500, frictionAngle=0.4, label='Mat'))
#O.materials.append(FrictMat(young=5e6, poisson=0.5, density=2500, 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
partCloud.makeCloud(minCorner=(0,leftLimitY,groundPosition+diameterPart+1e-4),maxCorner=(length,rightLimitY,groundPosition+depositionHeight),rRelFuzz=0., rMean=diameterPart/2.0, num = partNumber)
partCloud.toSimulation(material='Mat')

#########################
#### 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()',virtPeriod = 0.05,label = 'packing'),
 GlobalStiffnessTimeStepper(defaultDt = 1e-4, viscEl = True,timestepSafetyCoefficient = 0.7, label = 'GSTS'),
 NewtonIntegrator(damping=0.2, gravity= Vector3(0.,0.0,-9.81), label='newton')
 ]
#save the initial configuration to be able to recharge the simulation starting configuration easily
O.saveTmp()
#run
#O.run()

def gravityDeposition():
 if O.time<0.3 : return
 else:
  packing.virtPeriod = 0.2
  O.pause()
 return

Revision history for this message
Jérôme Duriez (jduriez) said :
#3

I gave another try with your new script (re-correcting the material definition ;-) ). With or without lateralWalls, I see:

- after a first "Play" (then automatic pause), the sample stops in both cases in a state that does not really look like in equilibrium: some spheres have very few contacts with others, such that it is hard to imagine they are balanced.

- after a 2d

Revision history for this message
Jérôme Duriez (jduriez) said :
#4

(sorry...)

- after a 2d play (then automatic pause) all the spheres lie on the floor, in a way that does not seem problematic to me. I ll send you some pictures directly.

Anyway, I won't be surprised that Yade uncomplete periodic simulations (with allowBiggerThanPeriod flag) face problems. See https://bugs.launchpad.net/yade/+bug/1112763.
(I can not be sure it is related to your case, though)

Revision history for this message
Raphaël Maurin (raphael-maurin) said :
#5

Hi Jerome !

Thank you for your tests ! On the picture you sent me by email, you can see very few particles at the bottom of the sample (one or two layers at most, which makes ~20 particles), but actually you introduced 120 particles ! And if you look carefully, some particles are completely overlapping...

So it seems you are able to reproduce the bug.
In my case, when I put the flag lateralWalls to 1 (line 11 of the script), it is completely different: you can see a deposit of about 5 layers of particles as expected, and no weird overlapping.

I saw the bug you mentioned and I believe that it is linked. Do you think the new version of the insertion collider you proposed would fix the problem ?

Raphael

Revision history for this message
Jérôme Duriez (jduriez) said :
#6

I do not know..

In my case, the major issue was at the very first step of the simulation. The insertion collider version I proposed solved at least the contact detection for this very first step.
For contact detection during the simulation, the fix was probably less efficient, as I mentionned in the bug report. (different parts of the collider code apply for very first step and next steps, if I remember well)

This is quite old for me now, but all my passed efforts are described on the bug report. It may be worth for you to give a try with what I proposed at that time, and it may give you some ideas to complete the fix, if necessary ?

Anyway, I would be very happy if this is finally useful (hoping that the collider did not change too much in the meantime..)

Revision history for this message
Launchpad Janitor (janitor) said :
#7

This question was expired because it remained in the 'Open' state without activity for the last 15 days.