Download this testcase.
Bidimensional body Sedimentation (Heat Transfer)
This example illustrates the sedimentation of hot dense solid particles in a viscous fluid under gravity in two dimensions using MigFlow. Particles are initialized in a dense hexagonal packing.
Keywords
DEM, FEM, Generation, heat transfer
Description
The test demonstrates: - How to generate a 2D rectangular fluid domain with named boundaries in Gmsh. - How to initialize a dense cloud of circular particles. - How to couple the particle solver(scontact) and the fluid solver(fluid). - How to set up and solve a thermal advection - diffusion problem(advdiff)
import os, sys, shutil
import numpy as np
import gmsh
from migflow import fluid, scontact, time_integration, advdiff, gmsh_io
gmsh.initialize()
Output Directory
Create a clean output directory for simulation results.
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("outputdir", nargs="?", default="output_2d_depot_hex_heat")
parser.add_argument("--tend", type=float, default=5,
help="final simulation time [s]")
args = parser.parse_args()
outputdir = args.outputdir
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
Geometrical parameters and mesh generation
A rectangular 2D mesh is generated using Gmsh with named boundaries.
height = 0.6 # domain height [m]
width = 0.4 # domain width [m]
mesh_size = 0.005 # element size [m]
h = 0.15 # particle bed height [m]
w = 0.4 # particle bed width [m]
origin = [-width / 2, -height / 2] # mesh origin (bottom-left)
eps = 1e-8 # small numerical tolerance
def gen_mesh(width, height, mesh_size, origin=np.array([0, 0])):
"""Generate a rectangular 2D mesh with physical boundaries."""
origin = np.asarray(origin)
gmsh.model.add("box")
gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height)
gmsh.model.occ.synchronize()
def get_line(x0, x1, eps=1e-6):
r = gmsh.model.get_entities_in_bounding_box(
x0[0] - eps, x0[1] - eps, -eps, x1[0] + eps, x1[1] + eps, eps, 1
)
return [tag for dim, tag in r]
h, w = height, width
gmsh.model.add_physical_group(
1, get_line(origin + [0, 0], origin + [w, 0]), name="Bottom"
)
gmsh.model.add_physical_group(
1, get_line(origin + [0, h], origin + [w, h]), name="Top"
)
gmsh.model.add_physical_group(
1, get_line(origin + [0, 0], origin + [0, h]), name="Left"
)
gmsh.model.add_physical_group(
1, get_line(origin + [w, 0], origin + [w, h]), name="Right"
)
gmsh.model.add_physical_group(2, [1], name="domain")
gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: mesh_size)
gmsh.model.mesh.generate(2)
gen_mesh(width, height, mesh_size, origin)
Physical Parameters
g = np.array([0, -9.81]) # gravity [m/s²]
r = 1.5e-3 # particle radius [m]
rhop = 8000 # steel density [kg/m³]
rho = 1000 # fluid density [kg/m³]
nu = 1e-6 # kinematic viscosity [m²/s]
mu = rho * nu # dynamic viscosity [Pa·s]
T0 = 20 # reference temperature [°C]
k = 0.6 # fluid thermal conductivity [W/m/K]
cp = 4.18e3 # fluid heat capacity [J/kg/K]
beta = 3e-4 # fluid thermal dilatation [/]
Tp = 50 # initial particle temperature [°C]
cpp = 1000 # particle heat capacity [J/kg/K]
kpp = 20 # particle thermal conductivity [W/m/K]
Particle Problem Initialization
Particles are placed in a hexagonal patch at the top of the domain.
p = scontact.ParticleProblem(2)
p.set_fixed_contact_geometry(0)
gmsh_io.load_msh_boundaries(p, None, ["Top", "Left", "Right", "Bottom"])
# Generate a hexagonal dense packing of particles
step = 2 * r + eps
lx, ly = w, h
y0 = height / 2 - r - eps
x = np.arange(-lx / 2, lx / 2 - 2 * r, step)
y = np.arange(y0, y0 - h, -np.sqrt(3) * step)
x2 = np.arange(-lx / 2 + r, lx / 2 - 2 * r, step)
y2 = np.arange(y0, y0 - h, -np.sqrt(3) * step) - np.sqrt(3) * step / 2
x, y = np.meshgrid(x, y)
x2, y2 = np.meshgrid(x2, y2)
x = np.concatenate([x.ravel(), x2.ravel()])
y = np.concatenate([y.ravel(), y2.ravel()])
x_off = w - ((x.max() + step / 2) - (x.min() - step / 2))
x += r + x_off / 2
# Sort particles by height for body ordering, for better visualization
order = np.argsort(y)
x, y = x[order], y[order]
for xi, yi in zip(x, y):
p.add_particle((xi, yi), r, np.pi * r**2 * rhop)
Fluid Problem Initialization
f = fluid.FluidProblem2(g, mu, rho, density_element="triangle_p1")
gmsh_io.load_msh(f, None)
for wall in ["Bottom", "Left", "Right", "Top"]:
f.set_wall_boundary(wall, velocity=[0, 0])
f.set_mean_pressure(0)
f.set_strong_boundary("Left", velocity=[0, 0])
f.set_strong_boundary("Right", velocity=[0, 0])
Advection - Diffusion Problem Initialization
c = advdiff.AdvDiffProblem2(
f=0,
k=k,
cp=cp,
rho=rho,
mu=mu,
velocity_element="triangle_p1",
density_element="triangle_p1",
)
gmsh_io.load_msh(c, None)
for wall in ["Bottom", "Left", "Right", "Top"]:
c.set_strong_boundary(wall, T0)
c.solution().fill(T0)
Output Fields Function
def get_fields(fluid, advdiff=None):
"""Return derived output fields for visualization."""
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),
"density": (fluid.density().get(), p1_element),
"temperature": (advdiff.solution().get().reshape(-1, 1), p1_element),
}
Simulation Parameters
outf = 10 # number of iterations between outputs
dt = 1e-3 # time step [s] -- the steel bed (rho_p = 8000) lands in ~1 ms; dt = 5e-3 cannot resolve the impact (six configurations failed at t ~ 1.2, dt = 1e-3 rides through it)
tEnd = args.tend # final time [s] (--tend)
t = 0
i = 0
Simulation Loop
Time integration of coupled fluid–particle motion.
mass = np.pi * p.r() ** 2 * rhop
Tp = np.full_like((p.volume()), Tp, dtype=np.float64) # initial particle temperature
cpp = np.full_like((p.volume()), cpp * rhop, dtype=np.float64) # steel heat capacity
kpp = np.full_like((p.volume()), kpp, dtype=np.float64) # steel thermal conductivity
while t < tEnd:
print(f"{i:4d}, {t:.6g}/{tEnd:.6g}", flush=True)
# thermal step on the lagged fluid state
c.set_particles(p.volume(), p.position(), p.velocity(), cpp * p.volume(), Tp)
c.velocity().set_from_host(f.velocity())
# Boussinesq with the temperature CLAMPED to the physical bounds [T0, Tp0]
# in the density formula only: the stabilized-P1 T field over/undershoots
# at fronts ([14, 79] from a [20, 50] problem at this mesh), and feeding
# those wiggles into the density sets the buoyancy forcing from a
# non-physical temperature. Clamping bounds the forcing amplitude; the
# solver and the transported T field are untouched.
T_for_rho = np.clip(c.solution().get().reshape(-1, 1), T0, 50.0)
f.density().set_from_host(rho * (1 - beta * (T_for_rho - T0)))
if i % outf == 0:
f.write_mig(outputdir, t, get_fields(f, c))
p.write_mig(outputdir, t, {"temperature": Tp.reshape(-1, 1)})
c.implicit_euler(dt)
# Coupled step: real contact forces in the prediction (library default),
# with the DEM settings the ~1 ms landing needs (tol 5e-4*r, max_nsub 200).
# Gate: at t = 20 the bed is settled and the two schemes agree (patankar
# 1.1e-3 / iqn 2.1e-2 m/s residual velocities).
# mu/rho_f are passed EXPLICITLY here because the density field is
# Boussinesq (set from T every step): left out, the coupling closure would
# read the field's mean, and the drag law is meant to see the reference
# density, not a temperature average of it.
time_integration.iterate_patankar(f, p, dt, 1, 5e-4 * r,
external_particles_forces=mass * g,
max_nsub=200, mu=mu, rho_f=rho)
nonzero = p.volume() > 0
Tp[nonzero] += (
np.asarray(c.get_flux_on_bodies()).ravel()[nonzero.ravel()] * dt / (cpp[nonzero] * p.volume()[nonzero])
)
t += dt
i += 1
Plot
python3 -m migflow.plot.migplot output_2d_depot_hex_heat --actors fluid particles --fluid-field pressure --element-type triangle_p1