Saithe Fish Simulation: ANSYS Fluent Dynamic Mesh Setup

     One of the most famous post on the blog can be read here. Worryingly😁, many fellow researchers and readers are interested in the aerodynamics of flexible robots 🤓. In this post, the dynamic mesh 🕸 settings used are shared 🥰. These settings are used to reproduce 🖨️ the results from [1], all those years ago. All in a hope that this post helps the readers in their scholarly work! 🎩

     Once the UDF 💻 has been acquired, the next step is to apply the UDF to the airfoil 🐠 geometry correctly ✔️. The airfoil geometry at the first time-step 🕰 i.e. at t = 0 for UDF 02 obtained from [1] is made available here. Once on the dynamic 🏎️ mesh page, select the options shown in Fig. 1. The options selected in Fig. 1 show the default parameters. Within Fig. 1, "wing" refers to the named selection that includes the only the airfoil geometry. Named selections can be created during the meshing process. The "wing" named selection is shown in Fig. 3.

Fig. 1, The dynamic mesh settings

     Before following the settings in Fig. 1, do remember to compile the UDF. To compile the UDF, please use the settings shown in Fig. 2. After selecting the UDF, select the options as shown in the Fig. 2 and then select Build and Load.

Fig. 2, Compile UDF

Fig. 3, Named selection for the dynamic mesh

     The maximum Lift ⬆️ force coefficient from the simulations performed using the method explained here is at 1.77 as compared to 1.68 [1]. The average Drag ⬅️ coefficient is at 0.097 as compared to 0.103 [1]. The obtained flow-field 🏟️ is shown in Fig. 4. Within Fig. 4, top row has v and u components of velocity while the bottom row shows pressure field❗

Fig. 4, The flow-field

     If you are still having trouble, switch to immersed boundary method. The immersed boundary method code your truly wrote, is available here. The validation of this code is available here, here, here and more generally here.

     If you want to hire me as your next shining post-doc or collaborate in research, please reach out! Thank you for reading!

References

[1] Shi, Fulong, Xin, Jianjian and Ou, Chuanzhong, Li, Zhiwei, Chang, Xing, Wan, Ling, "Effects of the Reynolds number and attack angle on wake dynamics of fish swimming in oblique flows", Physics of Fluids, 37(2), 025205, 2025 https://doi.org/10.1063/5.0252506 

A GPU Accelerated Simplified Immersed Boundary Method using Ray Casting

     Yours truly is an avid gamer, @fadoobaba (YouTube). The ray casting algorithm🛸; a fundamental 🧱 technique used in video game 🎮 development and computer graphics, has been implemented within the finite difference method 🏁 code yours truly has been developing. In abundant spare time, of course😼. In this post, this method is explained. For details about ray casting using matplotlib.path, refer to here. For validation of the code, refer to here, backwards-facing step, curved boundaries, here (moving cylinder) and here.

NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... 🙉

     If you plan to use these codes in your scholarly work, do cite this blog as:

     Fahad Butt (2025). S-IBM (https://fluiddynamicscomputer.blogspot.com/2025/07/a-simplified-immersed-boundary-method.html), Blogger. Retrieved Month Date, Year

     The first step is to setup the polygon 🐠. The polygon = cp.column_stack((x1, y1)) statement is used to combine x1 and y1 into a single array representing the polygon vertices. The following statements are used to store the x and y coordinates of the polygon vertices and nvert is the total number of vertices.

px = polygon[:, 0]

py = polygon[:, 1]

nvert = len(polygon)

     The boolean masks 👺 are then initialized. The following arrays track whether a grid point is inside the polygon based on horizontal and vertical ray intersections. False = outside ❌, True = inside ✅. Then a loop is implemented to for each edge 📏 of the polygon, defined by vertices i (current) and j (previous), with % nvert ensuring the polygon is closed (last vertex connects back to the first). (xi, yi) and (xj, yj) define the current edge and previous edge.

horizontal_inside = cp.zeros_like(test_x, dtype=bool) 

vertical_inside = cp.zeros_like(test_x, dtype=bool) 

for i in range(nvert):

    j = (i - 1) % nvert

xi, yi = px[i], py[i]

xj, yj = px[j], py[j]

     A check ⁉️ is performed to ensure a ray crosses the polygon edge once. cond1 statement checks if the test point lies between the y-values of the edge's endpoints i.e., a ray could cross it. intersect_x finds where the edge crosses a horizontal line at test_y. cond2 checks if the intersection lies to the right of the test point. ^= is XOR toggles the "inside" state each time the ray crosses an edge. The addition of small term ~1e-16 prevents division by zero for horizontal edges. Similar method is applied to verify the points using vertical ray. cond3 checks if the edge crosses a vertical ray (top to bottom) from (test_x, test_y). slope1 and intersect_y computes where the vertical ray at x=test_x intersects the edge. A point is inside only if both horizontal and vertical rays classify it as inside. A point is considered inside the foil only if both ray checks are true. curve is the 2D boolean mask of grid points inside the foil body.

cond1 = ((yi > test_y) != (yj > test_y))

slope = (xj - xi) / (yj - yi + 1e-16)

intersect_x = slope * (test_y - yi) + xi

cond2 = test_x < intersect_x 

horizontal_inside ^= cond1 & cond2

cond3 = ((xi > test_x) != (xj > test_x))

slope1 = (yj - yi) / (xj - xi + 1e-16)

intersect_y = slope1 * (test_x - xi) + yi

cond4 = test_y < intersect_y

vertical_inside ^= cond3 & cond4

inside = horizontal_inside & vertical_inside

curve = inside.reshape(X.T.shape)

     Using both horizontal and vertical rays reduces false positives (e.g., near sharp corners). The & ensures only points unambiguously inside are marked. Within Fig. 1, the points on the shape boundary, points inside 🪰 and outside 🐌 the shape boundary are shown. A point is a boundary if it is inside the body, but any of its neighbors is outside. Following statements are used to mark the boundary of the object. Boundary detects the "skin" of the body for applying no-slip, force, or stress conditions.

interior = curve[1:-1, 1:-1]

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)

     The statements boundary_indices = cp.where(boundary_mask) ... valid_right = (boundary_indices[0] + 1 < X.shape[0]) & (~curve[right_neighbors]) ... etc. extract boundary indices and their valid neighbors i.e. these lines extract the (i, j) grid indices of the surface and locate which neighbor cells are valid fluid neighbors (outside the body and within domain bounds). These are needed to compute normals or apply boundary conditions via interpolation or extrapolation from the fluid.

Fig. 1, Mesh cells

     For validation of the results from present simulations, the case of flow around a circular cylinder is selected. Fig. 1 shows the results from the code at Re 200. The drag coefficient obtained from this code is 1.396 while from the literature, the value is at 1.4 [1]. The lift coefficient from the code is at 0.000134. Within Fig. 2, top row shows u and v velocity components and bottom row shows pressure and vorticity.

Fig. 2, The flow-field

     For the second an more complex validation case, a swimming fish is simulated. The lift and drag coefficients obtained from the simulation are compared with experimental results. The drag coefficient from the current code is at 0.342 while from the published literature, the value is at 0.348. Maximum lift coefficient is at 7.243 and 8 from the present code as compared to the published literature [2]. Within Fig. 3 the u and v velocities, pressure and velocity streamlines are shown for St = 0.8 and Reynolds number of 500. The computational mesh near the fish is shown in Fig. 4. Within Fig.4, a zoomed in view towards the right shows the mesh at the trailing edge of the fish.

Fig.3, Post-processing of results

Fig. 4, The mesh

     This method allows handling arbitrary deforming / non-deforming shapes on a fixed Cartesian grid. In summary, the method has the following steps.

1. Generate polygon shape (e.g., airfoil, cylinder)

2. Flatten mesh for vectorized testing

3. Use ray casting to check if points are inside shape

4. Build a boolean mask of body region

5. Identify surface (boundary) points

6. Extract boundary indices for physics coupling

References

     [1] Braza M, Chassaing P, Minh HH. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics. 1986;165:79-130. doi:10.1017/S0022112086003014

     [2] Fulong Shi, Jianjian Xin, Chuanzhong Ou, Zhiwei Li, Xing Chang, Ling Wan; Effects of the Reynolds number and attack angle on wake dynamics of fish swimming in oblique flows. Physics of Fluids 1 February 2025; 37 (2): 025205 doi.org/10.1063/5.0252506

     Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!

A Simplified Immersed Boundary Method (S-IBM)

     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 = C_d and black = C_l

Fig. 4, Post processing from the simulations.

References

[1] Braza M, Chassaing P, Minh HH. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics. 1986;165:79-130. doi:10.1017/S0022112086003014

     Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!

18th Step of the 12 Steps to Navier-Stokes 😑 (Morphing Airfoil / Swimming Fish)

     This is probably the final step in the series. Unless, one fine morning yours truly randomly decides to write about moving and morphing airfoils 😑 i.e. a true virtual wind tunnel 🤓.

     This whole series is a continuation of a series created by Dr. Barbra 👩‍🏫 to teach about coding 🖳 of the Navier-Stokes 🌬 equations. The code made available this blog is an independent continuation of her open-source work. Yours truly has coded the scaled version of the equations with a modified pressure Poisson equation. The original series ends with the example of the channel flow 🔚. So far, yours truly has added flow around obstacles both curved and non-curved; flow around moving obstacles and species transport to the series. The Fluid-Structure Interaction (FSI) version of the code 💻 is available here.

     This post is about flow around shape changing 🪱 obstacles i.e. morphing wings. The development of these codes shared here, including the one mentioned in this post, started in the year 2020 and continued till the fall of 2023. All this effort was to create a digital CV so that yours truly can get himself accepted in to a funded PhD program. That happened in fall 2023 😌.

     The frequency of oscillations, the equations for amplitude and lateral oscillations 〰️ are all taken from the published literature about swimming fish 🐟. As a starting point, the NACA 0012 airfoil is used. A polygon is then created. Each grid index is then tested for presence inside or outside the polygon. An index is marked as the boundary of the polygon if one index has anyone neighboring index outside the polygon. The neighboring boundary indices are used to apply Neumann boundary conditions for pressure. The boundary indices are used to apply no-slip Dirichlet boundary conditions on the polygon boundary. More details can be read here.

     Because the camber of airfoil changes at each time-step, approximately 400 points are required on the airfoil boundary to properly capture the time-varying camber of the airfoil i.e. the fish 🐋. The output from the code is shown in Fig. 1. Within Fig. 1, top row has u and v components of velocity. Meanwhile, the bottom row shows the pressure and streamlines. The animation is shown in Fig. 2.

     If you plan to use these codes in your scholarly work, do cite this blog as:

     Fahad Butt (2025). 18th Step of the 12 Steps to Navier-Stokes (https://fluiddynamicscomputer.blogspot.com/2025/05/18th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year

Fig. 1, The result from code

 

Fig. 2, Animation

Code

#Copyright <2025> <FAHAD BUTT>

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#%% import libraries

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.path import Path

#%% define parameters

l_cr = 1 # characteristic length

h = l_cr / 400 # grid spacing

dt = 0.00001 # time step

L = 8 # domain length

D = 4 # domain depth

Nx = round(L / h) + 1 # grid points in x-axis

Ny = round(D / h) + 1 # grid points in y-axis

nu = 1 / 500 # kinematic viscosity

Uinf = 1 # free stream velocity / inlet velocity / lid velocity

cfl = dt * Uinf / h # cfl number

Re = round(l_cr * Uinf / nu) # Reynolds number

#%% intialization

u = Uinf * np.ones((Nx, Ny)) # x-velocity

v = np.zeros((Nx, Ny)) # y-velocity

p = np.zeros((Nx, Ny)) # pressure

X, Y = np.meshgrid(np.linspace(0, L, Nx), np.linspace(0, D, Ny)) # spatial grid

#%% pre calculate for speed

P1 = h / (16 * dt)

P2 = (2 / Re) * dt / h**2

P3 = dt / h

P4 = 1 - (4 * P2)

#%% motion parameters

freq = 0.8 # frequency of oscillation (from litrature)

factor = 2 * np.pi * freq # 2 * pi * f

TP = 1 / freq # time period

ns = int(6 * TP / dt) # number of time steps

t = 0.12 # airfoil thickness

c = 1 # airfoil chord length

num_points = 150 # points along chord

alpha = np.linspace(0, np.pi, num_points) # angle values

x_airfoil = 0.5 * (1 - np.cos(alpha)) * c # cosine spacing on x coordinates

x_frac = x_airfoil / c # normalized airfoil x coordinates

a0 = 0.02 # amplitude co efficient a0 (from litrature)

a1 = 0.01 # amplitude co efficient a1

a2 = 0.1 # amplitude co efficient a2

a_x = a0 + a1 * x_frac + a2 * x_frac**2 # amplitude of oscillation (from litrature)

yt = 5 * t * (0.2969 * np.sqrt(x_airfoil) - 0.1260 * x_airfoil - 0.3516 * x_airfoil**2 + 0.2843 * x_airfoil**3 - 0.1036 * x_airfoil**4) # airfoil co ordinates

foil_center_x = L / 3 # foil center x

foil_center_y = D / 2 # foil center y

x1 = foil_center_x + np.concatenate([x_airfoil, x_airfoil[::-1]]) # foil x

#%% solve 2D Navier-Stokes equations

for nt in range(ns):

    y_oscillation = a_x * np.cos(2 * np.pi * x_frac - factor * nt * dt) # lateral oscillation (from litrature)

    y_top = y_oscillation + yt # airfoil top

    y_bottom = y_oscillation - yt # airfoil bottom

    y1 = foil_center_y + np.concatenate([y_top, y_bottom[::-1]]) # foil y

    foil_path = Path(np.column_stack((x1, y1))) # foil region (polygon)

    points = np.vstack((X.ravel(), Y.ravel())).T # check if inside region

    custom_curve = foil_path.contains_points(points).reshape(X.shape) # test all grid points for inclusion in airfoil, contains_points returns True for each point inside the airfoil (polygon)

    curve = custom_curve.T # transpose to match grid orientation

    boundary_mask = np.zeros_like(curve, dtype = bool) # create boolean mask

    interior = curve[1:-1, 1:-1] # interior of foil

    right = curve[2:, 1:-1] # boolean arrays indicating if neighbors of each interior point are inside airfoil (polygon)

    left = curve[:-2, 1:-1]

    top = curve[1:-1, 2:]

    bottom = curve[1:-1, :-2]

    boundary = interior & ~(right & left & top & bottom) # mark if point is inside foil and at least one neighbor is outside foil (mark interface between solid and fluid)

    boundary_mask[1:-1, 1:-1] = boundary # fill the empty boundary

    custom_curve_boundary = boundary_mask.T # transpose back

    boundary_indices = np.where(custom_curve_boundary.T) # mark boundary indices

    right_neighbors = (boundary_indices[0] + 1, boundary_indices[1]) # right index (neighboring grid points in each direction)

    left_neighbors = (boundary_indices[0] - 1, boundary_indices[1]) # left index

    top_neighbors = (boundary_indices[0], boundary_indices[1] + 1) # top index

    bottom_neighbors = (boundary_indices[0], boundary_indices[1] - 1) # bottom index

    # identify valid neighbor indices, ensure the neighbor is not inside airfoil (polygon) i.e. is inside fluid)

    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])

    pn = p.copy()

    p[1:-1, 1:-1] = 0.25 * (pn[2:, 1:-1] + pn[:-2, 1:-1] + pn[1:-1, 2:] + pn[1:-1, :-2]) - P1 * (u[2:, 1:-1] - u[:-2, 1:-1] + v[1:-1, 2:] - v[1:-1, :-2]) # pressure

    p[0, :] = p[1, :] # dp/dx = 0 at x = 0

    p[-1, :] = p[-2, :] # dp/dx = 0 at x = L

    p[:, 0] = p[:, 1] # dp/dy = 0 at y = 0

    p[:, -1] = p[:, -2] # dp/dy = 0 at y = D

    p[curve] = 0 # p = 0 inside shape

    p[boundary_indices] = np.where(valid_right, p[right_neighbors], p[boundary_indices]) # right neighbor is fluid

    p[boundary_indices] = np.where(valid_left, p[left_neighbors], p[boundary_indices]) # left neighbor is fluid

    p[boundary_indices] = np.where(valid_top, p[top_neighbors], p[boundary_indices]) # top neighbor is fluid

    p[boundary_indices] = np.where(valid_bottom, p[bottom_neighbors], p[boundary_indices]) # bottom neighbor is fluid

    un = u.copy()

    vn = v.copy()

    u[1:-1, 1:-1] = un[1:-1, 1:-1] * P4 - P3 * (un[1:-1, 1:-1] * (un[2:, 1:-1] - un[:-2, 1:-1]) + vn[1:-1, 1:-1] * (un[1:-1, 2:] - un[1:-1, :-2]) + p[2:, 1:-1] - p[:-2, 1:-1]) + P2 * (un[2:, 1:-1] + un[:-2, 1:-1] + un[1:-1, 2:] + un[1:-1, :-2]) # x momentum

    u[0, :] = Uinf # u = Uinf at x = 0

    u[-1, :] = u[-2, :] # du/dx = 0 at x = L

    u[:, 0] = Uinf # u = Uinf at y = 0

    u[:, -1] = Uinf # u = Uinf at y = D

    u[curve] = 0 # u = 0 inside shape

    u[custom_curve_boundary.T] = 0 # no slip

    v[1:-1, 1:-1] = vn[1:-1, 1:-1] * P4 - P3 * (un[1:-1, 1:-1] * (vn[2:, 1:-1] - vn[:-2, 1:-1]) + vn[1:-1, 1:-1] * (vn[1:-1, 2:] - vn[1:-1, :-2]) + p[1:-1, 2:] - p[1:-1, :-2]) + P2 * (vn[2:, 1:-1] + vn[:-2, 1:-1] + vn[1:-1, 2:] + vn[1:-1, :-2]) # y momentum

    v[0, :] = 0 # v = 0 at x = 0

    v[-1, :] = v[-2, :] # dv/dx = 0 at x = L

    v[:, 0] = 0 # v = 0 at y = 0

    v[:, -1] = 0 # v = 0 at y = D

    v[curve] = 0 # u = 0 inside shape

    v[custom_curve_boundary.T] = 0 # no slip

    boundary_x = X[boundary_indices[1], boundary_indices[0]] # get actual x-coordinate values of all boundary points

    x_frac_boundary = (boundary_x - foil_center_x) / c # normalized streamwise position along the airfoil for each boundary point

    vinf_foil_boundary = (a0 + a1 * x_frac_boundary + a2 * x_frac_boundary**2) * factor * np.sin(2 * np.pi * x_frac_boundary - factor * nt * dt) # vertical velocity at boundaries

    v[top_neighbors] = vinf_foil_boundary # above the boundary

    v[bottom_neighbors] = vinf_foil_boundary # below the boundary

    v[right_neighbors] = vinf_foil_boundary # right of boundary

    v[left_neighbors] = vinf_foil_boundary # left of boundary

#%% post process

u1 = u.copy() # u-velocity for plotting with foil

v1 = v.copy() # v-velocity for plotting with foil

p1 = p.copy() # pressure for plotting with foil

u1[custom_curve.T] = np.nan # shape geometry for plotting

v1[custom_curve.T] = np.nan

p1[custom_curve.T] = np.nan

velocity_magnitude1 = np.sqrt(u1**2 + v1**2) # velocity magnitude with shape

# visualize velocity vectors and pressure contours

plt.figure(dpi = 500)

plt.contourf(X, Y, u1.T, 128, cmap = 'jet')

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 1)

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

plt.contourf(X, Y, v1.T, 128, cmap = 'jet')

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 1)

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

plt.contourf(X, Y, p1.T, 128, cmap = 'jet')

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 1)

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

plt.streamplot(X, Y, u1.T, v1.T, color = 'black', cmap = 'jet', density = 2, linewidth = 0.5, arrowstyle='->', arrowsize = 1)  # plot streamlines

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 1)

plt.gca().set_aspect('equal', adjustable = 'box')

plt.xlim([0, L])

plt.ylim([0, D])

plt.axis('off')

plt.show()

     Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!