Showing posts with label Simulation. Show all posts

Monday, 17 November 2025

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.

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[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[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!

Monday, 7 July 2025

A GPU Accelerated Simplified Immersed Boundary Method using Ray Casting

Yours truly is an avid gamer, @fadoobaba (YouTube). The ray casting algorithm🛸; a fundamental 🧱 technique used in video game 🎮 development and computer graphics, has been implemented within the finite difference method 🏁 code yours truly has been developing. In abundant spare time, of course😼. In this post, this method is explained. For details about ray casting using matplotlib.path, refer to here. For validation of the code, refer to here, backwards-facing step, curved boundaries, here (moving cylinder) and here.

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

If you plan to use these codes in your scholarly work, do cite this blog as:

Fahad Butt (2025). S-IBM (https://fluiddynamicscomputer.blogspot.com/2025/07/a-simplified-immersed-boundary-method.html), Blogger. Retrieved Month Date, Year

The first step is to setup the polygon 🐠. The polygon = cp.column_stack((x1, y1)) statement is used to combine x1 and y1 into a single array representing the polygon vertices. The following statements are used to store the x and y coordinates of the polygon vertices and nvert is the total number of vertices.

px = polygon[:, 0]

py = polygon[:, 1]

nvert = len(polygon)

The boolean masks 👺 are then initialized. The following arrays track whether a grid point is inside the polygon based on horizontal and vertical ray intersections. False = outside ❌, True = inside ✅. Then a loop is implemented to for each edge 📏 of the polygon, defined by vertices i (current) and j (previous), with % nvert ensuring the polygon is closed (last vertex connects back to the first). (xi, yi) and (xj, yj) define the current edge and previous edge.

horizontal_inside = cp.zeros_like(test_x, dtype=bool)

vertical_inside = cp.zeros_like(test_x, dtype=bool)

for i in range(nvert):

j = (i - 1) % nvert

xi, yi = px[i], py[i]

xj, yj = px[j], py[j]

A check ⁉️ is performed to ensure a ray crosses the polygon edge once. cond1 statement checks if the test point lies between the y-values of the edge's endpoints i.e., a ray could cross it. intersect_x finds where the edge crosses a horizontal line at test_y. cond2 checks if the intersection lies to the right of the test point. ^= is XOR toggles the "inside" state each time the ray crosses an edge. The addition of small term ~1e-16 prevents division by zero for horizontal edges. Similar method is applied to verify the points using vertical ray. cond3 checks if the edge crosses a vertical ray (top to bottom) from (test_x, test_y). slope1 and intersect_y computes where the vertical ray at x=test_x intersects the edge. A point is inside only if both horizontal and vertical rays classify it as inside. A point is considered inside the foil only if both ray checks are true. curve is the 2D boolean mask of grid points inside the foil body.

cond1 = ((yi > test_y) != (yj > test_y))

slope = (xj - xi) / (yj - yi + 1e-16)

intersect_x = slope * (test_y - yi) + xi

cond2 = test_x < intersect_x

horizontal_inside ^= cond1 & cond2

cond3 = ((xi > test_x) != (xj > test_x))

slope1 = (yj - yi) / (xj - xi + 1e-16)

intersect_y = slope1 * (test_x - xi) + yi

cond4 = test_y < intersect_y

vertical_inside ^= cond3 & cond4

inside = horizontal_inside & vertical_inside

curve = inside.reshape(X.T.shape)

Using both horizontal and vertical rays reduces false positives (e.g., near sharp corners). The & ensures only points unambiguously inside are marked. Within Fig. 1, the points on the shape boundary, points inside 🪰 and outside 🐌 the shape boundary are shown. A point is a boundary if it is inside the body, but any of its neighbors is outside. Following statements are used to mark the boundary of the object. Boundary detects the "skin" of the body for applying no-slip, force, or stress conditions.

interior = curve[1:-1, 1:-1]

right = curve[2:, 1:-1]

left = curve[:-2, 1:-1]

top = curve[1:-1, 2:]

bottom = curve[1:-1, :-2]

boundary = interior & ~(right & left & top & bottom)

The statements boundary_indices = cp.where(boundary_mask) ... valid_right = (boundary_indices[0] + 1 < X.shape[0]) & (~curve[right_neighbors]) ... etc. extract boundary indices and their valid neighbors i.e. these lines extract the (i, j) grid indices of the surface and locate which neighbor cells are valid fluid neighbors (outside the body and within domain bounds). These are needed to compute normals or apply boundary conditions via interpolation or extrapolation from the fluid.

Fig. 1, Mesh cells

For validation of the results from present simulations, the case of flow around a circular cylinder is selected. Fig. 1 shows the results from the code at Re 200. The drag coefficient obtained from this code is 1.396 while from the literature, the value is at 1.4 [1]. The lift coefficient from the code is at 0.000134. Within Fig. 2, top row shows u and v velocity components and bottom row shows pressure and vorticity.

Fig. 2, The flow-field

For the second an more complex validation case, a swimming fish is simulated. The lift and drag coefficients obtained from the simulation are compared with experimental results. The drag coefficient from the current code is at 0.342 while from the published literature, the value is at 0.348. Maximum lift coefficient is at 7.243 and 8 from the present code as compared to the published literature [2]. Within Fig. 3 the u and v velocities, pressure and velocity streamlines are shown for St = 0.8 and Reynolds number of 500. The computational mesh near the fish is shown in Fig. 4. Within Fig.4, a zoomed in view towards the right shows the mesh at the trailing edge of the fish.

Fig.3, Post-processing of results

Fig. 4, The mesh

This method allows handling arbitrary deforming / non-deforming shapes on a fixed Cartesian grid. In summary, the method has the following steps.

1. Generate polygon shape (e.g., airfoil, cylinder)

2. Flatten mesh for vectorized testing

3. Use ray casting to check if points are inside shape

4. Build a boolean mask of body region

5. Identify surface (boundary) points

6. Extract boundary indices for physics coupling

References

[1] Braza M, Chassaing P, Minh HH. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics. 1986;165:79-130. doi:10.1017/S0022112086003014

[2] Fulong Shi, Jianjian Xin, Chuanzhong Ou, Zhiwei Li, Xing Chang, Ling Wan; Effects of the Reynolds number and attack angle on wake dynamics of fish swimming in oblique flows. Physics of Fluids 1 February 2025; 37 (2): 025205 doi.org/10.1063/5.0252506

Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!

Thursday, 20 February 2025

Flow Simulation around Shapes (Includes Free Code)

This post is about the simulation of flow 🍃 around and through various objects 🪈 and obstacles ⭕. The simulated cases include flow through a partially blocked pipe, shown in Fig. 1. Flow around flat plates arranged in the shape of the letter "T", shown in Fig. 2. Flow around a wedge and flow around a triangle 🔺, which are shown in Fig. 3 and 4 respectively.

To create obstacles in Python 🐍, the Path statement is used. This is similar to the inpolygon statement in MATLAB 🧮. For example, the following code 🖳 is used to create an ellipse ⬭.

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

theta = np.linspace(0, np.pi, 100) # create equally spaced angles from 0 to pi

a = 0.05 # semi-major axis (along x)

b = 0.025 # semi-minor axis (along y)

shape_center_x = L / 2 # shape center x

shape_center_y = 0 # ellipse center x

x1 = shape_center_x + (a * np.cos(theta)) # x-coordinates

y1 = shape_center_y + (b * np.sin(theta)) # y-coordinates

shape_path = Path(np.column_stack((x1, y1))) # shape region

points = np.vstack((X.ravel(), Y.ravel())).T # mark inside region

custom_curve = shape_path.contains_points(points).reshape(X.shape) # create boolean mask

custom_curve_boundary = np.zeros_like(custom_curve, dtype=bool) # find boundary of shape

for i in range(1, Nx-1):

for j in range(1, Ny-1):

if custom_curve.T[i, j]: # inside shape

if (not custom_curve.T[i+1, j] or not custom_curve.T[i-1, j] or

not custom_curve.T[i, j+1] or not custom_curve.T[i, j-1]):

custom_curve_boundary.T[i, j] = True # mark as boundary

boundary_indices = np.where(custom_curve_boundary.T) # mark boundary

right_neighbors = (boundary_indices[0] + 1, boundary_indices[1]) # right index

left_neighbors = (boundary_indices[0] - 1, boundary_indices[1]) # left index

top_neighbors = (boundary_indices[0], boundary_indices[1] + 1) # top index

bottom_neighbors = (boundary_indices[0], boundary_indices[1] - 1) # bottom index

valid_right = ~custom_curve.T[right_neighbors] # check if points are on shape boundary

valid_left = ~custom_curve.T[left_neighbors]

valid_top = ~custom_curve.T[top_neighbors]

valid_bottom = ~custom_curve.T[bottom_neighbors]

The arrays x1 and y1 are created using the equation for ellipse with the required parameters. The ellipse region is marked using the Path statement. Indices inside the ellipse are marked using a Boolean mask. Nested loops are used to mark the ellipse boundary using a curve. The right, left, top and bottom neighbor points are identified to be used for the application of Neumann boundary conditions for pressure (no-slip). Within the time loop, "where" statement is used to identify and apply the Neumann boundary conditions on the ellipse wall. The same method applied on any shape, for example aero-foils, circles, nozzles etc. The complete code to reproduce Fig. 1 is made available 😇.

It should be noted that, that the mesh is still stairstep. It doesn't matter how small the mesh resolution 😲. This code is developed for educational and research purposes only as there is not much application to stairstep mesh in real world❗

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 libraries

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.path import Path

#%% define parameters

l_square = 1 # length of square

h = 0.025 / 25 # grid spacing

dt = 0.00001 # time step

L = 0.3 # domain length

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

#%% initialize variables

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

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

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

#%% create a shape

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

theta = np.linspace(0, np.pi, 100) # create equally spaced angles from 0 to pi

a = 0.05 # semi-major axis (along x)

b = 0.025 # semi-minor axis (along y)

shape_center_x = L / 2 # shape center x

shape_center_y = 0 # ellipse center x

x1 = shape_center_x + (a * np.cos(theta)) # x-coordinates

y1 = shape_center_y + (b * np.sin(theta)) # y-coordinates

shape_path = Path(np.column_stack((x1, y1))) # shape region

points = np.vstack((X.ravel(), Y.ravel())).T # mark inside region

custom_curve = shape_path.contains_points(points).reshape(X.shape) # create boolean mask

custom_curve_boundary = np.zeros_like(custom_curve, dtype=bool) # find boundary of shape

for i in range(1, Nx-1):

for j in range(1, Ny-1):

if custom_curve.T[i, j]: # inside shape

if (not custom_curve.T[i+1, j] or not custom_curve.T[i-1, j] or

not custom_curve.T[i, j+1] or not custom_curve.T[i, j-1]):

custom_curve_boundary.T[i, j] = True # mark as boundary

boundary_indices = np.where(custom_curve_boundary.T) # mark boundary

right_neighbors = (boundary_indices[0] + 1, boundary_indices[1]) # right index

left_neighbors = (boundary_indices[0] - 1, boundary_indices[1]) # left index

top_neighbors = (boundary_indices[0], boundary_indices[1] + 1) # top index

bottom_neighbors = (boundary_indices[0], boundary_indices[1] - 1) # bottom index

valid_right = ~custom_curve.T[right_neighbors] # check if points are on shape boundary

valid_left = ~custom_curve.T[left_neighbors]

valid_top = ~custom_curve.T[top_neighbors]

valid_bottom = ~custom_curve.T[bottom_neighbors]

#%% solve 2D Navier-Stokes equations

for nt in range(ns):

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

# apply pressure boundary conditions

p[0, :] = p[1, :] # dp/dx = 0 at x = 0

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

p[:, 0] = p[:, 1] # dp/dy = 0 at y = 0

p[:, -1] = p[:, -2] # dp/dy = 0 at y = D / 2

p[custom_curve.T] = 0 # p = 0 inside shape

# dp/dn = 0 at shape boundary (no slip)

p[boundary_indices] = np.where(valid_right, p[right_neighbors], p[boundary_indices]) # right neighbor is fluid

p[boundary_indices] = np.where(valid_left, p[left_neighbors], p[boundary_indices]) # left neighbor is fluid

p[boundary_indices] = np.where(valid_top, p[top_neighbors], p[boundary_indices]) # top neighbor is fluid

p[boundary_indices] = np.where(valid_bottom, p[bottom_neighbors], p[boundary_indices]) # bottom neighbor is fluid

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 boundary conditions

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

u[-1, :] = u[-2, :] # du/dx = 0 at x = L

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

u[:, -1] = u[:, -2] # du/dy = 0 at y = D / 2

u[custom_curve.T] = 0 # u = 0 inside shape

u[custom_curve_boundary.T] = 0 # no slip

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 boundary conditions

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 / 2

v[custom_curve.T] = 0 # u = 0 inside shape

v[custom_curve_boundary.T] = 0 # no slip

#%% post process

u1 = u.copy() # u-velocity for plotting with shape

v1 = v.copy() # v-velocity for plotting with shape

p1 = p.copy() # pressure for plotting with shape

# shape geometry for plotting

u1[custom_curve.T] = np.nan

v1[custom_curve.T] = np.nan

p1[custom_curve.T] = np.nan

velocity_magnitude1 = np.sqrt(u1**2 + v1**2) # velocity magnitude with shape

# visualize velocity vectors and pressure contours

plt.figure(dpi = 500)

plt.contourf(X, Y, u1.T, 128, cmap = 'jet')

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

plt.contourf(X, Y, v1.T, 128, cmap = 'jet')

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

plt.contourf(X, Y, p1.T, 128, cmap = 'jet')

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

plt.streamplot(X, Y, u1.T, v1.T, color = 'black', cmap = 'jet', density = 2, linewidth = 0.5, arrowstyle='->', arrowsize = 1) # plot streamlines

plt.plot(x1, y1, color='black', alpha = 1, linewidth = 2)

plt.gca().set_aspect('equal', adjustable = 'box')

plt.axis('off')

plt.show()

Fig. 1, Flow in a clogged pipe

Fig. 2, Flow around a blunt obstacle

Fig. 3, Flow around asymmetric wedge

Fig. 4, Flow around a triangle

Within Figs. 1 - 4, top row shows u and v components of velocity, bottom row shows pressure and streamlines 💫.

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

Wednesday, 8 January 2025

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!

The new version of the code is faster as the equations are simplified and many factors are precalculated, resulting in quicker execution times.

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[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[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 💫.

Artificial Compressibility

The artificial compressibility method is now implemented in the code. The output is flow inside the lid-driven cavity at Reynolds number 10,000. The smallest scale resolved is at 0.001 m and smallest time scale resolved is at 0.0001 s. This version of code seems to be more stable as compared to the one that uses pressure Poisson equation. The results are shown in Fig. 5.

Fig. 5, Look at all those secondary vortices 😚

#Copyright <2025> <FAHAD BUTT>

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#%% import libraries

import numpy as np

import matplotlib.pyplot as plt

#%% define parameters

l_cr = 1 # characteristic length

h = l_cr / 1000 # grid spacing

dt = 0.0001 # 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 / 10000 # kinematic viscosity

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

cfl = dt * Uinf / h # cfl number

travel = 20 # 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 = np.zeros((Nx, Ny)) # x-velocity

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

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

#%% pre calculate for speed

P1 = dt / h

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

P3 = 1 - (4 * P2)

#%% solve 2D Navier-Stokes equations

for nt in range(ns):

pn = p.copy()

p[1:-1, 1:-1] = pn[1:-1, 1:-1] - P1 * (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] * P3 - P1 * (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 = 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] * P3 - P1 * (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

#%% post processing

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

plt.figure(dpi = 200)

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

plt.colorbar()

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()

plt.figure(dpi = 200)

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

plt.colorbar()

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()

plt.figure(dpi = 200)

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

plt.colorbar()

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()

plt.figure(dpi = 200)

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()

velocity_magnitude = np.sqrt(u**2 + v**2) # calculate velocity magnitude

plt.figure(dpi = 200)

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

plt.streamplot(X, Y, u.T, v.T, color = 'black', cmap = 'jet', density = 2, linewidth = 0.1, arrowstyle='->', arrowsize = 0.5) # 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.axis("off")

plt.show()

Relaxation

The relaxation parameter along with equations is now added! The "omega" parameter can be adjusted to over / under relax the simulation according to requirements! At Reynolds number of 5000, the code provides up to 4x speed in convergence with over-relaxation.

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#%% import libraries

import numpy as np

import matplotlib.pyplot as plt

#%% define parameters

l_cr = 1 # characteristic length

h = l_cr / 600 # grid spacing

dt = 0.0001 # 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 / 5000 # 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

#%% intialization

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

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

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

u_new = u.copy() # x-velocity (relaxation)

v_new = v.copy() # y-velocity (relaxation)

p_new = p.copy() # pressure (relaxation)

#%% pre calculate for speed

P1 = dt / h

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

P3 = 1 - (4 * P2)

omega = 4 # relaxation parameter (omega < 5)

#%% solve 2D Navier-Stokes equations

for nt in range(ns):

pn = p.copy()

p_new[1:-1, 1:-1] = pn[1:-1, 1:-1] - P1 * (u[2:, 1:-1] - u[:-2, 1:-1] + v[1:-1, 2:] - v[1:-1, :-2]) # pressure

p[1:-1, 1:-1] = (1 - omega) * pn[1:-1, 1:-1] + omega * p_new[1:-1, 1:-1] # relaxation

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_new[1:-1, 1:-1] = un[1:-1, 1:-1] * P3 - P1 * (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[1:-1, 1:-1] = (1 - omega) * un[1:-1, 1:-1] + omega * u_new[1:-1, 1:-1] # relaxation

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_new[1:-1, 1:-1] = vn[1:-1, 1:-1] * P3 - P1 * (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[1:-1, 1:-1] = (1 - omega) * vn[1:-1, 1:-1] + omega * v_new[1:-1, 1:-1] # relaxation

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

#%% post processing

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

plt.figure(dpi = 200)

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

plt.colorbar()

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()

plt.figure(dpi = 200)

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

plt.colorbar()

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()

plt.figure(dpi = 200)

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

plt.colorbar()

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()

plt.figure(dpi = 200)

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()

velocity_magnitude = np.sqrt(u**2 + v**2) # calculate velocity magnitude

plt.figure(dpi = 200)

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

plt.streamplot(X, Y, u.T, v.T, color = 'black', cmap = 'jet', density = 2, linewidth = 0.1, arrowstyle='->', arrowsize = 0.5) # 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.axis("off")

plt.show()

Pre-Calculation

# 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_cr = 1 # characteristic length

h = l_cr / 120 # grid spacing

dt = 0.0001 # 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 = 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

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

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

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

P1 = h / (16 * dt)

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

P3 = dt / h

P4 = 1 - (4 * P2)

for nt 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]) # 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[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 = 0 at y = D

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 = 500)

plt.contourf(X, Y, u.T, 128, cmap = 'jet')

plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.xticks([0, L])

plt.yticks([0, D])

plt.xlabel('x [m]')

plt.ylabel('y [m]')

plt.axis('off')

plt.show()

plt.figure(dpi = 500)

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.xlim([0, L])

plt.ylim([0, D])

plt.axis('off')

plt.show()

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

Wednesday, 18 August 2021

Computational Fluid Dynamics Simulation of a Swimming Fish (Includes UDF)

This post is about the simulation of a swimming fish. The fish body is made of NACA 0020 and 0015 aero-foils (air-foils). The fluke is made of NACA 0025 aero-foil (air-foil), as shown in Fig. 1. the CAD files with computational domain modelled around the fish is available here. A tutorial is available here.

Fig. 1, The generic fish CAD model

The motion of the fish's body is achieved using a combination of two user-defined functions (UDF). The UDFs use DEFINE_GRID_MOTION script mentioned below, for the head/front portion. This is taken from the ANSYS Fluent software manual, available in its original form here. The original UDF is modified for present use as required. To move the mesh, dynamic mesh option within ANSYS Fluent is enabled; with smoothing and re-meshing options. The period of oscillation is kept at 2.0 s. The Reynolds number of flow is kept at 100,000; which is typical for a swimming fish.

/**********************************************************

node motion based on simple beam deflection equation

compiled UDF

**********************************************************/

#include "udf.h"

DEFINE_GRID_MOTION(undulating_head,domain,dt,time,dtime)

{

Thread *tf = DT_THREAD(dt);

face_t f;

Node *v;

real NV_VEC(omega), NV_VEC(axis), NV_VEC(dx);

real NV_VEC(origin), NV_VEC(rvec);

real sign;

int n;

/* set deforming flag on adjacent cell zone */

SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf));

sign = 0.15707963267948966192313216916398 * cos (3.1415926535897932384626433832795 * time);

Message ("time = %f, omega = %f\n", time, sign);

NV_S(omega, =, 0.0);

NV_D(axis, =, 0.0, 1.0, 0.0);

NV_D(origin, =, 0.7, 0.0, 0.0);

begin_f_loop(f,tf)

{

f_node_loop(f,tf,n)

{

v = F_NODE(f,tf,n);

/* update node if x position is greater than 0.02

and that the current node has not been previously

visited when looping through previous faces */

if (NODE_X(v) > 0.05 && NODE_X(v) < 0.7 && NODE_POS_NEED_UPDATE (v))

{

/* indicate that node position has been update

so that it's not updated more than once */

NODE_POS_UPDATED(v);

omega[1] = sign * pow (NODE_X(v), 0.5);

NV_VV(rvec, =, NODE_COORD(v), -, origin);

NV_CROSS(dx, omega, rvec);

NV_S(dx, *=, dtime);

NV_V(NODE_COORD(v), +=, dx);

}

end_f_loop(f,tf);

}

The computational mesh, as shown in Fig. 2, uses cut-cell method with inflation layers. The mesh has 2,633,133 cells. The near wall y+ is kept at 5. The Spalart-Allmaras turbulence model is used to model the turbulence. The second order upwind scheme is used to discretize the momentum and modified turbulent viscosity equations. The time-step for this study is kept at 100th/period of oscillation.

Fig. 2, The mesh and zoom in view of the trailing edge.

The animation showing fish motion is shown in Fig. 3. Within Fig. 3, the left side showcases the velocity iso-surfaces coloured by pressure and the vorticity iso-surfaces coloured by velocity magnitude is shown on the right.

Fig. 3, The animation.

Another animation showing the fish motion is shown in Fig.4. Within Fig. 4, the left side shows surface pressure while the right side shows pressure iso-surfaces coloured by vorticity.

Fig. 4, The animation.

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

Fluid Dynamics using the Computer

My Blogger Profile

Monday, 17 November 2025

A Simple Poisson's Equation for Pressure

Code 01

Code 2

Monday, 7 July 2025

A GPU Accelerated Simplified Immersed Boundary Method using Ray Casting

References

Thursday, 20 February 2025

Flow Simulation around Shapes (Includes Free Code)

Code

Wednesday, 8 January 2025

CFD Wizardry: A 50-Line Python Marvel

Code

Lid-Driven Cavity

Free-Jet

Heated Room

Backward - Facing Step (BFS)

Artificial Compressibility

Relaxation

Pre-Calculation

Wednesday, 18 August 2021

Computational Fluid Dynamics Simulation of a Swimming Fish (Includes UDF)