N = 51;
L = 1; W = 1; % Dimension(m)
uo = 1; % Initial Velocity of top wall (m/s)
Psio = 0; Re = 100; % Initial Stream Function & Reynolds Number
dx = 1/N; dy = 1/N; % Size of cell
dt = 1e-4; % Time Step
%% Intialization of real cell %%
w = zeros(N+2); % Omega
u = zeros(N+2); v = zeros(N+2); % Velocities in x & y direction
Psi = zeros(N+2);
%% B.C. for Fictitious cell values at kth time level %%
for i = 1:N
% Left boundary
w(i+1,1) = w(i+1,2) + 2*((-8*(Psi(i+1,2) - Psio)/dx^2) - w(i+1,2));
u(i+1,1) = u(i+1,2) + 2*(0 - u(i+1,2));
v(i+1,1) = v(i+1,2) + 2*(0 - v(i+1,2));
Psi(i+1,1) = Psi(i+1,2) + 2*(Psio - Psi(i+1,2));
% Right boundary
w(i+1,N+2) = w(i+1,N+1) + 2*((-8*(Psi(i+1,N+1) - Psio)/dx^2) - w(i+1,N+1));
u(i+1,N+2) = u(i+1,N+1) + 2*(0 - u(i+1,N+1));
v(i+1,N+2) = v(i+1,N+1) + 2*(0 - v(i+1,N+1));
Psi(i+1,N+2) = Psi(i+1,N+1) + 2*(Psio - Psi(i+1,N+1));
% Top boundary
w(1,i+1) = w(2,i+1) + 2*((-8*(Psi(2,i+1) - Psio)/dy^2) - (4*uo/dy) - w(2,i+1));
u(1,i+1) = u(2,i+1) + 2*(uo - u(2,i+1));
v(1,i+1) = v(2,i+1) + 2*(0 - v(2,i+1));
Psi(1,i+1) = Psi(2,i+1) + 2*(Psio - Psi(2,i+1));
% Bottom boundary
w(N+2,i+1) = w(N+1,i+1) + 2*((-8*(Psi(N+1,i+1) - Psio)/dy^2) - w(N+1,i+1));
u(N+2,i+1) = u(N+1,i+1) + 2*(0 - u(N+1,i+1));
v(N+2,i+1) = v(N+1,i+1) + 2*(0 - v(N+1,i+1));
Psi(N+2,i+1) = Psi(N+1,i+1) + 2*(Psio - Psi(N+1,i+1));
end
%% Solution for omega at (k+1)th time level %%
for time = 1:1e5
% Real cell
for i = 2:N+1
for j = 2:N+1
% Convective Term
% Upwind scheme is used for Real Boundary Cell
if i==2||i==N+1||j==2||j==N+1
Fe = ((u(i,j+1) + u(i,j))/2)*dy; % East cell
if Fe>=0
we = w(i,j);
else
we = w(i,j+1);
end
Fw = ((u(i,j-1) + u(i,j))/2)*(-dy); % West cell
if Fw>=0
ww = w(i,j);
else
ww = w(i,j-1);
end
Fn = ((v(i-1,j) + v(i,j))/2)*dx; % North cell
if Fn>=0
wn = w(i,j);
else
wn=w(i-1,j);
end
Fs = ((v(i+1,j) + v(i,j))/2)*(-dx); % South cell
if Fs>=0
ws = w(i,j);
else
ws = w(i+1,j);
end
else
% Quick scheme is used for Real Inner Cell
Fe = ((u(i,j+1) + u(i,j))/2)*dy;
if Fe>=0
we = (6*w(i,j)/8) + (3*w(i,j+1)/8) - w(i,j-1)/8;
else
we = (6*w(i,j+1)/8) + (3*w(i,j)/8) - w(i,j+2)/8;
end
Fw = ((u(i,j-1) + u(i,j))/2)*(-dy);
if Fw>=0
ww = (6*w(i,j)/8) + (3*w(i,j-1)/8) - w(i,j+1)/8;
else
ww = (6*w(i,j-1)/8) + (3*w(i,j)/8) - w(i,j-2)/8;
end
Fn = ((v(i-1,j) + v(i,j))/2)*dx;
if Fn>=0
wn = (6*w(i,j)/8) + (3*w(i-1,j)/8) - w(i+1,j)/8;
else
wn = (6*w(i-1,j)/8) + (3*w(i,j)/8) - w(i-2,j)/8;
end
Fs = ((v(i+1,j) + v(i,j))/2)*(-dx);
if Fs>=0
ws = (6*w(i,j)/8) + (3*w(i+1,j)/8) - w(i-1,j)/8;
else
ws = (6*w(i+1,j)/8) + (3*w(i,j)/8) - w(i+2,j)/8;
end
end
% Diffusion Term
Fde = (-1/Re)*(((w(i,j+1) - w(i,j))/dx)*dy); % East cell
Fdw = (-1/Re)*(((w(i,j) - w(i,j-1))/dx)*-dy); % West cell
Fdn = (-1/Re)*(((w(i-1,j) - w(i,j))/dy)*dx); % North cell
Fds = (-1/Re)*(((w(i,j) - w(i+1,j))/dy)*-dx); % South cell
% Final Discretised Form
wnew(i,j) = w(i,j) + (dt/(dx*dy))*(-Fe*we - Fw*ww - Fn*wn - Fs*ws - Fde - Fdw -
Fdn - Fds);
end
end
% B.C. for omega at(k+1)th time level %
for i = 1:N
wnew(1,i+1) = wnew(2,i+1) + 2*((-8*(Psi(2,i+1) - Psio)/dy^2)-(4*uo/dy) -
wnew(2,i+1));
% Top
wnew(N+2,i+1) = wnew(N+1,i+1) + 2*((-8*(Psi(N+1,i+1) - Psio)/dy^2) -
wnew(N+1,i+1));
% Bottom
wnew(i+1,1) = wnew(i+1,2) + 2*((-8*(Psi(i+1,2) - Psio)/dx^2) - wnew(i+1,2));
% Left
wnew(i+1,N+2) = wnew(i+1,N+1) + 2*((-8*(Psi(i+1,N+1) - Psio)/dx^2) -
wnew(i+1,N+1));
% Right
end
%% Solution for Stream Function at (k+1)th time level %%
Psirms = 1; Psinew = Psi;
while Psirms>1e-6
% Real cell
for i = 2:N+1
for j = 2:N+1
Psinew(i,j) = (Psinew(i+1,j) + Psinew(i-1,j) + Psinew(i,j+1) + Psinew(i,j-1) +
(dx^2)*wnew(i,j))/4;
end
end
% B.C. for Stream Function
for i = 1:N
Psinew(1,i+1)= Psinew(2,i+1) + 2*(Psio - Psinew(2,i+1)); % Top
Psinew(N+2,i+1) = Psinew(N+1,i+1) + 2*(Psio - Psinew(N+1,i+1)); % Bottom
Psinew(i+1,1) = Psinew(i+1,2) + 2*(Psio - Psinew(i+1,2)); % Left
Psinew(i+1,N+2)= Psinew(i+1,N+1) + 2*(Psio - Psinew(i+1,N+1)); % Right
end
Psirms = sqrt(abs((sum(sum((Psi - Psinew).^2)))/(N+2)^2));
Psi = Psinew;
end
%% Finding Velocities at (k+1)th time level %%
% Real Cell
for i = 2:N+1
for j = 2:N+1
u(i,j) = (Psi(i-1,j) - Psi(i+1,j))/(2*dy);
v(i,j) = (Psi(i,j-1) - Psi(i,j+1))/(2*dx);
end
end
% B.C. for velocities
for i = 1:N
u(1,i+1) = u(2,i+1) + 2*(uo - u(2,i+1)); % Top
v(1,i+1) = v(2,i+1) + 2*(0 - v(2,i+1)); % Top
u(N+2,i+1) = u(N+1,i+1) + 2*(0 - u(N+1,i+1)); % Bottom
v(N+2,i+1) = v(N+1,i+1) + 2*(0 - v(N+1,i+1)); % Bottom
u(i+1,1) = u(i+1,2) + 2*(0 - u(i+1,2)); % Left
v(i+1,1) = v(i+1,2) + 2*(0 - v(i+1,2)); % Left
u(i+1,N+2) = u(i+1,N+1) + 2*(0 - u(i+1,N+1)); % Right
v(i+1,N+2) = v(i+1,N+1) + 2*(0 - v(i+1,N+1)); % Right
end
% Checking for steady state convergence
wrms = (sqrt((sum(sum((wnew - w).^2)))/((N+2)^2)))/dt;
w = wnew;
end
%% Plotting and Comparison of Velocities %%
% Data of U Ghia et al.
Y = [1 0.9766 0.9688 0.9609 0.9531 0.8516 0.7344 0.6172 0.5000 0.4531 0.2813 0.1719
0.1016 0.0703 0.0625 0.0547 0.0000];
U = [1 0.84123 0.78871 0.73722 0.68717 0.23151 0.00332 -0.13641 -0.20581 -0.21090 -
0.15662 -0.10150 -0.06434 -0.04775 -0.04192 -0.03717 -0.0000];
X = [1 0.9688 0.9609 0.9531 0.9453 0.9063 0.8594 0.8047 0.5000 0.2344 0.2266 0.1563
0.0938 0.0781 0.0703 0.0625 0.000];
V = [0.0000 -0.05906 -0.07391 -0.08864 -0.10313 -0.16914 -0.22445 -0.24533 0.05454
0.17527 0.17507 0.16077 0.12317 0.10890 0.10091 0.09233 0.0000];
subplot(1,2,1)
plot(u(2:N+1,(N+3)/2),1-1/(2*N):-1/N:1/(2*N)',U,Y,'*')
title('Comparision of u vs y')
xlabel('u','fontsize',10)
ylabel('y','fontsize',10)
subplot(1,2,2)
plot(1/(2*N):1/N:1-1/(2*N),v((N+3)/2,2:N+1),X,V,'*')
title('Comparision of x vs v')
xlabel('x','fontsize',10)
ylabel('v','fontsize',10)
% End
