Size problem in yade simulation

Asked by Huang peilun on 2020-01-27

I did a gravity deposition test in a cylinder. The following is the original code and it works.

************************************************
from yade import pack

#material
SoilMat=FrictMat(young=25e6,poisson=0.2,density=2650,frictionAngle=0.28)
O.materials.append((SoilMat))

#define the chamber
O.bodies.append(geom.facetCylinder(center=(.5,.5,.5),radius=.5,height=1,wallMask=6,color=(1,0,0)))

#create foundation by making spheres in the cylinder.
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.05,rRelFuzz=.3)
cyl=pack.inCylinder((.5,.5,0),(.5,.5,1),.45)
sp=pack.filterSpherePack(cyl,sp,True)
sp.toSimulation(material=SoilMat)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 # update position using Newton's equations
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.05),
 # call the checkUnbalanced function (defined below) every 2 seconds
 PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.2*PWaveTimeStep()

def checkUnbalanced():
 print('%s: Unbalanced=%s' % (O.iter,unbalancedForce()))
 if O.iter<5000: return
 if unbalancedForce()<.05:
  O.pause()

O.saveTmp()
************************************************

However, when I magnify all units 1000 times, like the following code, some sphere just went through the bottom of the cylinder and I don't know why.

************************************************
from yade import pack

#material
SoilMat=FrictMat(young=25e6,poisson=0.2,density=2650,frictionAngle=0.28)
O.materials.append((SoilMat))

#define the chamber
O.bodies.append(geom.facetCylinder(center=(500,500,500),radius=500,height=1000,wallMask=6,color=(1,0,0)))

#create foundation by making spheres in the cylinder.
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1000,1000,1000),rMean=50,rRelFuzz=.3)
cyl=pack.inCylinder((500,500,0),(500,500,1000),450)
sp=pack.filterSpherePack(cyl,sp,True)
sp.toSimulation(material=SoilMat)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 # update position using Newton's equations
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.05),
 # call the checkUnbalanced function (defined below) every 2 seconds
 PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.2*PWaveTimeStep()

def checkUnbalanced():
 print('%s: Unbalanced=%s' % (O.iter,unbalancedForce()))
 if O.iter<5000: return
 if unbalancedForce()<.05:
  O.pause()

O.saveTmp()
************************************************

Thanks in advance!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-01-28
Last query:
2020-01-28
Last reply:
2020-01-27
Robert Caulk (rcaulk) said : #1

Have you tried increasing the stiffness of the facets relative to the particles?

Huang peilun (hpl16) said : #2

Thanks Robert.
I haven't tried increasing the stiffness of the facets because I don't know how. Can you show me how to do that?

Robert Caulk (rcaulk) said : #3

For example [1].

I will also give you the answer:

O.materials.append(FrictMat(young=desiredStiffness,label='walls'))
O.bodies.append(geom.facetCylinder(center=(.5,.5,.5),radius=.5,height=1,wallMask=6,color=(1,0,0)),material='walls')

However, I still recommend you read [2] carefully to ensure functionality.

Cheers.

Robert

[1]https://gitlab.com/yade-dev/trunk/blob/master/examples/triax-tutorial/script-session1.py#L54
[2]https://yade-dev.gitlab.io/trunk/tutorial-hands-on.html#singles

Huang peilun (hpl16) said : #4

Thanks!
I increased the stiffness but there still are spheres that go through the bottom of the cylinder. Here's my code.

********************************
from yade import pack

Mat=FrictMat(young=2.06e20)
SoilMat=FrictMat(young=25e6,poisson=0.2,density=2650,frictionAngle=0.28)
O.materials.append((SoilMat,Mat))

O.bodies.append(geom.facetCylinder(center=(500,500,500),radius=500,height=1000,wallMask=6,color=(1,0,0),material=Mat))

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1000,1000,1000),rMean=50,rRelFuzz=.3)
cyl=pack.inCylinder((500,500,0),(500,500,1000),450)
sp=pack.filterSpherePack(cyl,sp,True)
sp.toSimulation(material=SoilMat)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 # update position using Newton's equations
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.05),
 # call the checkUnbalanced function (defined below) every 2 seconds
 PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.2*PWaveTimeStep()

def checkUnbalanced():
 print('%s: Unbalanced=%s' % (O.iter,unbalancedForce()))
 if O.iter<5000: return
 if unbalancedForce()<.05:
  O.pause()

O.saveTmp()

Best Jan Stránský (honzik) said : #5

Hello,

welcome to the playing with physical dimensions and units :-)
I assume you want in both cases "same" results:
- you increase size by a factor f (here f=1000)
- cross sections A increases by factor f^2
- volume and mass and gravity force increases with factor f^3
- penetration depth (pd) increases with factor f
- you left stiffness E unchanged, increasing by factor 1
- repulsive force F=pd*E*A/L increases with f*1*f^2/f = f^2
now it becomes 1000x "smaller" then in the original case! and is not large enough to prevent gravity force pushing particles through the cylinder..
So to "scale" repulsive force the same way as gravity force, you have to scale the stiffness by the factor f, too. Then F=pd*E*A/L=f*f*f^2/f=f^3.

Using
SoilMat=FrictMat(young=25e9,...)
in the second script, I got very similar results as in the first one

cheers
Jan

Huang peilun (hpl16) said : #6

Thanks Jan!
I tried SoilMat=FrictMat(young=25e9,...) and it works!
Now I understand why these spheres go through the cylinder. It's very interesting!

Huang peilun (hpl16) said : #7

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