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, gmsh_io, volume_coupling as vc
np.random.seed(0)
Output Directory
Create an output directory for results and generate a 2D mesh using Gmsh.
outputdir = "output_2d_avalanche_unresolved" 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()
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()
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)
gmsh_io.load_msh_boundaries(p, 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)
first = p.n_bodies()
for xi, ri in zip(p1.position(), p1.r()):
p.add_particle(xi, ri, np.pi * ri**2 * rhop)
p.body_velocity()[first:] = 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.FluidProblem2(g, mu, rho)
gmsh_io.load_msh(f, 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)
# Tier A (body geometry on the mesh) and tier B (the legacy mixture closure)
# only: this primes the coupling so the first frame carries a porosity field,
# and time_integration.iterate re-sets all three tiers -- datum included -- at
# every step.
geo = vc.get_particles(f, p.position(), p.r(), density=vc.body_density(p))
clo = vc.get_particles_closure(f, geo, mu, rho, p.velocity(),
omega=p.omega(),
contact_forces=p.contact_forces())
vc.set_coupling_geometry(f, geo)
vc.set_coupling_closure(f, clo)
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_2d_avalanche_unresolved --actors fluid particles --fluid-field pressure