Download this testcase.

2D lid-driven cavity flow with PFEM

This example simulates a 2D lid-driven cavity flow using the Particle Finite Element Method (PFEM).

Keywords

PFEM, lid driven cavity, 2D, fluid

Description

The test demonstrates: - How to perform a PFEM simulation in a fixed domain (without free surfaces)

import os, sys, shutil
import numpy as np

Gmsh pfem branch

This test requires the alphashapes branch of GMSH to use the Particle Finite Element (PFEM) module. PFEM is a lagrangian method where the mesh is regenerated at each time step, making it suitable for problems with large deformations and free surfaces.

sys.path.insert(0, os.environ["GMSH_PFEM_DIR"])

import gmsh
from migflow import fluid, pfem

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

l_box = 1.0
h_box = 1.0

Mesh parameters

geo_mesh_size = l_box / 10
mesh_size = l_box / 40
alpha = 10.0

Initial fluid mesh

gmsh.model.add("ModelInit")
gmsh.model.occ.addRectangle(0, 0, 0, l_box, h_box)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
gmsh.model.mesh.generate(2)
nodeTags, coords, _ = gmsh.model.mesh.getNodes()

Solid domain

gmsh.model.add("ModelGeo")
gmsh.model.occ.addRectangle(0, 0, 0, l_box, h_box)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSizeCallback(lambda *args: mesh_size)
gmsh.model.mesh.generate(2)
nodeTags, coords, _ = gmsh.model.mesh.getNodes()
gmsh.model.addPhysicalGroup(1, [1], -1, "bottom")
gmsh.model.addPhysicalGroup(1, [2], -1, "right")
gmsh.model.addPhysicalGroup(1, [3], -1, "top")
gmsh.model.addPhysicalGroup(1, [4], -1, "left")
geoEntities = gmsh.model.getEntities(1)

PFEM mesh and size fields

gmsh.model.add("ModelFluid")
alphaDomainTag = gmsh.model.addDiscreteEntity(2, -1, [])
for dim, tag in geoEntities:
    gmsh.model.addDiscreteEntity(dim, tag, [])
alphaBoundaryTag = gmsh.model.addDiscreteEntity(1, -1, [])
gmsh.model.mesh.addNodes(2, alphaDomainTag, nodeTags, coords)
gmsh.model.addPhysicalGroup(1, [1], -1, "bottom")
gmsh.model.addPhysicalGroup(1, [2], -1, "right")
gmsh.model.addPhysicalGroup(1, [3], -1, "top")
gmsh.model.addPhysicalGroup(1, [4], -1, "left")
gmsh.model.addPhysicalGroup(1, [alphaBoundaryTag], -1, "freeSurface")
gmsh.model.addPhysicalGroup(2, [alphaDomainTag], -1, "domain")

Size fields

sizeFieldConstant = gmsh.model.mesh.field.add("Box")
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VIn", mesh_size)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "VOut", mesh_size)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMin", 0.0)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "XMax", l_box)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMin", 0.0)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMax", h_box)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "Thickness", 0.001)

Physical Parameters

g = np.array([0.0, 0.0])
rho = 1000
mu = 1e-3
Re = 1000
v_top = Re * mu / (rho * l_box)

Time parameters

cfl = 0.5
U = v_top
U_init = v_top
dt = mesh_size / U * cfl
t = 0
tEnd = 25000.0

Fluid problem

f = fluid.FluidProblem(2, g, mu, rho, advection=False)
f.set_wall_boundary("bottom")
f.set_wall_boundary("right")
f.set_wall_boundary("top", velocity=[v_top, 0])
f.set_wall_boundary("left")
f.set_strong_boundary("bottom", velocity=[0, 0])
f.set_strong_boundary("right", velocity=[0, 0])
f.set_strong_boundary("top", velocity=[v_top, 0])
f.set_strong_boundary("left", velocity=[0, 0])

Simulation Loop

Time integration of coupled fluid–particle motion.

i = 0
outf = 25
gmsh.option.setNumber("General.Verbosity", 0)

while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}, {dt:.6g}")
    # Update PFEM mesh
    nodetag, oelemtag, oparamcoord, _ = gmsh.model.mesh.computeAlphaShape(
        2,
        alphaDomainTag,
        alphaBoundaryTag,
        "ModelGeo",
        alpha,
        sizeFieldConstant,
        sizeFieldConstant,
        boundaryTolerance=0.01 * mesh_size,
        usePreviousMesh=False,
    )
    oparamcoord = oparamcoord.reshape((-1, 3))
    gmsh.write(outputdir + "/lastMesh.msh")
    ordered_node_tags = pfem.prepareMeshForMigflow(
        i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
    )

    if i % outf == 0:  # or t > 0.3:
        f.write_mig(outputdir, t)

    # nodes velocity, be aware that its dimension is (n_nodes, 3) not (n_nodes, 2)
    u_old = np.zeros_like(f.coordinates())
    u_old[:, :2] = f.velocity()

    f.implicit_euler(dt)

    dx = np.zeros((f.coordinates().shape[0], 3))
    u = np.zeros_like(dx)
    u[:, :2] = f.velocity()
    if i == 0:
        u_old = u
    dx = u * dt + 0.5 * (u - u_old) / dt * dt**2

    # Do not move the top boundary nodes
    top_nodes = np.unique(f.mesh_boundaries()[b"top"].flatten())
    dx[top_nodes, :] = 0.0

    # advect nodes and project if needed
    gmsh.model.mesh.advectMeshNodes(
        2,
        alphaDomainTag,
        alphaBoundaryTag,
        "ModelGeo",
        ordered_node_tags,
        dx.flatten(),
        0.01 * mesh_size,
    )

    # set new coordinates
    _, newCoords, _ = gmsh.model.mesh.getNodes(2, alphaDomainTag)
    f.set_coordinates(newCoords)

    t += dt
    i += 1

Plot

python3 -m migflow.plot.migplot output --actors fluid --fluid-field velocity --fluid-vmin 0 --fluid-vmax 0.001 --show-edges 1