Zou and He boundary problems

Hi All,

I’m currently trying to implement the Zou and He pressure boundary condition but the result that I get looks wrong. On the following pictures you see the z-velocity and the density in a channel that has uniform pressures defined at z=0 and z=z_max, the higher pressure beeing at z=0 (even though the z-velocity might suggest that it is the other way round):



The images correspond to the state after streaming. I’m trying to resolve this problem for three days now without success and hope that someone here can help me. The code I use looks like this (this code is for the boundary at z=z_max, D3Q19, numbering of f’s as in OpenLB):

                        T rhoIn=1.;
			T u0[3];
			T t1=u0[2]*u0[2]+u0[1]*u0[1]+u0[0]*u0[0];

			//bounce back of nonequilibrium part:

			//redistribution of Momentum:
			T partial_u_rho=-f[5]-f[1]-f[4]+f[10]+f[14]+f[13]; 
			T partial_v_rho=f[5]+f[11]+f[13]-f[4]-f[2]-f[14];

			for(int iPop=0;iPop<19;iPop++){
				feq=equilibrium(iPop, rhoIn, u0,  t1 );

I don’t understand the very beginning of the code: you successively attribute two different values to u0[2]. Which one is the right one? And, shouldn’t rhoIn contain the sum of the incoming f’s, instead of 1?

Since it is a pressure boundary condition, I prescribe the density, right? Two velocity components get set to zero. u0[2] get’s computed in two steps just to make it a little more readable.

It seems your equations are wrong. you can see correct equations in the fllowing papers:

1- “Laboratory validation of lattice Boltzmann method for modeling pore-scale ?ow in granular materials”
by: Muhammed E. Kutay , Ahmet H. Aydilek , Eyad Masad

2- “On pressure and velocity boundary conditions for the lattice Boltzmann BGK model”
by: Qisu Zou and Xiaoyi He

Oh, my mistake: I didn’t realize this was a pressure boundary condition. Then, my question is, are you using the Skordos trick to reduce round-off errors, or are you using plain BGK? In the latter case (plain BGK), I think there’s a 1 in excess in your formula. That means, the expression


should be replaced by


Also, I would have the tendency to say that this formula should be multiplied by a factor -1, but I am unsure of the orientation and sign of the wall in your problem.

Actually, we recently wrote a review paper on boundary conditions, in which all formulas are stated step-by-step. You may want to look at equation 27 in the preprint, and solve this equation for rho, to cross-verify the formula in your code.

Thanks a lot for your answers, I’ll check the papers you mentioned today.

I think that I have the same as you but with the D2Q9 model, in fact there a inreasing in pressure at vicinity of the outlet (boundary of low pressure).
have find a solution?
I noticed that the problem appears only when pressure gradient rises.