Download this testcase.
Darcy Flow Through a Bed of Circular Discs
This example simulates Darcy flow through a packed bed of circular disc-shaped particles in two dimensions. A pressure gradient is applied across the domain and the resulting flow field through the porous medium is computed.
Keywords
FEM, DEM, Darcy, Porous Media
import os
import gmsh
import shutil
import numpy as np
from migflow import fluid, scontact, gmsh_io
outputdir = "output"
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
Geometrical parameters
r = 5e-3 # grains radius
L = 25 * 2 * r # domain width
H = 25 * 2 * r # domain height
mesh_size = r / 4 # 0.1 # mesh size
# mesh_size = 0.1 # mesh size
# mesh_size = L/4 # mesh size
origin = np.array([-L / 2, -L / 2]) # mesh origin
Physical parameters
rho = 1000 # fluid density
rhop = 1500 # grains density
mu = 1 # dynamic viscosity
compacity = 0.5 # Compactness
porosity = 1 - compacity # porosity
p1 = 0
p0 = 10.0
g = np.array([-(p1 - p0) / (rho * L), 0]) # gravity
gmsh.initialize()
gmsh.model.add("mesh")
x0 = origin
x1 = origin + np.array([L, 0])
x2 = origin + np.array([L, H])
x3 = origin + np.array([0, H])
x = np.array([x0, x1, x2, x3])
edges = np.array([[0, 1], [1, 2], [2, 3], [3, 0]])
for xi in x:
gmsh.model.geo.add_point(xi[0], xi[1], 0, mesh_size)
for edge in edges:
gmsh.model.geo.add_line(edge[0] + 1, edge[1] + 1)
gmsh.model.geo.add_curve_loop([1, 2, 3, 4], 1)
gmsh.model.geo.add_plane_surface([1], 1)
gmsh.model.geo.synchronize()
def get_line(x0, x1, eps=1e-6):
bbentities = gmsh.model.get_entities_in_bounding_box(
x0[0] - eps, x0[1] - eps, -eps, x1[0] + eps, x1[1] + eps, eps, 1
)
return list(tag for dim, tag in bbentities)
# gmsh.model.add_physical_group(1,get_line([origin[0],origin[1]], [origin[0]+L,origin[1]]), name="Bottom")
# gmsh.model.add_physical_group(1,get_line([origin[0],origin[1]+L], [origin[0]+L,origin[1]+L]), name="Top")
gmsh.model.add_physical_group(
1, get_line([origin[0], origin[1]], [origin[0], origin[1] + L]), name="Left"
)
gmsh.model.add_physical_group(
1,
get_line([origin[0] + L, origin[1]], [origin[0] + L, origin[1] + L]),
name="Right",
)
gmsh.model.add_physical_group(2, [1], name="domain")
transform_y = np.array([[1, 0, 0, 0], [0, 1, 0, -H], [0, 0, 1, 0], [0, 0, 0, 1]])
transform_x = np.array([[1, 0, 0, +L], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
gmsh.model.mesh.set_periodic(1, [1], [3], transform_y.flatten().tolist())
gmsh.model.mesh.set_periodic(1, [2], [4], transform_x.flatten().tolist())
gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: mesh_size)
gmsh.model.mesh.generate(2)
p = scontact.ParticleProblem(2)
gmsh_io.load_msh_boundaries(p, None)
# Definition of the points where the grains are located
# The particles are first placed on a regular grid
# Then they are centered in the domain
a_g = np.pi * r**2
a_b = L**2
N = max(compacity * a_b / a_g, 1)
e = L / N**0.5
x = np.linspace(-L / 2, L / 2, int(N**0.5))
y = np.linspace(-H / 2, H / 2, int(N**0.5))
# x = np.arange(-L/2, L/2-r, e)
# y = np.arange(-L/2, L/2-r, e)
x, y = np.meshgrid(x, y)
for xi, yi in zip(x.flat, y.flat):
p.add_particle((xi, yi), r, r**2 * np.pi * rhop)
f = fluid.FluidProblem2(g, mu, rho)
gmsh_io.load_msh(f, None)
f.set_mean_pressure(0)
# f.interpolate(velocity_x=U)
Numerical parameters
dt = 1e-4 # time step
tEnd = 100 * dt # final time
outf = 10 # number of iterations between output files
t = 0.0
i = 0
Time Integration
def get_fields(fluid):
"""Return derived output fields for visualization."""
x = fluid.coordinates_fields()[fluid.field_indices(fluid.dimension())][:, 0]
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[0] * x, p1_element),
}
while t < tEnd:
print(f"t/tEnd={t:1.4f}/{tEnd:1.4f}, i={i:1d}")
if i % outf == 0:
p.write_mig(outputdir, t)
f.write_mig(outputdir, t, get_fields(f))
contacts = f.compute_node_force() if i > 0 else np.zeros_like(p.velocity())
f.set_particles(
p.delassus(), p.volume(), p.position(), p.velocity(), p.omega(), -contacts
)
f.implicit_euler(dt)
i += 1
t += dt
node_volume = f.node_volume().get()
U = np.sum(f.velocity()[:, 0] * node_volume) / np.sum(node_volume)
V = U / porosity
Re = V * rho * 2 * r / mu
print("Re : %2.g ---- t : %2.g ---- tEnd : %2.g" % (Re, dt, tEnd))
voidage = porosity ** (-1.8)
cross_section = 2 * r
Re = U * rho * cross_section / mu
Cd = voidage * (0.63 + 4.8 * (porosity * Re) ** (-0.5)) ** 2
gamma = Cd * 0.5 * rho * cross_section * U
vol = np.pi * r**2
N = compacity / vol
Dp_drag = -gamma * N * V
# The pressure gradient along the flow: fields_gradient is (node, field, dim).
dp = f.fields_gradient()[:, -1, 0]
dp = np.sum(dp * node_volume) / np.sum(node_volume) - rho * g[0]
error = np.sum(abs(dp - Dp_drag) / Dp_drag)
print(dp, Dp_drag, error)
print("Pressure gradient from simulation : ", dp)
print("Pressure gradient from drag model : ", Dp_drag)
print("Relative error : ", error)
Plot
python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field velocity