# Polyhedra Simulation - Numerical Parameters

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.
François

Here's my working code

###

from __future__ import print_function
from yade import plot, polyhedra_utils, export
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

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:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-09-15
Last query:
2020-09-15
2020-09-10 Jan Stránský (honzik) said on 2020-09-10: #1

Hello,

> 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