Yade

Problem with high value of material density

Asked by Ramin

I would like to reduce the computational cost (increase the time step) by up-scaling the particle density to 109 kg⁄m^3 but when I launch the simulation, the particles dissapeared after a few iterations. It looks like they doesn't contact with the bottom side of the box. My script is attached below. Please give me any suggestions how I should deal with it. Thanks

############################
from math import pi
import os, datetime
from yade import pack, plot, qt, geom

sector_angle = 2*pi
layers_count = 5

x, y, z = (1, 1, 1)
initialAngle = radians(10)
targetVoid = 1
density = 1e6

frictionAngle = initialAngle
# create the box
O.bodies.append(geom.facetBox((0, 0, 0), (x/2, y/2, z/2), wallMask=31))

# create spheres
sp = pack.SpherePack()
sp.makeCloud((-x/2, -y/2, -z/2), (x/2, y/2, z/2), rMean=.03, rRelFuzz=.5)
# add spheres to simulation
O.materials.append(CohFrictMat(label='mat',
                                young=1e8,
                                poisson=0.2,
                                frictionAngle=radians(17),
                                density=density,
                                isCohesive=True,
                                alphaKr=2,
                                alphaKtw=2,
                                momentRotationLaw=True,
                                etaRoll=0.15))
sp.toSimulation(material='mat')

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
            # handle sphere+sphere and facet+sphere collisions
            [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys()],
            [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    NewtonIntegrator(gravity=(0, 0, -9.81*10), damping=0.1),
    TranslationEngine(ids=[],
                      translationAxis=[0, 0, -1],
                      velocity=0.075,
                      dead=True,
                      label='cyl_engine'),
    PyRunner(command='checkPorosity()',
             iterPeriod=1000,
             label='porosity_checker',
             dead=False)
]

def checkPorosity():
    global frictionAngle
    u = utils.porosity()
    ev = u / (1 - u)
    if ev > targetVoid:
        frictionAngle *= 0.9
        utils.setContactFriction(frictionAngle)
    else:
        O.pause()

O.dt = .5 * PWaveTimeStep()

qt.View()

Question information

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

Hello,

TLTR: you cannot fake physics this way

> density to 109 kg⁄m^3
> density = 1e6

please be consistent

> It looks like they doesn't contact with the bottom side of the box

Instead of "it looks", always try to get (numerical) evidence.
I.e. print contacts and forces of the wall. You would see that the particles DO contact with the wall.
Even just visually, I can see some small particles clearly interacting with wall if I investigate your script step-by-step.

> I would like to reduce the computational cost (increase the time step) by up-scaling the particle density

You are caught in your own physics trap.
You significantly (unrealistically) increased density to increase time step (therefore you significantly changed the physics of your problem and the results would probably be anyway meaningless).
Now the gravity force applied to the particles is significantly higher. And the (unchanged) stiffness of particle-wall contact cannot hold such force.
To fix this, you would need to increase stiffness. Therefore you would decrease time step, ending at the very beginning.

What you can try is "selective density upscale". Instead of increasing density of all particles, you increase density only to the smallest particles - because the critical time step is determined by the smallest particles (for constant stiffness and density) [1].
"Faking" only the smallest particles, you would also not "fake" the whole system "significantly" (w.r.t. total mass of the whole system).

Also using GlobalStiffnessTimeStepper instead of using fixed time step might help in this cases.

Cheers
Jan

[1] https://yade-dem.org/doc/formulation.html#stability-considerations

Revision history for this message
Ramin (ginmaru) said : 		#2

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