%--------------------------------------------------------------------%
%----------ME4233-COMPUTATIONAL METHODS IN FLUID MECHANICS-----------%
%-------------------------ASSIGNMENT 2-------------------------------%
%-------------------------'Assign2.m'--------------------------------%
% Date Created: 22/10/2016
% Last Updated: 22/10/2016
% PROGRAM OUTLINE
% This program solves a steady state, 2D incompressible viscous flow problem
% using both the artificial compressibility method and the projection
% formulation. The flow prolem consists of a square both with an inflow in
% the bottom left corner, outflow in the central right section, with a
% moving wall at the vertical boundary.
% Two cases are solved in this method:
% Case 1 - Re=100, Uslide=Umax, dx=dy=1/128
% Case 2 - Re=200, Uslide=Umax, dx=dy=1/120
% Case 1
%% Variable Declaration
clear all
dx=1/150;
dy=1/150;
Re=200;
s=1/Re;
Umax=1;
Uslide=Umax;
Uin=-0.25*Umax;
h=((Re/2)*(dx^2)*(dy^2)/(dx^2 + dy^2));
uk=zeros(151,151);
vk=zeros(151,151);
psi=zeros(151,151);
vor=zeros(151,151);
vor_prev=zeros(151,151);
vor_final=zeros(151,151);
omega=1.95;
X = linspace(0,1,151);
Y = linspace(0,1,151);
%% Wall boundary Conditions
% Uniform Outflow
for i=1:31
uk(i,1)=(1/3)*Uin;
psi(i,1)=(1/3)*Uin*(i-1)*dy;
end
psi(1,1)=0;
% Parabolic Inflow
for i=61:76
uk(i,151)=-400*Uin*((i-1)*dy-0.4)*((i-1)*dy-0.5);
psi(i,151)=-400*Uin*((1/3)*((i-1)*dy)^3 - 0.45*((i-1)*dy)^2 +0.2*(i-1)*dy - 88/3300);
end
% Wall Sliding Velocity
for i=1:151
uk(151,i)=Uslide;
end
for k=1:5000
% Stream Function Boundary Conditions
for i=1:151
psi(151,i)=Uslide*dy + psi(150,i);
end
for i=1:151
psi(i,1)=psi(i,2);
psi(i,151)=psi(i,150);
end
if k == 1
else
for i=2:31
psi(i,1)=uk(i-1,1)*dy + psi(i-1,1);
end
for i=61:76
psi(i,151)=uk(i-1,151)*dy + psi(i-1,151);
end
for i=1:151
psi(1,i)= psi(2,i);
end
end
% Stream Function
for i=2:150
for j=2:150
psi(i,j)=0.25*omega*(psi(i,j+1)+psi(i+1,j)+psi(i,j-1)+psi(i-1,j)+vor(i,j)*(dx^2))+(1-omega)*psi(i,j);
end
end
% Velocities
for i=2:150
for j=2:150
uk(i,j)=(psi(i+1,j)-psi(i-1,j))*(1/(2*dy));
vk(i,j)=(psi(i,j+1)-psi(i,j-1))*(1/(2*dx));
end
end
% Vorticity Boundary Conditions
for i=1:151
vor_final(151,i)=(2/(dy^2))*(psi(151,i)-psi(150,i)-Uslide*dy);
vor_final(1,i)=(2/(dy^2))*(psi(1,i)-psi(2,i));
end
for i=32:151
vor_final(i,1)=(2/(dy^2))*(psi(i,1)-psi(i,2));
end
for i=1:60
vor_final(i,151)=(2/(dy^2))*(psi(i,151)-psi(i,150));
end
for i=77:151
vor_final(i,151)=(2/(dy^2))*(psi(i,151)-psi(i,150));
end
vor_final(76,151)=(2/(dx^2))*(2*vor_final(76,151)-psi(76,150)-psi(75,150));
vor_final(61,151)=(2/(dx^2))*(2*vor_final(61,151)-psi(62,151)-psi(61,150));
vor_final(31,1)=(2/(dx^2))*(2*vor_final(31,1)-psi(31,2)-psi(30,1));
% Parabolic Inflow
for i=62:75
vor_final(i,151)=(1/(2*dx))*(uk(i-1,151)-uk(i+1,151));
end
% Outflow Boundary Condition
for i=2:30
vor_final(i,1)=vor_final(i,2);
end
% Vorticity
vor=vor_final;
for i=2:150
for j=2:150
vor_final(i,j)=(vor(i,j+1)*(s-(dy*uk(i,j))/2)+vor(i,j-1)*(s+(dy*uk(i,j))/2)+vor(i+1,j)*(s-(dy*vk(i,j))/2)+vor(i-1,j)*(s+(dy*vk(i,j))/2)+vor_prev(i,j)*((dx*dx)/(2*h)-2*s))/((dx*dx)/(2*h)+2*s);
end
end
vor_prev=vor;
k=k+1
end
contourf(x,y,uk)
%contourf(x,y,vk)
%contour(x,y,sqrt(uk^2 + vk^2))
%surf(u)
%surf(v)
%plot(x,uk(76,:))
%plot(x,uk(20,:))
%plot(x,uk(120,:))