CFD Wizardry: A 50-Line Python Marvel

     In abundant spare time ⏳, yours truly has implemented the non-conservative and non-dimensional form of the discretized Navier-Stokes 🍃 equations. The code 🖳 in it's simplest form is less than 50 lines including importing libraries and plotting! 😲 For validation, refer here. More examples and free code is available here, here and here. Happy codding!

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 numpy as np

import matplotlib.pyplot as plt

l_square = 1 # length of square

h = l_square / 500 # grid spacing

dt = 0.00002 # 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 / 100 # kinematic viscosity

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

cfl = dt * Uinf / h # cfl number

travel = 200 # 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_square * Uinf / nu) # Reynolds number

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

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

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

for nt in range(ns): # solve 2D Navier-Stokes equations

    pn = p.copy()

    p[1:-1, 1:-1] = (pn[2:, 1:-1] + pn[:-2, 1:-1] + pn[1:-1, 2:] + pn[1:-1, :-2]) / 4 - h / (8 * dt) * (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] = 0 # p = 0 at y = D

    un = u.copy()

    vn = v.copy()

    u[1:-1, 1:-1] = (un[1:-1, 1:-1] - dt / (2 * h) * (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])) - dt / (2 * h) * (p[2:, 1:-1] - p[:-2, 1:-1]) + (1 / Re) * dt / h**2 * (un[2:, 1:-1] + un[:-2, 1:-1] + un[1:-1, 2:] + un[1:-1, :-2] - 4 * un[1:-1, 1:-1])) # 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] = Uinf # u = Uinf at y = D

    v[1:-1, 1:-1] = (vn[1:-1, 1:-1] - dt / (2 * h) * (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])) - dt / (2 * h) * (p[1:-1, 2:] - p[1:-1, :-2]) + (1 / Re) * dt / h**2 * (vn[2:, 1:-1] + vn[:-2, 1:-1] + vn[1:-1, 2:] + vn[1:-1, :-2] - 4 * vn[1:-1, 1:-1])) # 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

X, Y = np.meshgrid(np.linspace(0, L, Nx), np.linspace(0, D, Ny)) # spatial grid

plt.figure(dpi = 200)

plt.contourf(X, Y, v.T, 128, cmap = 'jet') # plot contours

plt.colorbar()

plt.streamplot(X, Y, u.T, v.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()

Lid-Driven Cavity

     The case of lid-driven cavity in the turbulent flow regime can now be solved in reasonable amount of time. The results are shown in Fig. 1. I stopped the code while the flow is still developing as you are reading a blog and not a Q1 journal. 😆 Within Fig. 1, streamlines, v and u component of velocity and pressure are shown going from left to right and top to bottom. At the center of Fig. 1, the velocity magnitude is superimposed. As this is DNS, the smallest spatial scale resolved is ~8e-3 m [8 mm]. While, the smallest time-scale ⌛ resolved is ~8e-4 s [0.8 ms].

Fig. 1, The results at Reynolds number of 10,000

Free-Jet

          The case of free jets in the turbulent flow regime can now be solved in reasonable amount of time. The results are shown in Fig. 2. I stopped the code while the flow is still developing. Once again, I remind you that you are reading a blog and not a Q1 journal. 😆 Within Fig. 2, streamlines and species are shown. As this is DNS, the smallest spatial scale resolved is ~0.02 m [2 cm]. While, the smallest time-scale ⌛ resolved is ~4e-4 s [0.4 ms]. The code for implementing species, in this case temperature using the energy equation is available on the previous post.

Fig. 2, Free jet at Reynolds number 10000

Heated Room

     The benchmark case of mixed convection in an open room in the turbulent flow regime can now be solved in reasonable amount of time as well. The results are shown in Fig. 3. I stopped the code while the flow field stopped showing any changes. 😆 As this is DNS, the smallest spatial scale resolved is ~0.0144 m [1.44 cm]. While, the smallest time-scale ⌛ resolved is ~1e-4 s [1 ms]. The code for implementing species, in this case temperature using the energy equation is available on the previous post. In the previous post, the momentum equation has no changes as the gravity vector is at 0 m/s2. For this example, Boussinesq assumption is used.

Fig. 3, Flow inside a heated room at Reynolds number of 5000

Backward - Facing Step (BFS)

     Another benchmark case of flow around a backwards facing step can now be solved in reasonable amount of time as well. The flow is fully turbulent. The results are shown in Fig. 4. I stopped the code while the flow field is still developing. 😆 As this is DNS, the smallest spatial scale resolved is ~0.01 m [1 cm]. While, the smallest time-scale ⌛ resolved is ~1e-4 s [1 ms]. As can be seen from Fig. 4, there are no abnormalities in the flow field.

Fig. 4, Flow around a backwards facing step at Reynolds number of 10000

PS: I fully understand, there is no such thing as 2D turbulence 🍃. Just don't kill the vibe please 💫.

     Thank you for reading! If you want to hire me as your next shinning post-doc, do let reach out!

Automotive Computational Fluid Dynamics (CFD) Analysis

     This post is about the numerical analysis of an Ahmed body. It was a new experience because of the area between the car's floor and the road which is different from most of the numerical analysis performed in open channel aeronautics and turbo-machinery.

     The numerical analysis was performed using the commercial software, SolidWorks Flow Simulation. The software 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.

     The software generates Cartesian mesh using immersed boundary method. The mesh had a cell size of 0.035 m in the far field regions within the computational domain. A fine mesh was need between the road the the car's floor to make sure the interaction of the car's floor with the road was captured accurately. Therefore the mesh between the car's floor and the road was refined to have a cell size of 0.00875 m. Another mesh control was applied around the body to refine the mesh with a cells size of 0.0175 m to capture the trailing vortices. The resulting mesh had 209,580 total cells, among those cells, 31,783 cells were at the solid fluid boundary. The computational domain size was ~1L x 1.12L x 3L where L being the vehicle's length. The computational domain along with the computational mesh is shown in Fig. 1.

Fig. 1 Mesh, computational domain and the boundary conditions.

     The red arrows within the Fig. 1 represents the inlet boundary condition of ambient (free-stream) velocity and the blue arrows represent the outlet boundary condition of the ambient pressure. The green arrows represents the co-ordinates axes direction.

     The results from the numerical analysis were compared with [1-3]. The results are within 10% of the experimental results. The velocity (superimposed by the velocity streamlines) and pressure profiles around the car body at various free-stream velocities is shown in Fig. 2.

Fig. 3 Velocity and pressure plots. From the top, Row 1, L-R; ambient velocity of 30 and 40 m/s. Row 2, L-R; ambient velocity of 60 and 80 m/s. Row 3, ambient velocity of 105 m/s.

     It was a good experience learning about automotive CFD after spending a long time in aeronautic/turbo-machinery CFD. Thank you for reading. Please share my work. If you would like to collaborate on a project please reach out.

[1] F.J.Bello-Millán, T.Mäkelä, L.Parras, C.delPino, C.Ferrera, "Experimental study on Ahmed's body drag coefficient for different yaw angles", Journal of Wind Engineering and Industrial Aerodynamics, Volume 157, October 2016, Pages 140-144.

[2] Guilmineau E., Deng G.B., Queutey P., Visonneau M. (2018) Assessment of Hybrid LES Formulations for Flow Simulation Around the Ahmed Body. In: Deville M. et al. (eds) Turbulence and Interactions. TI 2015. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol 135. Springer, Cham.

[3] A. Thacker, S.Aubrun, A.Leroy, P.Devinant, "Effects of suppressing the 3D separation on the rear slant on the flow structures around an Ahmed body", Journal of Wind Engineering and Industrial Aerodynamics, Volumes 107–108, August–September 2012, Pages 237-243.

Steady-State VS Transient Propeller Numerical Simulation Comparison

     This post is about the comparison between steady-state and transient computational fluid dynamics analysis of two different propellers. The propellers under investigation are 11x7 and 11x4.7 propellers. The first number in the propeller nomenclature is the propeller diameter and the second number represents the propeller pitch, both parameters are in inch. The transient analysis was carried out using the sliding mesh technique while the steady-state results were obtained by the local rotating region-averaging method. For details about 11x7 propeller click here, for the details about 11x4.7 propeller, click here.

 

     As expected, the propeller efficiencies of transient and steady-state analysis are within 0.9% of each other, as shown in Fig. 1-2. Therefore, it is advised to simulate propellers and horizontal axis wind turbines using the steady-state technique as long as no time-dependent boundary conditions are employed.

 

Fig. 1, Propeller efficiency plot.

  

 Fig. 2, Propeller efficiency plot.

 

     It can be seen from Fig. 3-4 that time taken by the steady-state simulation to converge is on average 42.37% less that the transient analysis.  The steady-state analysis takes considerably less time to give a solution then a transient analysis.

 

Fig. 3, Solution time.

 

Fig. 4, Solution time.

 

Thank you for reading. If you would like to collaborate on research projects, please reach out.