# 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):

http://img239.imageshack.us/img239/1600/zvelns2.jpg

http://img239.imageshack.us/img239/7995/densly8.jpg

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];
u0[0]=0;
u0[1]=0;
u0[2]=f[0]+f[1]+f[2]+f[4]+f[5]+f[10]+f[11]+f[13]+f[14]+2*(f[17]+f[12]+f[9]+f[15]+f[7]);
u0[2]=(rhoIn-1.-u0[2])/rhoIn;
T t1=u0[2]*u0[2]+u0[1]*u0[1]+u0[0]*u0[0];

//bounce back of nonequilibrium part:
f[16]=f[7]+rhoIn/6.*(u0[0]-u0[2]);
f[18]=f[9]+rhoIn/6.*(u0[1]-u0[2]);
f[3]=f[12]-rhoIn/3.*u0[2];
f[6]=f[15]-rhoIn/6.*(u0[0]+u0[2]);
f[8]=f[17]-rhoIn/6.*(u0[1]+u0[2]);

//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];
f[16]=f[16]-partial_u_rho/2.;
f[18]=f[18]-partial_v_rho/2.;
f[6]=f[6]+partial_u_rho/2.;
f[8]=f[8]+partial_v_rho/2.;

//Collision
for(int iPop=0;iPop<19;iPop++){
feq=equilibrium(iPop, rhoIn, u0,  t1 );
f[iPop]=f[iPop]+omega*(feq-f[iPop]);
}

``````

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.

hi,
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”

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

u0[2]=(rhoIn-1.-u0[2])/rhoIn;

should be replaced by

u0[2]=(rhoIn-u0[2])/rhoIn;

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.