Download this testcase.

Granular Column Collapse in Fluid¶

This test case simulates the collapse of a granular column in a viscous fluid. The column of grains can either be imported from a pre-deposited configuration (see ../depot/depot.py) or generated automatically using a regular hexagonal packing. The case demonstrates the coupling between a dense granular phase and a single-phase incompressible fluid.

Keywords¶

FEM, DEM, Unresolved, Friction

Description¶

The granular column initially rests in hydrostatic equilibrium inside a rectangular box filled with water. Upon release, the grains collapse and form a dense avalanche driven by gravity and resisted by viscous drag.

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

np.random.seed(0)

Output Directory¶

Create an output directory for results and generate a 2D mesh using Gmsh.

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

Hexagonal Packing Function¶

Define a helper function to create a regular hexagonal packing of grains when no pre-deposited configuration is available.

use_pre_deposit = False


def hexagonal_packing(p, x0, y0, lx, ly, rmin, rmax):
    """Create a hexagonal particle arrangement inside a rectangular region.

    Parameters
    ----------
    p : migflow.scontact.ParticleProblem
        Particle problem instance.
    x0, y0 : float
        Bottom-left corner of the region.
    lx, ly : float
        Width and height of the region.
    rmin, rmax : float
        Minimum and maximum particle radii.
    """
    for x in np.arange(x0 + rmax, x0 + lx - rmax, 2 * rmax):
        for y in np.arange(y0 + rmax, y0 + ly - rmax, 2 * (3**0.5) * rmax):
            z = np.random.random((1))
            r = rmin + (rmax - rmin) * z
            p.add_particle((x, y), r, np.pi * r**2 * rhop)
    for x in np.arange(2 * rmax + x0, x0 + lx - rmax, 2 * rmax):
        for y in np.arange(
            y0 + (3**0.5) * rmax + rmax, y0 + ly - rmax, 2 * (3**0.5) * rmax
        ):
            z = np.random.random((1))
            r = rmin + (rmax - rmin) * z
            p.add_particle((x, y), r, np.pi * r**2 * rhop)

Physical Parameters¶

The simulation parameters are defined below, including fluid viscosity, densities, gravity, and the typical particle radius.

g = np.array([0, -9.81])  # gravity [m/s²]
rhop = 1500  # particle density [kg/m³]
rho = 1000  # fluid density [kg/m³]
nu = 1e-6  # fluid kinematic viscosity [m²/s]
mu = rho * nu  # fluid dynamic viscosity [Pa·s]
r = 2e-3  # mean particle radius [m]
h = 0.2  # domain height [m]
lx = 0.1  # granular column width [m]
ly = 0.165  # granular column height [m]

Particle Problem¶

The particle problem is created and the boundary geometry is loaded. Depending on use_pre_deposit, the initial configuration is either imported from a previous depot testcase or generated on the fly.

p = scontact.ParticleProblem(2)
p.set_fixed_contact_geometry(0)
p.load_msh_boundaries(mesh_filename, ["Top", "Bottom", "Left", "Right"])

if use_pre_deposit:
    # Load compact column from depot testcase
    p1 = scontact.ParticleProblem(2)
    p1.read_mig("../depot/output", 10)
    p.add_particles(p1.position(), p1.r(), p1.mass(), v=p1.velocity())
else:
    # Generate hexagonal packing in the left side of the domain
    eps = 1e-4 * r
    hexagonal_packing(p, eps, eps, lx, ly, 0.9 * r, 1.1 * r)

# Set particle–particle and particle–wall friction
p.set_friction_coefficient(0.5)
p.write_mig(outputdir, 0)

Fluid Problem¶

The fluid problem is initialized using the same geometry and coupled to the particle phase. All walls are considered no-slip, and the mean pressure is set to zero.

f = fluid.FluidProblem(2, g, mu, rho)
f.load_msh(mesh_filename)
f.set_wall_boundary("Bottom")
f.set_wall_boundary("Left")
f.set_wall_boundary("Right")
f.set_wall_boundary("Top")
f.set_mean_pressure(0)
f.set_particles(
    p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), p.contact_forces()
)

Time Integration¶

The simulation runs for a short time interval with a fixed time step. At each iteration, the fluid phase is solved and the fluid-grain forces are applied to the particles, and then system evolves according to the contact dynamics solver. The simulation runs until the avalanche stabilizes.

i = 0
outf = 10  # number of iterations between outputs
dt = 1e-3  # time step [s]
t = 0  # initial time [s]
tEnd = 5  # final time [s]
mass = np.pi * p.r() ** 2 * rhop
tic = time.time()
while t < tEnd:
    print(f"{i = :4d} ----- {t = :g}")
    if i % outf == 0:
        p.write_mig(outputdir, t)
        f.write_mig(outputdir, t)
    time_integration.iterate(
        f, p, dt, min_nsub=1, external_particles_forces=g * mass, contact_tol=1e-4 * r
    )
    t += dt
    i += 1

Plot¶

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