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!

Datacenter Visualization (Verified and Validated)

     This simulation is done to create an aero-thermal digital twin of a datacenter using CFD. The details of datacenter are taken from [1]. The datacenter CAD model is shown in Fig. 1.

     The simulation employs κ − ε turbulence model with damping functions, SIMPLE-R (modified), as the numerical algorithm and second-order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The time derivatives are approximated with an implicit first-order Euler scheme. Flow simulation solves the Navier–Stokes equations, which are formulations of mass, momentum, and energy conservation laws for fluid flows. To predict turbulent flows, the Favre-averaged Navier–Stokes equations are used.

Fig. 1, Datacenter CAD

     A Cartesian mesh with octree refinement, cut-cell method and immersed boundary methods is used. Special mesh refinements are deployed in the areas of interest i.e. inlets and outlets and sharp edges of server racks and CRAH units to accurately capture aero-thermal gradients and vortices. The resulting computational mesh has 2,698,156 cells. The computational domain and mesh are shown in Fig. 2.

Fig. 2, Computational mesh and domain

     The results from the numerical analysis were compared with [1]. The results are in excellent agreement with previously published data. The animation in Fig. 3 shows thermal distribution inside datacenter using cut-plots. Fluid velocity distribution is also shown. The cut-plots are superimposed with streamlines and velocity vectors. These post processing features help identify hot-spots and recirculation zones. Design improvements can be made to reduce thee unwanted flow features. Within Fig. 3, Flow trajectories colored by air temperature are also shown. These show path the fluid takes between various inlets and outlets in the datacenter such as the CRAH system supply and return zones and inlets and outlets of servers. Fig. 4 shows various post processing tools available for diagnosing various issues from the aero-thermal perspective. These include iso-surfaces, cut-plots, flow trajectories and surface plots etc.

Fig. 3, The post processing animations

Fig. 4, The post processing images

     A comparison of the results from present simulations with previously published literature is shown in Fig. 5. it can be seen that out results are in close agreement with previously published numerical and experimental data [1]! The locations at which the data is extracted is shown in Fig. 6. Within Fig. 5, solid lines indicate present study, dashed lines indicate published numerical results and filled circles represent published experimental results.

Fig. 5, Comparison of results

Fig. 6, Location of data extraction

     If you want to collaborate on the research projects related to turbomachinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

[1] Wibron, Emelie, Anna-Lena Ljung, and T. Staffan Lundström. 2018. "Computational Fluid Dynamics Modeling and Validating Experiments of Airflow in a Data Center" Energies 11, no. 3: 644. https://doi.org/10.3390/en11030644