Download this testcase.

2D Hexagonal Polygon Deposition

This example demonstrates the deposition of hexagonal polygon-shaped particles in a rectangular box. It uses the oriented face contact algorithm with a fine contact detection distance to handle hexagonal grain interactions.

Keywords

DEM, Polygons, Hexagonal, Deposition

import os
import shutil

import numpy as np
from migflow import scontact, time_integration

Initialization

np.random.seed(0)

outputdir = "output"
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir, exist_ok=True)

Physical parameters

w = 1.0
r = 1e-2 * w
tol = 1e-4 * r

mass = 1.0
inertia = (mass * w**2) / 2.0

mu = 0.4
g = np.array([0.0, -9.81])

DEM problem

dem = scontact.ParticleProblem(2)
dem.set_oriented_faces(True)
dem.set_friction_coefficient(mu)
dem.set_separation_threshold(r / 2)
dem.contact_detection_d = 3 * r

Geometry utilities

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

Bodies

hexa = dem.add_body((0, 0), 1 / mass, 1 / inertia)

x_poly, edges_poly = add_polygon((0, 0), w, w, 12, 0)
x_smoothed, circle_centers = crop_polygon(x_poly, edges_poly, r)

for center in circle_centers:
    dem.add_particle_to_body(center, r, hexa)

for i in range(len(edges_poly)):
    dem.add_segment_to_body(
        x_smoothed[(2 * i + 1) % len(x_smoothed)],
        x_smoothed[(2 * i + 2) % len(x_smoothed)],
        hexa,
        "line",
    )

y_min = np.min(x_smoothed[:, 1])

octo = dem.add_body((0, -2 * y_min + w), 1 / mass, 1 / inertia)

x_poly, edges_poly = add_polygon((0, 0), w, w, 64, 0)
x_smoothed, circle_centers = crop_polygon(x_poly, edges_poly, r)

for center in circle_centers:
    dem.add_particle_to_body(center, r, octo)

for i in range(len(edges_poly)):
    dem.add_segment_to_body(
        x_smoothed[(2 * i + 1) % len(x_smoothed)],
        x_smoothed[(2 * i + 2) % len(x_smoothed)],
        octo,
        "line",
    )

Boundaries

l = 10 * w
bnd = dem.add_body((0, 0), 0, 0)

dem.add_segment_to_body(np.array([+l, +y_min]), np.array([-l, +y_min]), bnd, "line")
dem.add_segment_to_body(
    np.array([+l, +4 * y_min]), np.array([-l, +4 * y_min]), bnd, "line"
)
dem.add_segment_to_body(
    np.array([+l, -2 * y_min]), np.array([-l, -2 * y_min]), bnd, "line"
)

Disk

disk = dem.add_body((0, 4 * y_min + w), 1 / mass, 1 / inertia)
dem.add_particle_to_body((0, 0), w, disk)

Initial conditions

nonzero = dem.body_invert_mass().reshape(-1) != 0
dem.body_velocity()[nonzero] = np.array([1.0, 0.0])

dem.set_friction_coefficient(mu)

Time integration

i = 0
t = 0.0
dt = 1e-1
t_end = 5.0
outf = 1

particles_forces = np.zeros_like(dem.velocity())
body_forces = np.zeros_like(dem.body_velocity())
body_forces[nonzero] = mass * g

while t < t_end:
    print(f"Writing iteration {i} at time {t:.4f}")

    if i % outf == 0:
        dem.write_mig(outputdir, t)

    time_integration._advance_particles(
        dem,
        particles_forces,
        dt,
        1,
        tol,
        fbody=body_forces,
    )

    i += 1
    t += dt