Download this testcase.

Darcy Flow Through a Square Arrangement of Grains

This example simulates Darcy flow through a square lattice arrangement of circular grains in two dimensions. A fluid velocity is imposed across the periodic 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

np.random.seed(42)

outputdir = "output"
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Physical parameters

g = np.array([0, 0])  # gravity
rho = 1000  # fluid density
rhop = 1500  # grains density
mu = 1  # dynamic viscosity
compacity = 0.5  # Compactness
porosity = 1 - compacity  # porosity
V = 0.00001 / porosity  # fluid velocity

Geometrical parameters

r = 40e-3  # grains radius
L = 1  # domain width
h = r / 2  # mesh size
# h = 1/20                                    # mesh size
origin = np.array([-L / 2, -L / 2])  # mesh origin

gmsh.initialize()
gmsh.model.add("mesh")
gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, 2 * L, 2 * L)
gmsh.model.occ.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] + 2 * L, origin[1]]), name="Bottom"
)
gmsh.model.add_physical_group(
    1,
    get_line([origin[0], origin[1] + 2 * L], [origin[0] + 2 * L, origin[1] + 2 * L]),
    name="Top",
)
gmsh.model.add_physical_group(
    1, get_line([origin[0], origin[1]], [origin[0], origin[1] + 2 * L]), name="Left"
)
gmsh.model.add_physical_group(
    1,
    get_line([origin[0] + 2 * L, origin[1]], [origin[0] + 2 * L, origin[1] + 2 * L]),
    name="Right",
)
gmsh.model.add_physical_group(2, [1], name="domain")
gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: h)
gmsh.model.mesh.generate(2)

Fluid Problem

f = fluid.FluidProblem2(g, mu, rho)
gmsh_io.load_msh(f, None)
f.set_inflow("Left", velocity=[V, 0])
f.set_outflow("Right", 0.0, viscous_traction=True)
f.set_wall("Bottom", velocity=[0, 0])
f.set_wall("Top", velocity=[0, 0])
f.interpolate(velocity_x=V)

Particle Problem

p = scontact.ParticleProblem(2)

x_g = np.arange(0 + 2 * r, 2 * L - 2 * r, 3 * r) + origin[0]
y_g = np.arange(0 + 2 * r, 2 * L - 2 * r, 3 * r) + origin[1]
x_g, y_g = np.meshgrid(x_g, y_g)

x_g = x_g.flat
y_g = y_g.flat
polygons = np.array([x_g, y_g]).T.reshape(-1, 1, 2)
grains_size = np.random.uniform(0.6 * r, 1 * r, polygons.shape[0])
polygons = np.repeat(polygons, 4, axis=1)
polygons[:, 0, :] += np.array([-grains_size, -grains_size]).T
polygons[:, 1, :] += np.array([+grains_size, -grains_size]).T
polygons[:, 2, :] += np.array([+grains_size, +grains_size]).T
polygons[:, 3, :] += np.array([-grains_size, +grains_size]).T

volume = (2 * grains_size) ** 2
density = np.full_like(volume, rho)

for poly, ri in zip(polygons, grains_size):
    x_center = np.mean(poly, axis=0)
    p.add_particle(x_center, ri, ri**2 * 4 * rhop)

Numerical parameters

t = 0
i = 0
dt = 1e-3
tEnd = 1.0
outf = 10

Time Integration

def get_fields(fluid: fluid.FluidProblem2):
    """Return derived output fields for visualization."""
    p1_element = fluid.get_p1_element()
    grad_v = fluid.fields_gradient()[:, :2, :]
    vorticity = (grad_v[:, 1, 0] - grad_v[:, 0, 1]).reshape(-1, 1) * fluid.porosity().get()
    return {
        "pressure": (fluid.pressure(), p1_element),
        "velocity": (fluid.velocity(), p1_element),
        "porosity": (fluid.porosity().get(), p1_element),
        "vorticity": (vorticity, p1_element),
    }


gamma = 1e2 * np.ones(polygons.shape[0])
f.set_polygon_bodies(polygons, density, volume, gamma)


while t < tEnd:
    if i % outf == 0:
        print(f"t/tEnd={t:1.1f}/{tEnd:1.1f}, i={i:1d}")
        f.write_mig(outputdir, t, get_fields(f))
        p.write_mig(outputdir, t)
    force = f.get_forces_on_bodies()
    f.set_polygon_bodies(polygons, density, volume, gamma, contact_forces=-force)
    f.implicit_euler(dt)
    t += dt
    i += 1

Plot

python3 -m migflow.plot.migplot output --actors fluid particles --fluid-field porosity

Artifacts

output/animation.mp4