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_inletu_x2./3;
f5[0][y]=f7[0][y]-(f2[0][y]-f4[0][y])/2.+d_inletu_x/6;
f8[0][y]=f6[0][y]+(f2[0][y]-f4[0][y])/2.+d_inletu_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_outletu_x2./3;
f6[lx-1][y]=f8[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_outletu_x/6;
// 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];