yade-users team mailing list archive

[François Nader <question692853@xxxxxxxxxxxxxxxxxxxxx>]

New question #692853 on Yade:
https://answers.launchpad.net/yade/+question/692853

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()

###

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.