spheres_in_boxΒΆ

../_images/spheres_in_box.png
#!/usr/bin/env python-mpacts
from mpacts.contact.models.collision.hertz import Hertz
from mpacts.contact.detectors.multigrid import MultiGridContactDetector
from mpacts.contact.detectors.bruteforce import SmartBruteForce
import mpacts.geometrygenerators.pointgenerators as pointgen
import mpacts.geometrygenerators.quadgeometries as quadgen
import mpacts.commands.monitors.progress as pp
from mpacts.core.simulation import Simulation
from mpacts.core.valueproperties import Variable, VariableFunction
import mpacts.tools.random_seed as rs
from mpacts.io.datasave import DataSaveCommand
from mpacts.tools.setcontactdatastorage import SetContactDataStorage
from mpacts.core.new_units import registry as u
import mpacts.particles as prt

""" Simple example that puts some spherical particles into a box and let them settle.
    Output can be viewed by paraview.
"""

#We set the random seed only to ensure exactly equal results between different runs:
rs.set_random_seed( 1 )
mysim = Simulation("simulation", timestep=0.5*u('ms'))

Esphere  = Variable( "E_sphere"  , mysim('params'), value = 1*u("MPa")       )#Young's modulus of spheres
Ebox     = Variable( "E_box"     , mysim('params'), value = 1*u("MPa")       )#Young's modulus of box
Asphere  = Variable( "A_sphere"  , mysim('params'), value = 5*u("ms")        )#Normal damping spheres
Abox     = Variable( "A_box"     , mysim('params'), value = 5*u("ms")        )#Normal damping box
nusphere = Variable( "nu_sphere" , mysim('params'), value = 0.4              )#Poisson's ratio spheres
nubox    = Variable( "nu_box"    , mysim('params'), value = 0.4              )#Poisson's ratio box
rsphere  = Variable( "r_sphere"  , mysim('params'), value = 2*u("cm")        )#Sphere radius
rhosphere= Variable( "rho_sphere", mysim('params'), value = 1000*u("kg/m^3") )#Sphere density
out_int  = Variable( "out_int"   , mysim('params'), value = 50*u('ms')       )#Output interval vtk writer
cd_ue    = Variable( "cd_ue"     , mysim('params'), value = 22               )#Update every

msphere  = VariableFunction( "m_sphere", mysim('params')
                           , function = '4/3. * math.pi * $r_sphere$**3 * $rho_sphere$')

cd_kd    = VariableFunction( "cd_kd", mysim('params'), function = '1.5*$r_sphere$')

#-----------------------------------------------------------------------------------------------------
#Make a particle container for the spheres
spheres = prt.ParticleContainer("spheres", prt.Sphere, mysim)
#Add a data save command, that periodically saves the combined data of the whole mpacts simulation
DataSaveCommand("DataSaveCommand", mysim, executeInterval = out_int )

spheres.GravityCmd()
spheres.TimeIntegration_LeapFrog_Translation_Cmd()
spheres.TimeIntegration_LeapFrog_Rotation_Cmd()
spheres.VTKWriterCmd( executeInterval = out_int, select_all=True )

#-----------------------------------------------------------------------------------------------------
#Make a particle container for the box, that will consist out of RigidQuads
rigid = prt.ParticleContainer("box", prt.RigidBody.compose( prt.RigidQuad ), mysim)
rigid.VTKWriterCmd( executeOnce=True )#Only execute once, since the box does not change

# We calculate the normals of the quads (once) so that we don't need to do that everytime in the contact model:
rigid('RigidQuads').QuadNormalsCmd( x=rigid('controlPoints')['x']
                                  , executeOnce=True)

#-----------------------------------------------------------------------------------------------------
#this is the contact model between the box and the spheres.
cm_box_sphere = Hertz( "Hertz_box_sphere"
                     , pc1=rigid("RigidQuads")
                     , pc2= spheres
                     , A1=Abox, A2=Asphere
                     , E1=Ebox, E2=Esphere
                     , nu1=nubox, nu2=nusphere )

cm_sphere_sphere = Hertz("Hertz_sphere_sphere"
                        , pc1=spheres
                        , pc2=spheres
                        , A1=Asphere, A2=Asphere
                        , E1=Esphere, E2=Esphere
                        , nu1=nusphere, nu2=nusphere )

cd_box_sphere = SmartBruteForce( "CD_box_sphere", mysim
                               , cmodel = cm_box_sphere
                               , keep_distance = cd_kd
                               , update_every = cd_ue )

cm_sphere_sphere = MultiGridContactDetector( "CD_sphere_sphere", mysim
                                           , cmodel = cm_sphere_sphere
                                           , keep_distance = cd_kd
                                           , update_every = cd_ue )

pp.ProgressIndicator("PrintProgress", mysim, print_interval=1)

#-----------------------------------------------------------------------------------------------------
#Actually generate positions of the spheres that will be thrown in the box
xs = pointgen.create_evenly_spaced_in_box( (0.5, 0.5, 0.5)#Lx, Ly, Lz
                                         , (10 , 10 , 10 )#Nx, Ny, Nz
                                         , (0  , 0  , 0  )#Centroid
                                         , 0.01 )#randomness

spheres.add_and_set_particles( x = xs, r=rsphere.get_value(), m=msphere.get_value() )

box_size =  (0.8, 0.8, 2) 
(points, quads) = quadgen.create_box( box_size[0], box_size[1], box_size[2])

#Add a particle (the whole box)
box = rigid.add_particle()
#add the vertices of the quads to the controlpoints
box.controlPoints.add_and_set_particles( xBF = points)
#add the indices of the quads to the RigidQuads
box.RigidQuads.add_and_set_particles( vertexIndices = quads )
#There is only one position of the box, since all quads are rigid with respect to this point
#usually this point is the center of mass of the object.
box.x = (0,0,0)
#equally there is only one orientation of the box, given as a quaternion.
box.q = (1,0,0,0)

#-----------------------------------------------------------------------------------------------------
#print mysim.print_command_tree()
mysim.run_until( 1.5 * u.s )
#-----------------------------------------------------------------------------------------------------
#Verification that the output is identical to a reference value
from mpacts.tools.arrayhash import check_and_document_hash

check_and_document_hash( spheres['x'],  '0f8228dbba19d04b35e9107930e3706e', tolerance=1e-5*rsphere.get_value() )

#-----------------------------------------------------------------------------------------------------