%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%A CFD CODE FOR Stepped CAVITY%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Vorticity Stream Function %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Develoed By Sasan Tavakoli %%%%
clc
clear
RE=50;
D=1; %length of Problem
dx=0.01;
dy=0.01;
dt=5e-5;
l1=0.8; %length of step
u_tang=1; %Speed of moving wall
H1=0.8 ; %Height of Step
H2=1 ; %Height of Domainb
NX=D/dx;
NY=H2/dy;
Nu=u_tang*D/RE; % Kinematic Viscosity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% X-Locations:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:NX+1
X(i)=0+(i-1)*dx;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Y-Locations:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:NY+1
Y(j)=0+(j-1)*dy;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initial Conditions:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:NX
for j=2:NY
PSI(i,j)=0;
OMG(i,j)=0;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PSI Boundary Conditions:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%- in x-direction
%%Dawn wall befor step
for i=1:(l1/dx)+1
PSI(i,H1/dx+1)=0;
end
%%Dawn wall after step
for i=(l1/dx)+1:NX+1
PSI(i,1)=0;
end
%%Moving wall after step
for i=1:NX+1
PSI(i,NY+1)=0;
end
%%%- in y-direction
%%Left Wall
for j=H1/dx+1:NY+1
PSI(1,j)=0;
end
%Right wall
for j=H1/dx+1:NY+1
PSI(NX+1,j)=0;
end
max_err_OMG=1000; %Initial Error
itr=1;
time=0;
while max_err_OMG>= 1e-5 %begining of while
%%%%%%%%%%%%%%
% Boundary Conditions for OMG
%%%%%%%%%%%%%%
%%%%Dawn Wall Before Step
for i=1:(l1/dx)+1
OMG(i,H1/dx+1)=(2/dy^2)*(PSI(i,H1/dx+1)-PSI(i,H1/dx+2));
end
%%%%Dawn Wall After Step
for i=(l1/dx)+2:NX+1
OMG(i,1)=(2/dy^2)*(PSI(i,1)-PSI(i,2));
end
%%%%Moving Wall
for i=1:NX+1
OMG(i,NY+1)=(2/dy^2)*(PSI(i,NY+1)-PSI(i,NY))-2*u_tang/dy;
end
%%%%Left Wall
for j=H1/dx+1:NY+1
OMG(1,j)=(2/dx^2)*(PSI(1,j)-PSI(2,j));
end
%%%%Right Wall
for j=1:NY+1
OMG(NX+1,j)=(2/dx^2)*(PSI(NX+1,j)-PSI(NX,j));
end
%%%%%%%%%%%%%%%%%%%%%%%%%
%%Computational Procedure
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Block 1
for i=2:(l1/dx)+1
for j=H1/dx+2:NY
OMG_LI=OMG(i,j);
OMG(i,j)=OMG(i,j) + dt*Nu*( (OMG(i+1,j)-2*OMG(i,j)+OMG(i-1,j))/dx^2 + (OMG(i,j+1)-2*OMG(i,j)+OMG(i,j-1))/dy^2 ) + dt*(PSI(i+1,j)-PSI(i-1,j))/(2*dx)*(OMG(i,j+1)-OMG(i,j-1))/(2*dy) - dt*(PSI(i,j+1)-PSI(i,j-1))/(2*dy)*(OMG(i+1,j)-OMG(i-1,j))/(2*dx);
ERR_OMG(i,j)=abs(OMG(i,j)-OMG_LI);
end
end
%%%Block 2
for i=(l1/dx+2):NX
for j=2:NY
OMG_LI=OMG(i,j);
OMG(i,j)=OMG(i,j) + dt*Nu*( (OMG(i+1,j)-2*OMG(i,j)+OMG(i-1,j))/dx^2 + (OMG(i,j+1)-2*OMG(i,j)+OMG(i,j-1))/dy^2 ) + dt*(PSI(i+1,j)-PSI(i-1,j))/(2*dx)*(OMG(i,j+1)-OMG(i,j-1))/(2*dy) - dt*(PSI(i,j+1)-PSI(i,j-1))/(2*dy)*(OMG(i+1,j)-OMG(i-1,j))/(2*dx);
ERR_OMG(i,j)=abs(OMG(i,j)-OMG_LI);
end
end
%%%%%%finding the maximum error
max_err_OMG=ERR_OMG(1,1);
for i=2:NX
for j=2:NY
if ERR_OMG(i,j)>max_err_OMG
max_err_OMG=ERR_OMG(i,j);
end
end
end
%%%%%%%%%%%%%%%%%%%
%%Solution for psi
%%%%%%%%%%%%%%%%%%%
max_err_PSI=1000; % Initial Error for PSI
itpsi=1;
while max_err_PSI>1e-5
%Block 1
for i=2:(l1/dx)+1
for j=H1/dx+2:NY
PSI_LI=PSI(i,j);
PSI(i,j)=(1/4)*( PSI(i+1,j) + PSI(i-1,j) + PSI(i,j-1) + PSI (i,j+1) + OMG(i,j)*dx^2 );
ERR_PSI(i,j)=abs(PSI(i,j)-PSI_LI);
end
end
%Block 2
for i=(l1/dx+2):NX
for j=2:NY
PSI_LI=PSI(i,j);
PSI(i,j)=(1/4)*( PSI(i+1,j) + PSI(i-1,j) + PSI(i,j-1) + PSI (i,j+1) + OMG(i,j)*dx^2 );
ERR_PSI(i,j)=abs(PSI(i,j)-PSI_LI);
end
end
max_err_PSI=ERR_PSI(1,1);
for i=2:NX
for j=2:NY
if ERR_PSI(i,j)>max_err_PSI
max_err_PSI=ERR_PSI(i,j);
end
end
end
itpsi=itpsi+1 ;
end
time=time+dt;
itr=itr+1
% if itr>10000
% break
% end
end %finishing of while
%%%%%%%%%%%%%%%
%Velocity
%%%%%%%%%%%%%%%
%%%% u = d PSI / dx
%Block 1
for i=2:(l1/dx)+1
for j=H1/dx+2:NY
u(i,j)=(PSI(i,j+1)-PSI(i,j))/dy;
end
end
%%%% v = d PSI / dx
for i=2:(l1/dx)+1
for j=H1/dx+2:NY
v(i,j)=-(PSI(i+1,j)-PSI(i,j))/dx;
end
end
%%%Block 2
for i=(l1/dx+2):NX
for j=2:NY
u(i,j)=(PSI(i,j+1)-PSI(i,j))/dy;
end
end
%%%% v = d PSI / dx
for i=2:(l1/dx)+1
for j=2:NY
v(i,j)=-(PSI(i+1,j)-PSI(i,j))/dx;
end
end
%%%% u and v at Boundaries
for i=1:NX+1
u(i,1)=0;
v(i,1)=0;
u(i,NY+1)=u_tang;
v(i,NY+1)=0;
end
for j=1:NY+1
u(1,j)=0;
v(1,j)=0;
u(NX+1,j)=0;
v(NX+1,j)=0;
end
for i=1:NX+1
for j=1:NY+1
xx(i,j)=X(i);
yy(i,j)=Y(j);
end
end
figure
contour(xx,yy,PSI,30)
figure
contour(xx,yy,OMG,30)
for i=1:NX+1
for j=1:NY+1
xx(i,j)=X(i);
yy(i,j)=Y(j);
end
end
figure
contour(xx,yy,PSI,30)
daspect([1 1 1])
figure
contour(xx,yy,OMG,30)
daspect([1 1 1])