Strange behavior: sand pile settles in a squared shape ?

Asked by Luc Scholtès

Hi all,

When simulating the settlement of a relatively dense) packing on a plane surface, I obtain an equilibrium state which seems a bit awkward: from above, the particles tend to spread in a squared shape instead of a circular shape (as I would have expected). I tested different contact laws as well as different initial packing shapes and it seems to occur systematically (the squared shape is more or less evident depending on the case). The initial packing is built based on randomDensePack(): could this be the reason?

Please find below a MWE if you want to reproduce the "problem" (jut click and play for ~10 000 iterations when the GUI opens).

Any feedback would be appreciated.

Luc

---

# -*- coding: utf-8 -*-

from yade import pack

#### assembly
O.materials.append(FrictMat(density=2600,young=1e6,poisson=0.333,frictionAngle=radians(30),label='sphereMat'))

L=2.
D=L/2.
sphereRad=0.5*(min(L,D)/15.) # define sphere mean radius as a function of model dimensions / you can directly define a value
pred=pack.inCylinder((0,0,0),(0,0,L),D/2.)
#pred=pack.inSphere((0,0,0),D/2.)
O.bodies.append(pack.randomDensePack(pred,radius=sphereRad,rRelFuzz=0.333,spheresInCell=3000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,material='sphereMat'))

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

#### get spheres dimensions (pre-processing)
R=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   o.shape.color=(0.8,0.1,0.1)
   numSpheres+=1
   R+=o.shape.radius
Rmean=R/numSpheres

print('nb of spheres=',numSpheres,' | Rmean=',Rmean)

#### floor
O.materials.append(FrictMat(density=2600,young=1e6,poisson=0.333,frictionAngle=radians(30),label='floorMat'))
O.bodies.append(box(center=(0.,0.,zinf-Z/10.),extents=(X*100.0,Y*100.0,Z/10.),fixed=True,wire=False,color=(0,1,0),material='floorMat'))

##### engines
O.engines=[
    ForceResetter()
    ,InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()])

    ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
    )

    ,GlobalStiffnessTimeStepper(defaultDt=0.1*utils.PWaveTimeStep(),timestepSafetyCoefficient=0.5)
    ,NewtonIntegrator(damping=0.7,gravity=(0.,0.,-9.82),label='newton')
]

#### open yade GI
from yade import qt
v=qt.Controller()
v=qt.View()

#### Allows to reload the simulation
O.saveTmp()

print('press play now!)')

Question information

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

Hi Luc,

have a try with lower NewtonIntegrator.damping.
For zero damping, I would say the settlement is visually circular - although a more rigorous than "visual" metric should be defined for serious conclusion.

The numerical damping is computed not "isotropically", but component-wise and is therefore is "clearly non-physical, as it is not invariant with respect to coordinate system rotation" [1].
It could result in different behaviour along axes and "between" axes, resulting in "square-shaped" settlement.

Cheers
Jan

[1] https://yade-dem.org/doc/formulation.html#numerical-damping

Revision history for this message
Luc Scholtès (luc) said :
#2

Thank you Jan. Indeed, the damping coefficient plays a role! I guess I'd need to consider viscous damping at the contact scale to avoid such "non physical" behavior, right?

Luc

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

If the current non-physical behavior is a problem, then probably the only way is to use some physical one - either change the current NewtonIntegrator approach, or use interaction-level damping.
Jan

Revision history for this message
Luc Sibille (luc-sibille) said :
#5

Hello,

I don't know and use the function randomDensePack, but apparently it
makes a periodic packing (it is the information displayed in the console).
Then, is the packing really periodic? Is it in relation with the square
shape? Did you try to generate the initial assembly with makecloud for
instance, and does it change the final result?

Cheers,
Luc

Le 28/04/2021 à 10:01, Luc Scholtès a écrit :
> Question #696786 on Yade changed:
> https://answers.launchpad.net/yade/+question/696786
>
> Status: Answered => Solved
>
> Luc Scholtès confirmed that the question is solved:
> Thank you Jan. Indeed, the damping coefficient plays a role! I guess I'd
> need to consider viscous damping at the contact scale to avoid such "non
> physical" behavior, right?
>
> Luc
>

--
Luc Sibille
Université Grenoble Alpes / IUT1 de Grenoble
Laboratoire 3SR: Sols, Solides, Structures, Risques

Tel lab.: +33 (0)4 76 82 63 48
Tel IUT: +33 (0)4 76 82 53 36

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

> I don't know and use the function randomDensePack, but apparently it makes a periodic packing (it is the information displayed in the console).

If you use spheresInCell parameter, it:
- first "quickly" create a periodic cube ("quickly" w.r.t full domain, as the periodic packing is possibly much much smaller than the resulting domain)
- copy "roughly" the cube wherever needed to fully cover the domain (roughly = also outside bondaries of the domain)
- crop "finely" the copied cubes to form the required domain

> Then, is the packing really periodic?

Not necessarily (only by chance). The packing contains periodically repeated "cubes", but itself is periodic if and only if is is cropped in a very special way.

> Is it in relation with the square shape?

I would say no, or at least negligibly.

Cheers
Jan