Hi all,
I’m working on Poiseuille flow since a week and there are thing I don’t understand.
I read almost all the relative topics on the forum but there are still some dark points that I try to figure out.
I use incompressible BGK model and the flow is driven by a pressure drop (I use Zou/He BC). A half way Bounce back is applied on other walls.
The two issues that I would like understand are :
- why I obtain a vertical non-zero velocity ? I Zou and He paper they say that the error on the y-velocity is zero. (10^-15) and I think I have done things like they are descibed in the paper.
- why the code can not reach the imposed peak velocity (the pressure drop is construct in this way). I suppose the wo questions are linked.
I post the code I use, I started from the cylinder.m writen by J.Latt.
clear all
close all
format short g
% GENERAL FLOW CONSTANTS
lx = 45; % number of cells in x-direction
ly = 15; % number of cells in y-direction
Re = 10; % Reynolds number
maxT = 10000; % total number of iterations
tPlot = 100; % cycles
% 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];
%col = [2:(ly-1)];
col = [1:ly];
in = 1; % position of inlet
out = lx; % position of outlet
fIn=zeros(9,lx,ly);
fOut=zeros(9,lx,ly);
[y,x] = meshgrid(1:ly,1:lx); % get coordinate of matrix indices
obst(:,[1,ly]) = 1; % Location of top/bottom boundary
bbRegion = find(obst); % Boolean mask for bounce-back cells
% INITIAL CONDITION: u = 0
ux = zeros(lx,ly);
uy = zeros(lx,ly);
rho = 1;
% caracteristic length L=f(Bounce-Back)
% with HW BB, walls are in y=1-0.5 and y=ly+0.5
Length = ly;
omega = 1.5;
nu = (omega - 0.5)/3;
omega=1/omega;
uref_lb=nu*Re/Length;
rho0 = 1d0;
rho_out = rho0;
rho_in = 3*(8*abs(uref_lb)*nu*(lx-1)*rho0/Length^2)+rho_out;
for i=1:9
cu = 3*(cx(i)*ux + cy(i)*uy);
fIn(i,:,:) = t(i) .* ...
( rho + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) );
end
fid1 = fopen('hist.dat','w');
% MAIN LOOP (TIME CYCLES)
for cycle = 0:maxT
% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly);
uy = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly);
% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
% Inlet: Constant pressure
uy(:,in,col) = 0;
rho(:,in,col) = rho_in;
ux(:,in,col) = rho(:,in,col) - ( ...
sum(fIn([1,3,5],in,col)) + 2*sum(fIn([4,7,8],in,col)) );
% Outlet: Constant pressure
uy(:,out,col) = 0;
rho(:,out,col) = rho_out;
ux(:,out,col) = - rho(:,out,col) + ( ...
sum(fIn([1,3,5],out,col)) + 2*sum(fIn([2,6,9],out,col)) );
% MICROSCOPIC BOUNDARY CONDITIONS: INLET (Zou/He BC)
fIn(2,in,col) = fIn(4,in,col) + 2/3*ux(:,in,col);
fIn(6,in,col) = fIn(8,in,col) + 1/2*(fIn(5,in,col)-fIn(3,in,col)) ...
+ 1/2*uy(:,in,col) ...
+ 1/6*ux(:,in,col);
fIn(9,in,col) = fIn(7,in,col) + 1/2*(fIn(3,in,col)-fIn(5,in,col)) ...
- 1/2*uy(:,in,col) ...
+ 1/6*ux(:,in,col);
% MICROSCOPIC BOUNDARY CONDITIONS: OUTLET (Zou/He BC)
fIn(4,out,col) = fIn(2,out,col) - 2/3*ux(:,out,col);
fIn(8,out,col) = fIn(6,out,col) + 1/2*(fIn(3,out,col)-fIn(5,out,col)) ...
- 1/2*uy(:,out,col) ...
- 1/6*ux(:,out,col);
fIn(7,out,col) = fIn(9,out,col) + 1/2*(fIn(5,out,col)-fIn(3,out,col)) ...
+ 1/2*uy(:,out,col) ...
- 1/6*ux(:,out,col);
% COLLISION STEP
for i=1:9
cu = 3*(cx(i)*ux+cy(i)*uy);
fEq(i,:,:) = t(i) .* ...
( rho + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) );
fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:));
end
% OBSTACLE (BOUNCE-BACK)
% for i=1:9
% fOut(i,bbRegion) = fIn(opp(i),bbRegion);
% end
% STREAMING STEP
% for i=1:9
% fIn(i,:,:) = circshift(fOut(i,:,:), [0,cx(i),cy(i)]);
% end
for i=1:lx
for j=1:ly
for k=1:9
xnew = i + cx(k);
ynew = j + cy(k);
if (xnew>=1 && xnew <= lx) && (ynew>=1 && ynew <= ly)
fIn(k, xnew, ynew) = fOut(k, i, j);
end
end
end
end
% Bounce-Back
fIn(6,:,1) = fIn(opp(6),:,1);
fIn(3,:,1) = fIn(opp(3),:,1);
fIn(7,:,1) = fIn(opp(7),:,1);
fIn(8,:,ly) = fIn(opp(8),:,ly);
fIn(5,:,ly) = fIn(opp(5),:,ly);
fIn(9,:,ly) = fIn(opp(9),:,ly);
% VISUALIZATION
if (mod(cycle,tPlot)==0)
fprintf(fid1,'%12.8E %12.8E %12.8E %12.8E %12.8E \n',...
cycle,...
max(max(max(ux))),min(min(min(ux))),...
max(max(max(uy))),min(min(min(uy))));
fprintf('%12.4e %12.4e %12.4e %12.4e %12.4e \n',...
cycle,...
max(max(max(ux)))/uref_lb,min(min(min(ux)))/uref_lb,...
max(max(max(uy)))/uref_lb,min(min(min(uy)))/uref_lb);
%u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
u = reshape((uy),lx,ly);
%u(bbRegion) = nan;
figure(1);
imagesc(u');
axis equal off; drawnow
end
end
y_plot=[1:ly];
u_a(y_plot) = (rho_in-rho_out)/(6*nu*rho0*(lx-1))*(y_plot-(1-0.5)).*((ly+0.5)-y_plot);
u_x=reshape((ux),lx,ly);
figure(3); hold on; grid on;
plot(y_plot,u_x(10,y_plot)/uref_lb,'-c*');
plot(y_plot,u_x(30,y_plot)/uref_lb,'-g*');
plot(y_plot,u_x(40,y_plot)/uref_lb,'-b*');
plot(y_plot,u_a(y_plot)/uref_lb,'-ro');
xlabel('y');
ylabel('ux/u_{ref}');
legend('ux(x=10)','ux(x=30)','ux(x=50)','analytical ux');
Regards,
Ben