15th Step of the 12 steps to Navier-Stokes 😑

     FOANSS, i.e. Fadoobaba's Open Advanced Navier-Stokes 🌬 Solver is now capable to simulate flow around curved geometries. In this example, flow around a circular cylinder ⭕ is presented. The main challenge in handling the curved boundaries 🗻 in discrete world is application of the Neumann boundary condition for pressure. The method employed in FOANSS is to use polar ⭗ coordinate system for normal derivative calculation. The idea came while yours's truly was developing FOAMNE i.e. Fadoobaba's Open Advanced 🔨 Mechanics 🧠 Neural Engine. More details can be read here.

     In the CFD code, ∂p/dn is replaced by ∂p/dx * cos(θ) + ∂p/dy * sin(θ). After some trickery, the pressure on the curved surface can be represented by equation 1. Within equation 1, p is the pressure at current time-step. i and j are nodes in the x and y axis, respectively. cos_theta and sin_theta represent radial angles of the polar coordinate system (r θ).

((p[i + 1, j] + p[i - 1, j]) * cos_theta + (p[i, j + 1] + p[i, j - 1] * sin_theta)) / (2 * (cos_theta + sin_theta)) (1)

     At first, flow around a circular cylinder at Reynolds number 100 is simulated. The accuracy of the results is measured using the Strouhal number. For Reynolds number 100, Strouhal number for flow around a circular cylinder is 0.16. The Strouhal number obtained from the code on a coarse mesh, i.e. only 50 cells per cylinder diameter, is 0.17. This is within 8% of the published data. The Strouhal number will never reach 0.16 as the mesh is cartesian and the cylinder surface is represented as a stairstep. Still, very good for teaching purposes! The results are shown in Fig. 1. It can be seen that the vortex shedding phenomenon is calculated relatively accurately. A remainder that you are reading a blog, not a Q1 journal 😃. The vorticity plot is shown in Fig. 2.

Fig. 1, The results from the code

Code

#%% import libraries

import numpy as np

import matplotlib.pyplot as plt

#%% define spatial and temporal grids

l_square = 1 # length of square or diamater of circle

h = l_square / 50 # grid spacing

dt = 0.0001 # time step

L = 15 # domain length

D = 10 # 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 = 2 # 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

rho = 1 # fluid density

#%% initialize flowfield

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

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

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

x = np.linspace(0, L, Nx) # x-axis vector

y = np.linspace(0, D, Ny) # y-axis vector

XX, YY = np.meshgrid(x, y, indexing="ij") # grid for circle boundary

#%% create circle

x_c = 5 # circle center y

y_c = 5 # circle center y

radius = 0.5 # circle radius

mask = (XX - x_c)**2 + (YY - y_c)**2 <= radius**2 # circular region

tol = 0.125 * h # tolerance for boundary points

circle_boundary = np.abs(np.sqrt((XX - x_c)**2 + (YY - y_c)**2) - radius) <= tol # boundary mask

boundary_indices = np.argwhere(circle_boundary) # circle boundary

# compute radial distances and angles (polar coordinates)

r_values = []

cos_theta_values = []

sin_theta_values = []

for idx in boundary_indices:

    i, j = idx

    r = np.sqrt((XX[i, j] - x_c)**2 + (YY[i, j] - y_c)**2)

    cos_theta = (XX[i, j] - x_c) / r

    sin_theta = (YY[i, j] - y_c) / r

    r_values.append(r)

    cos_theta_values.append(cos_theta)

    sin_theta_values.append(sin_theta)

# define quadrant masks (circle is divided in 4 parts) and corner points

quadrant1 = (XX > x_c) & (YY > y_c) & circle_boundary

quadrant2 = (XX > x_c) & (YY < y_c) & circle_boundary

quadrant3 = (XX < x_c) & (YY < y_c) & circle_boundary

quadrant4 = (XX < x_c) & (YY > y_c) & circle_boundary

top_pt= (XX == x_c) & (YY > y_c) & circle_boundary

bottom_pt = (XX == x_c) & (YY < y_c) & circle_boundary

left_pt = (XX < x_c) & (YY == y_c) & circle_boundary

right_pt = (XX > x_c) & (YY == y_c) & circle_boundary

#%% 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 * rho / (8 * dt) * (u[2:, 1:-1] - u[:-2, 1:-1] + v[1:-1, 2:] - v[1:-1, :-2]) # pressure

    # boundary conditions for pressure

    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

    p[mask] = 0 # circle

    # apply Neumann boundary condition for pressure in polar coordinates

    p[quadrant1] = (((p[i + 1, j] + p[i - 1, j]) * cos_theta) + ((p[i, j + 1] + p[i, j - 1]) * sin_theta)) / (2 * (cos_theta + sin_theta))

    p[quadrant2] = (((p[i + 1, j] + p[i - 1, j]) * cos_theta) + ((p[i, j + 1] + p[i, j - 1]) * sin_theta)) / (2 * (cos_theta + sin_theta))

    p[quadrant3] = (((p[i + 1, j] + p[i - 1, j]) * cos_theta) + ((p[i, j + 1] + p[i, j - 1]) * sin_theta)) / (2 * (cos_theta + sin_theta))

    p[quadrant4] = (((p[i + 1, j] + p[i - 1, j]) * cos_theta) + ((p[i, j + 1] + p[i, j - 1]) * sin_theta)) / (2 * (cos_theta + sin_theta))

    p[top_pt] = p[i, j + 1]

    p[bottom_pt] = p[i, j - 1]

    p[left_pt] = p[i - 1, j]

    p[right_pt] = p[i + 1, j]

    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 * rho * h) * (p[2:, 1:-1] - p[:-2, 1:-1]) + nu * 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

    # boundary conditions for x-velocity

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

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

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

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

    u[mask] = 0 # circle

    u[circle_boundary] = 0 # circle

    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 * rho * h) * (p[1:-1, 2:] - p[1:-1, :-2]) + nu * 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

    # boundary conditions for y-velocity

    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

    v[mask] = 0 # circle

    v[circle_boundary] = 0 # circle

#%% post-processing

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

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

# visualize velocity vectors and pressure contours

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.colorbar(orientation='vertical')

plt.show()

plt.figure(dpi = 500)

plt.contourf(X, Y, v.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.colorbar(orientation='vertical')

plt.show()

plt.figure(dpi = 500)

plt.contourf(X, Y, p.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.colorbar(orientation='vertical')

plt.show()

plt.figure(dpi = 500)  # make a nice crisp image :)

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.axis('off')

plt.show()

Fig. 2, Vorticity with streamlines

     Thank you very much for reading! If you want to hire me as your new shinning post-doc or want to collaborate on the research, please reach out! The code might be available soon!

Airflow Simulation in Empty / Occupied Rooms and Environments

     One fine morning, I decided to code the Navier–Stokes equations [read if you are bored 🤣] . This post has the results of  flow simulation inside close environments of various aspect ratios. As is customary with all my CFD work using commercial and 🏡made CFD codes, this too is inspired by the free lectures of Dr. Lorena Barba.

The first case has an aspect ratio of 1:1 while for the second case, the aspect ratio is at 3:1. The airflow is at 0.4555 m/s for both cases. Both cases are isothermal at 293 K. There is no turbulence 🌬model [free code ] 🤑.

The smallest resolved scale (~4x smallest mesh size) for 1:1 case is ~ 0.0045 m and for the 3:1 case is at 0.0112 m. Time scales (~4x time-step size) are at 0.0004 s and 0.0004 s, respectively. Fig. 1 shows results for 1:1 aspect ratio while the results for 3:1 case are shown in Fig. 2. For both cases, flow enters from top-left and exists at bottom-right of the rooms. The boundary conditions are taken from [1]. I compared the results with Fluent simulations I ran at same boundary conditions and stopped my simulations when eye-balling didn't revealed any difference 😆. What you expect? This is a blog not a journal 😝.

Fig. 1, 1:1 aspect ratio

Fig. 2, 3:1 aspect ratio

 

Fig. 3, 3:1 aspect ratio with partition

Fig. 4, Flow inside ducts

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

     [1] Horikiri, Kana & Yao, Yufeng. Validation Study of Convective Airflow in an Empty Room, "Recent Researches in Energy, Environment, Devices, Systems, Communications and Computers", ISBN: 978-960-474-284-4

Lid-Driven Cavity MATLAB Code

     MATLAB code for 2D Lid-Driven Cavity. Includes labeled commands, plotting and is less than 100 lines of code. Resume possible, to resume comment close all, clear, clc, u, v, p and then run the code. Results are available here:

%% clear and close

close all

clear

clc

%% define spatial and temporal grids

h = 1/10; % grid spacing

cfl = h; % cfl number

L = 1; % cavity length

D = 1; % cavity depth

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

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

nu = 0.000015111; % kinematic viscosity

Uinf = 0.0015111; % free stream velocity / lid velocity

dt = h * cfl / Uinf; % time step

travel = 2; % times the disturbance travels entire length of computational domain

TT = travel * L / Uinf; % total time

ns = TT / dt; % number of time steps

l_square = 1; % length of square

Re = l_square * Uinf / nu; % Reynolds number

rho = 1.2047; % fluid density

%% initialize flowfield

u = zeros(Nx,Ny); % x-velocity

v = zeros(Nx,Ny); % y-velocity

p = zeros(Nx,Ny); % pressure

i = 2:Nx-1; % spatial interior nodes in x-axis

j = 2:Ny-1; % spatial interior nodes in y-axis

[X, Y] = meshgrid(0:h:L, 0:h:D); % spatial grid

maxNumCompThreads('automatic'); % select CPU cores

%% solve 2D Navier-Stokes equations

for nt = 1:ns

    pn = p;

    p(i, j) = ((pn(i+1, j) + pn(i-1, j)) * h^2 + (pn(i, j+1) + pn(i, j-1)) * h^2) ./ (2 * (h^2 + h^2)) ...

        - h^2 * h^2 / (2 * (h^2 + h^2)) * (rho * (1 / dt * ((u(i+1, j) - u(i-1, j)) / (2 * h) + (v(i, j+1) - v(i, j-1)) / (2 * h)))); % pressure poisson

    p(1, :) = p(2, :); % dp/dx = 0 at x = 0

    p(Nx, :) = p(Nx-1, :); % dp/dx = 0 at x = L

    p(:, 1) = p(:, 2); % dp/dy = 0 at y = 0

    p(:, Ny) = 0; % p = 0 at y = D 

    un = u;

    vn = v;

    u(i, j) = un(i, j) - un(i, j) * dt / (2 * h) .* (un(i+1, j) - un(i-1, j)) ...

        - vn(i, j) * dt / (2 * h) .* (un(i, j+1) - un(i, j-1)) - dt / (2 * rho * h) * (p(i+1, j) - p(i-1, j)) ...

        + nu * (dt / h^2 * (un(i+1, j) - 2 * un(i, j) + un(i-1, j)) + dt / h^2 * (un(i, j+1) - 2 * un(i, j) + un(i, j-1))); % x-momentum

    u(1, :) = 0; % u = 0 at x = 0

    u(Nx, :) = 0; % u = 0 at x = L

    u(:, 1) = 0; % u = 0 at y = 0

    u(:, Ny) = Uinf; % u = Uinf at y = D

    v(i, j) = vn(i, j) - un(i, j) * dt / (2 * h) .* (vn(i+1, j) - vn(i-1, j)) ...

        - vn(i, j) * dt / (2 * h) .* (vn(i, j+1) - vn(i, j-1)) - dt / (2  * rho * h) * (p(i, j+1) - p(i, j-1)) ...

        + nu * (dt / h^2 * (vn(i+1, j) - 2 * vn(i, j) + vn(i-1, j)) + dt / h^2 * (vn(i, j+1) - 2 * vn(i, j) + vn(i, j-1))); % y-momentum

    v(1, :) = 0; % v = 0 at x = 0

    v(Nx, :) = 0; % v = 0 at x = L

    v(:, 1) = 0; % v = 0 at y = 0

    v(:, Ny) = 0; % v = 0 at y = D

end

%% post-processing

velocity_magnitude = sqrt(u.^2 + v.^2); % velocity magnitude

% Visualize velocity vectors and pressure contours

hold on

contourf(X, Y, velocity_magnitude', 64, 'LineColor', 'none'); % contour plot

set(gca,'FontSize',40)

% skip = 20;

% quiver(X(1:skip:end, 1:skip:end), Y(1:skip:end, 1:skip:end),... % Velocity vectors

%     u1(1:skip:end, 1:skip:end)', v1(1:skip:end, 1:skip:end)', 1, 'k','LineWidth', 0.1);

hh = streamslice(X, Y, u', v', 5); % Streamlines

set(hh, 'Color', 'k','LineWidth', 0.1);

colorbar; % Add color bar

colormap hsv % Set color map

axis equal % Set true scale

xlim([0 L]); % Set axis limits

ylim([0 D]);

xticks([0 L]) % Set ticks

yticks([0 D]) % Set ticks

clim([0 0.95*max(velocity_magnitude(:))]) % Legend limits

title('Velocity [m/s]');

xlabel('x [m]');

ylabel('y [m]');

Cite as: Fahad Butt (2023). Lid-Driven Cavity (https://fluiddynamicscomputer.blogspot.com/2023/09/lid-driven-cavity-matlab-code.html), Blogger. Retrieved Month Date, Year

Heaving Flat Plate Computational Fluid Dynamics (CFD) Simulation (In-House CFD Code)

     The adventure 🏕 🚵 I started a while ago to make my own CFD 🌬 code / software 💻 for my digital CV and to make a shiny new turbulence mode (one-day perhaps) is going along nicely. This post is about a 2D 10% flat plate undergoing forced heaving motion. Heaving motion is achieved by Eqn. 1.

h_y = H_o*sin(2*π*f_h*t)                                              Eqn. 1

     w.r.t. Eqn. 1 reduced frequency is defined as (f_h*H_o/U∞)), f_h is the frequency of oscillations, while t is the instantaneous time. H_o is the heaving amplitude and U∞ is the free stream velocity. h_y is the position of the flat plate. The animation is shown in Fig. 1.

     The Strouhal number is 0.228 and Reynolds number is at 500. As we can see, the in-house CFD code works very well for this complex CFD simulation. Validation of this work will never be completed 😆. As soon, I will move on to the next project without completing this one. Anyway, discretized Navier-Stokes equations are available here, in both C++ and MATLAB formats if you want to validate this non sense yourself! Good Luck!

The animation from in-house CFD simulation

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

Home-Made Computational Fluid Dynamics (CFD), (Includes CFD code)

     One fine morning, I decided to code the 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 inspired by the free lectures of Dr. Lorena Barba.

Internal / External Fluid Dynamics - Lid-Driven Cavity

     I started using Lid-Driven Cavity example. Just because everyone else uses it to validate the code they write. The lid-driven cavity case is giving correct results up to ~Re 1,000 without any turbulence models or wall functions.

     The results shown in Fig. 1 correspond to at Re 1,000. It can be seen that the present code which is just vanilla Navier–Stokes; captures the vortices at both bottom edges well as compared to published data. But why would you want to use a vanilla CFD code which solved only ~Re 1,000?, anyways... Some validation... u-velocity at (0.5, 0.1719) is at -0.3869 m/s as compared to -0.38289 m/s [1]. Meanwhile, v-velocity at (0.2266, 0.5) is 0.3252 as compared to 0.30203. Furthermore, v-velocity at (0.8594, 0.5) is at -0.4128 m/s as compared to -0.44993 m/s. All in all, a good agreement with published data. An animation is uploaded here.  Fig. 1a shows a different variation of lid driven cavity.

Fig. 1, Lid driven cavity post-processing

Fig. 1a, Heated lid driven cavity post-processing

Internal / External Aerodynamics - Backward-Facing Step (BFS)

     The next example is the case of Backward Facing Step (BFS). This case is analogues to flow around a building or a pipe with different diameter or anything with sharp edge along the flow. Let's be honest, flow around a building 🏢 is at much higher Re than the measly ~1,000 this code solves. But for validation, here we are.

 

     The cases are self validating as the re-attachment length of the vortex just after the step is very well documented [4]. Fig. 2 compares the results from the vanilla code with published data. It can be seen that the results are in good agreement.

Fig. 2, Backward Facing Step (BFS) post processing at Re 200

Internal Fluid Dynamics - Aero-Thermal Pipe Flow

     This is the first case I ran as it is self validating. The product of area and velocity should be same at  the outlet as given at the inlet as a boundary condition, which it is! Furthermore, density calculated should be proportional to the temperature calculated, which it is! The results at Re ~3300 are shown in Fig. 3. This case is relevant to bio medical field as Reynolds number in human blood 🩸flow is around same flow regime. Moreover, flow through heated pipes etc. is also relevant.

Fig. 3, Flow through heated pipe / between flat plates

Clot Flow

     Fig. 4 shows an example of flow where a clot appears in a blood 🩸 vessel. Done on this same 🏡made CFD 🌬️ code used throughout this blog. Or in a pipe where oil 🛢️ flows. It can be seen how smooth /less chaotic the flow is without clot (less work for heart ❤️ ). Eat 🦐 🥑 🥙 🥕 🥒 🍲 healthy! Mechanical properties of fluid are taken as follows. 5 ltr/min flowrate. Blood vessel diameter 40 mm. Fluid density at 1,070 Kg/m³ and dynamic viscosity at 4.5 cP.

Fig. 4, Flow inside a blood vessel with and without clot

External Aerodynamics - Flow around a Square Cylinder

     Why not circular cylinder you ask? I am very lazy 😂. It is hard work to draw a circle using code. I am used to CAD. 😇. Again, low Reynolds number are very rare in practical applications but for the sake of completeness, I added this case as well. Fig. 5 shows flow around the cylinder at Re 100. Vorticity is shown in Fig. 5. The results are compared with experimental data. The Strouhal number from the home-made CFD code is at 0.158 as compared to a value of 0.148. The results are within 7%  of published literature [2-3].

Fig. 5, Post-processing

Thank you for reading! If you want to hire me as your most awesome PhD student, please reach out!

Code:

     I present to you the code-able discretized equations for solving fluid flow problems. At first, C++ version is presented followed by the MATLAB version. Equation 1-2 are Poisson equations for pressure. Equation 3 and 4 are x-momentum equations (without source). Equations 5 and 6 are y-momentum equations.

double p_ij = ((pn(i + 1, j) + pn(i - 1, j)) * dy * dy + (pn(i, j + 1) + pn(i, j - 1)) * dx * dx) / (2 * (dx * dx + dy * dy)) - dx * dx * dy * dy / (2 * (dx * dx + dy * dy)) * (rho * (1 / dt * ((u(i + 1, j) - u(i - 1, j)) / (2 * dx) + (v(i, j + 1) - v(i, j - 1)) / (2 * dy)) - ((u(i + 1, j) - u(i - 1, j)) / (2 * dx)) * ((u(i + 1, j) - u(i - 1, j)) / (2 * dx)) - 2 * ((u(i, j + 1) - u(i, j - 1)) / (2 * dy) * (v(i + 1, j) - v(i - 1, j)) / (2 * dx)) - ((v(i, j + 1) - v(i, j - 1)) / (2 * dy)) * ((v(i, j + 1) - v(i, j - 1)) / (2 * dy))));                (1)

p(i, j) = ((pn(i+1, j) + pn(i-1, j)) * dy^2 + (pn(i, j+1) + pn(i, j-1)) * dx^2) ./ (2 * (dx^2 + dy^2)) - dx^2 * dy^2 / (2 * (dx^2 + dy^2)) * (rho * (1/dt * ((u(i+1, j) - u(i-1, j)) / (2 * dx) + (v(i, j+1) - v(i, j-1)) / (2 * dy)) - ((u(i+1, j) - u(i-1, j)) / (2 * dx)).^2 - 2 * ((u(i, j+1) - u(i, j-1)) / (2 * dy) .* (v(i+1, j) - v(i-1, j)) / (2 * dx)) - ((v(i, j+1) - v(i, j-1)) / (2 * dy)).^2));                (2)

double u_ij = un(i, j) - un(i, j) * dt / dx * (un(i, j) - un(i - 1, j)) - vn(i, j) * dt / dy * (un(i, j) - un(i, j - 1)) - dt / (2 * rho * dx) * (p(i + 1, j) - p(i - 1, j)) + nu * (dt / (dx * dx) * (un(i + 1, j) - 2 * un(i, j) + un(i - 1, j)) + dt / (dy * dy) * (un(i, j + 1) - 2 * un(i, j) + un(i, j - 1)));                (3)

u(i, j) = un(i, j) - un(i, j) * dt/dx .* (un(i, j) - un(i-1, j)) - vn(i, j) * dt/dy .* (un(i, j) - un(i, j-1)) - dt / (2 * rho * dx) * (p(i+1, j) - p(i-1, j)) + nu * (dt/dx^2 * (un(i+1, j) - 2 * un(i, j) + un(i-1, j)) + (dt/dy^2 * (un(i, j+1) - 2 * un(i, j) + un(i, j-1))));                (4)

double v_ij = vn(i, j) - un(i, j) * dt / dx * (vn(i, j) - vn(i - 1, j)) - vn(i, j) * dt / dy * (vn(i, j) - vn(i, j - 1)) - dt / (2 * rho * dy) * (p(i, j + 1) - p(i, j - 1)) + nu * (dt / (dx * dx) * (vn(i + 1, j) - 2 * vn(i, j) + vn(i - 1, j)) + dt / (dy * dy) * (vn(i, j + 1) - 2 * vn(i, j) + vn(i, j - 1)));                (5)

v(i, j) = vn(i, j) - un(i, j) * dt/dx .* (vn(i, j) - vn(i-1, j)) - vn(i, j) * dt/dy .* (vn(i, j) - vn(i, j-1)) - dt / (2 * rho * dy) * (p(i, j+1) - p(i, j-1)) + nu * (dt/dx^2 * (vn(i+1, j) - 2 * vn(i, j) + vn(i-1, j)) + (dt/dy^2 * (vn(i, j+1) - 2 * vn(i, j) + vn(i, j-1))));                (6)
    

     Of course constants need to be defined, such as grid spacing in space and time, density, kinematic viscosity. These equations have been validated, as you might have read already! Happy coding!

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.

References

[1] U Ghia, K.N Ghia, C.T Shin, "High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method", Journal of Computational Physics, Volume 48, Issue 3, 1982, Pages 387-411, ISSN 0021-9991, https://doi.org/10.1016/0021-9991(82)90058-4
[2] Khademinejad, Taha & Talebizadeh Sardari, Pouyan & Rahimzadeh, Hassan. (2015). Numerical Study of Unsteady Flow around a Square Cylinder in Compare with Circular Cylinder
[3] Ávila, Ítalo & Santos, Gabriel & Ribeiro Neto, Hélio & Neto, Aristeu. (2019). Physical Mathematical and Computational Modeling of the Two-Dimensional Flow Over a Heated Porous Square Cylinder. 10.26678/ABCM.COBEM2019.COB2019-0854
[4] Irisarri, Diego & Hauke, Guillermo. (2019). Stabilized virtual element methods for the unsteady incompressible Navier–Stokes equations. Calcolo. 56. 10.1007/s10092-019-0332-5.