Download this testcase.
# MigFlow - Copyright (C) <2010-2026>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
Tridimensional Particle Sedimentation
This example illustrates the sedimentation of dense solid particles in a viscous fluid under gravity in three dimensions using MigFlow. Particles are initialized in a grid packing.
Keywords
DEM, FEM, Generation
Description
The test demonstrates: - How to generate a 3D rectangular fluid domain with named boundaries in Gmsh. - How to couple the particle solver (scontact) and the fluid solver (fluid). - How to export derived fields (pressure, velocity, porosity, dynamic pressure) for post-processing in Paraview or similar tools.
import os, time, shutil, sys
import gmsh
import numpy as np
from migflow import fluid, scontact, time_integration, gmsh_io
Output Directory
Create a clean output directory for simulation results.
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("outputdir", nargs="?", default="output_3d_depot_grid")
parser.add_argument("--tend", type=float, default=2,
help="final simulation time [s]")
args = parser.parse_args()
outputdir = args.outputdir
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
Physical parameters
g = np.array([0, -9.81, 0]) # gravity
r = 2e-3 # particles radius
rho = 1000 # fluid density
rhop = 1500 # particles density
nu = 1e-6 # kinematic viscosity
Atwood = (rhop - rho) / (rhop + rho)
print(f"Atwood number: {Atwood:1.4f}")
Geometrical parameters
w = 0.4 # domain width
h = 0.6 # domain height
d = 10 * r # domain depth
lx = w # Particles width
ly = h / 6 # Particles height
lz = d # Particles depth
lc = d / 2 # mesh size
Mesh generation
gmsh.initialize()
gmsh.model.add("box3d")
gmsh.model.occ.add_box(-w / 2, 0, -d / 2, w, h, d)
gmsh.model.occ.synchronize()
def get_surface(x0, x1, y0, y1, z0, z1, eps=1e-6):
r = gmsh.model.get_entities_in_bounding_box(
x0 - eps, y0 - eps, z0 - eps, x1 + eps, y1 + eps, z1 + eps, 2
)
return [tag for dim, tag in r]
gmsh.model.add_physical_group(
2, get_surface(-w / 2, +w / 2, h, h, -d / 2, +d / 2), name="Top"
)
gmsh.model.add_physical_group(
2, get_surface(-w / 2, +w / 2, 0, 0, -d / 2, +d / 2), name="Bottom"
)
gmsh.model.add_physical_group(
2, get_surface(-w / 2, -w / 2, 0, h, -d / 2, +d / 2), name="Left"
)
gmsh.model.add_physical_group(
2, get_surface(+w / 2, +w / 2, 0, h, -d / 2, +d / 2), name="Right"
)
gmsh.model.add_physical_group(
2, get_surface(-w / 2, +w / 2, 0, h, +d / 2, +d / 2), name="Front"
)
gmsh.model.add_physical_group(
2, get_surface(-w / 2, +w / 2, 0, h, -d / 2, -d / 2), name="Rear"
)
gmsh.model.add_physical_group(3, [1], name="domain")
# --- Uniform mesh control ---
gmsh.model.mesh.set_algorithm(3, 1, 10)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthFromPoints", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthExtendFromBoundary", 0)
gmsh.model.mesh.generate(3)
Function to generate a rectangular grid of particles
def gen_rect3d(origin, w, h, d, step):
"""Generate coordinates on a rectangular grid
Keyword arguments:
origin -- origin of top left corner
w : width
h : height
d : depth
step : maximal radius
"""
eps = 1e-8
x = np.arange(-w / 2 + step - eps, w / 2 - step + eps, 2 * step) + origin[0]
y = np.arange(origin[1] - step - eps, origin[1] - h + step + eps, -2 * step)
z = np.arange(-d / 2 + step - eps, d / 2 - step + eps, 2 * step) + origin[2]
x, y, z = np.meshgrid(x, y, z)
return x.reshape(-1), y.reshape(-1), z.reshape(-1)
Particle Problem Initialization
p = scontact.ParticleProblem(3)
gmsh_io.load_msh_boundaries(p, None)
x, y, z = gen_rect3d([0, h, 0], lx, ly, lz, r)
for xi, yi, zi in zip(x, y, z):
p.add_particle((xi, yi, zi), r, 4 / 3 * r**3 * np.pi * rhop)
Fluid Problem Initialization
f = fluid.FluidProblem3(g, nu * rho, rho)
gmsh_io.load_msh(f, None)
for tag in ["Left", "Right", "Front", "Rear", "Top", "Bottom"]:
f.set_wall_boundary(tag)
f.set_mean_pressure(0)
Function to get derived fields for output
def get_fields(fluid):
"""Return derived output fields for visualization."""
y = fluid.coordinates_fields()[fluid.field_indices(fluid.dimension())][:, 1]
p1_element = fluid.get_p1_element()
return {
"pressure": (fluid.pressure(), p1_element),
"velocity": (fluid.velocity(), p1_element),
"porosity": (fluid.porosity().get(), p1_element),
"u_solid": (fluid.u_solid(), p1_element),
"dynamic_pressure": (fluid.pressure() - rho * g[1] * (y - h / 2), p1_element),
}
Simulation Parameters
outf = 1 # number of iterations between output files
dt = 5e-3 # time step (halved from 1e-2 for the dt/2 rerun)
tEnd = args.tend # final time [s] (--tend)
t = 0
i = 0
Simulation Loop
Time integration of coupled fluid–particle motion.
mass = 4 / 3 * np.pi * p.r() ** 3 * rhop
# The coupled step is one library call: iterate_patankar carries the DEM
# solver's contact forces into the fluid prediction (the library default) and
# advances the grains. The stability requirement that comes with the real
# contact forces is the TIME STEP: dt*Fc/m must stay below the flow scale,
# which is what sets dt here.
while t < tEnd:
print(f"t/tEnd={t:1.4f}/{tEnd:1.4f}, i={i:1d}", flush=True)
if i % outf == 0:
p.write_mig(outputdir, t)
f.write_mig(outputdir, t, get_fields(f))
time_integration.iterate_patankar(f, p, dt, 1, 1e-3 * r,
external_particles_forces=mass * g)
i += 1
t += dt
Plot
python3 -m migflow.plot.migplot output_3d_depot_grid --actors fluid particles --fluid-field porosity