FOANSS, i.e. Fadoobaba's Open Advanced Navier-Stokes 🌬 Solver is now capable to simulate flow with moving objects 🤓. In this example, the objects are square and circular cylinders. These blog posts could have been a master's thesis or a peer-reviewed publication but yours truly decided to are write here for educational purposes! The
The simulated case is called the case of flow around an obstacle. For validation, read. This code and blog post is not endorsed or approved by Dr. Barba, I just continue the open-source work of her. The 13th, 14th, 15th and 16th Steps are available for reading along with code. The resulting motion of the square cylinder is shown in Fig. 1. The Fig. 1 is colored using the horizontal velocity. The Fig. 2 shows the example of moving circle. Within Fig. 2, the results of u and v velocities are shown on top while the pressure and streamlines are shown at the bottom. The code for stationery curved / arbitrary boundaries is made available here. A customary remainder that you are reading a blog, not a Q1 journal 😃.
If you plan to use these codes in your scholarly work, do cite this blog as:
Fahad Butt (2025). 17th Step of the 12 Steps to Navier-Stokes (https://fluiddynamicscomputer.blogspot.com/2025/04/17th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year
Fig. 1, The 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
#%% define parameters
l_cr = 1 # characteristic length
h = l_cr / 50 # grid spacing
dt = 0.0001 # time step
L = 10 # domain length
D = 5 # domain depth
Nx = round(L / h) + 1 # grid points in x-axis
Ny = round(D / h) + 1 # grid points in y-axis
nu = 1 / 100 # 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 = np.zeros((Nx, Ny)) # x-velocity
v = np.zeros((Nx, Ny)) # y-velocity
p = np.zeros((Nx, Ny)) # pressure
#%% pre calculate for speed
P1 = h / (16 * dt)
P2 = (2 / Re) * dt / h**2
P3 = dt / h
P4 = 1 - (4 * P2)
#%% add square
square_center_x = 0.75 # box starting position x-axis
square_center_y = D / 2 # box starting position y-axis
freq = 0.32 # f = 2 * St * Uinf / l_cr
factor = 2 * np.pi * freq # 2 * pi * f
a = 1 # amplitude of vibration
TP = 1 / freq # time period
ns = int(2 * TP / dt) # number of time steps
uinf_square = Uinf # horizontal velocity of box
#%% solve 2D Navier-Stokes equations
for nt in range(ns):
square_center_y_variable = square_center_y + a * np.sin(factor * dt * nt) # vertical position of box
square_center_x_variable = square_center_x + uinf_square * dt * nt # horizontal position of box
vinf_square = a * factor * np.cos(factor * dt * nt) # vertical velocity of box
square_left = round((square_center_x_variable - l_cr / 2) * Nx / L) # box left
square_right = round((square_center_x_variable + l_cr / 2) * Nx / L) # box right
square_bottom = round((square_center_y_variable - l_cr / 2) * Ny / D) # box bottom
square_top = round((square_center_y_variable + l_cr / 2)* Ny / D) # box top
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[square_left:square_right, square_bottom:square_top] = 0 # box geometry
p[square_left, square_bottom:square_top] = p[square_left - 1, square_bottom:square_top] # box left
p[square_right, square_bottom:square_top] = p[square_right + 1, square_bottom:square_top] # box right
p[square_left:square_right, square_bottom] = p[square_left:square_right, square_bottom - 1] # box bottom
p[square_left:square_right, square_top] = p[square_left:square_right, square_top + 1] # box top
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, :] = u[1, :] # du/dx = 0 at x = 0
u[-1, :] = u[-2, :] # du/dx = 0 at x = L
u[:, 0] = 0 # u = 0 at y = 0
u[:, -1] = 0 # u = 0 at y = D
u[square_left:square_right, square_bottom:square_top] = 0 # box geometry
u[square_left, square_bottom:square_top] = 0 # box left
u[square_right, square_bottom:square_top] = 0 # box right
u[square_left:square_right, square_bottom] = 0 # box bottom
u[square_left:square_right, square_top] = 0 # box top
u[square_left - 1, square_bottom:square_top] = uinf_square # box left
u[square_right + 1, square_bottom:square_top] = uinf_square # box right
u[square_left:square_right, square_bottom - 1] = uinf_square # box bottom
u[square_left:square_right, square_top + 1] = uinf_square # box top
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, :] = v[1, :] # dv/dx = 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[square_left:square_right, square_bottom:square_top] = 0 # box geometry
v[square_left, square_bottom:square_top] = 0 # box left
v[square_right, square_bottom:square_top] = 0 # box right
v[square_left:square_right, square_bottom] = 0 # box bottom
v[square_left:square_right, square_top] = 0 # box top
v[square_left - 1, square_bottom:square_top] = vinf_square # box left
v[square_right + 1, square_bottom:square_top] = vinf_square # box right
v[square_left:square_right, square_bottom - 1] = vinf_square # box bottom
v[square_left:square_right, square_top + 1] = vinf_square # box top
#%% post processing
u1 = u # copy for plotting
v1 = v
p1 = p
u1[square_left:square_right, square_bottom:square_top] = np.nan # create cavity for plotting
v1[square_left:square_right, square_bottom:square_top] = np.nan
p1[square_left:square_right, square_bottom:square_top] = np.nan
X, Y = np.meshgrid(np.linspace(0, L, Nx), np.linspace(0, D, Ny)) # spatial grid
plt.figure(dpi = 200)
plt.contourf(X, Y, u1.T, 128, cmap = 'jet') # plot contours
# 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.show()
plt.figure(dpi = 200)
plt.contourf(X, Y, v1.T, 128, cmap = 'jet') # plot contours
# 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.show()
plt.figure(dpi = 200)
plt.contourf(X, Y, p1.T, 128, cmap = 'jet') # plot contours
# 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.show()
plt.figure(dpi = 200)
plt.streamplot(X, Y, u1.T, v1.T, color = 'black', cmap = 'jet', density = 2, linewidth = 0.5, arrowstyle='->', arrowsize = 1) # plot streamlines
plt.gca().set_aspect('equal', adjustable='box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()
Code (Curved / Slanted Boundaries)
#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 / 20 # grid spacing
dt = 0.001 # time step
L = 24 # domain length
D = 5 # domain depth
Nx = round(L / h) + 1 # grid points in x-axis
Ny = round(D / h) + 1 # grid points in y-axis
nu = 1 / 100 # 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 = np.zeros((Nx, Ny)) # x-velocity
v = np.zeros((Nx, Ny)) # y-velocity
p = np.zeros((Nx, Ny)) # pressure
#%% pre calculate for speed
P1 = h / (16 * dt)
P2 = (2 / Re) * dt / h**2
P3 = dt / h
P4 = 1 - (4 * P2)
#%% create a shape
X, Y = np.meshgrid(np.linspace(0, L, Nx), np.linspace(0, D, Ny)) # spatial grid
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 = 0.75 # circle center x (starting point)
circle_center_y = D / 2 # circle center y
#%% motion parameters
freq = 0.32 # f = 2 * St * Uinf / l_cr
factor = 2 * np.pi * freq # 2 * pi * f
a = 1 # amplitude of vibration
TP = 1 / freq # time period
ns = int(6 * TP / dt) # number of time steps
uinf_circle = Uinf # horizontal velocity of box
#%% solve 2D Navier-Stokes equations
for nt in range(ns):
circle_center_x_variable = circle_center_x + uinf_circle * dt * nt # horizontal position of box
circle_center_y_variable = circle_center_y + a * np.sin(factor * dt * nt) # vertical position of circle
vinf_circle = a * factor * np.cos(factor * dt * nt) # vertical velocity of circle
x1 = circle_center_x_variable + (radius * np.sin(theta)) # x-coordinates
y1 = circle_center_y_variable + (radius * np.cos(theta)) # y-coordinates
cylinder_path = Path(np.column_stack((x1, y1))) # cylinder region
points = np.vstack((X.ravel(), Y.ravel())).T # check if inside region
custom_curve = cylinder_path.contains_points(points).reshape(X.shape) # create boolean mask
curve = custom_curve.T # transpose
boundary_mask = np.zeros_like(curve, dtype = bool) # find boundary of cylinder
interior = curve[1:-1, 1:-1] # interior of circle
right = curve[2:, 1:-1] # mark neighbours
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 cylinder and at least one neighbor is outside cylinder
boundary_mask[1:-1, 1:-1] = boundary # set 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
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
valid_right = ~curve[right_neighbors] # check if points are on shape boundary
valid_left = ~curve[left_neighbors]
valid_top = ~curve[top_neighbors]
valid_bottom = ~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, :] = 0 # p = 0 at x = 0
p[-1, :] = 0 # p = 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, :] = 0 # u = 0 at x = 0
u[-1, :] = 0 # u = 0 at x = L
u[:, 0] = 0 # u = 0 at y = 0
u[:, -1] = 0 # u = 0 at y = D
u[curve] = 0 # u = 0 inside shape
u[custom_curve_boundary.T] = 0 # no slip
u[top_neighbors] = uinf_circle # above the boundary
u[bottom_neighbors] = uinf_circle # below the boundary
u[right_neighbors] = uinf_circle # right of boundary
u[left_neighbors] = uinf_circle # left of boundary
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, :] = 0 # v = 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
v[top_neighbors] = vinf_circle # above the boundary
v[bottom_neighbors] = vinf_circle # below the boundary
v[right_neighbors] = vinf_circle # right of boundary
v[left_neighbors] = vinf_circle # left of boundary
#%% post process
u1 = u.copy() # u-velocity for plotting with circle
v1 = v.copy() # v-velocity for plotting with circle
p1 = p.copy() # pressure for plotting with circle
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', vmin = -2, vmax = 3)
plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)
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', vmin = -2, vmax = 2)
plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)
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', vmin = -10, vmax = 6)
plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)
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()
 |
| Fig. 2, Results from curved boundary simulation |
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!
