contact model for particle compaction

Asked by Ricardo Lorenzoni

Hi!!

I am conducting a study on compaction of particles. In my study, the particles must be independent, it is not interesting to know when a particular particle / structure would break due to the applied pressure.
It is common knowledge that, as the weight is increased over a layer of particles, its porosity reduces, however, this reduction does not happen in a linear way, but in a curve, where the greater the amount of weight added, the less the effect on compression.
However, I have already performed tests with several contact models for spherical particles and in none of the models tested so far, I obtained the expected behavior in compaction, always occurring in a linear way.
I would like to know if anyone would have any tips on contact model and physical model to use that will present the expected behavior.

Thank you for your attention and help.

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

> this reduction ... happen in ... a curve
> I would like to know if anyone would have any tips on contact model ... that will present the expected behavior.

please, provide what is THE expected behavior (specifically what is "a curve") before asking about suitable model

> However, I have already performed tests with several contact models

Please provide more info. What contact models AND with what settings. Sometimes a model must be properly adjusted to have desired response, with different parameters the model might behave seemingly not as it should.

> However, I have already performed tests with ... spherical particles

consider using non-spherical particles

cheers
Jan

Revision history for this message
Ricardo Lorenzoni (ricolorenzoni) said :
#2

Hi. Thanks for your answer.

figure 1 - https://pasteboard.co/JV5qiMz.png
In this first image, you can see the experimental results, where the curve shows that as more weight is added to the sample, the effect on the compacting of the grain mass will decrease.

figure 2 - https://pasteboard.co/JV5mPSf.png
In this second image, the results of the computer simulations are plotted, in which, we can see that the effect of increasing weight on the compaction of the particles is linear.

I have already carried out a series of tests with several contact and physical models, in none of them I was able to replicate the curve behavior, as shown in Figure 1.

Here is an example of the base code I am using for the simulations.
the variable "altura" indicates the height of the grain mass on the sample analyzed and, is set to 10 m.

sample code with cundallStrack contact model:

import sys, math, csv, os
from minieigen import *
from yade import pack, plot, utils, ymport, export

matWall=CohFrictMat(density=2700,young=7e10,poisson=0.334,label="parede")
idparede=O.materials.append(matWall)

r_mean=0.002963
altura=10

lx=0
hx=0.1
ly=0
hy=0.1
lz=0
hz=0.15

p1 = (lx, ly, lz)
p2 = (lx, hy, lz)
p3 = (hx, ly, lz)
p4 = (hx, hy, lz)

p5 = (lx, ly, hz)
p6 = (lx, hy, hz)
p7 = (hx, ly, hz)
p8 = (hx, hy, hz)

p9 = (lx, ly, hz+0.05)
p10= (lx, hy, hz+0.05)
p11= (hx, ly, hz+0.05)
p12= (hx, hy, hz+0.05)

# DEFINICAO DAS FACES

#base
base1=utils.facet([p1,p2,p3],mask=1,wire=True,material=idparede)
base2=utils.facet([p2,p3,p4],mask=1,wire=True,material=idparede)
O.bodies.append(base1)
O.bodies.append(base2)

#Wall1
parede11=utils.facet([p1,p2,p5],mask=1,wire=True,material=idparede)
parede12=utils.facet([p2,p5,p6],mask=1,wire=True,material=idparede)
O.bodies.append(parede11)
O.bodies.append(parede12)

#Wall2
parede21=utils.facet([p1,p3,p5],mask=1,wire=True,material=idparede)
parede22=utils.facet([p3,p5,p7],mask=1,wire=True,material=idparede)
O.bodies.append(parede21)
O.bodies.append(parede22)

#Wall3
parede31=utils.facet([p3,p4,p7],mask=1,wire=True,material=idparede)
parede32=utils.facet([p4,p7,p8],mask=1,wire=True,material=idparede)
O.bodies.append(parede31)
O.bodies.append(parede32)

#Wall4
parede41=utils.facet([p2,p4,p6],mask=1,wire=True,material=idparede)
parede42=utils.facet([p4,p6,p8],mask=1,wire=True,material=idparede)
O.bodies.append(parede41)
O.bodies.append(parede42)

#Wall5
parede51=utils.facet([p5,p6,p9],mask=1,wire=True,material=idparede)
parede52=utils.facet([p6,p9,p10],mask=1,wire=True,material=idparede)
O.bodies.append(parede51)
O.bodies.append(parede52)

#Wall6
parede61=utils.facet([p6,p8,p12],mask=1,wire=True,material=idparede)
parede62=utils.facet([p12,p10,p6],mask=1,wire=True,material=idparede)
O.bodies.append(parede61)
O.bodies.append(parede62)

#Wall7
parede71=utils.facet([p7,p8,p11],mask=1,wire=True,material=idparede)
parede72=utils.facet([p8,p11,p12],mask=1,wire=True,material=idparede)
O.bodies.append(parede71)
O.bodies.append(parede72)

#Wall8
parede81=utils.facet([p5,p7,p9],mask=1,wire=True,material=idparede)
parede82=utils.facet([p7,p9,p11],mask=1,wire=True,material=idparede)
O.bodies.append(parede81)
O.bodies.append(parede82)

particles=CohFrictMat(density=1163,young=2.6e6,poisson=0.25,alphaKr=0.05,frictionAngle=radians(22.6),label="particle")
idparticle=O.materials.append(particles)

O.bodies.append(sphere((0.01,0.01,0.01),radius=r_mean,material=idparticle))

# DEFINICAO DOS MOTORES
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],verletDist=.05*.29e-3),
   InteractionLoop(
      # handle sphere+sphere and facet+sphere collisions
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()] #modelo linear de contato
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.04),
   BoxFactory(maxParticles=-1, extents=(hx/2,hy/2,0.01), center=(hx/2,hy/2,hz+0.03), rMin=r_mean, rMax=r_mean, vMin=100, vMax=100, vAngle=0, massFlowRate=0.01*16.6, normal=(0,0,0), label='factory', materialId=idparticle),
   PyRunner(virtPeriod=0.02, command='verificaAltura()'),
]

global pesomassa
pesomassa=0
def calculaPeso():
 global pesomassa
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   pesomassa+=b.state.mass
 pesometro=(altura*((100*pesomassa)/15))
 adicionaPeso(pesometro)

global executa
executa=True

def verificaAltura():
 global executa
 if executa:
  count=0
  for b in O.bodies:
   if isinstance(b.shape, Sphere):
    if (b.state.vel[2]<0.05 and b.state.pos[2]>0.15 and b.state.pos[2]<0.16):
     count+=1
  if (count>16.6*4):
   executa=False
   O.engines[4].dead=True
   ajusta()

def ajusta():
 bodiesToBeDeleted=[]
 for b in O.bodies:
  if b.state.pos[2]>hz or b.state.pos[1]<ly or b.state.pos[1]>hy or b.state.pos[0]>hx or b.state.pos[0]<lx: # artificial condition for particles deletion,
       bodiesToBeDeleted.append(b)
 for b in bodiesToBeDeleted:
    O.bodies.erase(b.id)
 calculaPeso()

def calculaporosidade():
# array=[]
 if os.path.exists('analise_porosidade_variando_young_por_altura_'+str(altura)+'.csv'):
  conf='a+'
 else:
  conf='w'
 with open('analise_porosidade_variando_young_por_altura_'+str(altura)+'.csv', conf,newline='') as csvfile:
  writer=csv.writer(csvfile)
  if conf=='w':
   writer.writerow(["Porosidade","Tempo","eixo x/y","altura"])
   writer.writerow([str(voxelPorosity(800,((hx/5)*2,(hy/5)*2,0.01),((hx/5)*4,(hy/5)*4,0.04))),str(O.time),str(hx),str(altura)])
  else:
   writer.writerow([str(voxelPorosity(800,((hx/5)*2,(hy/5)*2,0.01),((hx/5)*4,(hy/5)*4,0.04))),str(O.time),str(hx),str(altura)])

global amax
amax=0
def calc_altura():
 global amax
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   if b.state.pos[2]>amax:
    amax=b.state.pos[2]

def adicionaPeso(peso):
 calc_altura()
 peso1=[lx+0.0001,ly+0.0001,amax+0.003]
 peso2=[lx+0.0001,hy-0.0001,amax+0.003]
 peso3=[hx-0.0001,ly+0.0001,amax+0.003]
 peso4=[hx-0.0001,hy-0.0001,amax+0.003]
 facepeso1=yade.utils.facet([peso1,peso2,peso3],dynamic=True,fixed=False )
 facepeso2=yade.utils.facet([peso2,peso3,peso4],dynamic=True,fixed=False)
 facepeso1.state.mass = peso/2 # (!)
 facepeso1.state.inertia = (0,0,0) # (!)
 facepeso2.state.mass = peso/2 # (!)
 facepeso2.state.inertia = (0,0,0) # (!)
 O.bodies.appendClumped([facepeso1,facepeso2])
 O.bodies[len(O.bodies)-1].state.blockedDOFs="xyXYZ"
 O.bodies[len(O.bodies)-2].state.blockedDOFs="xyXYZ"
 O.bodies[len(O.bodies)-3].state.blockedDOFs="xyXYZ"
 momento=O.time+0.01
 O.engines=O.engines+[PyRunner(command='O.pause()',virtPeriod=momento)]
 O.engines=O.engines+[PyRunner(virtPeriod=momento, command='calculaporosidade()')]

O.dt=0.9*utils.PWaveTimeStep()
O.saveTmp()
O.run()
O.wait()

for each physical and contact model, I changed the properties and materials of the particles and walls (when I thought it was necessary) and also, the contact and physical models, as described below.

elPerfPl Model:

wall -> FrictMat(density=2700,young=7e10,poisson=0.334,label="parede")
particles -> FrictMat(density=1163,young=2.6e6,poisson=0.25,frictionAngle=radians(22.6),label="particle")
InteractionLoop(
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
)

HertzMindlin Model:

wall -> ViscElMat(density=2700,young=7e10,poisson=0.334,label="parede")
particles -> ViscElMat(density=1163,young=2.6e6,poisson=0.25,frictionAngle=radians(22.6),label="particle")
InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_MindlinPhys()],
      [Law2_ScGeom_MindlinPhys_Mindlin()] #modelo linear de contato
)

ViscElPhys_Basic Model:

wall -> ViscElMat(density=2700,young=7e10,poisson=0.334,label="parede")
particles -> ViscElMat(density=1163,young=1.0e7,poisson=0.25,frictionAngle=radians(22.6),label="particle")
InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
      [Law2_ScGeom_ViscElPhys_Basic()] #modelo linear de contato
)

About using non-spherical particles like clumps, this option is in my plans but, for preliminarly tests, spherical particles have less computational costs to simulate, and, we hope that spherical particles could be used for the study with a good similarity to the real particles behavior.

Thanks for your help. :)

Revision history for this message
Jan Stránský (honzik) said :
#3

Hello,

the results may be affected by many factors, just some brainstorming.
- effect of walls (friction, stiffness, size w.r.t real sample...)
- particle size distribution (uniform size, real PSD, ...)
- porosity computation (voxel resolution in case of voxelPorosity)
- friction of particles themselves
- ...

Particles in different heights are subjected to pressure roughly proportional to height (?), so you can test smaller samples for different pressures and see the response. (instead of a big simulation).

Cundall-Strack model is linear, so some linearity in porosity might be expectable.
Hertz-Mindlin is, contrary, non-linear, producing higher force for increasing overlap, so I would expect the response as you want - lower decrease of porosity with increasing stress..
Same results for the basic linear model and Hertz-Mindlin model would be surprising..

So definitely for (some) different models you should get different trends.
Accurate fitting is a different story :-)

cheers
Jan

Can you help with this problem?

Provide an answer of your own, or ask Ricardo Lorenzoni for more information if necessary.

To post a message you must log in.