Zou & He pressure boundary problem

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_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_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_outlet
u_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];

Hello,

first about your bounceback boundaries. They seem problematic to me. You say that your numbering is like the one in OpenLB

1 8 7
\ | /
2—0---6
/ |
3 4 5

Then if you are doing fullway bounceback then the BC for any kind of walls should be

for i=1,…, 4
swap(f_i,f_(opposite(i)));
end for

and I’m not sure that’s what you are doing.

Then for your BC I didn’t check your code, but there are different possible sources of instabilities. One of them is the omega (1/tau, the inverse of the relaxation time). If it is too close to 2 (in fact we experienced instabilities for omega > 1.8) then the code becomes unstable. You can also reach too high velocities. Remember that we are in the low Mach number regime, therefore uMax << cs (cs is the speed of sound), and if your grad P becomes too high you may enter this regime.

Finally your corners seem a bit weird to me. You are not supposed to bounce back the whole f_i, but only the non equilibrium part of it.

I hope this helps a bit.

Correction: sorry numbering is not like the one in OpenLB but like this:
6 2 5
\ | /
3—0---1
/ | \
7 4 8

Anyway. Bounceback is

swap(f_1,f_3), swap(f_2,f_4), swap(f_5,f_7), swap(f_6,f_8). I think that you forgot swap(f_1,f_3). Except if you are doing halfway bounceback.

Thank you malaspin,
I tried with swaping (f_1,f_3) but no change in results.

Perhaps you haven’t understand me in your first message. In fact I haven’t a problem of stability but the velocity obtained from simulation is not the good compared with analytical solution; this only for high pressures. Also, for all my simulations the mach number remains low.

I repeat, it seems that I have the same problem as that of Mr Martin Thomas submited on PHORUM on July.

Can I send you by email a pdf report on which I explain the problem with curves? (I don’t know how to send curves here).

What do you mean by not good? What is the resolution? What happens to the error when if you increase the resolution?

Thank you dear malaspin for interest,
I will send you by email a brief pdf report in which is explained the problem. i’ll be at your disposal for any questions.

I answered to you by PM.

Thank you very much Orestis.
I’ll try to take account of your remarks.