Download this testcase.

2D Polygon Pile

This example simulates the formation of a pile from octagonal polygon-shaped particles falling into a rectangular box. The particles interact through frictional contacts, illustrating the packing behavior of non-spherical granular materials.

Keywords

DEM, Polygons, Pile, Non-spherical Particles

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

# n_particles = 5000
# n_particles = 2500
n_particles = 1000
# for n_particles in[5000]:
# Seed for the rng
np.random.seed(0)

output directory

outputdir = f"output_{n_particles}" 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 = 8
r_roundness = 0.05 * d_min
tol = 1e-3 * r_roundness

# 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]

Problem definition

dem = scontact.ParticleProblem(2)
dem.set_oriented_faces(True)
dem.set_friction_coefficient(0.4)
dem.set_separation_threshold(r_roundness)
dem.contact_detection_d = 3 * r_roundness
dem.set_friction_coefficient(friction_coefficient_particle_particle, "Sand", "Sand")
dem.set_friction_coefficient(friction_coefficient_particle_wall, "Sand", "Wall")

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

x_shift = domain_length
x_bounds = np.array([-domain_length * 20 + x_shift, domain_length * 20 + x_shift])
dem.add_segment_to_body([x_bounds[1], 0], [x_bounds[0], 0], floor, "Wall", "Wall")
dem.add_segment_to_body(
    [-domain_length / 2 + x_shift, 0],
    [-domain_length / 2 + x_shift, domain_height],
    wall_left,
    "Wall",
    "Wall",
)
dem.add_segment_to_body(
    [+domain_length / 2 + x_shift, domain_height],
    [+domain_length / 2 + x_shift, 0],
    wall_right,
    "Wall",
    "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)


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)
    for center in circle_centers:
        dem.add_particle_to_body(
            center / r_min * ri, r_roundness / r_min * ri, body, "Sand"
        )
    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", "Sand")

Simulation parameters

dt = 1e-4
tEnd = 0.1
outf = 10
iit = 0
t = 0
velocity_walls = 0.05

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)
    if t > 0.0:
        dem.body_velocity()[wall_right, 1] = velocity_walls
        dem.body_velocity()[wall_left, 1] = velocity_walls
    time_integration._advance_particles(
        dem, particles_forces, dt, 1, tol, fbody=body_forces
    )
    t += dt
    iit += 1