Download this testcase.

2D Polygon Particle Deposition

This example demonstrates the deposition of hexagonal polygon-shaped particles in a rectangular box. The particles interact through oriented contact faces, simulating the packing behavior of non-spherical granular materials.

Keywords

DEM, Polygons, Deposition, Non-spherical Particles

from migflow import scontact, time_integration
import numpy as np
import shutil, os, sys


import time

tic = time.time()

n_particles = 100
use_polygon = True
# use_polygon = False
# Seed for the rng
np.random.seed(0)

output directory

outputdir = f"output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Physical parameters of the particles

rho = 2650  # [kg.m**-3]
friction_coefficient_particle_particle = 0.4
friction_coefficient_particle_wall = 1.4

mm = 1e-3  # [m] reference length
d_min = 0.6 * mm  # [m]
r_min = d_min / 2
d_max = 2 * d_min
n_vertices = 6
r_roundness = 0.05 * d_min
r_roundness = 0.1 * d_min
tol = 1e-3 * 0.05 * d_min  # [m] tolerance for the contact solver

# Domain size in 2d of the cut of the cylinder
domain_length = 50 * mm  # [m]
domain_length = 25 * mm  # [m]
domain_height = 700 * mm  # [m]
x_shift = domain_length
x_bounds = np.array([-domain_length + x_shift, domain_length + x_shift])

Problem definition

dem = scontact.ParticleProblem(2)
dem.set_oriented_faces(True)
dem.set_predict_basis(False)
dem.set_separation_threshold(d_min / 3)
# dem.contact_detection_d = d_min
dem.contact_detection_d = d_min / 100

# Geometry = an immovable body
floor = dem.add_body(
    (0, 0), 0, 0
)  # [0, 0] = origin of the reference coordinates of the body
wall_left = dem.add_body((0, 0), 0, 0)
wall_right = dem.add_body((0, 0), 0, 0)

dem.add_segment_to_body([x_bounds[1], 0], [x_bounds[0], 0], floor, "Wall")
dem.add_segment_to_body(
    [-domain_length / 2 + x_shift, 0],
    [-domain_length / 2 + x_shift, domain_height],
    wall_left,
    "Wall",
)
dem.add_segment_to_body(
    [+domain_length / 2 + x_shift, domain_height],
    [+domain_length / 2 + x_shift, 0],
    wall_right,
    "Wall",
)

Particle initialisation

dem_particles = np.asarray([], int)
x0 = np.array([-domain_length / 2 + x_shift, 0.8 * domain_height])
x1 = np.array([-domain_length / 2 + x_shift, 0])
x2 = np.array([domain_length / 2 + x_shift, 0])
x3 = np.array([domain_length / 2 + x_shift, 0.8 * domain_height])
bnd = np.array([x0, x1, x2, x3])
radii = np.random.uniform(d_min / 2, d_max / 2, n_particles)
x, r = scontact.fastdeposit_2d(bnd, radii)
x[:, 1] += 2 * d_max
r = 0.9 * r


def add_polygon(x_center, w, h, n_vertices, theta_0=0.0):
    if n_vertices < 3:
        raise ValueError("n_vertices must be at least 3")
    theta = np.linspace(0, 2 * np.pi, n_vertices, endpoint=False) + theta_0
    x = w * np.cos(theta) + x_center[0]
    y = h * np.sin(theta) + x_center[1]
    vertices = np.array([x, y]).T
    nodes = np.arange(n_vertices)
    edges_0 = nodes.copy()
    edges_1 = np.append(nodes[1:], nodes[0])
    edges = np.array((edges_0, edges_1)).T
    return vertices, edges


def crop_polygon(vertices, edges, r):
    n_vertices = len(edges) * 2
    n_centers = len(vertices)
    new_vertices = np.zeros((n_vertices, 2))
    circle_centers = np.zeros((n_centers, 2))
    for i in range(len(edges)):
        v0 = vertices[edges[i][0]]
        v1 = vertices[edges[i][1]]
        v2 = vertices[edges[(i + 1) % len(edges)][1]]
        du = (v0 - v1) / np.linalg.norm(v0 - v1)
        dv = (v2 - v1) / np.linalg.norm(v2 - v1)
        cos_theta = np.dot(du, dv)
        alpha = np.arccos(cos_theta)
        tan_theta = np.tan(alpha / 2)
        a = r / tan_theta
        x0 = v1 + a * du
        x1 = v1 + a * dv
        t_dv = np.array([-dv[1], dv[0]])
        x_center = x1 + r * t_dv
        new_vertices[2 * i + 0] = x0
        new_vertices[2 * i + 1] = x1
        circle_centers[i] = x_center
    return new_vertices, circle_centers


x_ref, edges_ref = add_polygon((0, 0), r_min, r_min, n_vertices, 0)
x_smoothed, circle_centers = crop_polygon(x_ref, edges_ref, r_roundness)
for xi, ri in zip(x, r):
    polygon_mass = np.pi * ri**2 * rho  # [kg] mass of the particle
    polygon_inertia = 0.5 * polygon_mass * ri**2
    body = dem.add_body((xi[0], xi[1]), 1 / polygon_mass, 1 / polygon_inertia)
    if use_polygon:
        for center in circle_centers:
            dem.add_particle_to_body(
                center / r_min * ri, r_roundness / r_min * ri, body
            )
        for i in range(len(edges_ref)):
            x0 = x_smoothed[(2 * i + 1) % len(x_smoothed)] / r_min * ri
            x1 = x_smoothed[(2 * i + 2) % len(x_smoothed)] / r_min * ri
            dem.add_segment_to_body(x0, x1, body, "line")
    else:
        dem.add_particle_to_body((0, 0), ri, body)

Simulation parameters

dt = 1e-4
tEnd = 1000 * dt
outf = 10
iit = 0
t = 0

particles_forces = np.zeros_like(dem.velocity())
body_forces = np.divide(
    -9.81, dem.body_invert_mass(), where=dem.body_invert_mass() > 0
).reshape(-1, 1) * np.array([0, 1])
nonzero = dem.body_invert_mass().reshape(-1) != 0
dem.body_theta()[nonzero, :] = np.random.uniform(
    0, 2 * np.pi, (dem.body_theta()[nonzero, :].shape)
)

while t < tEnd:
    print(
        f"{iit=: 04d} -- t/tEnd: {t:.3f}/{tEnd: .3f} -- {t/tEnd*100:02.2f}% completed"
    )
    if iit % outf == 0:
        dem.write_mig(outputdir, t)
    dem.iterate(dt, particles_forces, body_forces, tol, 1)
    contacts = dem.contacts()
    reaction = contacts["reaction"]
    active_contacts = np.sum(np.linalg.norm(reaction, axis=1) > 0)
    print(f"  Active contacts: {active_contacts} / {len(reaction)}")
    # time_integration._advance_particles(dem, particles_forces, dt, 1, tol, fbody = body_forces)
    t += dt
    iit += 1
    dem.contacts()[:] = 0.0  # reset contact reactions for next iteration

toc = time.time()
print(f"Simulation time: {toc - tic:.2f} seconds")

Plot

python3 -m migflow.plot.migplot output --actors particles