Imported spheres go through the walls

Asked by Chien-Cheng Hung on 2020-12-17

Dear all,

I find a question that is weird to me. This question is followed by my previous one [1].
I've created and saved a stress-free dense packing with a script (MWE 1) and then I load it with another script (MWE 2) and run isotropic compression.
While I was running the isotropic compression, some of the imported spheres went through the walls. To fix this problem, I've tried to decrease the timestep by setting a lower safety factor (0.8 --> 0.3) and also increase the stiffness of both particles and walls (1e8 --> 1e10). Unfortunately, both methods don't work. I have tried ever higher
So I ran another simulation (MWE 3) to check if this also happens a non-imported loose packing where I generate particles with makeCloud with identical material properties (1e8) and timestep (factor 0.8). The process of isotropic compression is successful as I expected where no particles going through the walls. For this reason, I don't think the particle penetration issue is related to stiffness or timestep.
For the imported dense packing, I imagine it's also a "loose" packing since it's stress-free and there is no overlap between particles.
Under this situation, why would I get particle penetration issue with imported dense packing?
Does someone have any idea about it? I've been struggling with this for a while...
Thanks in advance!

The MWEs are provided below.

[1] https://answers.launchpad.net/yade/+question/694379

MWE 1
##########################################
######## Stress-free packing generation ############
#########################################

from yade import pack,plot,export
import numpy as np
sp=pack.SpherePack()
O.periodic=True
a=5
b=5
c=5
RADIUS_dimension=0.0250
OUT='mypacking'
length=a*(2*RADIUS_dimension)
height=b*(2*RADIUS_dimension)
width=c*(2*RADIUS_dimension)
thickness=RADIUS_dimension

### Particle size distribution
psdSizes=[.0232,.0254,.0276]
psdCumm=[0,0.5,1]

### friction angles
spFRIC=0.5

### boundary conditions
PI=1.e5

### material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)
sp.makeCloud((0,0,0),(length,3*height,width), psdSizes = psdSizes, psdCumm = psdCumm, periodic=True, seed =1)
O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,0),material='sphereMat') for s in sp])

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI,-PI,-PI),globUpdate=1,maxStrainRate=(10,10,10),doneHook='stressFree()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

### To achieve stress-free condition using stress control to unload the packing
def stressFree():
    print("start stress-free unloading")
    O.cell.trsf=Matrix3.Identity
    triax.stressMask=7
    triax.goal=(0,0,0)
    triax.maxStrainRate=(1,1,1)
    triax.doneHook='triaxDone()'
    triax.maxUnbalanced=10

def triaxDone():
 triax.dead=True
 O.pause()

O.run(10000000,1)

### Check total contact stress within the system
print ('Total stress (contacts) = ',getStress(O.cell.hSize[1,1]*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

### Reduce the radius of all particles to further make sure there is no overlap between particles
for b in range(len(O.bodies)):
    O.bodies[b].shape.radius = O.bodies[b].shape.radius * 0.99

O.step()
print ('Total stress (contacts) after size reduction = ',getStress(O.cell.hSize[1,1]*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

export.text(OUT+".isoCompression")

### Save the cell size
a = O.cell.size
np.save('cellSize',a)

MWE 2
#######################################################
######## Isotropic compression for imported dense packing ############
######################################################

from yade import pack,plot,export,ymport
import numpy as np

sp=pack.SpherePack()

O.periodic=True

RADIUS_dimension=0.0250
a = np.load('cellSize.npy')

length= a[0]
height= a[1]/3
width= a[2]
thickness=RADIUS_dimension

# friction angles
wallFRIC=90
spFRIC=0.5

# boundary conditions
PI=1.e5

### Material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(wallFRIC),normalCohesion=-1,shearCohesion=-1,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

upBox = utils.box(center=(length/2,2*height+1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

### Import isotropic-compressed packing and boundary
packing = ymport.text("mypacking.isoCompression",color=(1,1,0),material='sphereMat')

### Filter only desired packing
def sphereWanted(b):
    x,y,z = b.state.pos
    return x >= 0 and x <= length and y >= height and y <= 2*height and z >= 0 and z <= width

packing = [b for b in packing if sphereWanted(b)]

O.bodies.append(packing)

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=10,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 triax.dead=True
 O.pause()

O.step()
print ('Total stress (contacts) before isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

O.run(100000000,1)

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Total stress (contacts) after isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

MWE3
##########################################################
######## Isotropic compression for non-imported loose packing ############
#########################################################

from yade import pack,plot,export,ymport
import numpy as np

sp=pack.SpherePack()

O.periodic=True

RADIUS_dimension=0.0250

a = np.load('cellSize.npy')

length= a[0]
height= a[1]/3
width= a[2]
thickness=RADIUS_dimension

sp1=pack.SpherePack()

### Particle size distribution
psdSizes=[.0232,.0254,.0276]
psdCumm=[0,0.5,1]

# friction angles
wallFRIC=90
spFRIC=0.5

# boundary conditions
PI=1.e5

### Material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(wallFRIC),normalCohesion=-1,shearCohesion=-1,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

upBox = utils.box(center=(length/2,2*height+1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

sp1.makeCloud((0,height,width),(length,2*height,2*width), psdSizes = psdSizes, psdCumm = psdCumm, periodic=True, seed =1)

sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,0),material='sphereMat') for s in sp1])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=10,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 triax.dead=True
 O.pause()

print ('Total stress (contacts) before isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

O.run(100000000,1)

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Total stress (contacts) after isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

#############################################################

Best regards,
Chien-Cheng

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-12-17
Last query:
2020-12-17
Last reply:
2020-12-17
Best Jan Stránský (honzik) said : #1

Hello,

using periodic contact detection, particles can have arbitrary position in space, but for "ordinary" contacts, their positions are wrapped into the cell and it works well, no matter how "far" they are.
However, if you use allowBiggerThanPeriod [1]*, the interactions of the "large" body is handled differently, not always using the periodicity.

* probably, didn't search the source code in detail, but the docs and results suggest this and I remember something similar was discussed before

So you can try 2) or 3) below, optionally also with 1):

1) wrap all particles before saving (not sure how important it is for this case, but I consider it a good habit)
###
for b in O.bodies:
    b.state.pos = O.cell.wrap(b.state.pos)
export.text(OUT+".isoCompression")
###

2) enlarge the boxes to have "safe" widths
Using extents=(10*length,thickness/5.,10*width) in MWE2, all spheres remained within the walls

3) model the walls with multiple (e.g. 3x3 or 4x4?) boxes, not using allowBiggerThanPeriod

I personally would go 1)+2)

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.InsertionSortCollider.allowBiggerThanPeriod

PS: A few not important notes:

> for b in range(len(O.bodies)):
> O.bodies[b].shape.radius = O.bodies[b].shape.radius * 0.99

iterating bodies directly is more natural, more readable and more error prone:
for b in O.bodies:
    b.shape.radius = b.shape.radius * 0.99

Furthermore, you can use *= to reduce code duplication:
b.shape.radius *= 0.99

Chien-Cheng Hung (chiencheng) said : #2

Hi Jan,

Thanks for your detailed explanation, it works!
Appreciate your help!

Cheers,
Chien-Cheng

Chien-Cheng Hung (chiencheng) said : #3

Thanks Jan Stránský, that solved my question.