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!
No comments:
Post a Comment