%*****************Poiseuille flow using LBM***********************%%Based on procedures as explained in 'Lattice Gas Cellular Automata and Lattice  Boltzmann%Models'by Wolf Gladrow%code may have errors as it is my first experience with LBM.Readers are%suggested to check the code for errors.For feedback%vinuvargheseijk@gmail.comOMEGA=0.2;%Relaxation factorTAU=1/OMEGA;XMAX=20;%Mesh size in x-directionYMAX=20;%Mesh size in y-directionFORCING=(1.024)/YMAX^3;%Forcing term in Poiseuille equationKVISC = ((1/OMEGA)-0.5)/3;%lattice viscosityUX=(0.5*FORCING*YMAX*YMAX)/KVISC;%Velocity in x-direction cx=[1 0 -1 0 1 -1 -1 1 0];%components of lattice velocities in x-directioncy=[0 1 0 -1 1 1 -1 -1 0];%components of lattice velocities in y-directionjx=zeros(XMAX,YMAX);%creating an array for storing momentum in x-directionjy=zeros(XMAX,YMAX);%creating an array for storing momentum in y-directionrho=ones(XMAX,YMAX);%creating an array for storing densitiesf=zeros(XMAX,YMAX);%creating an array for storing distributionsfprop=zeros(XMAX,YMAX);%creating an array for storing propagating distributionsx=1:XMAX;y=1:YMAX;MAXIT=input('Number of iterations');%Initializataion of distributions        u=jx./rho;        v=jy./rho;      f(x,y,1) = (rho./9)*(1+3.*u+4.5.*u.*u-1.5.*(u.^2+v.^2));      f(x,y,2) = (rho./9)*(1+3.*v+4.5.*v.*v-1.5*(u.^2+v.^2));      f(x,y,3) = (rho./9)*(1-3.*u+4.5.*u.*u-1.5*(u.^2+v.^2));      f(x,y,4) = (rho./9)*(1-3.*v+4.5.*v.*v-1.5*(u.^2+v.^2));      f(x,y,5) = (rho./36)*(1+3.*(u+v)+4.5*(u+v).^2-1.5*(u.^2+v.^2));      f(x,y,6) = (rho./36)*(1+3.*(-u+v)+4.5*(-u+v).^2-1.5*(u.^2+v.^2));      f(x,y,7) = (rho./36)*(1-3.*(u+v)+4.5*(u+v).^2-1.5*(u.^2+v.^2));      f(x,y,8) = (rho./36)*(1+3.*(u-v)+4.5*(u-v).^2-1.5*(u.^2+v.^2));      f(x,y,9) = (4/9)*rho.*(1-1.5*(u.^2+v.^2));%Initialization completed%Assign first fprop to equilibrium distributions.      fprop=f;        feq=zeros(20,20,9);for iter=1:MAXIT      u=jx./rho;      v=jy./rho;                       feq(x,2:YMAX-1,1) = (rho(x,2:YMAX-1)./9).*(1+3.*u(x,2:YMAX-1)+4.5.*u(x,2:YMAX-1).*u(x,2:YMAX-1)-1.5.*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,2) = (rho(x,2:YMAX-1)./9).*(1+3.*v(x,2:YMAX-1)+4.5.*v(x,2:YMAX-1).*v(x,2:YMAX-1)-1.5.*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,3) = (rho(x,2:YMAX-1)./9).*(1-3.*u(x,2:YMAX-1)+4.5.*u(x,2:YMAX-1).*u(x,2:YMAX-1)-1.5.*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,4) = (rho(x,2:YMAX-1)./9).*(1-3.*v(x,2:YMAX-1)+4.5.*v(x,2:YMAX-1).*v(x,2:YMAX-1)-1.5.*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,5) = (rho(x,2:YMAX-1)./36).*(1+3.*(u(x,2:YMAX-1)+v(x,2:YMAX-1))+4.5*(u(x,2:YMAX-1)+v(x,2:YMAX-1)).^2-1.5*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,6) = (rho(x,2:YMAX-1)./36).*(1+3.*(-u(x,2:YMAX-1)+v(x,2:YMAX-1))+4.5*(-u(x,2:YMAX-1)+v(x,2:YMAX-1)).^2-1.5*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,7) = (rho(x,2:YMAX-1)./36).*(1-3.*(u(x,2:YMAX-1)+v(x,2:YMAX-1))+4.5*(u(x,2:YMAX-1)+v(x,2:YMAX-1)).^2-1.5*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,8) = (rho(x,2:YMAX-1)./36).*(1+3.*(u(x,2:YMAX-1)-v(x,2:YMAX-1))+4.5*(u(x,2:YMAX-1)-v(x,2:YMAX-1)).^2-1.5*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      feq(x,2:YMAX-1,9) = ((4/9)*rho(x,2:YMAX-1)).*(1-1.5*(u(x,2:YMAX-1).^2+v(x,2:YMAX-1).^2));      %propagating distributions fprop by applying kinetic equation.      fprop = (1.0-OMEGA).* f + OMEGA.* feq;      force =FORCING/6;      fprop(x,2:YMAX-1,1) = fprop(x,2:YMAX-1,1) + force;      fprop(x,2:YMAX-1,3) = fprop(x,2:YMAX-1,3) - force;      fprop(x,2:YMAX-1,5) = fprop(x,2:YMAX-1,5) + force;      fprop(x,2:YMAX-1,6) = fprop(x,2:YMAX-1,6) - force;      fprop(x,2:YMAX-1,7) = fprop(x,2:YMAX-1,7) - force;      fprop(x,2:YMAX-1,8) = fprop(x,2:YMAX-1,8) + force;       %Copying the boundary values.            fprop(x,1,:) = f(x,1,:);            fprop(x,YMAX,:) = f(x,YMAX,:);%        C6  C2  C5       ^ y%          \ | /          |%        C3-C0-C1         | %          / | \          | %        C7  C4  C8         -----> x       p=2:XMAX;       q=2:YMAX;       r=1:XMAX-1;       s=1:YMAX-1;       %Propagating C1              f(p,y,1) = fprop(p-1,y,1);      %Propagating C2       f(x,q,2) = fprop(x,q-1,2);      %Propagating C3       f(r,y,3) = fprop(r+1,y,3);      %Propagating C4          f(x,s,4) = fprop(x,s+1,4);       %Propagating C5        f(p,q,5) = fprop(p-1,q-1,5);      %Propagating C6       f(r,p,6) = fprop(r+1,p-1,6);      %Propagating C7      f(r,s,7) = fprop(r+1,s+1,7);      %Propagating C8      f(p,s,8) = fprop(p-1,s+1,8);            %Propagating C9      f(x,y,9) = fprop(x,y,9);            %Complete Bounce Back Boundary Conditions      %1.Implementing Periodic BC.f(x(1),y,1)=fprop(XMAX,y,1);f(XMAX,y,3)=fprop(x(1),y,3);f(x(1),2:YMAX,5)=fprop(XMAX,1:YMAX-1,5);f(XMAX,2:YMAX,6)=fprop(x(1),1:YMAX-1,6);f(XMAX,1:YMAX-1,7)=fprop(x(1),2:YMAX,7);f(x(1),1:YMAX-1,8)=fprop(XMAX,2:YMAX,8);   %2.Bounce Back Begins.temp=f(1:XMAX,1,2);f(1:XMAX,1,2)=f(1:XMAX,1,4);f(1:XMAX,1,4)=temp;temp=f(1:XMAX,YMAX,4);f(1:XMAX,YMAX,4)=f(1:XMAX,YMAX,2);f(1:XMAX,YMAX,4)=temp;temp=f(1:XMAX,1,5);f(1:XMAX,1,5)=f(1:XMAX,1,7);f(1:XMAX,1,7)=temp;temp=f(1:XMAX,YMAX,7);f(1:XMAX,YMAX,7)=f(1:XMAX,YMAX,5);f(1:XMAX,YMAX,5)=temp;temp=f(1:XMAX,1,6);f(1:XMAX,1,6)=f(1:XMAX,1,8);f(1:XMAX,1,8)=temp;temp=f(1:XMAX,YMAX,8);f(1:XMAX,YMAX,8)=f(1:XMAX,YMAX,6);f(1:XMAX,YMAX,6)=temp;rho(x,y) = f(x,y,1)+f(x,y,2)+f(x,y,3)+f(x,y,4)+f(x,y,5)+f(x,y,6)+f(x,y,7)+f(x,y,8)+f(x,y,9);jx(x,y)=f(x,y,1)-f(x,y,3)+f(x,y,5)-f(x,y,6)-f(x,y,7)+f(x,y,8);%Distributions multiplied by lattice velocities in x-directionsjy(x,y)=f(x,y,2)-f(x,y,4)+f(x,y,5)+f(x,y,6)-f(x,y,7)-f(x,y,8);%Distributions multiplied by lattice velocities in y-directionsuprofile(1:YMAX)=sum(jx(1:XMAX,1:YMAX)./rho(1:XMAX,1:YMAX))/XMAX;%Velocity profile in x-directionplot(uprofile)pause(.1)end