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];
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”
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
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.
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.