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
