In this post, details will be shared about a simplified version of the Immersed 🤿 Boundary Method for fluid dynamics 💨 simulations which yours truly has developed in abundant spare time ⏲, referred to as S-IBM 🤓. The purpose of this code is for teaching 🎓 purposes only. The mesh is Cartesian and the method used has no measures in place to resolve the stairstep 𓊍 mesh. The syntax share in this post is for Python 🐍, but this code can be translated to other languages such as C or MATLAB, easily.
The process starts with creating a polygon, a circular cylinder ⭕ for this post. A polygon is created using parametric equations for the circle. NOTE: Before writing these lines, several parameters such as domain size, mesh size, number of points in each coordinate direction etc. have been already defined and calculated. For more information, read here.
theta = np.linspace(0, 2 * np.pi, 100) # create equally spaced angles from 0 to 2 * pi
radius = 0.5; # radius of circle
circle_center_x = 3.3 # circle center x
circle_center_y = D / 2 # circle center x
x1 = circle_center_x + (radius * np.sin(theta)) # x-coordinates
y1 = circle_center_y + (radius * np.cos(theta)) # y-coordinates
Path from matplotlib is used to create a polygon approximation of a circle or an airfoil by using the x and y co-ordinates of the circle (x1, y1) or any shape. Next, coordinate pairs of the grid points are created using vstack and ravel commands. A check is performed to verify if a grid point falls inside the polygon. custom_curve is a boolean array. The values inside this array are True only for the points inside the polygon. Transpose is used to keep the curve orientation consistent with the Navier-Stokes equations.
circle_path = Path(np.column_stack((x1, y1)))
points = np.vstack((X.ravel(), Y.ravel())).T
custom_curve = circle_path.contains_points(points).reshape(X.shape)
curve = custom_curve.T
Another boolean mask is created for boundary points using the boundary_mask statement. Points that are inside the polygon (circle) are marked using using the interior statement.
boundary_mask = np.zeros_like(curve, dtype=bool)
interior = curve[1:-1, 1:-1]
Next a check is performed to mark each interior point's immediate neighbor. right means point to the right of each interior point, bottom means point to the bottom of each interior point etc. boundary statement is used to identify the boundary points (fluid-solid interface) according to the following logic. The point itself is inside the circle (interior = True) and at least one neighbor is outside. boundary_mask[1:-1, 1:-1] is used to fill the boolean with boundary. Consider a 5×5 a grid (zoomed near the polygon boundary). X marks solid points (inside circle), O marks fluid points, as shown in Fig. 1. For example, if interior = True, right = True, left = True, top = True and bottom = True → boundary = True & ~(True) = False i.e. an interior point. And for example, interior = True, right = False, left = True, top = True, bottom = True → boundary = True & ~(False) = True, i.e. a point on the right boundary. For points outside the polygon, interior = False → boundary = False (regardless of neighbors).
 |
| Fig. 1, A grid example |
right = curve[2:, 1:-1]
left = curve[:-2, 1:-1]
top = curve[1:-1, 2:]
bottom = curve[1:-1, :-2]
boundary = interior & ~(right & left & top & bottom)
boundary_mask[1:-1, 1:-1] = boundary
custom_curve_boundary = boundary_mask
In the next step, boundary_indices obtains indices of all boundary points. right_neighbors, left_neighbors etc. calculates indices of neighboring points in all directions. valid_right, valid_left etc. checks if the neighbor is within domain bounds (not at edge) and is outside the circle (fluid point). This allows proper application of no-slip boundary conditions. If for example, valid_right is True, pressure is taken from the right neighbor to be applied to the corresponding boundary index. A visualization of the points is shown in Fig. 2.
boundary_indices = np.where(custom_curve_boundary)
right_neighbors = (boundary_indices[0] + 1, boundary_indices[1])
left_neighbors = (boundary_indices[0] - 1, boundary_indices[1])
top_neighbors = (boundary_indices[0], boundary_indices[1] + 1)
bottom_neighbors = (boundary_indices[0], boundary_indices[1] - 1)
valid_right = (boundary_indices[0] + 1 < X.shape[0]) & (~curve[right_neighbors])
valid_left = (boundary_indices[0] - 1 >= 0) & (~curve[left_neighbors])
valid_top = (boundary_indices[1] + 1 < X.shape[1]) & (~curve[top_neighbors])
valid_bottom = (boundary_indices[1] - 1 >= 0) & (~curve[bottom_neighbors])
 |
| Fig. 2, Depiction of grid points in various regimes |
While plotting Fig. 2, several points are skipped. This is done to make plot clearer. This code can be easily implemented in other languages such as MATLAB (inpolygon), C (ray-casting) etc. At Re 200, this code predicts a Drag coefficient of 1.395 and a Lift coefficient of 0.053. These values are well within rage of the values reported by [1]. In this code, there is no interpolation of scheme or smoothing employed. The lift and drag force coefficients for one shedding cycle are shown in Fig. 3. Within Fig. 4, u and v velocities (top left and right), pressure and streamlines (bottom left and right) are shown.
 |
| Fig. 3, red = Cd and black = Cl |
 |
| Fig. 4, Post processing from the simulations. |
