Download this testcase.

Study of the drag on a fixed cylinder in a 3D granular flow

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

Keywords

DEM, 3D, Drag, 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 and the mesh is generated.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
geomesh_filename = "3d_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", "-3", geomesh_filename, "-o", mesh_filename])

Geometrical Parameters

r = 9e-3  # particle radius [m]
h, w, d = 1.0, 0.28, 0.05  # domain height, width, and depth [m]
h0 = -0.5  # initial bed offset
r_outer = 0.025  # fixed cylinder radius [m]

Physical Parameters

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

Particle Problem Initialization

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

Particle Initialization

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

    Parameters
    ----------
    origin : list(float)
        Origin of the bottom-left-front corner.
    w : float
        Width of the domain [m].
    h : float
        Height of the domain [m].
    d : float
        Depth of the domain [m].
    step : float
        Particle radius (spacing control).

    Returns
    -------
    tuple of np.ndarray
        Arrays of x, y, and z 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]
    z = np.arange(step, d - step, 2 * step) + origin[2]
    x, y, z = np.meshgrid(x, y, z)
    return x.ravel(), y.ravel(), z.ravel()


# Fill the domain around the fixed cylinder
x, y, z = gen_rect([0, h0, 0], w, 1.5 * h, d, r)
for xi, yi, zi in zip(x, y, z):
    if xi**2 + yi**2 > (r_outer + r) ** 2:
        p.add_particle((xi, yi, zi), r, (4 / 3) * np.pi * r**3 * rhop, "Sand")

# Add a reservoir of particles above the bed
x, y, z = gen_rect([0, -2 * h0, 0], 2 * w, 0.5 * h, d, r)
for xi, yi, zi in zip(x, y, z):
    p.add_particle((xi, yi, zi), r, (4 / 3) * np.pi * r**3 * rhop, "Sand")

p.write_mig(outputdir, 0)

Numerical Parameters

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

Helper Function to Accumulate Drag Forces

def accumulate(F, n_divide):
    """Accumulate total force on the inner fixed cylinder."""
    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 column height drops
    nonzero = p.r()[:, 0] != 0
    position = p.position()[nonzero, :]
    if np.amax(position[:, 1]) + r < 1:
        for xi, yi, zi in zip(x, y, z):
            p.add_particle((xi, yi, zi), r, (4 / 3) * np.pi * r**3 * 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 = (4 / 3) * np.pi * p.r() ** 3 * rhop
    F = np.zeros(3)
    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]};{F[2]};{t}\n")
    # Advance time
    t += dt
    i += 1