% This Program was structured to solve Flow passes via Channel
% Copyright@2011-amsinfo.org
%E-mail:sja210@lehigh.edu
% It was written in Matlab.
%:::::::::::::::::::::::::::::::::::

% Computing the velocity profile and Temperature revolution
%:::::::::::::::::::::::::::::::::
clear all
clc
L=0.5; % The length of channel
H=0.02; % The hieght of channel
AR=round(L/H); % Aspect ratio
m=40;   % Number of lattice nodes in y-direction
n=AR*m; % Number of lattice Nodes in x-direction
cx=[0,1,0,-1,0,1,-1,-1,1]; % Velocity vector in x-direction
cy=[0,0,-1,0,1,1,1,-1,-1]; % Velocity vector in y-direction
w=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]; % Weight Factor
% Macroscopic Parameters
Ui=0.0002; % Inlet velocity
nu=6e-7; % Kinematic viscosity in m^2/sec at the mean temperature
Re=(Ui*H)./nu
rhoo=1;
dx=1;          % Distance between two adjacent lattice nodes in x-direction
dy=dx;         % Distance between two adjacent lattice nodes in y-direction
alpha=1.58e-7; %Thermal diffusivity in at the mean temperature in m^2/sec
Tw=1;          %The wall temperature 
% MESOSCOPIC PARAMETERS
ui=0.1;  % Inlet lattice velocity
% initialaization the distribution function
g(1:9,1:m,1:n)=0;
f(1:9,1:m,1:n)=0;
geq(1:9)=0;
feq(1:9)=0;
th(1:m,1:n)=0;
visco=(ui*m)/Re;       % Kinematic lattice viscosity  
omega=1./(3*visco+0.5);% Relaxation frequencey
pr=nu/alpha;           % Prandtle number
alphal=visco./pr;
%alphal=1;
% omegat=1/(3*alphal+0.5);
omegat=omega;
% Initilaization of macroscopic
rho(1:m,1:n)=rhoo;
u(1:m,1:n)=0;
v(1:m,1:n)=0;
uy(1:m,1:n)=0;
% Inlet Boundary condition
u(2:m-1,1)=ui;
uy(2:m-1,1)=Ui;
v(2:m,1)=0;
Hbottom=0;
Htop=m-1;
Hplot=Hbottom:Htop;
Y=Hplot/Htop;
dY=1/(m);
tol=0.001; 
x(1)=0;
Xtop=n-1;
Xbottom=0;
Xplot=Xbottom:Xtop;
X=Xplot/Xtop;
error=1;     %MAXIMUIM INITIAL ERROR
d=0;
while error>tol
    d=d+1;
    for j=1:m
        for i=1:n
            sum=0;
            for k=1:9
                sum=sum+g(k,j,i);
            end
            th(j,i)=sum;
        end 
    end
    th1=th;
    for i=2:n
        for j=2:m-1
            usum1=0;
            vsum1=0;
            for k=1:9
                
                usum1=usum1+f(k,j,i)*cx(k);
                vsum1=vsum1+f(k,j,i)*cy(k);  
            end 
            u(j,i)=usum1./rho(j,i);
            v(j,i)=vsum1./rho(j,i);
        end 
    end
    u1=u;
    % EQUILIBURIM VELOCITY DISTRIBUTION FUNCTION
    for i=1:n
    for j=1:m
    t1=u(j,i)*u(j,i)+v(j,i)*v(j,i);
    for k=1:9
      t2=u(j,i)*cx(k)+v(j,i)*cy(k);
      feq(k)=rho(j,i)*w(k)*(1+3*t2+4.5*t2*t2-1.5*t1);
      f(k,j,i)=omega*feq(k)+(1-omega)*f(k,j,i);
    end 
    end
    end 
    
    % STREAMING STEP 
     for j=1:m
     for i=n:-1:2
     f(2,j,i)=f(2,j,i-1);          % Right to Left 
     end 
     for i=1:n-1
     f(4,j,i)=f(4,j,i+1);         % Left ot Right 
     end
     end 
     for j=m:-1:2                 % Top to Bottom
     for i=1:n
     f(5,j,i)=f(5,j-1,i);
     end 
     for i=n:-1:2
     f(6,j,i)=f(6,j-1,i-1);       % Top to Bottom
     end 
     for i=1:n-1
     f(7,j,i)=f(7,j-1,i+1);
     end 
     end 
     for j=1:m-1
     for i=1:n
     f(3,j,i)=f(3,j+1,i);
     end 
     for i=1:n-1
     f(8,j,i)=f(8,j+1,i+1);   % Bottom to top
     end 
     for i=n:-1:2
     f(9,j,i)=f(9,j+1,i-1);
     end 
     end 
     % Flow in on the Left boundary condition
     for j=1:m
     rhow=(f(1,j,1)+f(5,j,1)+f(3,j,1)+2*(f(4,j,1)+f(7,j,1)+f(8,j,1)))./(1-ui);
     f(2,j,1)=f(4,j,1)+(2*rhow*ui)./3;
     f(6,j,1)=f(8,j,1)+(rhow*ui)./6;
     f(9,j,1)=f(7,j,1)+(rhow*ui)./6;
     end 
%      for i=1:n
%      f(5,1,i)=f(3,1,i);
%      f(6,1,i)=f(8,1,i);   % Lower Boundary Condition-Bounce back
%      f(7,1,i)=f(9,1,i);
%      end 
%      for i=1:n
%      f(3,m,i)=f(5,m,i);
%      f(8,m,i)=f(6,m,i);   % Upper Boundary Condition-Bounce back
%      f(9,m,i)=f(7,m,i);
%      end
   % Lower Boundary Condition-Bounce back
     for i=1:n
     rhow1=(f(1,1,i)+f(5,1,i)+f(2,1,i)+2*(f(3,1,i)+f(9,1,i)+f(8,1,i)))./(1-v(1,i));
     f(5,1,i)=f(3,1,i)+(2*rhow1*v(1,i))./3;
     f(6,1,i)=f(8,1,i)+0.5*(f(4,1,i)-f(2,1,i))+(rhow1*v(1,i))./6+(rhow1*u(1,i))./2;;
     f(7,1,i)=f(9,1,i)-0.5*(f(4,1,i)-f(2,1,i))+(rhow1*v(1,i))./6-(rhow1*u(1,i))./2;;
     end 
     % Upper Boundary Condition-Bounce back
     for i=1:n
     rhow=(f(1,m,i)+f(2,m,i)+f(4,m,i)+2*(f(5,m,i)+f(7,m,i)+f(6,m,i)))./(1+v(m,i));
     f(3,m,i)=f(5,m,i)-(2*rhow*v(m,i))./3;
     f(8,m,i)=f(6,m,i)+0.5*(f(2,m,i)-f(4,m,i))-(rhow*v(m,i))./6-(rhow*u(m,i))./2;
     f(9,m,i)=f(7,m,i)+0.5*(f(4,m,i)-f(2,m,i))-(rhow*v(m,i))./6+(rhow*u(m,i))./2;
     end
     % Account for open Boundry condition
     for j=2:m
     f(2,j,n)=2*f(2,j,n-1)-f(2,j,n-2);
     f(6,j,n)=2*f(6,j,n-1)-f(6,j,n-2); % Outlet Boundary
     f(9,j,n)=2*f(9,j,n-1)-f(9,j,n-2);
     end
     %TEMPERATURE PROFILE 
    % EQUILIBURIM TEMPERATURE DISTRIBUTION FUNCTION
    for i=1:n
    for j=1:m
    for k=1:9
      geq(k)=th(j,i)*w(k)*(1+3*(u(j,i)*cx(k)+v(j,i)*cy(k)));
      g(k,j,i)=omegat*geq(k)+(1-omegat)*g(k,j,i);
    end 
    end
    end 
     % STREAMING STEP 
     for j=1:m
     for i=n:-1:2
     g(2,j,i)=g(2,j,i-1);          % Right to Left 
     end 
     for i=1:n-1
     g(4,j,i)=g(4,j,i+1);         % Left ot Right 
     end
     end 
     for j=m:-1:2                 % Top to Bottom
     for i=1:n
     g(5,j,i)=g(5,j-1,i);
     end 
     for i=n:-1:2
     g(6,j,i)=g(6,j-1,i-1);       % Top to Bottom
     end 
     for i=1:n-1
     g(7,j,i)=g(7,j-1,i+1);
     end 
     end 
     for j=1:m-1
     for i=1:n
     g(3,j,i)=g(3,j+1,i);
     end 
     for i=1:n-1
     g(8,j,i)=g(8,j+1,i+1);   % Bottom to top
     end 
     for i=n:-1:2
     g(9,j,i)=g(9,j+1,i-1);
     end 
     end
     % LEFT BOUNDARY CONDITION THE INLET TEMPERATURE IS GIVEN (ZERO INLET)
     for j=1:m
     g(2,j,1)=-g(4,j,1);
     g(6,j,1)=-g(8,j,1);
     g(9,j,1)=-g(7,j,1);
     end 
    % RIGHT BOUNDARY CONDITION THE INLET TEMPERATURE IS GIVEN (OPEN OUTLET)
    for j=2:m
     g(7,j,n)=2*g(7,j,n-1)-g(7,j,n-2);
     g(4,j,n)=2*g(4,j,n-1)-g(4,j,n-2);
     g(8,j,n)=2*g(8,j,n-1)-g(8,j,n-2);
     g(5,j,n)=2*g(5,j,n-1)-g(5,j,n-2);
     g(1,j,n)=2*g(1,j,n-1)-g(1,j,n-2);
     g(2,j,n)=2*g(2,j,n-1)-g(2,j,n-2);
     g(3,j,n)=2*g(3,j,n-1)-g(3,j,n-2);
     g(6,j,n)=2*g(6,j,n-1)-g(6,j,n-2);
     g(9,j,n)=2*g(9,j,n-1)-g(9,j,n-2);   
    end
   % TOP BOUNDARY CONDITION THE INLET TEMPERATURE IS GIVEN (T=Tw)
   for i=1:n
     g(9,m,i)=Tw*(w(9)+w(7))-g(7,m,i);
     g(8,m,i)=Tw*(w(8)+w(6))-g(6,m,i);
     g(3,m,i)=Tw*(w(3)+w(5))-g(5,m,i);
     g(2,m,i)=Tw*(w(2)+w(4))-g(4,m,i);
   end
   % BOTTOM BOUNDARY CONDITION THE INLET TEMPERATURE IS GIVEN (T=Tw)
   for i=1:n
     g(5,1,i)=Tw*(w(5)+w(3))-g(3,1,i);
     g(7,1,i)=Tw*(w(7)+w(9))-g(9,1,i);
     g(6,1,i)=Tw*(w(6)+w(8))-g(8,1,i);
   end 
% RESULTS-PROGRAM OUTPUT
% TEMPERATURE PROFILE 

% VELOCITY PROFILE 
for j=1:m
    for i=1:n
        ssum=0;
        for k=1:9
            ssum=ssum+f(k,j,i);
        end 
        rho(j,i)=ssum;
    end 
end 
for i=2:n
    for j=2:m-1
        usum=0;
        vsum=0;
        for k=1:9
            usum=usum+f(k,j,i)*cx(k);
            vsum=vsum+f(k,j,i)*cy(k);
        end 
        u(j,i)=usum./rho(j,i);
        v(j,i)=vsum./rho(j,i);
    end
end 
for j=1:m
    v(j,n)=0;
end 
sume=0;
     for i=2:n
          for j=2:m-1
           sume=sume+(u(j,i)-u1(j,i)).^2;
          end
     end
     error1=sqrt(sume/(n+1)*(m+1));
    
% MACROSCOPIC VELOCITY PROFILE 
      for i=2:n
        for j=2:m-1
            uy(j,i)=(u(j,i)*m*nu)./(H*visco);
        end 
      end
    for j=2:m-1
    for i=2:n-1
        ssumt=0;
        for k=1:9
            ssumt=ssumt+g(k,j,i);
        end 
        th(j,i)=ssumt;
    end 
    end 
    sumte=0;
     for i=2:n
          for j=2:m-1
           sumte=sumte+(th(j,i)-th1(j,i)).^2;
          end
     end
     error2=sqrt(sumte/(n+1)*(m+1));
     error=error1+error2
end % THE END OF THE MAIN PROGRAM
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
figure(1)
plot(Y,th(:,100));
hold on 
plot(Y,th(:,300),'*');
plot(Y,th(:,1000),'rs:');
hold off
legend ('x/H=100/50','x/H=300/50','x/H=1000/50')
xlabel('Y');ylabel('(T-Tin)/(Tw-Tin)'),title('TEMPERATURE PROFILES AT DIFFERENT CROSS-SECTION AREAS')
figure(2)
plot(Y,uy(:,100)/Ui);
hold on
plot (Y,uy(:,300)./Ui,'ok');
plot (Y,uy(:,1000)./Ui,'b:');
hold off
legend ('x/H=100/50','x/H=300/50','x/H=1000/50')
xlabel('Y');ylabel('uy/Ui'),title('VELOCITY PROFILES AT DIFFERENT CROSS-SECTION AREAS')