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

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.FluidProblem(2, g, mu, rho)
f.load_msh(None)
f.set_open_boundary("Left", velocity=[V, 0])
f.set_open_boundary("Right", pressure=0.0)
f.set_wall_boundary("Bottom", velocity=[0, 0])
f.set_wall_boundary("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)

csr_ptr, csr_element, csr_xi, csr_surface = f.get_mesh_polygons_overlap(polygons)

if False:
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from matplotlib.patches import Circle

    fig, ax = plt.subplots()
    ax.set_aspect("equal")
    triangulation = mpl.tri.Triangulation(
        f.coordinates()[:, 0], f.coordinates()[:, 1], f.elements()
    )
    ax.triplot(triangulation, linewidth=0.2, color="k")
    phi0 = 1 - csr_xi[..., 0] - csr_xi[..., 1]
    phi1 = csr_xi[..., 0]
    phi2 = csr_xi[..., 1]
    phi = np.array([phi0, phi1, phi2]).T  # qp, node
    xelt = f.coordinates()[f.elements()[csr_element], :2]  # qp, node, dim
    xqp = np.einsum("ijk, ij -> ik", xelt, phi)  # qp, dim
    for i, x_i in enumerate(xqp):
        p_area = csr_surface[i]
        r = np.sqrt(p_area / np.pi)
        c = Circle(x_i[:2], r / 2, facecolor="black")
        ax.add_patch(c)
    ax.set_xlim(origin[0], origin[0] + 2 * L)
    ax.set_ylim(origin[1], origin[1] + 2 * L)
    for poly in polygons:
        polygon = plt.Polygon(poly, fill=None, edgecolor="r")
        ax.add_patch(polygon)
    plt.show()

Numerical parameters

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

Time Integration

def get_fields(fluid: fluid.FluidProblem):
    """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()
    return {
        "pressure": (fluid.pressure(), p1_element),
        "velocity": (fluid.velocity(), p1_element),
        "porosity": (fluid.porosity(), p1_element),
        "vorticity": (vorticity, p1_element),
    }


csr_velocity = np.zeros_like(csr_xi)
csr_gamma = 1e2 * np.ones(csr_xi.shape[0])
csr_contact = np.zeros_like(csr_xi)
f.set_bodies(
    density,
    volume,
    csr_ptr,
    csr_element,
    csr_xi,
    csr_surface,
    csr_velocity,
    csr_gamma,
    csr_contact,
)


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)
    csr_force = f.get_bodies_csr_force().copy()
    f.set_bodies(
        density,
        volume,
        csr_ptr,
        csr_element,
        csr_xi,
        csr_surface,
        csr_velocity,
        csr_gamma,
        -csr_force,
    )
    f.implicit_euler(dt)
    t += dt
    i += 1

Plot

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