Indeed, the 13th step now exists! This case is called the case of the backward facing step (BFS)! ⬜ This is an unofficial continuation to this. If I get it approved by Dr. Barba, then it will be official. The original series is in Python but I coded this in MATLAB without using many MATLAB specific functions so the code can be translated to other programing languages 🖧 quite easily.
The Code
%% clear and close
close all
clear
clc
%% define spatial and temporal grids
l_square = 1; % length of square
h = l_square/50; % grid spacing
dt = 0.1; % time step
L = 21; % cavity length
D = 2; % 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.0060444; % free stream velocity / inlet velocity / lid velocity
cfl = dt*Uinf/h; % % cfl number
travel = 4; % times the disturbance travels entire length of computational domain
TT = travel*L/Uinf; % total time
ns = TT/dt; % number of time steps
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)+pn(i, j+1)+pn(i, j-1))/4 ...
-h*rho/(8*dt)*(u(i+1, j)-u(i-1, j)+v(i, j+1)-v(i, j-1)); % pressure poisson
p(1, :) = p(2, :); % dp/dx = 0 at x = 0
p(Nx, :) = 0; % p = 0 at x = L
p(:, 1) = p(:, 2); % dp/dy = 0 at y = 0
p(:, Ny) = p(:, Ny-1); % dp/dy = 0 at y = D
p(round(1:1*Nx/L), round(1:1*Ny/D)) = 0; % step geometry
p(round(1*Nx/L), round(1:1*Ny/D)) = p(round(1*Nx/L)+1, round(1:1*Ny/D)); % dp/dx = 0 at x = 1 and y = 0 to 1
p(1:round(1*Nx/L), round(1*Ny/D)) = p(1:round(1*Nx/L), round(1*Ny/D)+1); % dp/dy = 0 at x = 0 to 1 and y = 1
p(1:round(1*Nx/L), 1) = p(1:round(1*Nx/L), 2); % dp/dy = 0 at x = 0 to 1 and y = 1
un = u;
vn = v;
u(i, j) = un(i, j)-dt/(2 * h)*(un(i, j).*(un(i+1, j)-un(i-1, j))+vn(i, j).*(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)+un(i-1, j)+un(i, j+1)+un(i, j-1)-4*un(i, j)); % x-momentum
u(round(1:1*Nx/L), round(1:1*Ny/D)) = 0; % step geometry
u(1, round(1:1*Ny/D)) = 0; % u = 0 at x = 0 and y = 0 to 1
u(1, round(1*Ny/D:2*Ny/D)) = Uinf; % u = Uinf at x = 0 and y = 1 to 2
u(round(1*Nx/L), round(1:1*Ny/D)) = 0; % u = 0 at x = 1 and y = 0 to 1
u(1:round(1*Nx/L), round(1*Ny/D)) = 0; % u = 0 at x = 0 to 1 and y = 1
u(1:round(1*Nx/L), 1) = 0; % u = 0 at x = 0 to 1 and y = 1
u(Nx, :) = u(Nx-1, :); % du/dx = 0 at x = L
u(:, 1) = 0; % u = 0 at y = 0
u(:, Ny) = 0; % u = 0 at y = D
v(i, j) = vn(i, j)-dt/(2*h)*(un(i, j).*(vn(i+1, j)-vn(i-1, j))+vn(i, j).*(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)+vn(i-1, j)+vn(i, j+1)+vn(i, j-1)-4*vn(i, j)); % y-momentum
v(round(1:1*Nx/L), round(1:1*Ny/D)) = 0; % step geometry
v(1, round(1:1*Ny/D)) = 0; % v = 0 at x = 0 and y = 0 to 1
v(1, round(1*Ny/D:2*Ny/D)) = 0; % v = Uinf at x = 0 and y = 1 to 2
v(round(1*Nx/L), round(1:1*Ny/D)) = 0; % v = 0 at x = 1 and y = 0 to 1
v(1:round(1*Nx/L), round(1*Ny/D)) = 0; % v = 0 at x = 0 to 1 and y = 1
v(1:round(1*Nx/L), 1) = 0; % v = 0 at x = 0 to 1 and y = 1
v(Nx, :) = v(Nx-1, :); % dv/dx = 0 at x = L
v(:, 1) = 0; % u = 0 at y = 0
v(:, Ny) = 0; % u = 0 at y = D
end
%% post-processing
velocity_magnitude = sqrt(u.^2 + v.^2); % velocity magnitude
u1 = u; % u-velocity for plotting with box
v1 = v; % v-velocity for plotting with box
p1 = p; % p-velocity for plotting with box
u1(1:round(1*Nx/L) , 1:round(1*Ny/D)) = NaN; % step geometry
v1(1:round(1*Nx/L) , 1:round(1*Ny/D)) = NaN; % step geometry
p1(1:round(1*Nx/L) , 1:round(1*Ny/D)) = NaN; % step geometry
velocity_magnitude1 = sqrt(u1.^2 + v1.^2); % velocity magnitude with box
%% Visualize velocity vectors and pressure contours
hold on
contourf(X, Y, u1', 64, 'LineColor', 'none'); % contour plot
set(gca, 'FontSize', 20)
% skip = 20;
% quiver(X(1:skip:end, 1:skip:end), Y(1:skip:end, 1:skip:end),...
% u1(1:skip:end, 1:skip:end)', v1(1:skip:end, 1:skip:end)', 1, 'k','LineWidth', 0.1); % Velocity vectors
% hh = streamslice(X, Y, u1', v1',2); % Streamlines
% set(hh, 'Color', 'k','LineWidth', 01);
colorbar; % Add color bar
colormap parula % Set color map
axis equal % Set true scale
xlim([0 L]); % Set axis limits
ylim([0 D]);
xticks([0 9 L]) % Set ticks
yticks([0 D]) % Set ticks
xt = [0 21]; % draw top wall
yt = [2 2];
xb = [1 21]; % draw bottom wall
yb = [0 0];
xbox = [0 1 1 0 0]; % draw box
ybox = [0 0 1 1 0];
plot(xbox, ybox, 'k', 'LineWidth', 2)
plot(xt, yt, 'k', 'LineWidth', 2)
plot(xb, yb, 'k', 'LineWidth', 2)
% clim([0 max(velocity_magnitude(:))]) % legend limits
% title('Velocity [m/s]');
xlabel('x [m]');
ylabel('y [m]');
The results from this code at Re 400 are presented in Fig. 1. The re-attachment length is ~8 m from the trailing-edge of the box, which is same as previously published results, from example in [1].
 |
| Fig. 1, post-processes results |
Thank you for reading! If you want to hire me as your next PhD student, please do reach out!
References
[1] Irisarri, D., Hauke, G. Stabilized virtual element methods for the unsteady incompressible Navier–Stokes equations. Calcolo 56, 38 (2019). https://doi.org/10.1007/s10092-019-0332-5
