My code is about the uniform flow in the two parallel plates size 400x100 with the obstacle.
I hope that the fluid flow from the left hand side toward the right hand side and pass from the east boundary in right hand side.
But, when I run my code in matlab. The fluid flow from left to right.Until it near the east boundary,the fluid always reflect back to left side.
It is not flow pass from the east boundary.
What is wrong in my code?
Please help me. I’m new in LBM.
I have been trying to mend my code about 1 month. But now it is not better.
clc
clear all
m = 400;
n = 100;
uMax = 0.05;
Re = 100;
nu = uMax*n/Re;
omega = 1./(3*nu+1./2.);
maxT = 40000;
tPlot = 50;
%D2Q9 LATTICE CONSTANTS
w = [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9];
cx = [ 1, 0,-1, 0, 1,-1,-1, 1, 0];
cy = [ 0, 1, 0,-1, 1, 1,-1,-1, 0];
opp = [1, 2, 3, 4, 5, 6, 7, 8, 9];
col = [2:(n-1)];
in = 1;
out = m;
obst(:,1)=1;
obst(100:400,n)=1;
obst(1:100,61:100)=1;
bbRegion = find(obst);
[y,x] = meshgrid(1:n,1:m);
%INITIAL CONDITION
ux = zeros(m,n);
%L = n-2; y_phys = y-1.5;
%ux = 4*uMax/(L*L)*(y_phys.*L-y_phys.*y_phys);
uy = zeros(m,n);
rho =1;
for i =1:9
cu = 3*(cx(i)*ux+cy(i)*uy);
fln(i,:,:) = rho.*w(i).*...
(1+cu+1/2*(cu.*cu)-3/2*(ux.^2+uy.^2));
end
%MAIN LOOP
for cycle =1:maxT
%SUM
rho = sum(fln);
%ux = reshape((cx*reshape(fln,9,m*n)),1,m,n)./rho;
ux = abs(sum(fln([1,5,8],:,:)-fln([6,3,7],:,:))./rho)
uy = reshape((cy*reshape(fln,9,m*n)),1,m,n)./rho;
%COLLISION STEP
for i=1:9
cu = 3*(cx(i)*ux+cy(i)*uy);
fEq(i,:,:)= rho.*w(i).*...
(1+cu+1/2*(cu.*cu)-3/2*(ux.^2+uy.^2));
%fln(i,:,:)=fln(i,:,:)-omega.*(fln(i,:,:)-fEq(i,:,:));
fln(i,:,:)=fln(i,:,:).*(1-omega)+omega.*fEq(i,:,:);
end
%STREAMING STEP
for j = 1:n
%RIGHT TO LEFT
for i = m:-1:2
fln(1,i,j) = fln(1,i-1,j);
end
%LEFT TO RIGHT
for i = 1:m-1
fln(3,i,j) = fln(3,i+1,j);
end
end
for j = n:-1:2 %TOP TO BOTTOM
for i = 1:m
fln(4,i,j) = fln(4,i,j-1);
end
for i = m:-1:2
fln(8,i,j) = fln(8,i-1,j-1);
end
for i = 1:m-1
fln(7,i,j)=fln(7,i+1,j-1);
end
end
for j = 1:n-1 %BOTTOM TO TOP
for i = 1:m
fln(2,i,j) = fln(2,i,j+1);
end
for i = 1:m-1
fln(6,i,j) = fln(6,i+1,j+1);
end
for i = m:-1:2
fln(5,i,j) = fln(5,i-1,j+1);
end
end
%BOUNDARY CONDITION
%WEST B.C.
ux(:,in,2:60) = uMax;
uy(:,in,2:60) = 0;
rho(:,in,2:60) = 1./(1-ux(:,in,2:60)).*(...
sum(fln([9,2,4],in,2:60))+2*sum(fln([3,6,7],in,2:60)));
fln(1,in,2:60)=fln(3,in,2:60)+2/3*rho(:,in,2:60).*ux(:,in,2:60);
fln(5,in,2:60)=fln(7,in,2:60)-1/2*(fln(2,in,2:60)-fln(4,in,2:60))...
+1/2*rho(:,in,2:60).*uy(:,in,2:60)...
+1/6*rho(:,in,2:60).*ux(:,in,2:60);
fln(8,in,2:60)=fln(6,in,2:60)+1/2*(fln(2,in,2:60)-fln(4,in,2:60))...
-1/2*rho(:,in,2:60).*uy(:,in,2:60)...
+1/6*rho(:,in,2:60).*ux(:,in,2:60);
%NORTH B.C.
fln(4,:,1) = fln(2,:,1);
fln(7,:,1) = fln(5,:,1);
fln(8,:,1) = fln(6,:,1);
%SOUTH B.C.
fln(2,:,n) = fln(4,:,n);
fln(5,:,n) = fln(7,:,n);
fln(6,:,n) = fln(8,:,n);
%EAST B.C.
%TYPE1
%fln(1,out,col) = fln(1,out-1,col);
%fln(5,out,col) = fln(5,out-1,col);
%fln(8,out,col) = fln(8,out-1,col);
%fln(3,out,col) = fln(3,out-1,col);
%fln(6,out,col) = fln(6,out-1,col);
%fln(7,out,col) = fln(7,out-1,col);
%TYPE2
%rho(:,out,col) = 1;
%ux(:,out,col) = -1+1./(rho(:,out,col)).*(...
%sum(fln([9,2,4],out,col))+2*sum(fln([1,5,8],out,col)));
%uy(:,out,col) = 0;
%fln(3,out,col)=fln(1,out,col)-2/3*rho(:,out,col).*ux(:,out,col);
%fln(7,out,col)=fln(5,out,col)+1/2*(fln(2,out,col)-fln(4,out,col))...
%-1/2*rho(:,out,col).*uy(:,out,col)...
%-1/6*rho(:,out,col).*ux(:,out,col);
%fln(6,out,col)=fln(8,out,col)+1/2*(fln(4,out,col)-fln(2,out,col))...
%+1/2*rho(:,out,col).*uy(:,out,col)...
%-1/6*rho(:,out,col).*ux(:,out,col);
%TYPE3
fln(1,out,col) = 2*fln(1,out-1,col)-fln(1,out-2,col);
fln(5,out,col) = 2*fln(5,out-1,col)-fln(5,out-2,col);
fln(8,out,col) = 2*fln(8,out-1,col)-fln(8,out-2,col);
%type4
%ux(:,out,col)=ux(:,out-1,col);
%OBSTACLE
%TOP OF OBSTACLE
fln(2,1:100,61) = fln(4,1:100,61);
fln(5,1:100,61) = fln(7,1:100,61);
fln(6,1:100,61) = fln(8,1:100,61);
%EAST OF OBSTACLE
fln(1,100,61:100) = fln(3,100,61:100);
fln(5,100,61:100) = fln(7,100,61:100);
fln(8,100,61:100) = fln(6,100,61:100);
%WEST OF OBSTACLE
%fln(3,1,60:100) = fln(1,1,60:100);
%fln(7,1,60:100) = fln(5,1,60:100);
%fln(6,1,60:100) = fln(8,1,60:100);
%VISUALIZATION
if(mod(cycle,tPlot)==1)
u=reshape(sqrt(ux.^2+uy.^2),m,n);
u(bbRegion)=nan;
imagesc(u');
%contour(u');
axis equal off;drawnow
end
end