Basic Lid Driven Cavity problem solved using MATLAB

 

This is a basic code of computation. The method used here is not explained but I will try to upload some document later on. 

Formulation Method: Stream function - Vorticity Method

Discretization Method: Simple Second Order Central Difference Scheme

Iteration Method: Gauss Seidel (Point iteration) with no under/over relaxation

Convergence: RMS psi and omega <10 ^ -5

**N.B: For higher Reynolds number increase the grid density or use higher order upwind schemes/QUICK Scheme and for efficient iteration use line gauss seidel with little under relaxation/ADI Method.

CODE:

clear all;
clc;
format longG

m=101;n=101;
p=(m-2)*(n-2);
dx=1/(m-1);
dy=1/(n-1);
beta=dx/dy;
Re=500;
x=dx/2:dx:1-dx/2;  
y=0:dy:1;
psi=zeros(m,n); omega=zeros(m,n); u=zeros(m,n); v=zeros(m,n);
psi_prev=zeros(m,n); omega_prev=zeros(m,n);
error_psi=10; error_omega=10;
iteration=0;
%Initialisation of U
u(:,n)=1; u(1,1)=0; u(1,n)=0;
%Initialisation of omega
for j=1:n
    for i=1:m
        if (j == 1)
            omega(i,j)=(2/dy^2)*(psi(i,j)-psi(i,j+1));
        elseif (j == n)
            omega(i,j)=((2/dy^2)*(psi(i,j)-psi(i,j-1)))-((2/dy)*u(i,j));
        elseif (i == 1)
            omega(i,j)=(2/dx^2)*(psi(i,j)-psi(i+1,j));
        elseif (i == m)
            omega(i,j)=(2/dx^2)*(psi(i,j)-psi(i-1,j));
        else
            omega(i,j)=0;
                
        end
    end
end
%Gauss Seidel
% ||
while( error_psi>1e-5 || error_omega>1e-5)
    psi_prev=psi;
    omega_prev=omega;
    %solving stream function
for j=2:n-1
    for i=2:m-1
        psi(i,j)=(0.5/(1+beta^2))*(psi(i+1,j)+psi(i-1,j)+(beta^2*(psi(i,j+1)+psi(i,j-1)))+(dx^2*omega(i,j)));
    end
end
%solving vorticity
for j=2:n-1
    for i=2:m-1
        omega(i,j)=(0.5/(1+beta^2))*(((1-((psi(i,j+1)-psi(i,j-1))*(beta*Re/4)))*omega(i+1,j))+((1+((psi(i,j+1)-psi(i,j-1))*(beta*Re/4)))*omega(i-1,j))+ ((1+((psi(i+1,j)-psi(i-1,j))*(Re/(4*beta))))*((beta^2)*omega(i,j+1)))+((1-((psi(i+1,j)-psi(i-1,j))*(Re/(4*beta))))*((beta^2)*omega(i,j-1))));
            
    end
end
      
%Updation of boundary omega
for j=1:n
    for i=1:m
        if (j == 1)
            omega(i,j)=(2/dy^2)*(psi(i,j)-psi(i,j+1));
        elseif (j == n)
            omega(i,j)=((2/dy^2)*(psi(i,j)-psi(i,j-1)))-((2/dy)*u(i,j));
        elseif (i == 1)
            omega(i,j)=(2/dx^2)*(psi(i,j)-psi(i+1,j));
        elseif (i == m)
            omega(i,j)=(2/dx^2)*(psi(i,j)-psi(i-1,j));                 
        end
    end
end
error_psi=0;
error_omega=0;
for j=2:n-1
    for i=2:m-1
        error_psi=error_psi+(psi(i,j)-psi_prev(i,j))^2;
        error_omega=error_omega+(omega(i,j)-omega_prev(i,j))^2;
    end
end
error_psi=sqrt(error_psi/p);
error_omega=sqrt(error_omega/p);

fprintf("Iteration=%d ",iteration);
fprintf(" error PSI=%f  error Omega=%f\n",error_psi,error_omega);
iteration=iteration+1;
end

%update velocities
for j=2:n-1
    for i=2:m-1
        u(i,j)=(0.5/dy)*(psi(i,j+1)-psi(i,j-1));
        v(i,j)=(-0.5/dy)*(psi(i+1,j)-psi(i-1,j));
    end
end

figure
contour(x,y,psi(2:m,:)',50);
colormap jet
colorbar;
figure
contour(x,y,omega(2:m,:)',50);
colormap jet
colorbar;
figure
contour(x,y,u(2:m,:)',50);
colormap jet
colorbar;

 

OUTPUT 

PSI Contour:

 

 

 Omega Contour

 

 

 u-vel Contour