Polyhedra Simulation - Numerical Parameters

Asked by François Nader on 2020-09-10

Hello,

I am trying to model a deposit of polyhedral particles in YADE.
I modified a pre-existing example in YADE's documentation (textExport.py).
However, I am not sure about the choice of the timestep.
And I would like to know if it possible to artificially increase the inertia tensor or damp the particles' rotation, because I have a very high residual rotational energy.
Can you please inform me about the choice of numerical parameters and whether the ones I use in the following example make sense from a numerical point of view.
Thanks in advance!
François

Here's my working code

###

from __future__ import print_function
from yade import plot, polyhedra_utils, export
from yade import qt
import numpy as np
import sys, time, os

if not(os.path.isdir('VTK')):
    os.mkdir('VTK')

m = PolyhedraMat()
m.density = 2700 #kg/m^3
m.young = 1E6 #Pa
m.poisson = 20000/1E6
m.frictionAngle = 0.67474 #rad

freq = 100 #output frequency, in iterations

dx = 0.07
nx = 8
dy = 0.07
ny = 8
dz = 0.07
nz = 8

O.bodies.append(utils.wall(0,axis=0,sense=1, material = m))
O.bodies.append(utils.wall((nx+1)*dx,axis=0,sense=-1, material = m))
O.bodies.append(utils.wall(0,axis=1,sense=1, material = m))
O.bodies.append(utils.wall((ny+1)*dy,axis=1,sense=-1, material = m))
O.bodies.append(utils.wall(0,axis=2,sense=1, material = m))
O.bodies.append(utils.wall((nz+1)*dz,axis=2,sense=-1, material = m))

for k in range(0,nz):
    for j in range(0,ny):
        for i in range(0,nx):
            t = polyhedra_utils.polyhedra(m,size = (0.05,0.05,0.05),seed = 1)
            t.state.pos = (dx*(1+i),dy*(1+j),dz*(1+k))
            u = np.random.random()
            v = np.random.random()
            w = np.random.random()
            t.state.ori = Quaternion((math.sqrt(1-u)*math.sin(2*math.pi*v), math.sqrt(1-u)*math.cos(2*math.pi*v), math.sqrt(u)*math.sin(2*math.pi*w)), math.sqrt(u)*math.cos(2*math.pi*w))
            O.bodies.append(t)

vtkExporter = export.VTKExporter('VTK/',0)
def vtkOutput():
    vtkExporter.exportPolyhedra(what=[('n','b.id'),('vel','b.state.vel')])
    vtkExporter.exportInteractions(what=[('kn','i.phys.kn'),('Rn','i.phys.normalForce'),('Rt','i.phys.shearForce')])
    vtkExporter.exportFacets(what= [('pos','b.state.pos')])
    vtkExporter.exportFacetsAsMesh(what=[('pos','b.state.pos')])

vtkOutput() #initial state VTK

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
   ),
   NewtonIntegrator(damping=0.99,gravity=( 0,0,-9.81)),
   PyRunner(command='print_state()',iterPeriod=freq),
   PyRunner(command='vtkOutput()',iterPeriod=freq),

]

O.dt=0.025*polyhedra_utils.PWaveTimeStep()
#O.dt=0.00025

qt.Controller()
V = qt.View()

O.saveTmp()
O.run()

###

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-09-15
Last query:
2020-09-15
Last reply:
2020-09-10
Best Jan Stránský (honzik) said : #1

Hello,

below are simple answers to your questions:

> However, I am not sure about the choice of the timestep.

chose such time step that leave the simulation stable.
Unfortunately I don't know if there is something like PWaveTimeStep for the polyhedra law..
trial-and-error time step determination can always be used

> And I would like to know if it possible to artificially increase the inertia tensor or damp the particles' rotation, because I have a very high residual rotational energy.

Yes, it is possible to artificially increase inertia tensor (b.state.inertia = ...).
You can also block degrees of freedom totally (b.state.blockedDOFs).
If such results better correspond to expected results, why not.

> Can you please inform me about the choice of numerical parameters and whether the ones I use in the following example make sense from a numerical point of view.

any set of numerical parameters may be good, it depends on your problem

> m.young = 1E6 #Pa

note that in the case of Law2_PolyhedraGeom_PolyhedraPhys_Volumetric, the meaning is not Pa, but has different meaning

Recently I had a discussion of friend of mine using a few polyhedrons and their collisions observing the same, too much rotation. They tried PotentialBlocks with better results, so you also may try those.

cheers
Jan

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