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.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.FluidProblem(2, 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