Download this testcase.
2D droplet simulation with PFEM
This example simulates a drop of fluid falling into a bulk of the same fluid. We use the setup from ‘On the effect of standard PFEM remeshing on volume conservation in free-surface fluid flow problems’ by Franci and Cremonesi, 2016.
Keywords
PFEM, droplet, 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
H = 0.07
B = 0.3
R = 0.025
drop_center = np.array([B / 2, 2 * H + R, 0])
Mesh parameters
geo_mesh_size = B / 10
mesh_size = R / 5
smin = 0.4 * 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, B, H)
disk = gmsh.model.occ.addDisk(drop_center[0], drop_center[1], drop_center[2], R, R)
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, B, 5 * H)
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", B)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMin", 0.0)
gmsh.model.mesh.field.setNumber(sizeFieldConstant, "YMax", 5 * H)
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,
)
# Put the nodes on the bubble boundary back onto the exact circle (this is necessary due to the refinement)
_, elnbubble = gmsh.model.mesh.getElementsByType(1, alphaBoundaryTag)
nt_bubble = np.unique(elnbubble)
for n in nt_bubble:
c, _, _, _ = gmsh.model.mesh.getNode(n)
if c[1] > 2 * H:
rad = np.linalg.norm(c - drop_center)
c[0] = drop_center[0] + R * (c[0] - drop_center[0]) / rad
c[1] = drop_center[1] + R * (c[1] - drop_center[1]) / rad
gmsh.model.mesh.setNode(n, c, [])
Physical Parameters
g = np.array([0.0, -9.81])
rho = 1000
rho_bubble = 1
mu = 1e-1
sigma = 0.0
Time parameters
cfl = 0.3
U = 0.1
U_init = U
dt = mesh_size / U * cfl
t = 0
tEnd = 1.05
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)
max_U = np.max([U_init, np.max(np.linalg.norm(f.velocity(), axis=1))])
dt = smin / max_U * cfl
t += dt
i += 1
Plot
python3 -m migflow.plot.migplot output --actors fluid --fluid-field velocity --fluid-vmin 0.0 --fluid-vmax 0.5 --show-edges 1