Download this testcase.
2D Bubbling Fluidized Bed¶
This example simulates a two-dimensional bubbling fluidized bed in which solid granular particles are suspended by an upward viscous fluid flow. Particles are injected through a narrow nozzle at the base of the domain and interact with the fluid through unresolved drag coupling.
Keywords¶
FEM, DEM, Fluidized Bed, Unresolved
import os, sys, shutil
import gmsh
import numpy as np
from migflow import advdiff, 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)
Geometrical parameters and mesh generation¶
A rectangular 2D mesh is generated using Gmsh with named boundaries.
cm = 0.01
np.random.seed(0)
def gen_mesh(h, w, n, lc, o):
o = np.asarray(o)
gmsh.model.add("box")
p0 = gmsh.model.occ.add_point(o[0], o[1], 0, lc)
p1 = gmsh.model.occ.add_point(o[0] + (w - n) / 2, o[1], 0, lc)
p2 = gmsh.model.occ.add_point(o[0] + (w + n) / 2, o[1], 0, lc)
p3 = gmsh.model.occ.add_point(o[0] + w, o[1], 0, lc)
p4 = gmsh.model.occ.add_point(o[0] + w, o[1] + h, 0, lc)
p5 = gmsh.model.occ.add_point(o[0], o[1] + h, 0, lc)
l0 = gmsh.model.occ.add_line(p0, p1)
l1 = gmsh.model.occ.add_line(p1, p2)
l2 = gmsh.model.occ.add_line(p2, p3)
l3 = gmsh.model.occ.add_line(p3, p4)
l4 = gmsh.model.occ.add_line(p4, p5)
l5 = gmsh.model.occ.add_line(p5, p0)
cloop = gmsh.model.occ.add_curve_loop([l0, l1, l2, l3, l4, l5])
gmsh.model.occ.add_plane_surface([cloop])
gmsh.model.occ.synchronize()
gmsh.model.add_physical_group(1, [l0, l2], name="Bottom")
gmsh.model.add_physical_group(1, [l4], name="Top")
gmsh.model.add_physical_group(1, [l3, l5], name="Wall")
gmsh.model.add_physical_group(1, [l1], name="Nozzle")
gmsh.option.setNumber("Mesh.Algorithm", 2)
gmsh.model.mesh.generate(2)
index_v = 2
index_m = 0
# ========= Input parameters =========
v_b = [1.2, 1.54, 1.71][index_v] # bottom velocity [m/s]
mass = [75e-3, 125e-3][index_m] # mass of the particles [kg]
friction = 0.3 # friction coefficient [/]
print(f"{v_b = }, {mass = }, {friction = }")
# ========= Geometry =========
height = 25 * cm # domain height [m]
width = 8 * cm # domain width [m]
depth = 1.5 * cm # domain depth [m]
nozzle = 1.3 * cm # nozzle width [m]
origin = [-width / 2, 0] # leftmost bottom corner [m]
# lc = 0.3*cm # mesh size [m]
lc = 0.15 * cm # mesh size [m]
gen_mesh(height, width, nozzle, lc, origin)
# ========= Fluid properties =========
g = np.array([0, -9.81]) # gravity [m/s**2]
rho = 1.204 # nitrogen density [kg/m**3]
mu = 2.0e-5 # nitrogen dynamic viscosity [kg/(ms)]
nu = mu / rho # nitrogen kinematic viscosity [m**2/s]
k = 2e-2 # nitrogen thermal conductivity [W/(m K)]
cp = 1010 # nitrogen heat capacity [J/(Kg K)]
beta = 3e-3 # nitrogen thermal dilatation [/]
u_b, v_b = 0, v_b # bottom velocity [m/s]
# ========= Grains properties =========
d = 0.1 * cm # glass diameter [m]
r = d / 2 # glass radii [m]
rhog = 2500 # glass density [kg/m**3]
cpg = 840 # glass heat capacity [J/(kg K)]
kg = 1.4 # glass thermal conductivity [W/(m K)]
young_modulus = 6e7 # glass young modulus [kg/(m s**2)]
poisson_ratio = 0.22 # glass poisson ratio [/]
T0 = -20 # offset temperature [K]
Tf0 = T0 + 20 # initial fluid temperature [K]
Tp0 = T0 + 90 # initial particle temperature [K]
hw = 350 # wall convective heat transfer [W/(m**2 K)]
volume_ratio = 0.82
# ========= Particle Problem =========
p = scontact.ParticleProblem(2)
volume = mass / rhog
n_p3d = int(volume / (4 / 3 * np.pi * r**3))
n_p1d = int(depth / d)
n_p = int(n_p3d / n_p1d)
h = 1.0
read_input = False
eps = r
# Generate a hexagonal dense packing of particles
step = 2 * r + eps
lx, ly = width, height
x = np.arange(-lx / 2, lx / 2 - 2 * r, step)
y = np.arange(r + eps, h - r - eps, np.sqrt(3) * step)
x2 = np.arange(-lx / 2 + r, lx / 2 - 2 * r, step)
y2 = np.arange(r + eps, h - r - eps, 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()])
dx = x.max() - x.min()
offset_x = (width - dx) / 2
offset_y = y.min() - origin[1] - r
x += offset_x
y -= offset_y
# Sort particles by height for body ordering, for better visualization
order = np.argsort(y)
x, y = x[order], y[order]
x = x[:n_p]
y = y[:n_p]
for xi, yi in zip(x, y):
p.add_particle((xi, yi), r, np.pi * r**2 * rhog, material="particle")
bl = np.array([origin[0], origin[1]])
br = np.array([origin[0] + width, origin[1]])
tl = np.array([origin[0], origin[1] + height])
tr = np.array([origin[0] + width, origin[1] + height])
p.add_boundary_segment(bl, br, "bottom", "wall")
p.add_boundary_segment(bl, tl, "left", "wall")
p.add_boundary_segment(tl, tr, "top", "wall")
p.add_boundary_segment(br, tr, "right", "wall")
p.set_friction_coefficient(1.0, "particle", "wall")
p.set_friction_coefficient(friction, "particle", "particle")
mcp = cpg / p.delassus()[:, 0, 0]
ks = np.full(p.n_particles(), kg)
e = np.full(p.n_particles(), young_modulus)
v = np.full(p.n_particles(), poisson_ratio)
# ========= Fluid Problem =========
f = fluid.FluidProblem(
2, g, mu, rho, volume_drag=12 * mu / (depth**2), density_element=b"triangle_p1"
)
f.load_msh(None)
f.set_open_boundary("Bottom", velocity=[u_b, v_b])
f.set_open_boundary("Top", pressure=0)
f.set_wall_boundary("Nozzle")
f.set_wall_boundary("Wall")
f.set_strong_boundary("Bottom", velocity=[u_b, v_b])
f.set_strong_boundary("Nozzle", velocity=[0, 0])
Advection-Diffusion Problem Initialization¶
hp_mean = 2.0 * hw / depth
c = advdiff.AdvDiffProblem(
2, 0, k, cp * rho, rho, mu, velocity_element=b"triangle_p1", volumetric_drag=hp_mean
)
c.load_msh(None)
c.set_open_boundary("Top")
c.set_strong_boundary("Bottom", solution=Tf0)
c.set_neumann_boundary("Nozzle", solution=Tf0, h_w=hw)
c.set_neumann_boundary("Wall", solution=Tf0, h_w=hw)
c.solution()[:] = Tf0
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(), p1_element),
"u_solid": (fluid.u_solid(), p1_element),
"density": (fluid.density(), p1_element),
"temperature": (advdiff.solution().reshape(-1, 1), p1_element),
}
print(f"Start simulation with input parameters : {v_b=} ---- {mass=} ---- {friction=}")
Simulation Loop¶
Time integration of coupled fluid–particle motion. ========= Numerical parameters =========
t = 0
i = 0
dt_fluid = 1e-3
dt_heat = 1e-3
n_substeps = int(dt_heat / dt_fluid)
dt = dt_heat
tEnd = 1.0
outf = 10
# ========= Time integration =========
mass = np.pi * p.r() ** 2 * rhog
Tp = np.full_like((p.volume()), Tp0, dtype=np.float64) # initial particle temperature
cpp = np.full_like(
(p.volume()), cpg * rhog, dtype=np.float64
) # heat capacity times density
kpp = np.full_like((p.volume()), kg, dtype=np.float64) # thermal conductivity (unused)
nonzero = p.volume() > 0
while t < tEnd:
print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
c.set_particles(
volume_ratio * p.volume(), p.position(), p.velocity(), cpp * p.volume(), kpp, Tp
)
f.set_particles(
p.delassus(),
volume_ratio * p.volume(),
p.position(),
p.velocity(),
p.omega(),
p.contact_forces(),
)
if i % outf == 0:
f.write_mig(outputdir, t, get_fields(f, c))
p.write_mig(outputdir, t, {"temperature": Tp.reshape(-1, 1)})
c.velocity()[:] = f.velocity().copy()
c.implicit_euler(dt_heat)
Tp[nonzero] += (
c.particles_heat()[nonzero] * dt / (cpp[nonzero] * p.volume()[nonzero])
)
dT = c.solution().reshape(-1, 1) - Tf0
f.density()[:] = rho * (1 - np.fmin(beta * dT, 0.3))
# Solve fluid-particle coupling
f.implicit_euler(dt_fluid)
fext = g * mass + f.compute_node_force()
time_integration._advance_particles(p, fext, dt_fluid, 1, 1e-3 * r)
t += dt
i += 1
Plot¶
python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field density --element-type triangle_p1