% Define problem parameters
L = 1; % size of domain
Ul = 1; % velocity of upper lid
Re = 1000; % Reynolds number
nu = Ul*L/Re; % viscosity
nx = 40; % number of cells in x direction
ny = 40; % number of cells in y direction
dx = L/nx; % cell size in x direction
dy = L/ny; % cell size in y direction
dt = 0.001; % time step
nt = 100; % number of time steps
% Create grid
x = linspace(0, L, nx+1);
y = linspace(0, L, ny+1);
[X, Y] = meshgrid(x, y);
% Set initial conditions
u = zeros(nx+1, ny+2);
v = zeros(nx+2, ny+1);
p = zeros(nx+2, ny+2);
% Define coefficients for finite volume discretization
ae = repmat(1/nu*dt/dx^2, nx+2, ny+2);
aw = repmat(1/nu*dt/dx^2, nx+2, ny+2);
an = repmat(1/nu*dt/dy^2, nx+2, ny+2);
as = repmat(1/nu*dt/dy^2, nx+2, ny+2);
ap = ae+aw+an+as+repmat(1/dt, nx+2, ny+2);
b = zeros(nx+2, ny+2);
% Main loop for time integration
for n = 1:nt
% Calculate intermediate velocities
ue = u(2:end-1, 2:end-1) + dt*(-(p(2:end-1, 3:end)-p(2:end-1, 2:end-1))/dy + nu*(u(3:end, 2:end-1)-2*u(2:end-1, 2:end-1)+u(1:end-2, 2:end-1))/dx^2);
uw = u(1:end-2, 2:end-1) + dt*(-(p(2:end-1, 2:end-1)-p(2:end-1, 1:end-2))/dy + nu*(u(2:end-1, 2:end-1)-2*u(1:end-2, 2:end-1)+u(1:end-3, 2:end-1))/dx^2);
vn = v(2:end-1, 2:end-1) + dt*(-(p(3:end, 2:end-1)-p(2:end-1, 2:end-1))/dx + nu*(v(2:end-1, 3:end)-2*v(2:end-1, 2:end-1)+v(2:end-1, 1:end-2))/dy^2);
vs = v(2:end-1, 1:end-2) + dt*(-(p(2:end-1, 2:end-1)-p(1:end-2, 2:end-1))/dx + nu*(v(2:end-1, 2:end-1)-2*v(2:end-1, 1:end-2)+v(2:end-1, 1:end-3))/dy^2);
end
% Impose boundary conditions
% Velocity boundary conditions
u(1,:) = Ul; % Top wall
v(1,:) = 0;
u(end,:) = 0; % Bottom wall
v(end,:) = 0;
u(:,1) = 0; % Left wall
v(:,1) = 0;
u(:,end) = 0; % Right wall
v(:,end) = 0;
% Pressure boundary conditions
% Top wall (Neumann boundary condition)
for i = 2:Nx-1
Ap =zeros (1, i);
Ap(i,end) = 0; % diagonal element
Aw =zeros (1, i);
Aw(i,end) = 0; % west coefficient
Ae = zeros (1, i);
Ae(i,end) = 0; % east coefficient
An = zeros (1, i);
An(i,end) = 1/dy; % north coefficient
As = zeros (1, i);
As(i,end) = -1/dy; % south coefficient
b(i,end) = 0; % right-hand side
end
% Bottom wall (Neumann boundary condition)
for i = 2:Nx-1
Ap(i,1) = 0; % diagonal element
Aw(i,1) = 0; % west coefficient
Ae(i,1) = 0; % east coefficient
An(i,1) = -1/dy; % north coefficient
As(i,1) = 1/dy; % south coefficient
b(i,1) = 0; % right-hand side
end
% Left wall (Neumann boundary condition)
for j = 2:Ny-1
Ap(1,j) = 1; % diagonal element
Aw(1,j) = 0; % west coefficient
Ae(1,j) = -1/dx; % east coefficient
An(1,j) = 0; % north coefficient
As(1,j) = 0; % south coefficient
b(1,j) = 0; % right-hand side
end
% Right wall (Dirichlet boundary condition)
for j = 2:Ny-1
Ap(end,j) = 1; % diagonal element
Aw(end,j) = -1/dx; % west coefficient
Ae(end,j) = 0; % east coefficient
An(end,j) = 0; % north coefficient
As(end,j) = 0; % south coefficient
b(end,j) = 0; % right-hand side
end
