An immersed boundary method for simulating fluid-solid interaction and particle-laden multiphase flows.
Requirement: g++
with C++11
support and OpenFOAM v6. Other compilers weren't tested but should work.
Step 1, follow the official guide to install OpenFOAM v6. Test afterwards to ensure a working installation (e.g., by running some serial and parallel tutorial cases).
Step 2, whether on your laptop or HPC cluster, simply type three commands.
git clone https://github.com/ChenguangZhang/sdfibm.git
cd sdfibm/src
make
The solver binary is ./src/sdfibm
. It is recommended that you soft-link it to a system-wide path (e.g., sudo ln -s ~/sdfibm/src/sdfibm /usr/local/bin/sdfibm
), then you can simply type sdfibm
to run a simulation without specifying the whole binary path.
-
The convention of coordinate system is that we face the
$x$ -$y$ plane, with$z$ -axis pointing towards us. -
A 2D simulation uses a one-cell-thick 3D mesh,
sdfibm
mandates that the "thickness direction" is in$z$ and the$z$ -span of the mesh is$[-0.5,0.5]$ . In other words, the 2D dynamics occurs exactly on the$z=0$ , aka the$x$ -$y$ plane. In addition, due to the different ways of handling the geometry, theon_twod
option insolidDict
must be set to 1.
Cases under the example folder are some good starting points of using the code. That include
- 2D "flow past cylinder" at
$Re=200$ (animation below)
-
the Taylor-Couette flow with a rotational core. A MATLAB finite-difference script is also there, which solves a reduced-order PDE that validates sdfibm
-
an ellipse that wobbles and settles in a rectangular container
-
./examples\sedimentation
is a case where 100 circular particles settle in a square container (a video version you can pause is here https://www.youtube.com/watch?v=C6U9X9zatz8).The particle trajectories below clearly show the circulation caused by the counter flow.
The plane definition (e.g., in the sedimentation case) can be a bit confusing, especially its orientation. We know a plane is specified by an on-plane point pos
below) and a direction
If we want to make a left wall---meaning the fluid is to its right---then pos
variable to pos
is that it is on the plane, so
Not directly related with immersed boundary method, a tool is created to smoothly initialize the phase field of VOF simulations. The word "smooth" here is essential when simulating low Capillary number flows. The simplest method (as used by setFields
) that simply sets the phase fraction at mesh cells to 0/1 creates a staircase phase interface, which generates capillary waves that pollute the simulation easily. My tool initializes the VOF phase field accurately and smoothly. The image below is the initialized field of the (2D) ./tool_vof/example
case, which shows
- the various shapes can be used
- the shapes can overlap each other (see lower right)
The case is ready for interFoam
and can be run right after initialization. The result is below. As expected, the droplet oscillations are smooth without anything spurious.
The animation below compares smooth (left) vs. staircase (right) initializations, using the oscillation of inviscid 2D elliptical droplets (setFields
of OpenFOAM---there the spurious capillary waves caused by the staircase boundary are clearly visible. Even those waves appear diminished later on, they are hard to control or separate from the physical results. Thus, the left simulation is strongly preferred.
Lastly, the tool works with arbitrary polyhedron meshes, see figure below (left: current too, right: setFields
, contoured value: 0.5).
Initially I thought about writing locations from sdfibm, but actually ParaView contour can do the same thing. It is also more flexible, because you can specify the contour value to slightly offset from the real solid surface.
Here are the steps for 2D cases.
- First, read only the front surface. In case you don't have a separate front patch, you can create a slice of the 3D volume.
- Second, apply the "Contour" filter. Certainly, it should be contoured by As and a reasonable (single) contour value is 0.5.
- Third, apply the "ExtractSurface" filter.
- Fourth, from the menu, File -> Save Data ..., choose csv format. ParaView will ask you to choose the data to save. The final output will contain the xy location of the contour and the data (such as U, p, etc.) at those contour locations.
Cheers.
-
libshape
contains shape definitions using SDF-
ishape.h
the virtual base class all shapes inherit from. It also contains many SDF transformation and Boolean utilities -
shapefactory.h
implement the factory design pattern, it holds a library of shapes that are used in theSolidCloud
class -
template.h
a heavily annotated template for creating new shapes, place to change are marked byCHANGE
-
circle.h
a circle of radius$r$ -
ellipse.h
an ellipse of radius$ra$ and$rb$ -
kernel.h
legacy code for volume fraction calculation in square/cubic Cartesian mesh
-
-
libmotion
contains motion constraints. A rigid body has six degrees of freedom (linear velocity in$x$ ,$y$ ,$z$ , and angular velocity in$x$ ,$y$ ,$z$ ), correspondingly many of the file names contains six digits of[0,1,2]
. The convention is: 0 if the DOF is killed, 1 if it is free (i.e., unconstrained), and 2 if the DOF is prescribed.-
imotion.h
the virtual base class all motions inherit from. -
motion01mask.h
generic 0-1 mask of the degrees of freedom. It largely supersedes other files whose name contains 01 strings.
-
-
libcollision
for collision detection, currently only sphere-sphere and sphere-plane -
solidcloud.h (.cpp)
the main class. It holds record of all solids and is responsible for reading the input, creating the shape/motion/material library, interacting with the fluid, and storing data into thecloud.out
output file -
solid.h
corresponds to the individual solid. It keeps and updates the kinematic and dynamic information of a solid. It also holds pointers to shape/motion/material instances obtained from the corresponding libraries in theSolidCloud
class. -
main.cpp
implements the projection method and the direct forcing IBM
Information of the solid bodies is written to the cloud.out
file, which has the simple plain text format: if there are sdfibm
writes
names = ["t", "x", "y", "z", "vx", "vy", "vz", "fx", "fy", "fz", "EulerAx", "EulerAy", "EulerAz", "wx", "wy", "wz", "Tx", "Ty", "Tz"]
In words, they are:
Time, position-x, position-x, position-z, velocity-x, velocity-y, velocity-z, force-x, force-y, force-z, Euler_angle-x, Euler_angle-y, Euler_angle-z, angular_velocity-x, angular_velocity-y, angular_velocity-z, torque-x, torque-y, torque-z
Cite the following papers if you use this code.
-
Chenguang Zhang, Chunliang Wu, and Krishnaswamy Nandakumar. Effective geometric algorithms for immersed boundary method using signed distance field. Journal of Fluids Engineering, 141(6):061401, 2019.
-
Chenguang Zhang. sdfibm: a Signed Distance Field Based Discrete Forcing Immersed Boundary Method in OpenFOAM. Computer Physics Communications (accepted 2020).
sdfibm
is licensed under GNU GENERAL PUBLIC LICENSE (GPLv3). You can use, copy and modify sdfibm
freely. If you distribute your software based on (any part of) sdfibm
, you are obliged to distribute the source code of your software under GPL, in addition to any modifications made to sdfibm
.
In case the constraints prevent you from using sdfibm
, you can obtain a commercial license by contacting me.
I prefer 🍔 to coffee if you want to brighten my day :)
To check the version of your sdfibm
repository, use git log -1
. The minus one means the latest commit.