This is a sample script that simulates the compaction of a powder using a PID controller. The system consists of the particles (spheres), 4 bounding walls defining a rectangular box, and a moving mesh (PID controller) along the z direction that is used to compress the powder at a pressure of 50 bar. Note that the behaviour of PID controllers is defined at the onset of a simulation, and there is no need to specifically call the 'move' method since the velocity of the mesh is set by the pressure it 'feels'.
The input script and mesh files are found here.
First, a python dictionary (params) defining the simulation parameters is created as shown below:
from pygran import simulation
from pygran.params import organic
# Setup material parameters
organic['cohesionEnergyDensity'] = 1e5
organic['youngsModulus'] = 6e7
# Assume the wall is non-cohesive
wall = organic.copy()
wall['cohesionEnergyDensity'] = 0.0
# Assume no friction between the bounding walls and the particles
PID = wall.copy()
wall['coefficientFriction'] = 0
# Setup arguments for the PID mesh
args_pid = {'scale':9.8e-4, 'move':(0, 0, 6e-3), 'com':(0, 0, 0), 'axis':(0, 0, -1), 'vel_max':10.0, \
'ctrlPV':'force', 'target_va':5.0}
# Setup arguments for the bounding meshes
args_wallX1 = {'scale':1e-3, 'rotate':('axis', 0, 1, 0, 'angle', 90), 'move':(0, 0, 2e-3)}
args_wallX2 = {'scale':1e-3, 'rotate':('axis', 0, 1, 0, 'angle', -90), 'move':(0, 0, 2e-3)}
args_wallY1 = {'scale':1e-3, 'rotate':('axis', 1, 0, 0, 'angle', 90, 'rotate', 'axis', 0, 1, 0, 'angle', 90), \
'move':(0, 0, 2e-3)}
args_wallY2 = {'scale':1e-3, 'rotate':('axis', 1, 0, 0, 'angle', -90, 'rotate', 'axis', 0, 1, 0, 'angle', 90), \
'move':(0, 0, 2e-3)}
# Define the length of the bounding box
xdim = 1e-3
# Create a dictionary of physical parameters
params = {
# Define the system
'boundary': ('p','f','f'), # fixed BCs
'box': (-xdim, xdim, -2*xdim, 2*xdim, 0, 10*xdim),
'model': simulation.models.HertzMindlin,
# Define component(s)
'species': ({'material': organic, 'radius': ('constant', 2e-4)},
),
# Timestep
'dt': 1e-6,
'nns_skin': 1e-4,
# Apply gravitional force in the negative direction along the z-axis
'gravity': (9.81, 0, 0, -1),
# Number of simulation steps (non-pygran variable)
'run': 1e5,
'traj': {'mfile': 'mesh*.vtk', 'freq':1000, 'pfile': 'traj.dump'},
'output': 'output',
# Import surface mesh
'mesh': {
'PID': {'file': 'mesh/wall-2D.stl', 'mtype': 'mesh/surface/stress/servo', 'material': PID, \
'args': args_pid},
'wallX1': {'file': 'mesh/wall-2D-side.stl', 'mtype': 'mesh/surface/stress', 'material': wall, \
'args': args_wallX1},
'wallX2': {'file': 'mesh/wall-2D-side.stl', 'mtype': 'mesh/surface/stress', 'material': wall, \
'args': args_wallX2},
'wallY1': {'file': 'mesh/wall-2D-side.stl', 'mtype': 'mesh/surface/stress', 'material': wall, \
'args': args_wallY1},
'wallY2': {'file': 'mesh/wall-2D-side.stl', 'mtype': 'mesh/surface/stress', 'material': wall, \
'args': args_wallY2},
}
}
The params dictionary is then used to create a pygran.DEM class and run the simulation.
# Create an instance of the DEM class
sim = simulation.DEM(**params)
# Setup a primitive wall along the xoy plane at z=0
sim.setupWall(species=1, wtype='primitive', plane = 'zplane', peq = 0.0)
# Create an insertion box of height 4 times its length/width
box = ('block',-xdim, xdim, -xdim, xdim, 0, 8*xdim)
# Insert particles by volume fraction
sim.insert(species=1, value=1.0, mech='volumefraction_region', region=box)
# Run the simulation
sim.run(params['run'], params['dt'])
# Play video
import io
video = io.open('movie/PID.mp4', 'r+b').read()