spheres_between_cylindersΒΆ

../_images/spheres_between_cylinders.png
#!/usr/bin/env python-mpacts
import mpacts.particles as prt
from mpacts.core.simulation import Simulation
from mpacts.core.valueproperties import Variable, VariableFunction
from mpacts.commands.monitors.progress import ProgressIndicator
from mpacts.contact.detectors.multigrid import MultiGridContactDetector
from mpacts.contact.models.collision.hertz.hertz_basic import HertzCoulomb
from mpacts.geometrygenerators.pointgenerators import create_evenly_spaced_in_box
from mpacts.tools.setcontactdatastorage import SetContactDataStorage
from mpacts.io.datasave import DataSaveCommand
import mpacts.geometrygenerators.quadgeometries as quadgen
from mpacts.tools.random_seed import set_random_seed
#-------------------------------------------------------------------------

#To ensure that the random numbers are always chosen the same.
#This enables hash checking of the final result (see below).
#(Randomness comes from the initial positioning here.)
set_random_seed(1)

# Specifying some simulation parameters and constructing the simulation tree:

mysim = Simulation("simulation", timestep=2.5e-5)
SetContactDataStorage("Vector")
p = mysim('params')

Variable( 'endtime',p, value=3.0 )
Variable( 'rho', p, value= 2000)
Variable( 'r', p, value=0.04)
Variable( "E", p, value=2.5e6)
Variable( "nu", p, value = 0.4)
Variable( "A", p, value = 0.001)
Variable( "mu", p, value = 0.25)
Variable( "ct", p, value = 5.0)
Variable( "cdrag", p, value=0.1)
Variable( "cd_ue", p, value = 25)
Variable( "out_int", p, value = 0.05)

VariableFunction( 'm', p, function = '4/3.*math.pi*$rho$*$r$**3' )
VariableFunction( "cd_kd", p, function = '0.5*$r$')

#-------------------------------------------------------------------------
# Constructing a box particle pc in which the spheres will drop.
# This is designed as a rectangular box, with transversal bars composed of
# rigid cylinders:

pc = prt.ParticleContainer( "pc", prt.RigidBody.compose(prt.RigidQuad, prt.RigidCylinder), mysim)

pc('RigidQuads').QuadNormalsCmd( x=pc('controlPoints')['x'], executeOnce=True)

cylx = [(0, 0, -1), (0, 0, 1)]

# Positions of the bars in the 'xy' plane. The bars are parallel to the
# 'z' direction
bar_pos = [ (-0.666, 0, 0), (-0.333, 0, 0), (0, 0, 0)
          , (0.333, 0, 0), (0.666, 0, 0), (-0.83333, -0.5, 0)
          , (-0.5, -0.5, 0), (-0.16666, -0.5, 0), (0.16666, -0.5, 0)
          , (0.5, -0.5, 0), (0.83333, -0.5, 0)]

# Positions and connectivity of the box.
(cub_ctrl, cub_vil) = quadgen.create_box(2, 2, 2, False)

# The cylinder controlPoints and vertexIndicis are computed.
# We have to be careful to 'offset' the vertexIndices with the size of the
# box controlPoints
bar_ctrl = []
bar_vil = []
for i, bar in enumerate(bar_pos):
    bar_vil.append((2 * i + len(cub_ctrl), 2 * i + 1 + len(cub_ctrl)))
    for endpnt in cylx:
        bar_ctrl.append((bar[0], bar[1], endpnt[2]))

# Collect all control points for '1' box together in one list
box_ctrl = cub_ctrl + bar_ctrl

# Create a box particle.
# Because the whole box plus bars is ONE particle, it would be trivial to add more 'boxes'
# or e.g. translate / rotate a complete box in a simulation
box = pc.add_particle()
box.x = (0, 0, 0)
box.q = (1, 0, 0, 0)
box.v = (0, 0, 0)
box.w = (0, 0, 0)
box.controlPoints.add_and_set_particles(xBF=box_ctrl)
box.RigidCylinders.add_and_set_particles(vertexIndices=bar_vil, r=0.09)
box.RigidQuads.add_and_set_particles(vertexIndices=cub_vil)
#-------------------------------------------------------------------------

# Sphere positions are distributed initially in a uniform cartesian grid
# above the (open) box:
xs = create_evenly_spaced_in_box( (1.8, 1.8, 1.8)
                                , (8, 8, 8)
                                , (0, 1.8, 0)
                                , 0.001)

# Create a particle pc for the spheres:
spheres = prt.ParticleContainer("Spheres", prt.Sphere, mysim)

# Spheres are a simple geometry, so we can directly add the particles in the pc
spheres.add_and_set_particles(x=xs, r=p('r').get_value(), m=p('m').get_value())

# Adding commands for spheres: gravity, drag force and leap frog time
# integration:
spheres.GravityCmd()
spheres.DragForceCmd(c=p('cdrag'))
spheres.TimeIntegration_LeapFrog_Translation_Cmd()
spheres.TimeIntegration_LeapFrog_Rotation_Cmd()
DataSaveCommand("DataSaveCommand", mysim, executeInterval = p('out_int') )

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

# Adding Hertz Coulomb contact models between spheres, quads and cylinders:

spherebarsCM = HertzCoulomb( pc1=pc("RigidCylinders"),
                             pc2=spheres,
                             E1=p('E'), E2=p('E'),
                             nu1=p('nu'), nu2=p('nu'),
                             A1=p('A'), A2=p('A'),
                             c_t=p('ct'), mu=p('mu'))

sphereboxCM = HertzCoulomb( pc1=pc("RigidQuads"),
                            pc2=spheres,
                            E1=p('E'), E2=p('E'),
                            nu1=p('nu'), nu2=p('nu'),
                            A1=p('A'), A2=p('A'),
                            c_t=p('ct'), mu=p('mu'))

spheresphereCM = HertzCoulomb( pc1=spheres,
                               pc2=spheres,
                               E1=p('E'), E2=p('E'),
                               nu1=p('nu'), nu2=p('nu'),
                               A1=p('A'), A2=p('A'),
                               c_t=p('ct'), mu=p('mu'))

# Adding multigrid contact detectors for all the contact models:

spherebarsCD = MultiGridContactDetector( "sphere-bar-cd", mysim,
                                         cmodel=spherebarsCM,
                                         update_every=p('cd_ue'),
                                         keep_distance=p('cd_kd') )

sphereboxCD = MultiGridContactDetector( "sphere-box-cd", mysim,
                                        cmodel=sphereboxCM,
                                        update_every=p('cd_ue'),
                                        keep_distance=p('cd_kd') )

spheresphereCD = MultiGridContactDetector( "Sphere-sphere-cd", mysim,
                                           cmodel=spheresphereCM,
                                           update_every=p('cd_ue'),
                                           keep_distance=p('cd_kd') )

#-------------------------------------------------------------------------
# Finally, adding output writer commands, a progress printer and starting
# the simulation:

spheres.VTKWriterCmd(executeInterval=p('out_int'), select_all = True)

pc.VTKWriterCmd(executeOnce=True, select_all=True)

ProgressIndicator("PrintProgress", mysim, print_interval=2)

mysim.run_until(p('endtime').get_value())

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

#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'], 'fa374f4211a842d17434d4663d846a12', tolerance=1e-5) #reason unknown!
check_and_document_hash( spheres['x'], 'ca33ba9b4a81d7d43d033761e2f65ee4', tolerance=1e-5)


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