How to freeze the particles in the imported dense granular packing

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

Hi all,

This question follows my previous one [1].
After successfully imported a dense packing with a specific range of size, I would like to further isotropically compress the cropped dense packing with additional two boundary walls.
When the isotropic compression is running, I find that all the particles also start to move around in any direction. I don't want this to happen since I want to maintain the initial thickness of the packing and all the particle position. I can imagine the reason behind this is that there are forces stored at the grain contacts whenever it is compacted.
I've tried to use "blockDOFs" to fix the position of all the particles but in this way, the desired state of isotropic compression cannot be reached.
So my question is how do I freeze the position of all the particles without removing their original energy states and still reach the desired state of isotropic compression?

Thanks in advance!

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

Below are my two MWEs for packing generation and imported dense packing.

####################### Packing generation MWE ###########################

from yade import pack,plot,export
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
#####Mair and Hazzard (2007)'s microproperties
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='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

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

O.run(10000000,1)
export.text(OUT+".isoCompression")

#################### Isotropic compression MWE #######################

from yade import pack,plot,export,ymport

sp=pack.SpherePack()

O.periodic=True

a=2.5
b=2.5
c=2.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

# Old packing size for spheres generation
a_old=5
b_old=5
c_old=5
length_old=a_old*(2*RADIUS_dimension)
height_old=b_old*(2*RADIUS_dimension)
width_old=c_old*(2*RADIUS_dimension)
thickness_old=RADIUS_dimension

### Packing ratio
packingRatio = a / a_old

### 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.5*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.5*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 <= packingRatio*length_old and y >= packingRatio*height_old and y <= 2*packingRatio*height_old and z >= 0 and z <= packingRatio*width_old

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

O.bodies.append(packing)

### Fix particle position
#for b in O.bodies:
# if isinstance(b.shape,Sphere):
# b.state.blockedDOFs='xyzXYZ'
# b.state.dynamics = False
# b.state.vel = Vector3(0,0,0)

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=1e-3,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.run(100000000,1)

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

#################### END ###############################

Cheers,
Chien-Cheng

Question information

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

Hello,

> I would like to further isotropically compress the cropped dense packing with additional two boundary walls.
> I want to maintain the initial thickness of the packing and all the particle position.

Either you "maintain the initial thickness of the packing and all the particle position" or "isotropically compress the packing", not both (as you already realized).

Basic dualism in mechanics, either you know/prescribe displacement (rotation,strain,...) or force (moment, stress,...), never both.

> I can imagine the reason behind this is that there are forces stored at the grain contacts whenever it is compacted.

Yes, you saved not-stress-free (overlapping) packing, so after loading, repulsive forces are computed (from the overlaps).

> So my question is how do I freeze the position of all the particles without removing their original energy states and still reach the desired state of isotropic compression?

no way (see above).

What is "their original energy states"?

In the second simulation, I see large gap between spherical particles and walls. Why are there the gaps? They are reason that "all the particles start to move around in any direction", so maybe reducing the gap would make the simulation "nicer"?

cheers
Jan

Chien-Cheng Hung (chiencheng) said : #2

Hi Jan,

Thanks for your reply.

> Basic dualism in mechanics, either you know/prescribe displacement (rotation,strain,...) or force (moment, stress,...), never both.

I see the reason why... thanks for the explanation.

> Yes, you saved not-stress-free (overlapping) packing, so after loading, repulsive forces are computed (from the overlaps).

Is it possible that I can save a stress-free packing? Or by any chance that I could temporarily remove the repulsive forces, do the isotropic compression, and then recover the initial repulsive forces (just a thought)?

> What is "their original energy states"?

I just meant the initial repulsive forces. I will try to be careful with the word I use next time...

> In the second simulation, I see large gap between spherical particles and walls. Why are there the gaps? They are reason that "all the particles start to move around in any direction", so maybe reducing the gap would make the simulation "nicer"?

I've tried that but the simulation does not really get better since the real problem is that the moving particles would always go through the boundary walls. I've tried to increase the stiffness of the walls but still get the same. Do you have any idea about this?

Thanks, again.

Cheers,
Chien-Cheng

Jan Stránský (honzik) said : #3

> Is it possible that I can save a stress-free packing?

yes.
But saving is the least part of the job, first you have to create such packing :-)
One option is to decrease radius of all particles such that the resulting packing has no overlap.
Or move particles somehow (e.g. central affinity) such that particles have no overlap.
Or ...

In any case, you have to somehow modify the not-stress-free packing, but some of the modifications may be considered negligible (I expect overlaps much much smaller than particle size, so changing particle's radii by some factor of overlap may be OK)

> Or by any chance that I could temporarily remove the repulsive forces,

this can be somehow tricked (blocking DOFs, increasing mass, decreasing stiffness...), but you should know why you do such modification.
I expect there are more "natural" solutions to the problem.

> do the isotropic compression, and then recover the initial repulsive forces (just a thought)?

Why do you need to recover the repulsive forces after compression?
The compressed packing has their natural forces and positions, "unrelated" to the original input.

> I've tried that but the simulation does not really get better since the real problem is that the moving particles would always go through the boundary walls. I've tried to increase the stiffness of the walls but still get the same. Do you have any idea about this?

please provide the code you use for this case, otherwise it would be mere guessing..

cheers
Jan

Can you help with this problem?

Provide an answer of your own, or ask Chien-Cheng Hung for more information if necessary.

To post a message you must log in.