Dear all,

I am new user of LBM, I’ve coded the D2Q9 model to smulate the Pouseuille flow with Zou & He pressure boundary conditions.

For low pressure gradient it seems that velocity obtained by simulation agrees very well with analytical solution. but when I increase pressure gradient results began to diverge.

I noticed that for high pressures, when I sketch the lengthwise change of pressure, there is an abnormal pressure increase in vicinity of channel outlet (boundary of low pressure) !!!. (Sorry I dont’t know how to atach graphs)

I don’t know if there is a problem in my code, or it is normal for these coditions.

I think that Mr MartinThomas had the same problem.

Thanks for help

The code I use looks like this (plane flow in x direction, code is for the all boundaries, D2Q9, numbering of f’s as in OpenLB):

// Horizontal boundaries “bounce-back boundaries”

for (x=1; x<lx-1;x++)

{

//bottom boundary

f2[x][0]=f4[x][0];

f5[x][0]=f7[x][0];//-(f1[x][0]-f3[x][0])/2;

f6[x][0]=f8[x][0];//+(f1[x][0]-f3[x][0])/2;

//upper boundary

f4[x][ly-1]=f2[x][ly-1];

f7[x][ly-1]=f5[x][ly-1];//+(f1[x][ly-1]-f3[x][ly-1])/2;

f8[x][ly-1]=f6[x][ly-1];//-(f1[x][ly-1]-f3[x][ly-1])/2;

}

// Vertical boundaries “Inlet & outlet pressure boundaries”

for (y=1; y<ly-1;y++)

{

//left boundary; inlet pressure

u_x=1.-(f0[0][y]+f2[0][y]+f4[0][y]+2*(f3[0][y]+f6[0][y]+f7[0][y]))/d_inlet;

f1[0][y]=f3[0][y]+d_inlet*u_x*2./3;

f5[0][y]=f7[0][y]-(f2[0][y]-f4[0][y])/2.+d_inlet*u_x/6;
f8[0][y]=f6[0][y]+(f2[0][y]-f4[0][y])/2.+d_inlet*u_x/6;

//right boundary; outlet pressure

u_x=(f0[lx-1][y]+f2[lx-1][y]+f4[lx-1][y]+2*(f1[lx-1][y]+f5[lx-1][y]+f8[lx-1][y]))/d_outlet-1.;

f3[lx-1][y]=f1[lx-1][y]-d_outlet

*u_x*2./3;

f6[lx-1][y]=f8[lx-1][y]-(f2[lx-1][y]-f4[lx-1][y])/2-d_outlet

*u_x/6;*

f7[lx-1][y]=f5[lx-1][y]+(f2[lx-1][y]-f4[lx-1][y])/2-d_outletu_x/6;

f7[lx-1][y]=f5[lx-1][y]+(f2[lx-1][y]-f4[lx-1][y])/2-d_outlet

// printf(“Pressure boundary condition\n”);

}

// Corner nodes “pressure & bounce_back boundaries”

// left bottom node x=0, y=0

f1[0][0]=f3[0][0];

f2[0][0]=f4[0][0];

f5[0][0]=f7[0][0];

f6[0][0]=f8[0][0]=(d_inlet-f0[0][0])/2-f3[0][0]-f4[0][0]-f7[0][0];

// left upper node x=0, y=ly-1

f1[0][ly-1]=f3[0][ly-1];

f4[0][ly-1]=f2[0][ly-1];

f8[0][ly-1]=f6[0][ly-1];

f5[0][ly-1]=f7[0][ly-1]=(d_inlet-f0[0][ly-1])/2-f2[0][ly-1]-f3[0][ly-1]-f6[0][ly-1];

// right bottom node x=lx-1, y=0

f3[lx-1][0]=f1[lx-1][0];

f2[lx-1][0]=f4[lx-1][0];

f6[lx-1][0]=f8[lx-1][0];

f5[lx-1][0]=f7[lx-1][0]=(d_outlet-f0[lx-1][0])/2-f1[lx-1][0]-f4[lx-1][0]-f8[lx-1][0];

// right upper node x=lx-1, y=ly-1

f3[lx-1][ly-1]=f1[lx-1][ly-1];

f4[lx-1][ly-1]=f2[lx-1][ly-1];

f7[lx-1][ly-1]=f5[lx-1][ly-1];

f6[lx-1][ly-1]=f8[lx-1][ly-1]=(d_outlet-f0[lx-1][ly-1])/2-f1[lx-1][ly-1]-f2[lx-1][ly-1]-f5[lx-1][ly-1];