Download this testcase.

Study of the drag on a fixed disk in a 2D granular flow

This example illustrates how to compute the drag force exerted on a fixed circular obstacle immersed in a granular flow. The particles fall under gravity between two lateral walls and impact the fixed inner cylinder.

Keywords

DEM, Generation, Post-processing

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

Output Directory

The output directory is (re)created

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
geomesh_filename = "2d_mesh.geo"
mesh_filename = f"{outputdir}/mesh.msh"
csv_file = f"{outputdir}/Drag.csv"
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", geomesh_filename, "-o", mesh_filename])

Geometrical parameters

r = 6e-3  # particle radius [m]
h, w = 1.0, 0.28  # particle region height & width [m]
h0 = -0.5  # particle bed offset
r_outer = 0.025  # radius of the fixed disk [m]

Physical Parameters

rhop = 2500  # particle density [kg/m³]
g = np.array([0, -9.81])  # gravity [m/s²]
friction = 0.2  # friction coefficient

Particle Problem Initialization

Particle problem setup

Load mesh boundaries (requires a mesh.msh file)

p = scontact.ParticleProblem(2)
p.load_msh_boundaries(mesh_filename, ["Inner"], material="Steel")
p.load_msh_boundaries(
    mesh_filename,
    ["Right", "Left", "RightUp", "RightDown", "LeftUp", "LeftDown"],
    material="Plexi",
)
p.set_friction_coefficient(friction, "Sand", "Sand")
p.set_friction_coefficient(friction, "Sand", "Steel")

Particle initialization

Add granular particles around the fixed disk.

def gen_rect(origin, w, h, step):
    """Generate coordinates on a rectangular grid.

    Parameters
    ----------
    origin : list(float)
        Origin of the top-left corner.
    w : float
        Width of the rectangular domain.
    h : float
        Height of the rectangular domain.
    step : float
        Particle radius (spacing control).

    Returns
    -------
    tuple of np.ndarray
        Arrays of x and y coordinates.
    """
    eps = 1e-8
    x = np.arange(-w / 2 + step - eps, w / 2 - step + eps, 2 * step) + origin[0]
    y = np.arange(step, h - step, 2 * step) + origin[1]
    x, y = np.meshgrid(x, y)
    return x.ravel(), y.ravel()


# Fill the domain around the disk
x, y = gen_rect([0, h0], w, 1.5 * h, r)
for xi, yi in zip(x, y):
    if xi**2 + yi**2 > (r_outer + r) ** 2:
        p.add_particle((xi, yi), r, np.pi * r**2 * rhop, "Sand")
# Add a reservoir of particles above
x, y = gen_rect([0, -2 * h0], 2 * w, 0.5 * h, r)
for xi, yi in zip(x, y):
    p.add_particle((xi, yi), r, np.pi * r**2 * rhop, "Sand")

Numerical Parameters

dt = 1e-3  # time step [s]
tEnd = 2.0  # final time [s]
outf = 10  # output frequency
t, i = 0.0, 0

Helper function to accumulate the drag forces

def accumulate(F, n_divide):
    """Accumulate the total force on the inner fixed disk."""
    F += np.sum(p.get_boundary_forces("Inner"), axis=0) / n_divide

Time integration loop

while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
    # Add new particles if the column becomes too low
    nonzero = p.r()[:, 0] != 0
    position = p.position()[nonzero, :]
    if np.amax(position[:, 1]) + r < 1:
        for xi, yi in zip(x, y):
            p.add_particle((xi, yi), r, np.pi * r**2 * rhop)
    # Remove particles below the lower limit
    p.remove_bodies_flag(p.body_position()[:, 1] > -0.6)
    # Write outputs
    if i % outf == 0:
        p.write_mig(outputdir, t)
    # Compute contact and body forces
    mass = np.pi * p.r() ** 2 * rhop
    F = np.zeros(2)
    time_integration.iterate(
        None,
        p,
        dt,
        1,
        contact_tol=1e-3 * r,
        external_particles_forces=g * mass,
        after_sub_iter=lambda n_divide: accumulate(F, n_divide),
    )
    # Write drag results
    with open(csv_file, "a") as file1:
        file1.write(f"{F[0]};{F[1]};{t}\n")
    # Advance time
    t += dt
    i += 1

Plot

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