Download this testcase.
2D dam break simulation with PFEM
This example simulates a 2D dam break flow using the Particle Finite Element Method (PFEM). We use the setup from Cremonesi et al., “A State of the Art Review of the Particle Finite Element Method (PFEM)”, 2020.
# Keywords
# --------
# PFEM, dam break, 2D, fluid
#
# Description
# -----------
# The test demonstrates:
# - How to perform a PFEM simulation with free surface dynamics.
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 = 0.146
b = 0.175
l_box = 4 * L
h_box = 8 * L
l_fluid = L
h_fluid = 2 * L
Mesh parameters
geo_mesh_size = l_fluid / 10
mesh_size = l_fluid / 10
smin = 0.3 * mesh_size
smax = mesh_size
dmax = 5 * smax
alpha = 1.3
Initial fluid mesh
gmsh.initialize()
gmsh.model.add("ModelInit")
rect = gmsh.model.occ.addRectangle(0, 0, 0, l_fluid, h_fluid)
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: geo_mesh_size)
gmsh.model.mesh.generate(2)
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)
sizeFieldDistFS = gmsh.model.mesh.field.add("AlphaShapeDistance")
gmsh.model.mesh.field.setNumber(sizeFieldDistFS, "Tag", alphaBoundaryTag)
gmsh.model.mesh.field.setNumber(sizeFieldDistFS, "SamplingLength", 0.1 * smin)
sizeFieldRefine = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(sizeFieldRefine, "InField", sizeFieldDistFS)
gmsh.model.mesh.field.setNumber(sizeFieldRefine, "SizeMin", smin)
gmsh.model.mesh.field.setNumber(sizeFieldRefine, "SizeMax", smax)
gmsh.model.mesh.field.setNumber(sizeFieldRefine, "DistMin", 0.0)
gmsh.model.mesh.field.setNumber(sizeFieldRefine, "DistMax", dmax)
gmsh.model.mesh.computeAlphaShape(
2,
alphaDomainTag,
alphaBoundaryTag,
"ModelGeo",
alpha,
sizeFieldConstant,
sizeFieldRefine,
False,
)
Physical Parameters
g = np.array([0.0, -9.81])
rho = 1000
rho_bubble = 1
mu = 1e-3
sigma = 0.0
Time parameters
cfl = 0.3
U = 10.0
U_init = U
dt = mesh_size / U * cfl
t = 0
tEnd = 1
Fluid problem
f = fluid.FluidProblem2(g, mu, rho, advection=False)
f.set_wall_boundary("bottom")
f.set_wall_boundary("right")
f.set_wall_boundary("top")
f.set_wall_boundary("left")
f.set_strong_boundary("bottom", velocity_y=0)
f.set_strong_boundary("right", velocity_x=0)
f.set_strong_boundary("top", velocity_y=0)
f.set_strong_boundary("left", velocity_x=0)
f.set_open_boundary("freeSurface", pressure=0, viscous_flag=False)
Simulation Loop
Time integration of coupled fluid–particle motion.
i = 0
outf = 10
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,
sizeFieldRefine,
sizeFieldRefine,
boundaryTolerance=0.01 * mesh_size,
usePreviousMesh=True,
)
oparamcoord = oparamcoord.reshape((-1, 3))
gmsh.write(outputdir + "/lastMesh.msh")
ordered_node_tags = pfem.prepareMeshForMigflow(
i, alphaDomainTag, f, nodetag, oelemtag, oparamcoord
)
pfem.applySurfaceTension_v0(
f, sigma, bubbleCondition=True, g=g, rho_bubble=rho_bubble
)
if i % outf == 0:
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
# 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.0 --fluid-vmax 3.0 --show-edges 1 --bounds 0 0.584 0 1.168 --width 500 --height 1000