Lid Driven Cavity Code using SIMPLE in MATLAB
 

Hello, I am trying to code for Lid Driven Cavity problem using finite volume discretization as mentioned in Versteeg. The result is very near to the exact solution but I can find some discrepancies. I've only attached functions for x, y momentum and pressure correction. The other functions are only for post processing.

I am using central differencing throughout. Upwinding does not improve my solution.

Although the results might look good but as I decrease the grid spacing, the result deteriorates.

Thanks in advance.

%% Main Program

tic
clc;
clear all;

format longG;

imax=10;
jmax=10;



Re=1;

dx=1/(imax-1);
dy=1/(jmax-1);

x=0:dx:1;
y=0:dy:1;

u=zeros(imax+1,jmax);
u_star=zeros(imax+1,jmax);
u_old=zeros(imax+1,jmax);

v_star=zeros(imax,jmax+1);
v=zeros(imax,jmax+1);
v_old=zeros(imax,jmax+1);

p=zeros(imax+1,jmax+1);
p_star=zeros(imax+1,jmax+1);

velocity=1;

u(1,2:jmax-1)=2*velocity-u(2,2:jmax-1);
u_star(1,2:jmax-1)=2*velocity-u_star(2,2:jmax-1);

alpha_p=0.1;
alpha=0.7;
max_iteration=10000;

for iteration=1:max_iteration

[u,u_P,A]=x_momentum(dx,dy,imax,jmax,u_star,v_star,p,Re,vel ocity,alpha);
[v,v_P,A]=y_momentum(dx,dy,imax,jmax,u_star,v_star,p,Re,vel ocity,alpha);
u_old=u;
v_old=v;
[p_star,u_star,v_star,A]=P_correction(dx,dy,imax,jmax,u,v,p,velocity,u_P,v _P,alpha_p);

p=p_star;

[div]=divergence(dx,dy,imax,jmax,u_star,v_star);
[velocity_residue]=converge(u_star,v_star,u_old,v_old);
%disp(iteration);

fprintf(' %d %d \n',iteration,velocity_residue);
if velocity_residue<1e-5
break;
elseif velocity_residue>1e2
break;
end



end

u=u_star;
v=v_star;

for i=1:imax
for j=1:jmax

up(i,j)=0.5*(u(i+1,j)+u(i,j));
vp(i,j)=0.5*(v(i,j+1)+v(i,j));
end
end

[psi]=stream_function(imax,jmax,dx,dy,u,v);
[omega]=vorticity(imax,jmax,dx,dy,u,v);
my_plot(x,y,up,vp,p,imax,jmax,omega,psi,velocity)
toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[u,u_P,A]=x_momentum(dx,dy,imax,jmax,u_star,v_star,p,Re,vel ocity,alpha)


n=1;
u=zeros(imax+1,jmax);
u_P=zeros(imax+1,jmax);
A=zeros((imax-1)*(jmax-2),(imax-1)*(jmax-2));
b=zeros((imax-1)*(jmax-2),1);


for i=2:imax
for j=2:jmax-1

F_e=0.5*(u_star(i,j)+u_star(i,j+1))*dy;
F_w=0.5*(u_star(i,j)+u_star(i,j-1))*dy;
F_n=0.5*(v_star(i-1,j)+v_star(i-1,j+1))*dx;
F_s=0.5*(v_star(i,j)+v_star(i,j+1))*dx;

De=dy/(Re*dx);
Dw=dy/(Re*dx);
Dn=dx/(Re*dy);
Ds=dx/(Re*dy);


a_E=0.5*F_e-De;
a_W=-0.5*F_w-Dw;
a_N=0.5*F_n-Dn;
a_S=-0.5*F_s-Ds;

a_P=0.5*(F_e-F_w+F_n-F_s)+De+Dw+Dn+Ds;
u_P(i,j)=a_P;



if i==2 && j==2
%disp('1')
A(n,n)=a_P -0.5*F_n + Dn;
A(n,n+1)=a_E;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i,j)-p(i,j+1))*dy - velocity*F_n + 2*Dn*velocity ;

elseif i==2 && j~=jmax-1
%disp('2')
A(n,n)=a_P -0.5*F_n + Dn;
A(n,n-1)=a_W;
A(n,n+1)=a_E;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i,j)-p(i,j+1))*dy - velocity*F_n + 2*Dn*velocity ;

elseif i==2 && j==jmax-1
%disp('3')
A(n,n)=a_P -0.5*F_n + Dn;
A(n,n-1)=a_W;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i,j)-p(i,j+1))*dy - velocity*F_n + 2*Dn*velocity;

elseif j==2 && i~=imax
%disp('4')
A(n,n)=a_P;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i,j)-p(i,j+1))*dy ;

elseif j==jmax-1 && i~=imax
%disp('5')
A(n,n)=a_P;
A(n,n-1)=a_W;
A(n,n-(imax-2))=a_N;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i,j)-p(i,j+1))*dy ;

elseif i==imax && j==2
%disp('6');
A(n,n)=a_P +0.5*F_s + Ds;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
b(n,1)=(p(i,j)-p(i,j+1))*dy ;

elseif i==imax && j~=jmax-1
%disp('7')
A(n,n)=a_P +0.5*F_s + Ds;
A(n,n-1)=a_W;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
b(n,1)=(p(i,j)-p(i,j+1))*dy ;

elseif i==imax && j==jmax-1
%disp('8')
A(n,n)=a_P +0.5*F_s + Ds;
A(n,n-1)=a_W;
A(n,n-(imax-2))=a_N;
b(n,1)=(p(i,j)-p(i,j+1))*dy ;

else
%disp('9')
A(n,n)=a_P;
A(n,n-1)=a_W;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i,j)-p(i,j+1))*dy ;




end

n=n+1;

end
end


u1=inv(A)*b;
n=1;

for i=2:imax
for j=2:jmax-1

u(i,j)=u1(n,1);
n=n+1;
end
end




u(2:imax,1)=0; %left
u(2:imax,jmax)=0; %right
u(1,:)=2*velocity-u(2,:); %top
u(imax+1,:)=-u(imax,:); %bottom



end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[v,v_P,A]=y_momentum(dx,dy,imax,jmax,u_star,v_star,p,Re,vel ocity,alpha)



v=zeros(imax,jmax+1);
n=1;
v_P=zeros(imax,jmax+1);
A=zeros(imax-1*jmax-2,imax-1*jmax-2);
b=zeros((imax-1)*(jmax-2),1);


for i=2:imax-1
for j=2:jmax

F_e=0.5*(u_star(i,j)+u_star(i+1,j))*dy;
F_w=0.5*(u_star(i,j-1)+u_star(i+1,j-1))*dy;
F_n=0.5*(v_star(i-1,j)+v_star(i,j))*dx;
F_s=0.5*(v_star(i,j)+v_star(i+1,j))*dx;

De=dy/(Re*dx);
Dw=dy/(Re*dx);
Dn=dx/(Re*dy);
Ds=dx/(Re*dy);


a_E=0.5*F_e-De;
a_W=-0.5*F_w-Dw;
a_N=0.5*F_n-Dn;
a_S=-0.5*F_s-Ds;
a_P=0.5*(F_n-F_s+F_e-F_w)+De+Dw+Dn+Ds;
v_P(i,j)=a_P;
if i==2 && j==2
%disp('1')
A(n,n)=a_P +0.5*F_w + Dw;
A(n,n+1)=a_E;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif i==2 && j~=jmax
%disp('2')
A(n,n)=a_P;
A(n,n-1)=a_W;
A(n,n+1)=a_E;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif i==2 && j==jmax
%disp('3')
A(n,n)=a_P -0.5*F_e + De;
A(n,n-1)=a_W;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif j==2 && i~=imax-1
%disp('4')
A(n,n)=a_P +0.5*F_w + Dw;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif j==jmax && i~=imax-1
%disp('5')
A(n,n)=a_P -0.5*F_e + De;
A(n,n-1)=a_W;
A(n,n-(imax-2))=a_N;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif i==imax-1 && j==2
%disp('6');
A(n,n)=a_P +0.5*F_w + Dw;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif i==imax-1 && j~=jmax
%disp('7')
A(n,n)=a_P;
A(n,n-1)=a_W;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

elseif i==imax-1 && j==jmax
%disp('8')
A(n,n)=a_P -0.5*F_e + De;
A(n,n-1)=a_W;
A(n,n-(imax-2))=a_N;
b(n,1)=(p(i+1,j)-p(i,j))*dx ;

else

%disp('9')
A(n,n)=a_P;
A(n,n-1)=a_W;
A(n,n+1)=a_E;
A(n,n-(imax-2))=a_N;
A(n,n+(imax-2))=a_S;
b(n,1)=(p(i+1,j)-p(i,j))*dx;


end





n=n+1;







end
end


v1=inv(A)*b;
n=1;
for i=2:imax-1
for j=2:jmax

v(i,j)=v1(n,1);
n=n+1;
end
end




v(1,2:jmax)=0; %top
v(imax,2:jmax)=0; %bottom
v(1:imax,1)=-v(1:imax,2); %left
v(1:imax,jmax+1)=-v(1:imax,jmax); %right


end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function[p_star,u_star,v_star,A]=P_correction(dx,dy,imax,jmax,u,v,p,velocity,u_P,v _P,alpha_p)

p_star=zeros(imax+1,jmax+1);
n=1;
A=zeros((imax-1)^2,(jmax-1)^2);
b=zeros((imax-1)^2,1);
for i=2:imax
for j=2:jmax




P_e=-(dy^2)/u_P(i,j);
P_w=-(dy^2)/u_P(i,j-1);
P_n=-(dx^2)/v_P(i-1,j);
P_s=-(dx^2)/v_P(i,j);


if i==2 && j==2
%disp('1')
A(n,n)=-(P_e+P_s);
A(n,n+1)=P_e;
A(n,n+(imax-1))=P_s;

elseif i==2 && j~=jmax
%disp('2')
A(n,n)=-(P_w+P_e+P_s);
A(n,n-1)=P_w;
A(n,n+1)=P_e;
A(n,n+(imax-1))=P_s;

elseif i==2 && j==jmax
%disp('3')
A(n,n)= -(P_w+P_s);
A(n,n-1)=P_w;
A(n,n+(imax-1))=P_s;

elseif j==2 && i~=imax
%disp('4')

A(n,n)=-(P_e+P_n+P_s);
A(n,n+1)=P_e;
A(n,n-(imax-1))=P_n;
A(n,n+(imax-1))=P_s;

elseif j==jmax && i~=imax
%disp('5')
A(n,n)=-(P_w+P_n+P_s) ;
A(n,n-1)=P_w;
A(n,n-(imax-1))=P_n;
A(n,n+(imax-1))=P_s;

elseif i==imax && j==2
%disp('6');
A(n,n)=-(P_e+P_n);
A(n,n+1)=P_e;
A(n,n-(imax-1))=P_n;

elseif i==imax && j~=jmax
%disp('7')
A(n,n)=-(P_e+P_n+P_w);
A(n,n+1)=P_e;
A(n,n-1)=P_w;
A(n,n-(imax-1))=P_n;

elseif i==imax && j==jmax
%disp('8')
A(n,n)=-(P_w+P_n);
A(n,n-1)=P_w;
A(n,n-(imax-1))=P_n;

else
%disp('9')
A(n,n)=-(P_e+P_n+P_s+P_n);
A(n,n-1)=P_w;
A(n,n+1)=P_e;
A(n,n-(imax-1))=P_n;
A(n,n+(imax-1))=P_s;

end


b(n,1)=-u(i,j)*dy + u(i,j-1)*dy - v(i-1,j)*dx + v(i,j)*dx;
n=n+1;

end
end


A(1,:)=0;
A(1,1)=1;
%b(1,1)=0;
p_prime=inv(A)*b;
n=1;

dp=zeros(imax+1,jmax+1);

for i=2:imax
for j=2:jmax

dp(i,j)=p_prime(n,1);
n=n+1;

end
end

p_star=p+0.1*dp;


p_star(imax+1,2:jmax)=p_star(imax,2:jmax); %bottom
p_star(1,2:jmax)=p_star(2,2:jmax); %top
p_star(2:imax,1)=p_star(2:imax,2); %left
p_star(2:imax,jmax+1)=p_star(2:imax,jmax); %right


du=zeros(imax+1,jmax);
for i=2:imax
for j=2:jmax-1

du(i,j)=(dp(i,j)-dp(i,j+1))*dy/u_P(i,j);


end
end


u_star=u+du;


u(2:imax,1)=0; %left
u(2:imax,jmax)=0; %right
u(1,:)=2*velocity-u(2,:); %top
u(imax+1,:)=-u(imax,:); %bottom

dv=zeros(imax,jmax+1);
for i=2:imax-1
for j=2:jmax


dv(i,j)=(dp(i+1,j)-dp(i,j))*dx/v_P(i,j);



end
end

v_star=v+dv;


v(1,2:jmax)=0; %top
v(imax,2:jmax)=0; %bottom
v(1:imax,1)=-v(1:imax,2); %left
v(1:imax,jmax+1)=-v(1:imax,jmax); %right


end