Download this testcase.

Drag on a Fixed Disk in a 2D Immersed Granular Flow

This example studies the drag force acting on a fixed circular obstacle (a “disk”) immersed in a granular flow saturated by a viscous fluid.

Keywords

Drag, DEM, FEM, Boundary Force, 2D-3D Ratio

Description

The test demonstrates: - How to set up a mixed granular–fluid simulation in a 2D geometry. - How to impose an inflow boundary condition to simulate a stream. - How to insert, remove, and constrain particles dynamically. - How to compute and record the total drag force acting on a solid obstacle. A 2D–3D scaling factor is applied to account for the reduced dimensionality.

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

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

h_p = 0.3  # particle bed height [m]
w_p = 0.3  # particle bed width [m]
w = 0.3  # particle bed width [m]
h = 1.4  # total domain height [m]
r_inner = 0.025  # radius of the fixed disk [m]
ratio_2d_3d = 0.7  # 2D-3D correction

Physical Parameters

v_imposed = -0.05  # stream velocity
rhop = 2500  # particle density
rho = 1000  # fluid density
r = 3e-3  # particle radius
g = np.array([0, -9.81])  # gravity
nu = 1e-6  # fluid kinematic viscosity
mu = rho * nu
friction = 0.2

Particle Problem Initialization

Load mesh boundaries and define contact properties.

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

Particle Generation

The following function generates particle coordinates on a regular grid over a rectangular region.

def gen_rect(origin, w, h, step):
    """Generate coordinates on a rectangular grid
    Keyword arguments:
        origin -- origin of top left corner
        w : width
        h : height
        step : maximal radius
    """
    eps = 1e-4 * step
    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.reshape(-1), y.reshape(-1)


x, y = gen_rect([0, -h / 2], w, 4 * h_p, r)
for xi, yi in zip(x, y):
    if xi**2 + yi**2 > (r_inner + r) ** 2:
        p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand")
x, y = gen_rect([0, 0.5 * h / 1.4], w, 0.5 * h_p, 1.2 * r)
for xi, yi in zip(x, y):
    p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand")

Fluid Problem Initialization

The fluid problem is loaded from the same mesh and coupled to the particle phase.

f = fluid.FluidProblem(2, g, nu * rho, rho)
f.load_msh(mesh_filename)
f.set_wall_boundary("Inner")
f.set_wall_boundary("Left", velocity=[0, v_imposed])
f.set_wall_boundary("Right", velocity=[0, v_imposed])
f.set_open_boundary("Bottom", velocity=[0, v_imposed])
f.set_strong_boundary("Left", velocity=[0, v_imposed])
f.set_strong_boundary("Right", velocity=[0, v_imposed])
f.set_strong_boundary("Bottom", velocity=[0, v_imposed])
f.set_open_boundary("Top", pressure=0)

Numerical Parameters

dt = 5e-3  # time step
t = 0  # initial time
tEnd = 1.5  # final time
i = 0  # iteration number
outf = 2  # iterations between data frames

Computation Loop

Function to accumulate the force applied to the inner disk by the contact networks.

def accumulate(F, n_divide):
    F += np.sum(p.get_boundary_forces("Inner"), axis=0) / n_divide


while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
    # Remove particles that have fallen outside the domain
    alert = 4 * r
    keep = (p.body_position()[:, 1] > -h / 1.4 * 0.55) & (
        p.body_position()[:, 1] < h / 2 - alert
    )
    p.remove_bodies_flag(keep)

    # Constrain particles at the bottom
    nonzero = p.r()[:, 0] != 0
    position = p.position()[nonzero, :]
    p.body_invert_mass()[p.body_position()[:, 1] < -0.5] = 0
    p.body_invert_inertia()[p.body_position()[:, 1] < -0.5] = 0
    a = p.body_velocity()[p.n_bodies() - position.shape[0] :, :]
    b = p.body_omega()[p.n_bodies() - position.shape[0] :, :]
    c = p.body_invert_mass()[p.n_bodies() - position.shape[0] :, :]
    a[(c == 0)[:, 0], :] = [0, v_imposed]
    b[(c == 0)[:, 0], 0] = 0

    # Insert new particles at the top when needed
    if np.amax(position[:, 1]) + r < 0.5:
        for xi, yi in zip(x, y):
            p.add_particle((xi, yi), r, r**2 * np.pi * rhop, "Sand")

    # Set particles to the FEM module
    f.set_particles(
        p.delassus(),
        p.volume() * ratio_2d_3d,
        p.position(),
        p.velocity(),
        p.omega(),
        p.contact_forces(),
    )

    # Write output files
    if i % outf == 0:
        p.write_mig(outputdir, t)
        f.write_mig(outputdir, t)

    # Iterate
    F = np.zeros(2)  # to store inner disk forces
    mass = np.pi * p.r() ** 2 * rhop
    time_integration.iterate(
        f,
        p,
        dt,
        10,
        contact_tol=1e-3 * r,
        external_particles_forces=g * mass,
        after_sub_iter=lambda n_divide: accumulate(F, n_divide),
        use_predictor_corrector=False,
    )

    with open(csv_file, "a") as file1:
        F += f.boundary_force()[0, :]
        file1.write(str(F[0]) + ";" + str(F[1]) + ";" + str(t) + "\n")
    t += dt
    i += 1

Plot

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