Hi jmazo,
I’ve looked at your code, I found a flaw it
To fix the problem yourself, you change the location of the boundary conditions imposed
1 - Put the boundary conditions at the wall between the COLLISION STEP and Streeming
2 - boundary conditions at the inlet and outlet should be placed after Streeming
%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%
% program %
%
% c
% c C7 C3 C6 ^ y
% c \ | / |
% c C4-C1-C2 |
% c / | \ |
% c C8 C5 C9 -----> x
% c
% c
% C
% c the lattice nodes are numbered as follows.
% c
% c
% c 1 ----------–…------------
% c | | | | | |
% c | | | | | |
% c 2 ----------–…------------
% c | | | | | |
% c : : : : : :
% c : : : : : :
% c | | | | | |
% c ly-1 ----------–…------------
% c | | | | | |
% c | | | | | |
% c ly ----------–…------------
% c 1 2 lx-1 lx
%% %%%%%%%%%%%%%%%
clc
clear all
close all
tic
% GENERAL FLOW CONSTANTS
ly=20;
lx=10*ly;
r=ones(1,lx);
pOut=1/3;
pIn=2*pOut;
Kno=.0194;
%Kni=.011;
fluid=zeros(ly,lx);
fluid((2:ly-1),:)=1;
bbRegion = find(fluid==0);
maxT = 40000; % total number of iterations
tPlot = 10; % cycles for graphical output
% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
[y,x] = meshgrid(1:ly,1:lx);
% INITIAL CONDITION: (rho=1, u=0) ==> fIn(i) = t(i)
fEq = zeros(9,ly,lx);
fOut = zeros(9,ly,lx);
ro=ones(1,lxly);
ro(:,:)=3pOut;
fIn = reshape( t’ * ro, 9, ly, lx);
rho=zeros(ly,lx);
ux =zeros(ly,lx);
uy =zeros(ly,lx);
err=1;
uu=ones(ly,lx);
cycle=0;
pBar=zeros(1,lx);
KnBar=zeros(1,lx);
TauBar=zeros(1,lx);
while cycle<=300000
cycle=cycle+1;
% MACROSCOPIC VARIABLES
for j=1:ly
for i=1:lx
%if fluid(j,i)==1
rho(j,i)=fIn(1,j,i)+fIn(2,j,i)+fIn(3,j,i)+fIn(4,j,i)…
+fIn(5,j,i)+fIn(6,j,i)+fIn(7,j,i)+fIn(8,j,i)+fIn(9,j,i);
end
end
sumrhoIn =(sum(rho(:,1 )))/ly;
sumrhoOut=(sum(rho(:,lx)))/ly;
rho(:,1) =rho(:,1 )*(3*pIn/sumrhoIn );
rho(:,lx)=rho(:,lx)*(3*pOut/sumrhoOut);
for j=1:ly
for i=1:lx
ux (j,i)=((fIn(2,j,i)+fIn(6,j,i)+fIn(9,j,i))...
-(fIn(4,j,i)+fIn(7,j,i)+fIn(8,j,i)))./rho(j,i);
uy (j,i)=((fIn(3,j,i)+fIn(6,j,i)+fIn(7,j,i))...
-(fIn(5,j,i)+fIn(8,j,i)+fIn(9,j,i)))./rho(j,i);
%else
rho([1,ly],i)=1;
ux ([1,ly],i)=0;
uy ([1,ly],i)=0;
%end
end
end
for i=1:lx
pBar(1,i)=sum(rho(1:ly,i))/(3*ly);
KnBar(1,i)=(Kno*pOut)/pBar(1,i);
%KnBar(1,i)=(Kni*pIn)/pBar(1,i);
TauBar(1,i)=.5+(sqrt(6/pi))*(KnBar(1,i)*ly);
%TauBar(1,i)=.5+(sqrt(6/pi))*(KnBar(1,i)*lx);
end
%%%%%%%%%%%%%%
% COLLISION STEP
%%%%%%%%%%%%%%
for i=1:lx
for j=1:ly
if fluid(j,i)==1
for k=1:9
cu = 3*(cx(k)ux(j,i)+cy(k)uy(j,i));
fEq(k,j,i)=rho(j,i)t(k)(1+cu+(.5cucu)-…
((3/2)((ux(j,i)^2)+(uy(j,i)^2))));
fOut(k,j,i)=fIn(k,j,i) - (1/TauBar(1,i))(fIn(k,j,i)-fEq(k,j,i));
end
end
end
end
A1=0.8183;
A2=0.6531;
for i=1:lx
r(i)=1;%1/(1+( (sqrt(pi/6))* ((1/(4KnBar(1,i)(lx^2)))+A1+((2*A2-(8/pi))*KnBar(1,i)) ) ));
end
%%%%%%%%%%%%%%%%%%%
% LOWER WALL
%%%%%%%%%%%%%%%%%%%
for i=1:lx
fOut(3,ly,i)=fIn(5,ly,i);
fOut(6,ly,i)=r(i)*fIn(8,ly,i)+(1-r(i))*fIn(9,ly,i);
fOut(7,ly,i)=r(i)*fIn(9,ly,i)+(1-r(i))*fIn(8,ly,i);
end
%%%%%%%%%%%%%%%%%%%
% UPPER WALL
%%%%%%%%%%%%%%%%%%%
for i=1:lx
fOut(5,1,i)=fIn(3,1,i);
fOut(8,1,i)=r(i)*fIn(6,1,i)+(1-r(i))*fIn(7,1,i);
fOut(9,1,i)=r(i)*fIn(7,1,i)+(1-r(i))*fIn(6,1,i);
end
% for i=1:9
% fOut(i,bbRegion) = fIn(opp(i),bbRegion);
% end
% for i=1:9
% fIn(i,:,
= circshift(fOut(i,:,:), [0,cy(i),cx(i)]);
% end
%%%%%%%%%%%%%%%%
%%%% Streeming
%%%%%%%%%%%%%%%%
for kk=1
%%% EAST
k=1+1;
for j=1:ly
for i=lx:-1:2
fIn(k,j,i)=fOut(k,j,i-1);
end
end
%%% NORTH
k=2+1;
for j=1:ly-1
for i=1:lx
fIn(k,j,i)=fOut(k,j+1,i);
end
end
%%% WEST
k=3+1;
for j=1:ly
for i=1:lx-1
fIn(k,j,i)=fOut(k,j,i+1);
end
end
%%% SOUTH
k=4+1;
for j=ly:-1:2
for i=1:lx
fIn(k,j,i)=fOut(k,j-1,i);
end
end
%%% NORTH-EAST
k=5+1;
for j=1:ly-1
for i=lx:-1:2
fIn(k,j,i)=fOut(k,j+1,i-1);
end
end
%%% NORTH-WEST
k=6+1;
for j=1:ly-1
for i=1:lx-1
fIn(k,j,i)=fOut(k,j+1,i+1);
end
end
%%% SOUTH-WEST
k=7+1;
for j=2:ly
for i=1:lx-1
fIn(k,j,i)=fOut(k,j-1,i+1);
end
end
%%% SOUTH-EAST
k=8+1;
for j=2:ly
for i=2:lx
fIn(k,j,i)=fOut(k,j-1,i-1);
end
end
%%% CENTER
k=0+1;
for j=1:ly
for i=1:lx
fIn(k,j,i)=fOut(k,j,i);
end
end
end
%%%%%%%%%%%%%%%%%%%
% INLET BOUNDARY0
% &&&&&&&&
% OUTLET BOUNDARY
%%%%%%%%%%%%%%%%%%%
for j=2:ly-1
i=1;
%%%%%%%%%%
% Zou He
%%%%%%%%%%
RhoIn=3*pIn;
Uin=1-(sum(fIn([3,5,1],j,i))+2*sum(fIn([7,4,8],j,i)))/RhoIn;
fIn(2,j,i) = fIn(4,j,i) + 2/3*RhoIn*Uin;
fIn(9,j,i) = fIn(7,j,1) + 1/2.*(fIn(3,j,1)-fIn(5,j,1)) + 1/6*RhoIn*Uin;
fIn(6,j,i) = fIn(8,j,1) - 1/2.*(fIn(3,j,1)-fIn(5,j,1)) + 1/6*RhoIn*Uin;
end
% % Bottom left (inlet)
x=1; y=ly;
ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pIn;
fIn(2,y,x)=fIn(4,y,x);
fIn(3,y,x)=fIn(5,y,x);
fIn(6,y,x)=fIn(8,y,x);
fIn(7,y,x)=1/2(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(6,y,x)+fIn(8,y,x)));
fIn(9,y,x)=fIn(7,y,x);
% Top left (inlet)
x=1; y=1;
ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pIn;
fIn(2,y,x)=fIn(4,y,x);
fIn(5,y,x)=fIn(3,y,x);
fIn(9,y,x)=fIn(7,y,x);
fIn(8,y,x)=1/2(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(7,y,x)+fIn(9,y,x)));
fIn(6,y,x)=fIn(8,y,x);
for j=2:ly-1
i=lx;
%%%%%%%%%%
% Zou He
%%%%%%%%%%
% RhoOut=3pOut;
% Uout=-1+(sum(fIn([3,5,1],j,i))+2sum(fIn([2,6,9],j,i)))/RhoOut;
% fIn(4,j,i) = fIn(2,j,i)- 2/3RhoOutUout;
% fIn(8,j,i) = fIn(6,j,i)+ 1/2.(fIn(3,j,i)-fIn(5,j,i))- 1/6RhoOutUout;
% fIn(7,j,i) = fIn(9,j,i)- 1/2.(fIn(3,j,i)-fIn(5,j,i))- 1/6RhoOutUout;
%%%%%%%%%%%%%%
% Open Buondry
%%%%%%%%%%%%%%
fIn(4,j,i)= 2*fIn(4,j,i-1) - fIn(4,j,i-2);
fIn(8,j,i)= 2*fIn(8,j,i-1) - fIn(8,j,i-2);
fIn(7,j,i)= 2*fIn(7,j,i-1) - fIn(7,j,i-2);
end
%Top right (outlet)
x=lx; y=1;
ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pOut;
fIn(4,y,x)=fIn(2,y,x);
fIn(5,y,x)=fIn(3,y,x);
fIn(8,y,x)=fIn(6,y,x);
fIn(7,y,x)=1/2(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(6,y,x)+fIn(8,y,x)));
fIn(9,y,x)=fIn(7,y,x);
% Bottom right (outlet)
x=lx; y=ly;
ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pOut;
fIn(4,y,x)=fIn(2,y,x);
fIn(3,y,x)=fIn(5,y,x);
fIn(7,y,x)=fIn(9,y,x);
fIn(8,y,x)=1/2(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(7,y,x)+fIn(9,y,x)));
fIn(6,y,x)=fIn(8,y,x);
% VISUALIZATION
if (mod(cycle,tPlot)==0)
disp(cycle)
u =(sqrt(ux.^2+uy.^2));
%u(bbRegion) = nan;
subplot(3,2,1)
%figure(8)
uIn=sum(sqrt(ux(2:ly-1,1).^2+uy(2:ly-1,1).^2))/(ly-2);
imagesc(u(1:ly,1:lx)./uIn);
%colorbar
axis equal off; drawnow
yy=linspace (0,1,ly-2);
%figure(3);
subplot(3,2,2)
plot(yy, ((sqrt(ux(2:ly-1,lx).^2+uy(2:ly-1,lx).^2)))./ (max (sqrt(ux(2:ly-1,lx).^2+uy(2:ly-1,lx).^2))) );%( sum((sqrt(ux(2:ly-1,lx).^2+uy(2:ly-1,lx).^2)))/(ly-2) ) ); %
xlabel('y=0: y=L')
ylabel('U/umax @ x=L ')
xL=linspace (0,1,lx);
%figure(5);
subplot(3,2,3)
plot(xL, ((sqrt(ux(ly/2,1:lx).^2+uy(ly/2,1:lx).^2)))./( sum((sqrt(ux(2:ly-1,1).^2+uy(2:ly-1,1).^2)))/(ly-2) ) );%./max (sqrt(ux(a+1:3*a,lx).^2+uy(a+1:3*a,lx).^2))) )
xlabel('X=0: X=L')
ylabel('U @ Y/2 ')
Px=(rho(ly/2,1:lx))./3;
P1x=zeros(1,lx);
dp=zeros(1,lx);
for i=1:lx
P1x(1,i)=pIn+(((i-1)/(lx-1))*(pOut-pIn));
dp(1,i)=(Px(i)-P1x(i))/pOut;
end
xx=linspace (0,1,lx);
%figure(11);
subplot(3,2,4)
plot(xx,dp)
maxdp=max(dp);
if maxdp<=0
maxdp=1;
end
axis([0 1 0 (maxdp+.1*maxdp)])
xlabel('X=0: X=L')
ylabel(' dP ')
%figure(10);
subplot(3,2,5)
plot(xx,Px./pOut,xx,P1x./pOut,'r')
axis([0 1 pOut/pOut pIn/pOut])
% subplot(3,2,5)
% plot(xx,P1x./pOut,‘r’)
xlabel(‘X=0: X=L’)
ylabel(’ Px ')
% Pbarx=(sum(rho(2:ly-1,1:lx)))./(3*(ly-2));
% xx=linspace (0,1,lx);
% figure(11);
% plot(xx,Pbarx./pOut)
end
% Px=(rho(ly/2,1:lx))./3;
% Knx=zeros(1,lx);
% for i= 1:lx
% Knx(i)=(pOutKno)./Px(i);
% omega(i)=1./(.5+(((sqrt(6/pi))((L))).*Knx(i)));
% end
u =(sqrt(ux.^2+uy.^2));
err=sum(sum(abs(uu-u)));
uu=u;
end
I know that some people simply copy f(LX-1,