20th Step of the 12 Steps to Navier-Stokes 😑 (Heated Room)

     In abundant spare time ⏳, yours truly has updated the code for 3D incompressible Navier–Stokes 🍃 equations using the finite-difference method with the ability to simulate species transport, for example heat or temperature. This post is a continuation of this post.

As before, this code 🖳 is fully vectorized 😲 with the only change being an addition of coupled energy equation. Steps 13 to 19 are available here. 13, 14, 15, 16, 17, 18, 19. 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 weird reason, you plan to use this code in your scholarly work, do cite this blog as:

     Fahad Butt (2026). S-IBM (https://fluiddynamicscomputer.blogspot.com/2026/08/20th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year

     Flow inside an empty room 📦 provides a complex 💢 flow 🌬 with very simple boundary conditions. This case has been widely used by researchers to validate the HVAC code. The case is solved at ~Re 4,000 without any turbulence models or wall functions. The room used in the simulation is 1 x 1 x 1 m. The inlet and outlet vents are 0.02 and 0.023 m wide. The floor of the room is at a higher temperature and remaining walls are at a lower temperature. The code is a simple 3D extension of this 2D code. The code to reproduce plots shown within Fig. 1 is available here ❗

#%% import libraries

import cupy as cp

import matplotlib.pyplot as plt

#%% define parameters

l_cr = 1 # characteristic length

h = 0.02 / 2 # grid spacing

dt = 0.001 # time step

L = 1 # domain length

D = 1 # domain depth

W = 1 # domain width

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 / 4000 # kinematic viscosity

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

cfl = dt * Uinf / h # cfl number

travel = 5 # 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

Pr = 0.71465 # Prandtl number

alpha = 0.025969 # thermal conductivity

g = -9.81 # earth gravity

Tinf = 1 # free stream temperature

#%% 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

T = (20/35) * cp.ones((Nx, Ny, Nz)) # temperature

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

P10 = dt / (Re * Pr * h**2)

P11 = 1 - (6 * P10)

P12 = dt / (2 * h)

P13 = dt * alpha * g

P14 = round(0.023 * Ny / D)

P15 = round(0.98 * Ny / D)

#%% 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]))) # divergence

    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 = 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, :, 0:P15] = 0 # u at x = 0

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

    u[-1, :, 0:P14] = u[-2, :, 0:P14] # u = 0 at # x = L

    u[-1, :, P14:] = 0 # u = 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 = Uinf at z = W

    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, :, 0:P14] = v[-2, :, 0:P14] # v = 0 at # x = L

    v[-1, :, P14:] = 0 # v = 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 = 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]) - P13 * T[1:-1, 1:-1, 1:-1] # z momentum

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

    w[-1, :, 0:P14] = w[-2, :, 0:P14] # w = 0 at # x = L

    w[-1, :, P14:] = 0 # w = 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 = P8 * w + P9 * wn

    Tn = T.copy()

    T[1:-1, 1:-1, 1:-1] = Tn[1:-1, 1:-1, 1:-1] * P11 - P12 * (un[1:-1, 1:-1, 1:-1] * (Tn[2:, 1:-1, 1:-1] - Tn[:-2, 1:-1, 1:-1]) + vn[1:-1, 1:-1, 1:-1] * (Tn[1:-1, 2:, 1:-1] - Tn[1:-1, :-2, 1:-1]) + wn[1:-1, 1:-1, 1:-1] * (Tn[1:-1, 1:-1, 2:] - Tn[1:-1, 1:-1, :-2])) + P10 * (Tn[2:, 1:-1, 1:-1] + Tn[:-2, 1:-1, 1:-1] + Tn[1:-1, 2:, 1:-1] + Tn[1:-1, :-2, 1:-1] + Tn[1:-1, 1:-1, 2:] + Tn[1:-1, 1:-1, :-2]) # energy

    T[0, :, :] = 0.43 # x = 0

    T[-1, :, :] = 0.43 # x = L

    T[:, 0, :] = 0.43 # y = 0

    T[:, -1, :] = 0.43 # y = D

    T[:, :, 0] = Tinf # z = 0

    T[:, :, -1] = 0.43 # z = W

    T = P8 * T + P9 * Tn

    if nt % 1000 == 0:

        print("% Complete", 100 * nt / ns)

#%% post process

fig = plt.figure(dpi = 500)

ax = fig.add_subplot(111, projection = '3d')

ax.contourf(X[:, Ny // 2, :].get(), u[:, 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', pad=-2)

ax.tick_params(axis='y', pad=-2)

ax.tick_params(axis='z', pad=-2)

ax.set_xlabel('x [m]', labelpad = -15)

ax.set_ylabel('y [m]', labelpad = -15)

ax.set_zlabel('z [m]', labelpad = -15)

plt.gca().set_aspect('equal')

ax.view_init(elev = 22.5, azim = -45)

plt.show()

     The results from post processing are shown in Fig. 1. Within Fig. 1, the u, v and w components of velocities and room temperature are shown along the plane of flow.

Fig. 1, post processing

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

19th Step of the 12 Steps to Navier-Stokes 😑 (3D-Lid Driven Cavity)

     One fine morning (in abundant spare time ⏳, of course), yours truly decided to code the 3D incompressible Navier–Stokes 🍃 equations using the finite-difference method. This post has the results of this adventure 🏞️ (so-far). As is customary with all my CFD work using commercial and home-made CFD codes, this too is an unofficial continuation of the series by Dr. Lorena Barba.

The code 🖳 shared in this post is fully vectorized 😲. Steps 13 to 18 are available here. 13, 14, 15, 16, 17, 18. 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, 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/19th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year

     Lid-Driven Cavity 🕳 provides a complex 💢 flow 🌬 with very simple boundary conditions. Literally, everyone else uses this case to validate the code they write. The lid-driven cavity case is solved at ~Re 1,000 without any turbulence models or wall functions. The cavity used in the simulation is 1 x 1 x 1 m. This code is a simple 3D extension of this code. Well, the code is double in length as the 3D version requires under-relaxation and has one more momentum equation. Still, it is less than 100 lines ❗ 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 = 1 # domain length

D = 1 # domain depth

W = 1 # domain width

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 / 1000 # kinematic viscosity

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

cfl = dt * Uinf / h # cfl number

travel = 10 # 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

#%% 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]))) # divergence

    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 = 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, :, :] = 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[:, :, 0] = 0 # u = 0 at z = 0

    u[:, :, -1] = Uinf # u = Uinf at z = W

    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, :, :] = 0 # v = 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 = 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, :, :] = 0 # w = 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 = P8 * w + P9 * wn

    if nt % 1000 == 0:

        print("% Complete", 100 * nt / ns)

#%% post process

fig = plt.figure(dpi = 500)

ax = fig.add_subplot(111, projection = '3d')

ax.contourf(X[:, Ny // 2, :].get(), u[:, 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', pad = -2)

ax.tick_params(axis='y', pad = -2)

ax.tick_params(axis='z', pad = -2)

ax.set_xlabel('x [m]', labelpad = -15)

ax.set_ylabel('y [m]', labelpad = -15)

ax.set_zlabel('z [m]', labelpad = -15)

plt.gca().set_aspect('equal')

ax.view_init(elev = 22.5, azim = -45)

plt.show()

     The results from post processing are shown in Fig. 1. Within Fig. 1, the u, v and w components of velocities are shown along the plane of flow. The resulting pressure field and the velocity streamlines are shown within Fig. 2.

Fig. 1, Velocity components

Fig. 2, The pressure field and streamlines

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

A Simple Poisson's Equation for Pressure

     In abundant spare time 🕰️, yours truly has removed ❌ the non-linear 〰️ terms from the Pressure Poisson Equation (PPE) which is derived 🎡 by summing the divergence of momentum equations and then applying conservation of mass ⚖️. The non-linear terms create a problem in convergence. Therefore, inspired by the style of work of the "Skipper"🐧, yours truly removed the problem causing non-linear terms 😀. The derived PPE is mentioned by equation 1. The PPE used in the code yours truly has developed is mention in equation 2. Within equations 1 and 2, the u and v are components of velocity along x and y-axis. The pressure is represented by p and while, t represents the time.

∂²p/∂x² + ∂²p/∂y² = 1/∆t * (∂u/∂x + ∂v/∂y - (∂u/∂x)² - (∂v/∂y)² - 2*∂v/∂y*∂u/∂x) [1]

∂²p/∂x² + ∂²p/∂y² = 1/(2 * ∆t) * (∂u/∂x + ∂v/∂y) [2]

     For the same grid and for solving the same problem, the time-step supported by equation 2 is 40x more as compared to equation 1. Therefore, the resultant compute per watt is 40x less for equation 2 as compared to equation 1 🤯. The code for equation 1 is shown first, followed by the code for equation 2. The results are compared via streamlines and pressure contours, with in Fig. 1. The benchmark case solved is of the lid-driven cavity which offer simple implementation and very complex flow physics. For validation of the code, refer to here, here and here.

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.

     If you are mad 👨‍🔬 enough to use this code in your scholarly 👨‍🏫 work, then please remember to cite: Fahad Butt (2025). S-PPE (https://fluiddynamicscomputer.blogspot.com/2025/11/a-simple-poissons-equation-for-pressure.html), Blogger. Retrieved Month Date, Year

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

Code 01

#%% 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.00005 # time step

L = 1 # domain length

D = 1 # domain depth

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

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

nu = 1 / 400 # kinematic viscosity

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

cfl = dt * Uinf / h # cfl number

travel = 10 # 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) # Osborne Reynolds and his number :)

#%% intialization (t = 0)

u = cp.zeros((Nx, Ny)) # x-velocity

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

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

X, Y = cp.meshgrid(cp.linspace(0, L, Nx), cp.linspace(0, D, Ny), indexing = 'ij') # spatial grid

#%% pre calculate for speed

P1 = h / (8 * dt)

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

P3 = dt / h

P4 = 1 - (4 * P2)

#%% solve 2D incompressible Navier-Stokes equations

for _ in range(ns):

    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]) - (u[2:, 1:-1] - u[:-2, 1:-1])**2 - (v[1:-1, 2:] - v[1:-1, :-2])**2 - (2 * (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

    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 = Uinf at x = 0

    u[-1, :] = 0 # u = 0 at x = L

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

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

    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

Code 2

#%% 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.002 # time step

L = 1 # domain length

D = 1 # domain depth

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

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

nu = 1 / 400 # kinematic viscosity

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

cfl = dt * Uinf / h # cfl number

travel = 10 # 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) # Osborne Reynolds and his number :)

#%% intialization (t = 0)

u = cp.zeros((Nx, Ny)) # x-velocity

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

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

X, Y = cp.meshgrid(cp.linspace(0, L, Nx), cp.linspace(0, D, Ny), indexing = 'ij') # spatial grid

#%% pre calculate for speed

P1 = h / (16 * dt)

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

P3 = dt / h

P4 = 1 - (4 * P2)

#%% solve 2D incompressible Navier-Stokes equations

for _ in range(ns):

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

    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 = Uinf at x = 0

    u[-1, :] = 0 # u = 0 at x = L

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

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

    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

Fig. 1, pressure range for both cases is 0 (blue) till 1 (red).

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

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 yours truly wrote 🤓, is available here. Of course, this was done in abundant spare time 🕰️. 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