% Assignment 6 MMV042A 2003
% Boundary layer eqn. 
% Tangential flow along a flat plate
% Marching procedure in x-direction
% using TDMA in y-direction
clear all;
x=[];y=[];P=[];Q=[];U=[];V=[];T=[];
viscosity = 18.1e-6; % viscosity
density = 1.19 ;
kvisc = viscosity/density; % kinematic viscosity
Pr = 0.72; % Prandtl number
Umax = 5; % free stream velocity
maxres = 1e-5; % max residual
maxit = 3;
nx = input('Number of increment in x-dir= ')
ny = input('Number of increment in y-dir= ')
ni = nx +1;
nj = ny +1;
% Increment length in y-dir 0.05 - 0.2 mm
% Increment length in x-dir 0.1 - 1.0 mm
dx = input('Increment length in x-dir (mm)= ')
dy = input('Increment length in y-dir (mm)=  ')
dx = dx/1000;
dy = dy/1000;
x=0:dx:((ni-1)*dx)
y=0:dy:(nj-1)*dy
% Boundary and initial values
for i=1:ni
   for j=1:nj
      U(i,j)=Umax;
      if j==1; U(i,j)=0; end;
      V(i,j) = 0; 
      T(i,j) = 0;
      if j==1; T(i,j)=1.0; end;
   end;
end;  
counter = 0;
for i = 2:ni 
   counter = 0;
   sumres = 1;
   while ((sumres > maxres)&(counter < maxit) )
      P(1)=0;Q(1)=0;
      sumres=0;
%% momentuum eqn      
      for j = 2:nj-1
         a = U(i,j)/dx+2*kvisc/dy^2;
         b = kvisc/dy^2-V(i,j)/(2*dy);
         c = kvisc/dy^2+V(i,j)/(2*dy);
         d = U(i,j)*U(i-1,j)/dx;
         P(j)=b/(a-c*P(j-1));
         Q(j)=(c*Q(j-1)+d)/(a-c*P(j-1));
         sumres=sumres+abs(a*U(i,j)-(b*U(i,j+1)+c*U(i,j-1)+d)); 
      end;
      for j=nj-1:-1:2
         U(i,j) = P(j)*U(i,j+1)+Q(j);
      end;
      P(1)=0;Q(1)=0;
%%continuity eqn.
      for j = 2:nj-1
         a = 1/dy;
         b = 0;
         c = 1/dy;
         d = -((U(i,j)-U(i-1,j))+(U(i,j-1)-U(i-1,j-1)))/(2*dx);
         P(j)=b/(a-c*P(j-1));
         Q(j)=(c*Q(j-1)+d)/(a-c*P(j-1));
         sumres=sumres+abs(a*V(i,j)-(b*V(i,j+1)+c*V(i,j-1)+d)); 
      end;
      for j=nj-1:-1:2
         V(i,j) = P(j)*V(i,j+1)+Q(j);
      end;
      counter = counter +1; % counts number of iterations    
   end;
   TotalResidual = sumres
end;
for i = 2:ni % Temperature field
   P(1)=0;Q(1)=T(i,1);  % TDMA
%   for j = 2:nj-1
%      a =???
%      b =???
%      c =???
%      d =???
%      P(j)=b/(a-c*P(j-1));
%      Q(j)=(c*Q(j-1)+d)/(a-c*P(j-1));
%    end;
%    for j=nj-1:-1:2
%         T(i,j) = P(j)*T(i,j+1)+Q(j);
%    end;
end;
for i=2:ni
   % calc. Nusselt number and friction factor
end;   
subplot(2,1,1);pcolor(x,y,U');shading interp;title('U-velocity');colorbar;
% subplot(2,1,2);pcolor(x,y,T');shading interp;title('Temperature');colorbar;