# corner nodes 2d channel boundary condition (zou/he)

Hi there
I’ve been struggling to set the corner nodes for a 2d channel

here is the setup

i = inlet
o = outlet
Q = problematic nodes
X = WALL
``` QXXXXXXXXXXXXXXXXXXXXXXXQ i o i o i o i o QXXXXXXXXXXXXXXXXXXXXXXXQ ```
for an inlet FLUX boundary I ended up this:
-----------INLET--------
South:
f6 = 1/2 (rho - f0) + (f3 + f4 + f7)
f8 = f6
f5 = f7 + rho/2 (ux0)
f1 = f3
f2 = f4

## North: f1 = f3 f2 = f4 f6 = f8 f5 = 1/2 (rho - f0) - (f3 + f2 + f6) f7 = f5

on the outlet I have a zero gradient boundary!

anyway it does not work,
I want to check with anyone if my formulas are derived correctly

also if anyone has the pressure boundaries for the outlet corners and could share them it would be great.

the outlet is on the right side, but the phorum does not like my spaces

Dear jmazo

The explicit expressions for Zou and He BC’s can be found, for instance, in this reference:
http://arxiv.org/PS_cache/comp-gas/pdf/9611/9611001v1.pdf

Nevertheless, and for the incompressible model, here it goes an example on how implement Zou and He for the D2Q9 model:

``````% Inlet pressure boundary condition
x=1;
for y=2:LY-1
rho(x,y)=rho_inlet;
u(x,y)=rho(x,y)/rho0-(f(x,y,2)+f(x,y,4)+f(x,y,9)+2*(f(x,y,3)+f(x,y,6)+f(x,y,7)))./rho0;
v(x,y)=0;

f(x,y,1)=f(x,y,3)+2/3*rho0*u(x,y);
f(x,y,5)=f(x,y,7)+1/6*rho0*u(x,y)+1/2*(f(x,y,4)-f(x,y,2))+(1/2)*rho0*v(x,y);
f(x,y,8)=f(x,y,6)+1/6*rho0*u(x,y)+1/2*(f(x,y,2)-f(x,y,4))-(1/2)*rho0*v(x,y);
end

% Outlet pressure boundary condition
x=LX;
for y=2:LY-1
rho(x,y)=rho_outlet;
u(x,y)=-rho(x,y)/rho0+((f(x,y,9)+f(x,y,2)+f(x,y,4))+2*(f(x,y,5)+f(x,y,1)+f(x,y,8)))./rho0;
v(x,y)=0;

f(x,y,3)=f(x,y,1)-(2/3)*rho0*u(x,y);
f(x,y,7)=f(x,y,5)-(1/6)*rho0*u(x,y)-(1/2)*rho0*v(x,y)...
+(1/2)*(f(x,y,2)-f(x,y,4));
f(x,y,6)=f(x,y,8)-(1/6)*rho0*u(x,y)+(1/2)*rho0*v(x,y)...
+(1/2)*(f(x,y,4)-f(x,y,2));

end
``````

etc…

``````% Corner nodes

% Bottom left (inlet)
x=1; y=1;

u(x,y)=0;
v(x,y)=0;
rho(x,y)=rho_inlet;
f(x,y,1)=f(x,y,3);
f(x,y,2)=f(x,y,4);
f(x,y,5)=f(x,y,7);
f(x,y,6)=1/2*(rho(x,y)-(f(x,y,9)+f(x,y,1)+f(x,y,2)+f(x,y,3)+f(x,y,4)...
+f(x,y,5)+f(x,y,7)));
f(x,y,8)=f(x,y,6);

% Top left (inlet)
x=1; y=LY;

u(x,y)=0;
v(x,y)=0;
rho(x,y)=rho_inlet;
f(x,y,1)=f(x,y,3);
f(x,y,4)=f(x,y,2);
f(x,y,8)=f(x,y,6);
f(x,y,7)=1/2*(rho(x,y)-(f(x,y,9)+f(x,y,1)+f(x,y,2)+f(x,y,3)+f(x,y,4)...
+f(x,y,6)+f(x,y,8)));
f(x,y,5)=f(x,y,7));
``````

I hope that based on the previous examples you can generalize for the other planes and corners…

Regards and good work

Goncalo

Hi there

2 things

I’m guessing the rho0 = rho_Inlet?

and f9 = f0?

thanks

hello there

I have implemented everthing but I keep getting the same error

here is my example code.

``````
% mycavity
clear
iMax = 100000000000;

%grid size

lx = 400;
ly = 100;

ux  = 0.04; % horizontal lid velocity
uy  = 0;    % vertical lid velocity
rhoOut = 1.0;  %outlet pressure
obst_r = ly/10+1;
Re    = 100;  % Reynolds number

uMax   = 0.1;      % maximum velocity of Poiseuille inflow
nu     = uMax * 2.*obst_r / Re;  % kinematic viscosity
tau  = 1. / (3*nu+1./2.);
%tau = 1.0;

% D2Q9 LATTICE CONSTANTS

w0 = 4/9;
wm = 1/9;
wd = 1/36;

fc1 = 3.0;
fc2 = 9/2;
fc3 = 3/2;

%Arrays
% rho
rho = zeros(ly,lx);
u_x = zeros(ly,lx);
u_y = zeros(ly,lx);
vor = zeros(ly,lx);

% Equilibrium functions arrays
f0 = zeros(ly,lx);
f1 = zeros(ly,lx);
f2 = zeros(ly,lx);
f3 = zeros(ly,lx);
f4 = zeros(ly,lx);
f5 = zeros(ly,lx);
f6 = zeros(ly,lx);
f7 = zeros(ly,lx);
f8 = zeros(ly,lx);

% temporal storage
tmpf0 = zeros(ly,lx);
tmpf1 = zeros(ly,lx);
tmpf2 = zeros(ly,lx);
tmpf3 = zeros(ly,lx);
tmpf4 = zeros(ly,lx);
tmpf5 = zeros(ly,lx);
tmpf6 = zeros(ly,lx);
tmpf7 = zeros(ly,lx);
tmpf8 = zeros(ly,lx);

%---------------Create domain--------------------
% 1 = fluid
% 0 = Wall

D = ones(ly,lx);

for x=1:lx;
D(1,x) = 0;  %top wall
D(ly,x)= 0;  %bottom wall
end

roout = rhoOut;

% set the fluid to rest (rho=0, u=0) ==> fi(i) = t(i)
for y=1:ly;
L = ly-2;
y_phys = y-1.5;
for x=1:lx;

%if D(y,x) == 1 ;
%vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);

vxin = 0.04;

%           f0(y,x)  = w0 * rhoOut * (1                             - (3/2)*vxin*vxin);
%           f1(y,x)  = wm * rhoOut * (1 + 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
%           f2(y,x)  = wm * rhoOut * (1                             - (3/2)*vxin*vxin);
%           f3(y,x)  = wm * rhoOut * (1 - 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
%           f4(y,x)  = wm * rhoOut * (1                             - (3/2)*vxin*vxin);
%           f5(y,x)  = wd * rhoOut * (1 + 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
%           f6(y,x)  = wd * rhoOut * (1 - 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
%           f7(y,x)  = wd * rhoOut * (1 - 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
%           f8(y,x)  = wd * rhoOut * (1 + 3*vxin + (9/2)*vxin*vxin  - (3/2)*vxin*vxin);
f0(y,x)  = w0;
f1(y,x)  = wm;
f2(y,x)  = wm;
f3(y,x)  = wm;
f4(y,x)  = wm;
f5(y,x)  = wd;
f6(y,x)  = wd;
f7(y,x)  = wd;
f8(y,x)  = wd;
%end;

end
end

for i=1:ly;
for j=1:lx;

% if D(i,j) ~= 0;
rom =  f0(i,j) + f1(i,j)+ f2(i,j)+ f3(i,j)+ f4(i,j)+ f5(i,j)+ f6(i,j)+ f7(i,j)+ f8(i,j);

uxm =            f1(i,j)         - f3(i,j)         + f5(i,j)- f6(i,j)- f7(i,j)+ f8(i,j);
uym =                   + f2(i,j)         - f4(i,j)+ f5(i,j)+ f6(i,j)- f7(i,j)- f8(i,j);

rho(i,j) = rom;
u_x(i,j) = uxm / rom;
u_y(i,j) = uym / rom;

% end
end
end

for k=1:iMax% MAIN LOOP

%calculate macroscopic variables
%rho, u_x, u_y
for i=1:ly;
for j=1:lx;

%if D(i,j) ~= 0;
rom =  f0(i,j) + f1(i,j)+ f2(i,j)+ f3(i,j)+ f4(i,j)+ f5(i,j)+ f6(i,j)+ f7(i,j)+ f8(i,j);

uxm =            f1(i,j)         - f3(i,j)         + f5(i,j)- f6(i,j)- f7(i,j)+ f8(i,j);
uym =                   + f2(i,j)         - f4(i,j)+ f5(i,j)+ f6(i,j)- f7(i,j)- f8(i,j);

rho(i,j) = rom;
u_x(i,j) = uxm / rom;
u_y(i,j) = uym / rom;

%end
end
end

%INLET

x=1;
rho_inlet = 1.02;
rho0 = rho_inlet;
for y=2:ly-1
rho(y,x)=rho_inlet;
u_x(y,x)=rho(y,x)/rho0-(f2(y,x)+f4(y,x)+f0(y,x)+2*(f3(y,x)+f6(y,x)+f7(y,x)))./rho0;
u_y(y,x)=0;

f1(y,x)=f3(y,x)+2/3*rho0*u_x(y,x);
f5(y,x)=f7(y,x)+1/6*rho0*u_x(y,x)+1/2*(f4(y,x)-f2(y,x))+(1/2)*rho0*u_y(y,x);
f8(y,x)=f6(y,x)+1/6*rho0*u_x(y,x)+1/2*(f2(y,x)-f4(y,x))-(1/2)*rho0*u_y(y,x);
end

% Bottom left (inlet)

x=1; y=ly;

u_x(y,x)=0;
u_y(y,x)=0;
rho(y,x)=rho_inlet;
f1(y,x)=f3(y,x);
f2(y,x)=f4(y,x);
f5(y,x)=f7(y,x);

f6(y,x)=1/2*(rho(y,x)-(f0(y,x)+f1(y,x)+f2(y,x)+f3(y,x)+f4(y,x)...
+f5(y,x)+f7(y,x)));

f8(y,x)=f6(y,x);

% Top left (inlet)
x=1; y=1;

u_x(y,x)=0;
u_y(y,x)=0;
rho(y,x)=rho_inlet;
f1(y,x)=f3(y,x);
f4(y,x)=f2(y,x);
f8(y,x)=f6(y,x);

f7(y,x)=1/2*(rho(y,x)-(f0(y,x)+f1(y,x)+f2(y,x)+f3(y,x)+f4(y,x)...
+f6(y,x)+f8(y,x)));

f5(y,x)=f7(y,x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outlet pressure boundary condition %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=lx;

rho_outlet = 1;
rho0 = rho_outlet;
for y=2:ly-1
rho(y,x)=rho_outlet;
u_x(y,x)=-rho(y,x)/rho0+((f0(y,x)+f2(y,x)+f4(y,x))+2*(f5(y,x)+f1(y,x)+f8(y,x)))./rho0;
u_y(y,x)=0;

f3(y,x)=f1(y,x)-(2/3)*rho0*u_y(y,x);

f7(y,x)=f5(y,x)-(1/6)*rho0*u_x(y,x)-(1/2)*rho0*u_y(y,x)...
+(1/2)*(f2(y,x)-f4(y,x));

f6(y,x)=f8(y,x)-(1/6)*rho0*u_x(y,x)+(1/2)*rho0*u_y(y,x)...
+(1/2)*(f4(y,x)-f2(y,x));

end

% for i=1:ly
%          f3(i,lx) = f1(i,lx-1);
%          f7(i,lx) = f5(i,lx-1);
%          f6(i,lx) = f8(i,lx-1);
% end

% Top right (outlet)

x=lx; y=1;

u_x(y,x)=0;
u_y(y,x)=0;
rho(y,x)=rho_outlet;
f1(y,x)=f3(y,x);
f2(y,x)=f4(y,x);
f5(y,x)=f7(y,x);

f6(y,x)=1/2*(rho(y,x)-(f0(y,x)+f1(y,x)+f2(y,x)+f3(y,x)+f4(y,x)...
+f5(y,x)+f7(y,x)));

f8(y,x)=f6(y,x);

% Bottom right (outlet)
x=lx; y=ly;

u_x(y,x)=0;
u_y(y,x)=0;
rho(y,x)=rho_outlet;
f1(y,x)=f3(y,x);
f4(y,x)=f2(y,x);
f8(y,x)=f6(y,x);

f7(y,x)=1/2*(rho(y,x)-(f0(y,x)+f1(y,x)+f2(y,x)+f3(y,x)+f4(y,x)...
+f6(y,x)+f8(y,x)));

f5(y,x)=f7(y,x);

%=====================================================================
%============== COLLISION ============================================
%=====================================================================

%Compute equilibrium functions
for i=1:ly;
for j=1:lx;
% if D(i,j) ~= 0;
%get needed values
rq = rho(i,j);
vx = u_x(i,j);
vy = u_y(i,j);

rt0 = w0 * rq;
rt1 = wm * rq;
rt2 = wd * rq;

v_sq_term = fc3*(vx*vx + vy*vy);

% Evaluate the local equilibrium f values in all directions
f0eq = rt0 * (1 - v_sq_term);
f1eq = rt1 * (1 + fc1*vx +        fc2*vx*vx -                v_sq_term);
f2eq = rt1 * (1 + fc1*vy +        fc2*vy*vy -                v_sq_term);
f3eq = rt1 * (1 - fc1*vx +        fc2*vx*vx -                v_sq_term);
f4eq = rt1 * (1 - fc1*vy +        fc2*vy*vy -                v_sq_term);
f5eq = rt2 * (1 + fc1*(vx + vy) + fc2*(vx + vy)*(vx + vy) -  v_sq_term);
f6eq = rt2 * (1 + fc1*(-vx + vy)+ fc2*(-vx + vy)*(-vx + vy)- v_sq_term);
f7eq = rt2 * (1 + fc1*(-vx - vy)+ fc2*(-vx - vy)*(-vx - vy)- v_sq_term);
f8eq = rt2 * (1 + fc1*(vx - vy) + fc2*(vx - vy)*(vx - vy) -  v_sq_term);

tmpf0(i,j) = f0(i,j) - (f0(i,j) - f0eq)*tau;
tmpf1(i,j) = f1(i,j) - (f1(i,j) - f1eq)*tau;
tmpf2(i,j) = f2(i,j) - (f2(i,j) - f2eq)*tau;
tmpf3(i,j) = f3(i,j) - (f3(i,j) - f3eq)*tau;
tmpf4(i,j) = f4(i,j) - (f4(i,j) - f4eq)*tau;
tmpf5(i,j) = f5(i,j) - (f5(i,j) - f5eq)*tau;
tmpf6(i,j) = f6(i,j) - (f6(i,j) - f6eq)*tau;
tmpf7(i,j) = f7(i,j) - (f7(i,j) - f7eq)*tau;
tmpf8(i,j) = f8(i,j) - (f8(i,j) - f8eq)*tau;
%end;

end
end

% Bounce-Back
for i=1:ly;
for j=1:lx;
if D(i,j) == 0;
f1old = f1(i,j);
f2old = f2(i,j);
f3old = f3(i,j);
f4old = f4(i,j);
f5old = f5(i,j);
f6old = f6(i,j);
f7old = f7(i,j);
f8old = f8(i,j);

tmpf1(i,j) = f3old;
tmpf2(i,j) = f4old;
tmpf3(i,j) = f1old;
tmpf4(i,j) = f2old;
tmpf5(i,j) = f7old;
tmpf6(i,j) = f8old;
tmpf7(i,j) = f5old;
tmpf8(i,j) = f6old;
end
end
end

%=====================================================================
%=====================================================================

% Streaming
for i=1:ly
in = i - 1;
is = i + 1;

if i == 1
in = 1;
end

if i == ly
is = ly;
end

for j=1:lx
je = j + 1;
jw = j - 1;

if j == 1
jw = 1;
end

if j == lx
je = lx;
end

f0(i,j) = tmpf0(i,j);
f1(i,je) = tmpf1(i,j);
f2(in,j) = tmpf2(i,j);
f3(i,jw) = tmpf3(i,j);
f4(is,j) = tmpf4(i,j);
f5(in,je) = tmpf5(i,j);
f6(in,jw) = tmpf6(i,j);
f7(is,jw) = tmpf7(i,j);
f8(is,je) = tmpf8(i,j);

end

end

% VISUALIZATION
if (mod(k,1)==0)
u = sqrt(u_x.^2+u_y.^2);

imagesc(u);
%contour(u);

axis equal off; drawnow
end

end

``````

Hello

Since I use rho0 because the expressions that I sent you are for the incompressible model. Nevertheless, doing what you did, i.e saying that rho0=rho (the density fluctuation) you are using the standard model and that is also OK.

Regarding your code. I did a quick check and notice two things.

First, you use bounce-back at the walls so it is not problematic to extended the bounceback to to corners between the walls and the inlets/outlets… In fact Zou/He at a corner with zero velocity is equivalent to the expression of bounceback at the node.

Second, there are distribution functions at the corners that are ill-defined! Check the bottom right and top right corners and compare them to mine (see below). Apart from that everything appears to be correct in the BC model. If after you have put the corners right your problem still subsists I would say that your error source comes from some other part of your code…

% Bottom right (outlet)
x=LX; y=1;

``````u(x,y)=0;
v(x,y)=0;
rho(x,y)=rho_outlet;
f(x,y,3)=f(x,y,1);
f(x,y,2)=f(x,y,4);
f(x,y,6)=f(x,y,8);
f(x,y,7)=1/2*(rho(x,y)-(f(x,y,9)+f(x,y,1)+f(x,y,2)+f(x,y,3)+f(x,y,4)...
+f(x,y,6)+f(x,y,8)));
f(x,y,5)=f(x,y,7);

% Top right (outlet)
x=LX; y=LY;

u(x,y)=0;
v(x,y)=0;
rho(x,y)=rho_outlet;
f(x,y,3)=f(x,y,1);
f(x,y,4)=f(x,y,2);
f(x,y,7)=f(x,y,5);
f(x,y,6)=1/2*(rho(x,y)-(f(x,y,9)+f(x,y,1)+f(x,y,2)+f(x,y,3)+f(x,y,4)...
+f(x,y,5)+f(x,y,7)));
f(x,y,8)=f(x,y,6);
``````

Good luck

Regards

Goncalo

Hi Goncalo,

I have sent you a message about some questions about Zou/He boundary conditon in D3Q19 model. Please have a check. Many thanks.

Beryl

Dear jmazo,

I have gone through your code, are you able to simulation the problem correctly?, i am beginner, i just wondering if you could share your correct code, this will through some light on my problem. Thanks

Winter

Hi Goncalo,

I have a question for you in this regard. Can we use Zou-He wall velocity boundary condition for no-slip by setting the velocites to zero. Please help.

I am using the original compressible model and I did exactly what you mentioned and used Zou-He for no slip conditions. I get the steady state solution but the y-velocity is not zero and is the order of 10e-6.

Also I am not able to get the analytical maximum velocity. I am getting lesser that analytical max velocity.

``````
% Zou-He No-Slip Boundary, except for corner nodes.
if y==1   % Bottom wall

f(x,y,2) = f(x,y,4);
f(x,y,5) = f(x,y,7) - 0.5*(f(x,y,1) - f(x,y,3));
f(x,y,6) = f(x,y,8) + 0.5*(f(x,y,1) - f(x,y,3));

elseif y==YMAX  % Top Wall

f(x,y,4) = f(x,y,2);
f(x,y,7) = f(x,y,5) + 0.5*(f(x,y,1) - f(x,y,3));
f(x,y,8) = f(x,y,6) - 0.5*(f(x,y,1) - f(x,y,3));

end

``````

Thanks,
Narender

Hi Narender

Let me just recapitulate. You are simulating a plane Poiseuille flow driven by a given pressure gradient. You are using pressure BCs at the inlet/outlet and zero velocity BCs at the walls. All these are modeled using the Zou/He BC model. Furthermore, you ensure the continuity between walls and inlet/outlet through a proper corner treatment as that posted here before.

If you did do that, and you are using the incompressible LB model you will recover the exact (up to the machine accuracy) hydrodynamic solution.

If you are using the standard (weakly-compressible) LB model you can never obtain the exact solution and this is explained as follows.
Doing Chapman-Enskog analysis in the standard LB model you recover the compressible Navier-Stokes equations. However, the Poiseuille flow is a solution of the incompressible Navier-Stokes equations. Hence, as long as there is a pressure gradient , as in this simple problem, there will always be a density gradient implying that the strict incompressibility solution is violated. So don’t expect reaching this solution! This so-called compressibility error is always finite, though it decreases quadratically with the mesh size.

Another exercise that we can do to prove that the flow you get isn’t exactly the Poiseuille solution is checking the velocity streamwise invariance, i.e. the so-called fully-developed condition, which is characteristic of this flow.
Mass conservation at steady-state establishes that (rhou)in=(rhou)out. Moreover, we are using pressure BCs, which means (rho)in>(rho)out. Therefore, we will have (u)in<(u)out. This means that at steady-state the standard LB model can’t even respect the fully developed condition in this simple flow!

In summary, the problem you are reporting does not come from the boundary conditions. If you expect machine accurate results you have to blame the standard model.

I hope I have helped

Regards

Goncalo

Hi Goncalo,

Thanks a lot for your explanation of the problem. It was very clear and helpful. Yes, I observe that when I fine tune the mesh I go closer to the accurate solution. I get some y-component velocity (with compressible model and zou-he bc models) just after the wall nodes. Is this common? Please give some insight into this.

Thanks,
Narender.

Hello Narender,

Yes, I expect the y-velocity component to be different from zero as a result of the fact that the flow is not exactly fully developed.
I understand that you force uy=0 at the inlet/outlet boundaries but inside the domain due to compressibility ux(0)>ux(LX) and then because of continuity (i.e. mass conservation at steady-state) you will have to have uy /= 0! A nice way to see this is to plot the flow streamlines. You will see that as a result of the deceleration in ux the streamlines won’t be parallel as in an exact fully developed flow (with uy=0) but approach each other very slowly, the slope predicted by the (uy/ux) ratio!

Regards,

Goncalo

Hi Goncalo,

Thank you for the reply. It was helpful for me in understanding the problem better. My next step is to benchmark velocity boundary condition on inlet and outflow boundary condition on the outlet.

Can you please let me know if there is a way of specifying a corner boundary condition for a flow velocity boundary condition. (We get a single point discontinuity). Also what do we do for an outflow boundary condition.

Regards,
Narender.

Hi Narender,

Please give a look at some other older posts in this forum. I think you might find some answers. Regarding corner treatment check for instance:

Regarding the outlflow boundary condititon I have no experience I know that some people simply copy f(LX-1, to f(LX,:), i.e. the last flow region to the outlet region. Yet, I am not sure if this procedure is the most accurate. Please take a look at the following paper, perhaps it will shed some light over this problem in particular.

Regards,

Goncalo

Hi Goncalo,

Thanks for the reply. My results for the poiseulle flow (Incompressible LBGK, SRT) are as follows. Please see if they look good to you. Also my results vary from the analytical solution as I change XDim, YDim, Tau. It seems like I need to get a perfect combination of XDim, YDim, Tau to match it exactly with analytical solution.

``````
The value of TAU           :0.87
The value of NU            :0.12333
The value of XMAX          :17
The value of YMAX          :5
The value of RHO_0         :3
The value of P Inlet       :1.1
The value of P Outlet      :1
The value of U Max         :0.033784
Enter the number of iterations: 3000
>> u

u =

0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0
0    0.0253    0.0338    0.0253         0

>> v

v =

1.0e-015 *

0         0         0         0         0
0   -0.0463    0.0093    0.0324   -0.0046
0   -0.0278   -0.0324   -0.0185   -0.0046
-0.0046   -0.0463   -0.0509   -0.0278         0
-0.0046   -0.1573   -0.0925   -0.0093         0
0   -0.1018   -0.1064   -0.0278   -0.0046
-0.0046   -0.0370   -0.0833   -0.0463         0
0   -0.0601   -0.0046         0    0.0046
0   -0.0416   -0.0231   -0.0046         0
0.0046   -0.0648   -0.0046    0.0694   -0.0046
0   -0.0093    0.0139    0.0278         0
0    0.0231    0.0555    0.0555         0
-0.0046   -0.0324    0.0185    0.0648   -0.0046
-0.0046   -0.0601    0.0046    0.0046    0.0046
-0.0046   -0.0185    0.0046    0.0046         0
0         0    0.0093    0.0093    0.0046
0   -0.0046   -0.0046    0.0046    0.0046

>> rho

rho =

3.3000    3.3000    3.3000    3.3000    3.3000
3.2812    3.2812    3.2812    3.2813    3.2813
3.2625    3.2625    3.2625    3.2625    3.2625
3.2437    3.2437    3.2437    3.2437    3.2437
3.2250    3.2250    3.2250    3.2250    3.2250
3.2062    3.2062    3.2062    3.2062    3.2062
3.1875    3.1875    3.1875    3.1875    3.1875
3.1687    3.1687    3.1687    3.1687    3.1687
3.1500    3.1500    3.1500    3.1500    3.1500
3.1312    3.1312    3.1312    3.1312    3.1312
3.1125    3.1125    3.1125    3.1125    3.1125
3.0937    3.0937    3.0937    3.0937    3.0937
3.0750    3.0750    3.0750    3.0750    3.0750
3.0562    3.0562    3.0562    3.0562    3.0562
3.0375    3.0375    3.0375    3.0375    3.0375
3.0187    3.0187    3.0187    3.0187    3.0187
3.0000    3.0000    3.0000    3.0000    3.0000

``````

In the above run, I choose tau=0.87 and any value below that doesn’t work for this grid. Is this normal. Please clarify.

Thank you very much for your support.

Regards,
Narender

Hi again,

The results seem correct. As far as the v velocity concerns, the Hagen-Poiseuille solution is reproduced up to the machine accuracy. The u and rho also appear correct, up to the digits you display.

You have to keep in mind that when doing mesh refinement tests the Reynolds number, Re=LYUmax/kvisc, must be fixed. (Note kvisc = 1/3(tau-1/2) )

Generally, tau is also a priori fixed, as its value influences stability and also accuracy in some boundary schemes like the bounceback method.

Having the Re and the tau values defined I am still free to choose a value for LY or U max. (Note: LY = nº of nodes in width * delta_x).

I usually choose to vary LY, by halving it. This way I will obtain the U max, which tells me the Mach number!

I think you get my point… in your case what is establishing the flow Re is the P inlet - P outlet. Hence, you have to find a relation here!

However, and since this flow is an analytical solution of the LB method, you are lucky here! In theory, as long as P inlet - P outlet is small, i.e. still within the low Mach number range then results will always match the hydrodynamic solution exactly.

Hence, what I think is happening is that you increase the mesh size but you don’t change the number of iteration steps. Hence, you are not allowing your flow to reach the time-independent hydrodynamic solution!!! Am I right?

Regards

Goncalo

Hi Goncalo,

Thank you for checking my results. I allow sufficient time but the code blows up (NaN).

You are right, I am doing the mesh refinement and I dared to forget that my flow is actually changing with this. I mean, the values of pressure, rho0 etc., should also be changed to simulate the same flow (same Re).

Regards,
Narender

Hi jmazo,
I’ve looked at your code, I found a flaw it
To fix the problem yourself, you change the location of the boundary conditions imposed
1 - Put the boundary conditions at the wall between the COLLISION STEP and Streeming
2 - boundary conditions at the inlet and outlet should be placed after Streeming

%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%
% program %
%
% c
% c C7 C3 C6 ^ y
% c \ | / |
% c C4-C1-C2 |
% c / | \ |
% c C8 C5 C9 -----> x
% c
% c
% C
% c the lattice nodes are numbered as follows.
% c
% c
% c 1 ----------–…------------
% c | | | | | |
% c | | | | | |
% c 2 ----------–…------------
% c | | | | | |
% c : : : : : :
% c : : : : : :
% c | | | | | |
% c ly-1 ----------–…------------
% c | | | | | |
% c | | | | | |
% c ly ----------–…------------
% c 1 2 lx-1 lx
%% %%%%%%%%%%%%%%%
clc
clear all
close all
tic
% GENERAL FLOW CONSTANTS

ly=20;
lx=10*ly;

r=ones(1,lx);

pOut=1/3;

pIn=2*pOut;

Kno=.0194;
%Kni=.011;

fluid=zeros(ly,lx);
fluid((2:ly-1),:)=1;
bbRegion = find(fluid==0);

maxT = 40000; % total number of iterations
tPlot = 10; % cycles for graphical output

% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];

[y,x] = meshgrid(1:ly,1:lx);
% INITIAL CONDITION: (rho=1, u=0) ==> fIn(i) = t(i)
fEq = zeros(9,ly,lx);
fOut = zeros(9,ly,lx);
ro=ones(1,lxly);
ro(:,:)=3
pOut;
fIn = reshape( t’ * ro, 9, ly, lx);

rho=zeros(ly,lx);
ux =zeros(ly,lx);
uy =zeros(ly,lx);

err=1;
uu=ones(ly,lx);
cycle=0;
pBar=zeros(1,lx);
KnBar=zeros(1,lx);
TauBar=zeros(1,lx);

while cycle<=300000
cycle=cycle+1;
% MACROSCOPIC VARIABLES
for j=1:ly
for i=1:lx
%if fluid(j,i)==1
rho(j,i)=fIn(1,j,i)+fIn(2,j,i)+fIn(3,j,i)+fIn(4,j,i)…
+fIn(5,j,i)+fIn(6,j,i)+fIn(7,j,i)+fIn(8,j,i)+fIn(9,j,i);
end
end
sumrhoIn =(sum(rho(:,1 )))/ly;
sumrhoOut=(sum(rho(:,lx)))/ly;

``````rho(:,1) =rho(:,1 )*(3*pIn/sumrhoIn  );
rho(:,lx)=rho(:,lx)*(3*pOut/sumrhoOut);

for j=1:ly
for i=1:lx
ux (j,i)=((fIn(2,j,i)+fIn(6,j,i)+fIn(9,j,i))...
-(fIn(4,j,i)+fIn(7,j,i)+fIn(8,j,i)))./rho(j,i);
uy (j,i)=((fIn(3,j,i)+fIn(6,j,i)+fIn(7,j,i))...
-(fIn(5,j,i)+fIn(8,j,i)+fIn(9,j,i)))./rho(j,i);
%else
rho([1,ly],i)=1;
ux ([1,ly],i)=0;
uy ([1,ly],i)=0;
%end
end

end

for i=1:lx
pBar(1,i)=sum(rho(1:ly,i))/(3*ly);
KnBar(1,i)=(Kno*pOut)/pBar(1,i);
%KnBar(1,i)=(Kni*pIn)/pBar(1,i);
TauBar(1,i)=.5+(sqrt(6/pi))*(KnBar(1,i)*ly);
%TauBar(1,i)=.5+(sqrt(6/pi))*(KnBar(1,i)*lx);
end
``````

%%%%%%%%%%%%%%
% COLLISION STEP
%%%%%%%%%%%%%%
for i=1:lx
for j=1:ly
if fluid(j,i)==1
for k=1:9
cu = 3*(cx(k)ux(j,i)+cy(k)uy(j,i));
fEq(k,j,i)=rho(j,i)t(k)(1+cu+(.5
cu
cu)-…
((3/2)((ux(j,i)^2)+(uy(j,i)^2))));
fOut(k,j,i)=fIn(k,j,i) - (1/TauBar(1,i))
(fIn(k,j,i)-fEq(k,j,i));
end
end
end
end

``````A1=0.8183;
A2=0.6531;
``````

for i=1:lx
r(i)=1;%1/(1+( (sqrt(pi/6))* ((1/(4KnBar(1,i)(lx^2)))+A1+((2*A2-(8/pi))*KnBar(1,i)) ) ));
end

%%%%%%%%%%%%%%%%%%%
% LOWER WALL
%%%%%%%%%%%%%%%%%%%
for i=1:lx
fOut(3,ly,i)=fIn(5,ly,i);
fOut(6,ly,i)=r(i)*fIn(8,ly,i)+(1-r(i))*fIn(9,ly,i);
fOut(7,ly,i)=r(i)*fIn(9,ly,i)+(1-r(i))*fIn(8,ly,i);
end
%%%%%%%%%%%%%%%%%%%
% UPPER WALL
%%%%%%%%%%%%%%%%%%%
for i=1:lx
fOut(5,1,i)=fIn(3,1,i);
fOut(8,1,i)=r(i)*fIn(6,1,i)+(1-r(i))*fIn(7,1,i);
fOut(9,1,i)=r(i)*fIn(7,1,i)+(1-r(i))*fIn(6,1,i);
end

% for i=1:9
% fOut(i,bbRegion) = fIn(opp(i),bbRegion);
% end

% for i=1:9
% fIn(i,:, = circshift(fOut(i,:,:), [0,cy(i),cx(i)]);
% end

%%%%%%%%%%%%%%%%
%%%% Streeming
%%%%%%%%%%%%%%%%
for kk=1
%%% EAST
k=1+1;
for j=1:ly
for i=lx:-1:2
fIn(k,j,i)=fOut(k,j,i-1);
end
end
%%% NORTH
k=2+1;
for j=1:ly-1
for i=1:lx
fIn(k,j,i)=fOut(k,j+1,i);
end
end
%%% WEST
k=3+1;
for j=1:ly
for i=1:lx-1
fIn(k,j,i)=fOut(k,j,i+1);
end
end
%%% SOUTH
k=4+1;
for j=ly:-1:2
for i=1:lx
fIn(k,j,i)=fOut(k,j-1,i);
end
end
%%% NORTH-EAST
k=5+1;
for j=1:ly-1
for i=lx:-1:2
fIn(k,j,i)=fOut(k,j+1,i-1);
end
end
%%% NORTH-WEST
k=6+1;
for j=1:ly-1
for i=1:lx-1
fIn(k,j,i)=fOut(k,j+1,i+1);
end
end
%%% SOUTH-WEST
k=7+1;
for j=2:ly
for i=1:lx-1
fIn(k,j,i)=fOut(k,j-1,i+1);
end
end
%%% SOUTH-EAST
k=8+1;
for j=2:ly
for i=2:lx
fIn(k,j,i)=fOut(k,j-1,i-1);
end
end
%%% CENTER
k=0+1;
for j=1:ly
for i=1:lx
fIn(k,j,i)=fOut(k,j,i);
end
end

end

%%%%%%%%%%%%%%%%%%%
% INLET BOUNDARY0
% &&&&&&&&
% OUTLET BOUNDARY
%%%%%%%%%%%%%%%%%%%

for j=2:ly-1
i=1;

``````%%%%%%%%%%
% Zou He
%%%%%%%%%%
RhoIn=3*pIn;
Uin=1-(sum(fIn([3,5,1],j,i))+2*sum(fIn([7,4,8],j,i)))/RhoIn;
fIn(2,j,i) = fIn(4,j,i) + 2/3*RhoIn*Uin;
fIn(9,j,i) = fIn(7,j,1) + 1/2.*(fIn(3,j,1)-fIn(5,j,1)) + 1/6*RhoIn*Uin;
fIn(6,j,i) = fIn(8,j,1) - 1/2.*(fIn(3,j,1)-fIn(5,j,1)) + 1/6*RhoIn*Uin;
``````

end

% % Bottom left (inlet)

x=1; y=ly;

ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pIn;
fIn(2,y,x)=fIn(4,y,x);
fIn(3,y,x)=fIn(5,y,x);
fIn(6,y,x)=fIn(8,y,x);
fIn(7,y,x)=1/2
(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(6,y,x)+fIn(8,y,x)));
fIn(9,y,x)=fIn(7,y,x);

% Top left (inlet)
x=1; y=1;

ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pIn;
fIn(2,y,x)=fIn(4,y,x);
fIn(5,y,x)=fIn(3,y,x);
fIn(9,y,x)=fIn(7,y,x);
fIn(8,y,x)=1/2
(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(7,y,x)+fIn(9,y,x)));
fIn(6,y,x)=fIn(8,y,x);

for j=2:ly-1
i=lx;

``````%%%%%%%%%%
% Zou He
%%%%%%%%%%
``````

% RhoOut=3pOut;
% Uout=-1+(sum(fIn([3,5,1],j,i))+2
sum(fIn([2,6,9],j,i)))/RhoOut;
% fIn(4,j,i) = fIn(2,j,i)- 2/3RhoOutUout;
% fIn(8,j,i) = fIn(6,j,i)+ 1/2.(fIn(3,j,i)-fIn(5,j,i))- 1/6RhoOutUout;
% fIn(7,j,i) = fIn(9,j,i)- 1/2.
(fIn(3,j,i)-fIn(5,j,i))- 1/6RhoOutUout;

``````%%%%%%%%%%%%%%
% Open Buondry
%%%%%%%%%%%%%%
fIn(4,j,i)= 2*fIn(4,j,i-1) - fIn(4,j,i-2);
fIn(8,j,i)= 2*fIn(8,j,i-1) - fIn(8,j,i-2);
fIn(7,j,i)= 2*fIn(7,j,i-1) - fIn(7,j,i-2);
``````

end

%Top right (outlet)

x=lx; y=1;

ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pOut;
fIn(4,y,x)=fIn(2,y,x);
fIn(5,y,x)=fIn(3,y,x);
fIn(8,y,x)=fIn(6,y,x);
fIn(7,y,x)=1/2
(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(6,y,x)+fIn(8,y,x)));
fIn(9,y,x)=fIn(7,y,x);

% Bottom right (outlet)
x=lx; y=ly;

ux(y,x)=0;
uy(y,x)=0;
rho(y,x)=3pOut;
fIn(4,y,x)=fIn(2,y,x);
fIn(3,y,x)=fIn(5,y,x);
fIn(7,y,x)=fIn(9,y,x);
fIn(8,y,x)=1/2
(rho(y,x)-(fIn(1,y,x)+fIn(2,y,x)+fIn(3,y,x)+fIn(4,y,x)+fIn(5,y,x)+fIn(7,y,x)+fIn(9,y,x)));
fIn(6,y,x)=fIn(8,y,x);

``````% VISUALIZATION
if (mod(cycle,tPlot)==0)
disp(cycle)
u =(sqrt(ux.^2+uy.^2));
%u(bbRegion) = nan;
subplot(3,2,1)
%figure(8)
uIn=sum(sqrt(ux(2:ly-1,1).^2+uy(2:ly-1,1).^2))/(ly-2);
imagesc(u(1:ly,1:lx)./uIn);
%colorbar
axis equal off; drawnow

yy=linspace (0,1,ly-2);
%figure(3);
subplot(3,2,2)
plot(yy,    ((sqrt(ux(2:ly-1,lx).^2+uy(2:ly-1,lx).^2)))./ (max (sqrt(ux(2:ly-1,lx).^2+uy(2:ly-1,lx).^2))) );%(  sum((sqrt(ux(2:ly-1,lx).^2+uy(2:ly-1,lx).^2)))/(ly-2)  ) ); %
xlabel('y=0: y=L')
ylabel('U/umax @ x=L ')

xL=linspace (0,1,lx);
%figure(5);
subplot(3,2,3)
plot(xL, ((sqrt(ux(ly/2,1:lx).^2+uy(ly/2,1:lx).^2)))./(  sum((sqrt(ux(2:ly-1,1).^2+uy(2:ly-1,1).^2)))/(ly-2)  ) );%./max (sqrt(ux(a+1:3*a,lx).^2+uy(a+1:3*a,lx).^2)))   )
xlabel('X=0: X=L')
ylabel('U  @ Y/2 ')

Px=(rho(ly/2,1:lx))./3;
P1x=zeros(1,lx);
dp=zeros(1,lx);
for i=1:lx
P1x(1,i)=pIn+(((i-1)/(lx-1))*(pOut-pIn));
dp(1,i)=(Px(i)-P1x(i))/pOut;
end
xx=linspace (0,1,lx);
%figure(11);
subplot(3,2,4)
plot(xx,dp)
maxdp=max(dp);
if maxdp<=0
maxdp=1;
end
axis([0 1 0 (maxdp+.1*maxdp)])
xlabel('X=0: X=L')
ylabel(' dP ')

%figure(10);
subplot(3,2,5)
plot(xx,Px./pOut,xx,P1x./pOut,'r')
axis([0 1 pOut/pOut pIn/pOut])
``````

% subplot(3,2,5)
% plot(xx,P1x./pOut,‘r’)
xlabel(‘X=0: X=L’)
ylabel(’ Px ')

% Pbarx=(sum(rho(2:ly-1,1:lx)))./(3*(ly-2));
% xx=linspace (0,1,lx);
% figure(11);
% plot(xx,Pbarx./pOut)

``````end
``````

% Px=(rho(ly/2,1:lx))./3;
% Knx=zeros(1,lx);
% for i= 1:lx
% Knx(i)=(pOutKno)./Px(i);
% omega(i)=1./(.5+(((sqrt(6/pi))
((L))).*Knx(i)));
% end

``````u =(sqrt(ux.^2+uy.^2));
err=sum(sum(abs(uu-u)));
uu=u;
``````

end

Hi everyone,

Is this corner treatment applicable for slip velocity cases, too? I am trying to simulate a problem with a flexible vortex generator.

Regards,
Aysan