# %% # 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 # ---- # .. code-block:: shell # # 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 # %% # Artifacts # --------- # - output/animation.mp4