# Still a poiseuille flow problem

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 :

1. 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.
2. 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

Hello!
I meet over your same problem. How do you settle down it?

Have you tried
a) using a body force instead of a inlet/outlet density?
b) using either Zou and He on the walls or using on-node bounce back?

My guess (and it is only a guess) is something ugly is happening at the corners. Or that due to a density variation (which is imposed by the boundary conditions) you are seeing waves which are creating a non-zero vertical velocity. This could be related to the initialisation (Junk has a paper about initialisation the LBE consistently).

Hello!
I use the period boundary condtion boundary and bounce back boundary on the wall, adding the body force into the collision term ,can gain good results compared with the analytic results.
but, I use Zou/He pressure boundary on the inlet and outlet, Zou/He velocity boundary on top and bottom wall, results were smaller than the analytical results. If adding the body force into the collision term, results were bigger than the analytical results. Above two condtions, there was the same pheomenon which the velocity of x-direction gradually increased, had the maximum values on the middle channel.
Can you resovle it ? How do you resolve it/? Can i share your code? thank you very much.

This may be a really stupid question (so I apologise if I sound like an idiot), but why do you need a body force and Zou and He in/out boundaries? Just one of these is sufficient for Poiseuille (i.e., an inlet and outlet density will give you a pressure gradient -so you can have in/out boundaries rather than periodic - or you can impose a driving force in a periodic domain. Perhaps you are doing something more fancy involving gravity or other forces but if you want to check your code then it is a good idea to do the simplest problem thing first.

Hi,
Sorry, I “solve” the problem few month ago.
Like pleb01 said “something ugly is happening at the corners”. With a proper corner treatment I have a zero (1d-15) y-velocity.
Of course you have to choice a problem description (ie with inlet/outlet or periodic) and in both cases you should have a zero y-velocity.
Regards,
Ben

Hello!
I used zou/he pressure boundary on the inlet and outlet, hence there was pressure gradient at x-direction flow. at the same time, used bounce back scheme on the top and bottom wall, four corner nodes were specially treated. i comply with zou/he’s “on pressure and velocity boundary conditions for the lattice boltzmann BGK model”,1997.I gain results smaller than the analytical results.

Hello bent
The following is my codes with c++,including four corner nodes special treatment.

#include
#include
#include
#include
#include
#include
#include

using namespace std;
const int Q = 9; //D2Q9 model
const int NX = 20; //x-direction
const int NY = 10; //y-direction
const double U = 0.1; //inlet velocity
double cr=25;

int e[Q][2] = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
double w[Q] = {4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36};
// Macroscopic variables
double rho[NX+1][NY+1],u[NX+1][NY+1][2],u0[NX+1][NX+1][2];
//Velocity Space
double f[NX+1][NY+1][Q],feqq[NX+1][NY+1][Q],F[NX+1][NY+1][Q]
,gx[NX+1][NY-1],gy[NX+1][NY-1],ftemp[NX+1][NY+1][Q],temp;
// 矩空间密度分布函数
double m[NX+1][NY+1][Q],meqq[NX+1][NY+1][Q],M[NX+1][NY+1][Q];

int opposite[Q]={0,3,4,1,2,7,6,5,6};
int image[NX+1][NY+1];
int i,j,k,ip,jp,n;
double Re,rho0,tau,nu,error,omega,cx,cy;
double s[Q]; //D2Q9松弛系数矩阵
//D2Q0转换矩阵
double MM[Q][Q] = {{1,1,1,1,1,1,1,1,1},{-1,-1,-1,-1,2,2,2,2,-4},{ -2,-2,-2,-2,1,1,1,1,4}
,{1,0,-1,0,1,-1,-1,1,0},{-2,0,2,0,1,-1,-1,1,0},{0,1,0,-1,1,1,-1,-1,0},{0,-2,0,2,1,1,-1,-1,0}
,{1,-1,1,-1,0,0,0,0,0},{ 0,0,0,0,1,-1,1,-1,0}};
double INM[Q][Q];//D2Q9转换矩阵的逆矩阵
//double jx,jy,jz;
//
double inrho = 1.1;
double outrho = 1.0;
double deltaP;
double umax;
double uc[NX+1][NY+1];
//
void init();
void imagec();
double feqc(int k,double rho,double u[2]);
double meqc(int k,double rho,double u[2]);
double force(int k,double u[2],double gx,double gy);

void evolution();
void outputdata();
void Error();

int main()
{
using namespace std;
init();
imagec();
for(n = 0;n<=40000 ;n++)
{
evolution();
if(n%100 == 0)
{
Error();
cout<<“it times”<<n<<" error="<<error<<endl;

if(n >= 100)
{
if(n%100== 0)
outputdata();
if (error<1.0e-6)
break;
}
}
}
return 0;
}
//
//初始子程序
void init()
{

rho0 = 1.0;
Re = 10;
// nu = 0.161; //U2.0cr/Re;
// tau = 3.0nu + 0.5;
omega= 0.9; //1.0/tau;
tau = 1.0/omega;
nu = 1.0/3.0
(tau-0.5);
cx=0.2NX;
cy=0.5
NY;

std::cout<<“tau=”<<tau<<endl;
std::cout<<“nu=”<<nu<<endl;
// system(“pause”);

for(i = 0; i <=NX ; i++) //初始化
for(j = 0; j <=NY ; j++)
{
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = rho0;

for(k = 0;k < Q;k++)
{
f[i][j][k]=feqc(k,rho[i][j],u[i][j]);
// feqq[i][j][k]=feqc(k,rho[i][j],u[i][j]);
// m[i][j][k]=meqc(k,rho[i][j],u[i][j]);
}
}
// u[0][0][0] = 0;
// u[0][NY][0] = 0;
//
deltaP =(inrho-outrho)/NX;
std::cout<<“deltaP=”<<deltaP<<endl;
// system(“pause”);
umax = deltaP/(2nu);
std::cout<<“umax=”<<umax<<endl;
system(“pause”);
}
//
// 图像处理子程序
void imagec()
{
for(i = 0; i <= NX; i++)
for(j = 0; j <= NY; j++)
{
// if (((i-cx)
(i-cx)+(j-cy)(j-cy))<=crcr)
// if(i==1&&i==NX)
// {
// image[i][j]=1;
// }
// else
// {
image[i][j]=0;
// }
}
for(i=0; i<=NX; i++)
{
image[i][0]=1;
}
for(i=0; i<=NX; i++)
{
image[i][NY]=1;
}

}
//
//计算平衡态函数-速度空间
double feqc(int k,double rho,double u[2])
{
double eu,uv,feq;
eu = (e[k][0]u[0] + e[k][1]u[1]);
uv = (u[0]u[0] + u[1]u[1]);
feq = w[k]rho(1.0 + 3.0
eu +4.5
eu
eu - 1.5
uv);
return feq;
}

//
//计算平衡态函数-速度矩空间
double meqc(int k,double rho,double u[2])
{
double jx,jy,meq;
jx = u[1];
jy = u[2];
switch(k)
{
case 0: meq = rho;
break;
case 1: meq = -2.0rho + 3.0(jxjx + jyjy);
break;
case 2: meq = rho -3.0*(jxjx + jyjy); //energy square
break;
case 3: meq = jx; //momentum in x-dir
break;
case 4: meq = -jx; //energy flux in x-dir—/-jx
break;
case 5: meq = jy; //momentum in y-dir
break;
case 6: meq = -jy; //energy flux in y-dir—/-jy
break;
case 7: meq = (jxjx - jyjy); //diagonal comp of stress tensor
break;
case 8: meq = jxjy;
break;
}
return meq;
}
//
//计算外力
double force(int k,double u[2],double gx,double gy)
{
double fx,fy,eu;
eu = (e[k][0]u[0] + e[k][1]u[1]);
fx=w[k]
(1.0-0.5/tau)
(3.0
(e[k][0]-u[0])+9.0eue[k][0])gx;
fy=w[k]
(1.0-0.5/tau)(3.0(e[k][1]-u[1])+9.0eue[k][1])*gy;
return (fx+fy);
}

//
void evolution()
{

// 计算平衡态函数-速度矩空间
// for(i = 0; i <=NX ;i++)
// for(j = 0; j <=NY ; j++)
// for(k = 0;k < Q;k++)
// {

// meqq[i][j][k]=meqc(k,rho[i][j],u[i][j]);
// }
//
//

// computing the macroscopic variables
for(i = 0;i <= NX;i++)
for(j = 0;j <= NY;j++)
{
u0[i][j][0] = u[i][j][0];
u0[i][j][1] = u[i][j][1];
rho[i][j] = 0;
u[i][j][0] = 0;
u[i][j][1] = 0;
// if(image[i][j]==0)
{
for(k = 0;k < Q;k++)
{

``````rho[i][j] +=f[i][j][k];
u[i][j][0] += e[k][0]*f[i][j][k];
u[i][j][1] += e[k][1]*f[i][j][k];
``````

}
u[i][j][0]= u[i][j][0]/rho[i][j];//+0.5gx[i][j]/rho[i][j];
u[i][j][1]= u[i][j][1]/rho[i][j];//+0.5
gy[i][j]/rho[i][j];
}
}

//computing the equilbrium distribution functions
for(i = 0; i <= NX ;i++)
for(j = 0; j <= NY ; j++)
for(k = 0;k < Q;k++)
{
if(image[i][j]==0)
{
feqq[i][j][k]=feqc(k,rho[i][j],u[i][j]);
}
}

//The standard bounceback
for(i = 1; i <= NX-1; i++)
for(j = 0; j <= NY; j++)
// for(k = 0; k< 9; k++)
{

// if (((i-cx)(i-cx)+(j-cy)(j-cy))<=cr*cr)
// if(image[i][j]==1)
// {
// u[i][j][0]=0;
// u[i][j][1]=0;
if(j==0)
{
f[i][j][5]=f[i][j][7];
f[i][j][2]=f[i][j][4];
f[i][j][6]=f[i][j][8];
}
if(j==NY)
{
f[i][j][7]=f[i][j][5];
f[i][j][4]=f[i][j][2];
f[i][j][8]=f[i][j][6];
}

// temp = f[i][j][k];
// f[i][j][k] = f[i][j][opposite[k]];
// f[i][j][opposite[k]]=temp;

//

// temp=f[i][j][1];
// f[i][j][1]=f[i][j][3];
// f[i][j][3]=temp;
//
// temp=f[i][j][2];
// f[i][j][2]=f[i][j][4];
// f[i][j][4]=temp;
//
// temp=f[i][j][5];
// f[i][j][5]=f[i][j][7];
// f[i][j][7]=temp;
//
// temp=f[i][j][6];
// f[i][j][6]=f[i][j][8];
// f[i][j][8]=temp;
// }
}

//collision
for(i=0;i<=NX;i++)
for(j=0;j<=NY;j++)
for(k=0;k<Q;k++)
{
// gx[i][j] = deltaP; //(inrho-outrho)/NX;

// gy[i][j] = 0.0;
if(image[i][j]==0)
{//if
// F[i][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau;
f[i][j][k] = f[i][j][k] + omega*(feqq[i][j][k]-f[i][j][k]);// +3.0w[k](e[k][0]deltaP+e[k][1]0.0)rho[i][j];
// F[i][j][k] = inverM
M[i][j][k];
//f[i][j][k] = f[i][j][k] + omega
(feqq[i][j][k]-f[i][j][k])+force(k,u[i][j],gx[i][j],gy[i][j]);
// M[i][j][k] = m[ip][jp][k] + s(k)
(meqq[ip][jp][k]-m[ip][jp][k]);
}//if

}

//

//streaming/propagation
for(i=0;i<=NX;i++)
for(j=0;j<=NY;j++)
for(k=0;k<Q;k++)
{
//
ip = i + e[k][0];
jp = j + e[k][1];

`````` if(ip>=0 && jp>=0 && ip<=NX && jp<=NY)
{
``````

//
ftemp[ip][jp][k]= f[i][j][k];
}

}
//
for(i=0;i<=NX;i++)
for(j=0;j<=NY;j++)
for(k=0;k<Q;k++)
{
f[i][j][k]=ftemp[i][j][k];
}

// streaming end
//Zou and He pressure boundary on inlet/west side
for(i = 0;i <=NX;i++)
for(j = 1;j <=NY-1;j++)
{
if(i==0)
{
rho[i][j] = inrho;
u[i][j][0] = 1.0 - (f[i][j][0]+f[i][j][2]+f[i][j][4]+2.0*(f[i][j][3]+f[i][j][6]+f[i][j][7]))/rho[i][j];
u[i][j][1] = 0.0;
//
f[i][j][1] = f[i][j][3] + 2.0/3.0 * rho[i][j]*u[i][j][0];
f[i][j][5] = f[i][j][7] + 1.0/6.0 * rho[i][j]*u[i][j][0] + 1.0/2.0 * (f[i][j][4] - f[i][j][2]);
f[i][j][8] = f[i][j][6] + 1.0/6.0 * rho[i][j]*u[i][j][0] + 1.0/2.0 * (f[i][j][2] - f[i][j][4]);
}

// Zou and He pressure boundary on outlet/east side
if(i==NX)
{
rho[i][j] = outrho;
u[i][j][0] = -1.0 + (f[i][j][0]+f[i][j][2]+f[i][j][4]+2.0*(f[i][j][1]+f[i][j][5]+f[i][j][8]))/rho[i][j];
u[i][j][1] = 0.0;

``````		  f[i][j][3] = f[i][j][1] -2.0/3.0 * rho[i][j]*u[i][j][0];
f[i][j][7] = f[i][j][5] - 1.0/6.0 * rho[i][j]*u[i][j][0] + 1.0/2.0 * (f[i][j][2] - f[i][j][4]);
f[i][j][6] = f[i][j][8] - 1.0/6.0 * rho[i][j]*u[i][j][0] + 1.0/2.0 * (f[i][j][4] - f[i][j][2]);
}
}
``````

//
//Zou and He velocity boundary on south/bottom side
for(i = 1;i <=NX-1;i++)
for(j = 0;j <=NY;j++)
{
if(j==0)
{
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = (f[i][j][0]+f[i][j][1]+f[i][j][3]+2.0*(f[i][j][4]+f[i][j][7]+f[i][j][8]))/1.0;
//
f[i][j][2] = f[i][j][4];
f[i][j][5] = f[i][j][7] + 1.0/2.0 * (f[i][j][1] - f[i][j][3]);
f[i][j][6] = f[i][j][8] + 1.0/2.0 * (f[i][j][3] - f[i][j][1]);
}

// Zou and He velocity boundary on north/top side
if(j==NY)
{
// rho[i][j] = outrho;
rho[i][j]= (f[i][j][0]+f[i][j][1]+f[i][j][3]+2.0*(f[i][j][2]+f[i][j][5]+f[i][j][6]))/1.0;
u[i][j][0] = 0;
u[i][j][1] = 0;
f[i][j][4] = f[i][j][2];
f[i][j][7] = f[i][j][5] + 1.0/2.0 * (f[i][j][1] - f[i][j][3]);
f[i][j][8] = f[i][j][6] + 1.0/2.0 * (f[i][j][3] - f[i][j][1]);
}
}

// Four corner nodes specific treatment
// Bottom left corner
i=0;
j=0;
// {
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = rho[i][j+1];
f[i][j][1] = f[i][j][3];
f[i][j][2] = f[i][j][4];
f[i][j][5] = f[i][j][7];
f[i][j][6] = 1.0/2.0*(rho[i][j]-(f[i][j][0]+2.0*(f[i][j][3]+f[i][j][4]+f[i][j][7])));
f[i][j][8] = f[i][j][6];
// }
// Top left corner
i=0;
j=NY;
// {
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = rho[i][j-1];
f[i][j][1] = f[i][j][3];
f[i][j][4] = f[i][j][2];
f[i][j][8] = f[i][j][6];
f[i][j][5] = 1.0/2.0*(rho[i][j]-(f[i][j][0]+2.0*(f[i][j][2]+f[i][j][3]+f[i][j][6])));
f[i][j][7] = f[i][j][5];
// }
// Bottom right corner
i=NX;
j=0;
// {
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = rho[i][j+1];
f[i][j][3] = f[i][j][1];
f[i][j][2] = f[i][j][4];
f[i][j][6] = f[i][j][8];
f[i][j][5] = 1.0/2.0*(rho[i][j]-(f[i][j][0]+2.0*(f[i][j][1]+f[i][j][4]+f[i][j][8])));
f[i][j][7] = f[i][j][5];
// }
// Top right corner
i=NX;
j=NY;
// {
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = rho[i][j-1];
f[i][j][3] = f[i][j][1];
f[i][j][4] = f[i][j][2];
f[i][j][7] = f[i][j][5];
f[i][j][6] = 1.0/2.0*(rho[i][j]-(f[i][j][0]+2.0*(f[i][j][1]+f[i][j][2]+f[i][j][5])));
f[i][j][8] = f[i][j][6];
// }

//ZouHepressureBC

} //evolve
//

void Error()
{
double temp1,temp2;
temp1 = 0;
temp2 = 0;
for(i = 0;i <=NX;i++)
for(j =0;j <=NY;j++)
{
temp1 += ((u[i][j][0] - u0[i][j][0])(u[i][j][0] - u0[i][j][0]) +
(u[i][j][1] - u0[i][j][1])
(u[i][j][1] - u0[i][j][1]));
temp2 += (u[i][j][0]*u[i][j][0] + u[i][j][1]*u[i][j][1]);
}
temp1 = sqrt(temp1);
temp2 = sqrt(temp2);
error = temp1/(temp2 + 1e-30);
}

void outputdata()
{
if(n%100==0)
{
ofstream file;
char temp[10];
string str;
itoa(n,temp,10);
str = string(temp);
str += “u.dat”;
file.open(str.c_str());
file<<“Title=“LBM Lid Driven Flow”\n”<<“VARIABLES=“X”,“Y”,“U”,“V”,“rho”\n”<<“ZONE T=“BOX”,I=”
<<NX+1<<",J="<<NY+1<<",F=POINT"<<endl;
if (file.is_open())
{

`````` for ( int i = 0; i <= NX; i++ )
for ( int j = 0; j <= NY; j++ )
{
``````

// file<<double(i)<<" “<<double(j)<<” “<<u[i][j][0]<<” “<<u[i][j][1]<<” “<<rho[i][j]<<endl;
// file<<setiosflags(ios::left);
file<<setiosflags(ios::right)<<double(i)<<” “<<double(j)<<” “<<setiosflags(ios::scientific)<<setprecision(6)<<u[i][j][0]<<” “<<u[i][j][1]<<” "<<rho[i][j]<<resetiosflags(ios::scientific)<<endl;
}
}
}

}

Sorry ydt, I don’t read a piece of c …
you have matlab solutions of the poiseuille flow here
http://lbmworkshop.com/?page_id=24

Im really late to the game but i believe you missed a rho in first line for the inlet Zou He condition.Inlet density is 1.3477. this is missing in the inlet term for fIn(2,in,col) = fIn(4,in,col) + 2/3*ux(:,in,col);