Download this testcase.

Rotating Drum

This example illustrates a circular shear flow generated by a rotating inner cylinder (drum) surrounded by an outer stationary wall. Grains are placed in the annular gap and interact with the fluid flow through drag and pressure force.

Keywords

FEM, DEM, LMGC

Description

The test demonstrates: - How to define moving boundary conditions as functions of position. - How to generate an initial particle packing in a circular domain. - How to couple the particle and fluid solvers in MigFlow. - Optional use of lmgc90 instead of scontact for contact resolution.

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

use_lmgc90 = False
if use_lmgc90:
    from migflow import lmgc90Interface

Output Directory and Mesh Generation

The output directory is created, and a 2D annular mesh is generated using Gmsh.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
initdir = "init"
geomesh_filename = "2d_mesh.geo"
mesh_filename = f"{outputdir}/mesh.msh"

shutil.rmtree(outputdir, ignore_errors=True)
shutil.rmtree(initdir, ignore_errors=True)
os.makedirs(outputdir)
os.makedirs(initdir)

subprocess.call(["gmsh", "-2", geomesh_filename, "-o", mesh_filename])

Initial Particle Generation

Function to create an initial packing of grains within the annular domain.

def genInitialPosition(initdir, r, rout, rin, rhop, use_lmgc90):
    """Generate a particle packing and write it to an initial output file.

    Parameters
    ----------
    initdir : str
        Directory where the initial state is written.
    r : float
        Maximum particle radius.
    rout : float
        Outer radius of the drum.
    rin : float
        Inner radius of the drum.
    rhop : float
        Particle density.
    use_lmgc90 : bool
        If True, uses LMGC90 instead of scontact for contacts.
    """
    # Create particle problem
    p = scontact.ParticleProblem(2)
    print(mesh_filename)
    p.load_msh_boundaries(mesh_filename, ["Outer", "Inner"], material="Steel")

    # Define grid of possible particle centers
    x = np.arange(rout, -rout, 2.5 * -r)
    x, y = np.meshgrid(x, x)
    R2 = x**2 + y**2

    # Keep only points inside the annular region
    keep = (R2 < (rout - r) ** 2) & (R2 > (rin + r) ** 2)
    x = x[keep]
    y = y[keep]

    # Add particles with slight density variation
    for xi, yi in zip(x, y):
        if yi < rout:
            rhop1 = random.choice([rhop * 0.9, 1.1 * rhop, rhop])
            p.add_particle((xi, yi), r, r**2 * np.pi * rhop1, "Sand")

    p.write_mig(initdir, 0)

Physical and Numerical Parameters

g = np.array([0, 0])  # no gravity
rho = 1.253e3  # fluid density [kg/m³]
rhop = 1000  # particle density [kg/m³]
nu = 1e-3  # kinematic viscosity [m²/s]
mu = rho * nu  # dynamic viscosity [Pa·s]
rout = 0.0254  # outer radius [m]
rin = 0.0064  # inner radius [m]
r = 397e-6 / 2  # grain radius [m]

Particle Problem Initialization

genInitialPosition(initdir, r, rout, rin, rhop, use_lmgc90)

if use_lmgc90:
    friction = 0.1
    lmgc90Interface.scontactTolmgc90(initdir, 2, 0, friction)
    p = lmgc90Interface.ParticleProblem(2)
else:
    p = scontact.ParticleProblem(2)
    p.read_mig(initdir, 0)
    p.set_friction_coefficient(0.1, "Sand", "Sand")
    p.set_friction_coefficient(0.1, "Sand", "Steel")
p.set_fixed_contact_geometry(0)

Fluid Problem Initialization

The fluid solver is initialized, with a rotating inner boundary and fixed outer wall.

f = fluid.FluidProblem(2, g, nu * rho, rho)
f.load_msh(mesh_filename)
f.set_mean_pressure(0)
f.set_wall_boundary("Outer", velocity=[0, 0])
f.set_strong_boundary(
    "Inner",
    velocity=[lambda x: (x[:, 1] / rout) * 0.1, lambda x: (-x[:, 0] / rout) * 0.1],
)

# Coupling with particles
f.set_particles(
    p.delassus(),
    p.volume(),
    p.position(),
    p.velocity(),
    p.omega(),
    p.contact_forces(),
)

Simulation Parameters

dt = 5e-3  # time step [s]
tEnd = 10  # total simulation time [s]
outf = 10  # number of iterations between outputs
t = 0  # time [s]
i = 0  # iteration counter

Computational Loop

mass = np.pi * p.r() ** 2 * rhop  # particle masses
while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
    if i % outf == 0:
        p.write_mig(outputdir, t)
        f.write_mig(outputdir, t)
    time_integration.iterate(
        f,
        p,
        dt,
        use_predictor_corrector=False,
        external_particles_forces=g * mass,
        contact_tol=1e-3 * r,
    )
    t += dt
    i += 1

Plot

python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field velocity