Backwards facing step is the simplest case in external aerodynamics. In this post, the code is configured to simulate the case of flow around a 3D backwards facing step. The code 🖳 shared in this post is fully vectorized 😲. Steps 13 to 20 are available here. 13, 14, 15, 16, 17, 18, 19, 20. A thorough validation of the 2D code is available here and here. 🤓
NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... 🙉
If for some reason (why? 😁), you plan to use these codes in your scholarly work, do cite this blog as:
Fahad Butt (2026). S-IBM (https://fluiddynamicscomputer.blogspot.com/2026/08/21st-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year
Literally, everyone else uses this case to validate the code they write. The case is solved at Re 200 without any turbulence models or wall functions. This code is a simple 3D extension of this code. The code is almost double in length as the 3D version requires under-relaxation and has one more momentum equation. In less than 150 lines, flow around the backwards facing step can be simulated with excellent accuracy. The code is available here.
#%% import libraries
import cupy as cp
import matplotlib.pyplot as plt
#%% define parameters
l_cr = 1 # characteristic length
h = l_cr / 100 # grid spacing
dt = 0.001 # time step
L = 9 # domain length
D = 1 # domain width
W = 2 # domain depth
Nx = round(L / h) # grid points in x-axis
Ny = round(D / h) # grid points in y-axis
Nz = round(W / h) # grid points in z-axis
nu = 1 / 200 # kinematic viscosity
Uinf = 1 # free stream velocity / inlet velocity / lid velocity
cfl = dt * Uinf / h # cfl number
travel = 1 # times the disturbance travels entire length of computational domain
TT = travel * L / Uinf # total time
ns = int(TT / dt) # number of time steps
Re = round(l_cr * Uinf / nu) # Reynolds number
#%% intialization
u = cp.zeros((Nx, Ny, Nz)) # x-velocity
v = cp.zeros((Nx, Ny, Nz)) # y-velocity
w = cp.zeros((Nx, Ny, Nz)) # z-velocity
p = cp.zeros((Nx, Ny, Nz)) # pressure
X, Y, Z = cp.meshgrid(cp.linspace(0, L, Nx), cp.linspace(0, D, Ny), cp.linspace(0, W, Nz), indexing = 'ij') # spatial grid
#%% pre calculate for speed
P1 = 1 / (2 * h * dt)
P2 = 1 / (4 * h * h)
P3 = h**2
P4 = 1 / 6
P5 = (2 / Re) * dt / h**2
P6 = dt / h
P7 = 1 - (6 * P5)
P8 = 0.75 # under relaxation
P9 = 1 - P8
#%% create cube
cube_length = 1
cube_width = D
cube_depth = 1
center_x = cube_length / 2
center_y = cube_width / 2
center_z = cube_depth / 2
left_face = round((center_x - (cube_length / 2)) * Nx / L)
right_face = round((center_x + (cube_length / 2)) * Nx / L)
back_face = round((center_y - (cube_width / 2)) * Ny / D)
front_face = round((center_y + (cube_width / 2)) * Ny / D) - 1
bottom_face = round((center_z - (cube_depth / 2)) * Nz / W)
top_face = round((center_z + (cube_depth / 2)) * Nz / W)
#%% solve 3D Navier-Stokes equations
for nt in range(ns):
dudx = u[2:, 1:-1, 1:-1] - u[:-2, 1:-1, 1:-1]
dvdy = v[1:-1, 2:, 1:-1] - v[1:-1, :-2, 1:-1]
dwdz = w[1:-1, 1:-1, 2:] - w[1:-1, 1:-1, :-2]
pn = p.copy()
b = P1 * (dudx + dvdy + dwdz) - P2 * (dudx**2 + dvdy**2 + dwdz**2 + 2 * ((u[1:-1, 2:, 1:-1] - u[1:-1, :-2, 1:-1]) * (v[2:, 1:-1, 1:-1] - v[:-2, 1:-1, 1:-1]) + (u[1:-1, 1:-1, 2:] - u[1:-1, 1:-1, :-2]) * (w[2:, 1:-1, 1:-1] - w[:-2, 1:-1, 1:-1]) + (v[1:-1, 1:-1, 2:] - v[1:-1, 1:-1, :-2]) * (w[1:-1, 2:, 1:-1] - w[1:-1, :-2, 1:-1])))
p[1:-1, 1:-1, 1:-1] = P4 * (pn[2:, 1:-1, 1:-1] + pn[:-2, 1:-1, 1:-1] + pn[1:-1, 2:, 1:-1] + pn[1:-1, :-2, 1:-1] + pn[1:-1, 1:-1, 2:] + pn[1:-1, 1:-1, :-2] - P3 * b) # mass
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[:, :, 0] = p[:, :, 1] # dp/dz = 0 at z = 0
p[:, :, -1] = p[:, :, -2] # dp/dz = 0 at z = H
p[left_face:right_face, back_face:front_face, bottom_face:top_face] = 0
p[left_face, back_face:front_face, bottom_face:top_face] = p[left_face, back_face:front_face, bottom_face:top_face]
p[right_face, back_face:front_face, bottom_face:top_face] = p[right_face + 1, back_face:front_face, bottom_face:top_face]
p[left_face:right_face, front_face, bottom_face:top_face] = p[left_face:right_face, front_face, bottom_face:top_face]
p[left_face:right_face, back_face, bottom_face:top_face] = p[left_face:right_face, back_face, bottom_face:top_face]
p[left_face:right_face, back_face:front_face, bottom_face] = p[left_face:right_face, back_face:front_face, bottom_face]
p[left_face:right_face, back_face:front_face, top_face] = p[left_face:right_face, back_face:front_face, top_face + 1]
p = P8 * p + P9 * pn
un = u.copy()
vn = v.copy()
wn = w.copy()
u[1:-1, 1:-1, 1:-1] = un[1:-1, 1:-1, 1:-1] * P7 - P6 * (un[1:-1, 1:-1, 1:-1] * (un[2:, 1:-1, 1:-1] - un[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (un[1:-1, 2:, 1:-1] - un[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (un[1:-1, 1:-1, 2:] - un[1:-1, 1:-1, :-2]) + p[2:, 1:-1, 1:-1] - p[:-2, 1:-1, 1:-1]) + P5 * (un[2:, 1:-1, 1:-1] + un[:-2, 1:-1, 1:-1] + un[1:-1, 2:, 1:-1] + un[1:-1, :-2, 1:-1] + un[1:-1, 1:-1, 2:] + un[1:-1, 1:-1, :-2]) # x momentum
u[0, :, :] = Uinf # u = 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[:, :, 0] = 0 # u = 0 at z = 0
u[:, :, -1] = 0 # u = 0 at z = W
u[left_face:right_face, back_face:front_face, bottom_face:top_face] = 0
u[left_face, back_face:front_face, bottom_face:top_face] = 0
u[right_face, back_face:front_face, bottom_face:top_face] = 0
u[left_face:right_face, front_face, bottom_face:top_face] = 0
u[left_face:right_face, back_face, bottom_face:top_face] = 0
u[left_face:right_face, back_face:front_face, bottom_face] = 0
u[left_face:right_face, back_face:front_face, top_face] = 0
u = P8 * u + P9 * un
v[1:-1, 1:-1, 1:-1] = vn[1:-1, 1:-1, 1:-1] * P7 - P6 * (un[1:-1, 1:-1, 1:-1] * (vn[2:, 1:-1, 1:-1] - vn[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (vn[1:-1, 2:, 1:-1] - vn[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (vn[1:-1, 1:-1, 2:] - vn[1:-1, 1:-1, :-2]) + p[1:-1, 2:, 1:-1] - p[1:-1, :-2, 1:-1]) + P5 * (vn[2:, 1:-1, 1:-1] + vn[:-2, 1:-1, 1:-1] + vn[1:-1, 2:, 1:-1] + vn[1:-1, :-2, 1:-1] + vn[1:-1, 1:-1, 2:] + vn[1:-1, 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[:, :, 0] = 0 # v = 0 at z = 0
v[:, :, -1] = 0 # v = 0 at z = W
v[left_face:right_face, back_face:front_face, bottom_face:top_face] = 0
v[left_face, back_face:front_face, bottom_face:top_face] = 0
v[right_face, back_face:front_face, bottom_face:top_face] = 0
v[left_face:right_face, front_face, bottom_face:top_face] = 0
v[left_face:right_face, back_face, bottom_face:top_face] = 0
v[left_face:right_face, back_face:front_face, bottom_face] = 0
v[left_face:right_face, back_face:front_face, top_face] = 0
v = P8 * v + P9 * vn
w[1:-1, 1:-1, 1:-1] = wn[1:-1, 1:-1, 1:-1] * P7 - P6 * (un[1:-1, 1:-1, 1:-1] * (wn[2:, 1:-1, 1:-1] - wn[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (wn[1:-1, 2:, 1:-1] - wn[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (wn[1:-1, 1:-1, 2:] - wn[1:-1, 1:-1, :-2]) + p[1:-1, 1:-1, 2:] - p[1:-1, 1:-1, :-2]) + P5 * (wn[2:, 1:-1, 1:-1] + wn[:-2, 1:-1, 1:-1] + wn[1:-1, 2:, 1:-1] + wn[1:-1, :-2, 1:-1] + wn[1:-1, 1:-1, 2:] + wn[1:-1, 1:-1, :-2]) # z momentum
w[0, :, :] = 0 # w = 0 at x = 0
w[-1, :, :] = w[-2, :, :] # dw/dx = 0 at # x = L
w[:, 0, :] = 0 # w = 0 at y = 0
w[:, -1, :] = 0 # w = 0 at y = D
w[:, :, 0] = 0 # w = 0 at z = 0
w[:, :, -1] = 0 # w = 0 at z = W
w[left_face:right_face, back_face:front_face, bottom_face:top_face] = 0
w[left_face, back_face:front_face, bottom_face:top_face] = 0
w[right_face, back_face:front_face, bottom_face:top_face] = 0
w[left_face:right_face, front_face, bottom_face:top_face] = 0
w[left_face:right_face, back_face, bottom_face:top_face] = 0
w[left_face:right_face, back_face:front_face, bottom_face] = 0
w[left_face:right_face, back_face:front_face, top_face] = 0
w = P8 * w + P9 * wn
if nt % 1000 == 0:
print("% Complete", 100 * nt / ns)
#%% post process
u1 = u.copy() # u-velocity for plotting with box
v1 = v.copy() # v-velocity for plotting with box
w1 = w.copy() # v-velocity for plotting with box
p1 = p.copy() # pressure for plotting with box
u1[left_face:right_face, back_face:front_face, bottom_face:top_face] = cp.nan
v1[left_face:right_face, back_face:front_face, bottom_face:top_face] = cp.nan
w1[left_face:right_face, back_face:front_face, bottom_face:top_face] = cp.nan
p1[left_face:right_face, back_face:front_face, bottom_face:top_face] = cp.nan
fig = plt.figure(dpi = 500)
ax = fig.add_subplot(111, projection = '3d')
ax.contourf(X[:, Ny // 2, :].get(), u1[:, Ny // 2, :].get(), Z[:, Ny // 2, :].get(), zdir = 'y', offset = D / 2, levels = 128, cmap='jet', alpha = 0.5) # vertical plane
ax.set_xlim(0, L)
ax.set_ylim(0, D)
ax.set_zlim(0, W)
ax.set_xticks([0, L])
ax.set_yticks([0, D])
ax.set_zticks([0, W])
ax.tick_params(axis='x', labelsize = 4, pad = -2)
ax.tick_params(axis='y', labelsize = 4, pad = -2)
ax.tick_params(axis='z', labelsize = 4, pad = -2)
ax.set_xlabel('x [m]', fontsize = 4, labelpad = -15)
ax.set_ylabel('y [m]', fontsize = 4, labelpad = -15)
ax.set_zlabel('z [m]', fontsize = 4, labelpad = -15)
plt.gca().set_aspect('equal')
ax.view_init(elev = 22.5, azim = -45)
plt.show()
The resultant velocities and pressure field is shown in Fig. 1. The velocity iso-surfaces and streamlines along the plane of flow are shown with Fig. 2.
 |
| Fig. 1, Post-processing |
 |
| Fig. 2, More post-processing |
Thank you very much for reading! If you want to hire me as your next shining post-doc or collaborate in research, please reach out!
