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
Output Directory¶
Create a clean output directory for simulation results.
outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
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)
p.load_msh_boundaries(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.FluidProblem(3, g, nu * rho, rho)
f.load_msh(None)
f.set_wall_boundary(["Left", "Right", "Front", "Rear", "Top", "Bottom"])
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.pressure_index()][:, 1].reshape(-1, 1)
p1_element = fluid.get_p1_element()
return {
"pressure": (fluid.pressure(), p1_element),
"velocity": (fluid.velocity(), p1_element),
"porosity": (fluid.porosity(), 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 = 1e-2 # time step
tEnd = 2 # final time
t = 0
i = 0
Simulation Loop¶
Time integration of coupled fluid–particle motion.
mass = 4 / 3 * np.pi * p.r() ** 3 * rhop
while t < tEnd:
print(f"t/tEnd={t:1.4f}/{tEnd:1.4f}, i={i:1d}")
f.set_particles(
p.delassus(),
p.volume(),
p.position(),
p.velocity(),
p.omega(),
p.contact_forces(),
)
if i % outf == 0:
p.write_mig(outputdir, t)
f.write_mig(outputdir, t, get_fields(f))
f.implicit_euler(dt)
fext = g * mass + f.compute_node_force()
time_integration._advance_particles(p, fext, dt, 1, 1e-3 * r)
i += 1
t += dt