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

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)
p.load_msh_boundaries(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.FluidProblem(2, g, mu, rho)
f.load_msh(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.pressure_index()][:, 0].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[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()

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

dp = f.fields_gradient()[:, -1]
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