Outflow Boundary Conditions (Help please)

nkoosukuntla · March 15, 2011, 6:15pm

Hi All,

I am working on a Poiseulle flow with the following conditions

Uniform velocity Inlet (Zou He)
Zero Velocity gradient outlet (Outflow boundary)
Zou-He No Slip on Walls (With corner node implementation)

I am using a Incompressible BGK LBM. Parameters are as follows
Tau = 0.85
U = 0.05

When I use caorse grids (eg., 5x17) I see rapid increase in density and when I use fine grids (eg., 80x400) I see that I am getting closer to the analytical solution. How can we explain this phenomena? Please clarify. I am stuck here and needs some help.

Thanks,
Narender

stanabng2501 · March 16, 2011, 5:17pm

I am sorry i couldnot answer your question but i have a query regarding the same issue

can some1 t tell me is it necessary to add zero velocity zuo and he boundary condition at the wall alongwith bounceback conditions??
also please let me know wht is this corner boundary implementation…

nkoosukuntla · March 16, 2011, 7:26pm

Hi,

If you are using Bounce Back for no-slip that should be enough. Corner node implementation comes into picture only when you are using the Zou-He on all sides of the boundaries.

Eg: You use Zou-He inlet velocity of a channel and Zou-he no-slip on walls of a cannel then you need to use special treatment for corner nodes at inlet

Thanks,
Narender

stanabng2501 · March 17, 2011, 6:34am

so if i have understood correctly
then if you use pressure boundary condition at inlet and outlet and full bounceback conditions then i dun need any other boundary conditions isnit??

Thanks
Stanley
India

nkoosukuntla · March 17, 2011, 6:24pm

Corner node implementation is not needed if you use the following boundary conditions… But since you are using the Bounce back be aware that the no-slip you get is halfway between the last node and next last node…

X- Bounce back
I - Inlet Pressure
o- Outlet Pressure

xxxxxxxxxxxxxxxxxxxxxxxxxxx
I------------------------------------o
I------------------------------------o
I------------------------------------o
I------------------------------------o
I------------------------------------o
xxxxxxxxxxxxxxxxxxxxxxxxxxx

Corner Node implementation is needed if you use the following boundary conditions… You get the no slip exactly on the last node…

x- Zou-He No-slip
I - Inlet Pressure Zou-He
Q- Corner Node Zou-He
o- Outlet Zou-He

QxxxxxxxxxxxxxxxxxxxxxxxxxQ
I-------------------------------------o
I-------------------------------------o
I-------------------------------------o
I-------------------------------------o
I-------------------------------------o
QxxxxxxxxxxxxxxxxxxxxxxxxxQ

Hope I am clear…

Best,
Narender

stanabng2501 · March 18, 2011, 6:26am

yeh i think its clear now…
Now i understand why i was getting highest velocity at the boundary after using full bounceback condiions !!

Thanks for your effort
Stanley
India

mohammad3564 · October 15, 2011, 3:02pm

Hello Dear
Actually I have done my matlab code for Channel with Multi relaxation time but i think it has a problem with boundary condition and for example after 130000 iteration my code is crashed ,I was wondering If you could check my code and then say me what kind of boundary condition I should apply for Channel . I’m looking forward to hearing from you. Thank you for your consideration and time.

Yours faithfully

Mohammad Pourtousi

my matlab code:

clear;clc;

%%%% definitions %%%%%%%

lx=400;

ly=30;

hight1 =10;

hight2 =10;

lenght1=160;

lenght2=100;

lyInlet = hight1+1:ly;

lxInlet = 1:lenght1;

lxOutlet = lenght1+lenght2 : lx ;

lyOutlet = hight2 + 1 : ly ;

lyLeftCavity = 1:hight1;

lxBottomCavity = lenght1:lenght1+lenght2;

lyRightCavity = 1:hight2;

rho(1:lx+1,1:ly)=1;

Uo=0.05;

Re=100;

Nu = Uo*(ly-hight1)/Re;

tau = 3*Nu+0.5;

u=zeros(lx,ly);

v=zeros(lx,ly);

mEq=ones(9,lx,ly);

f=ones(9,lx,ly);

w=[4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36];

ex=[0,1,0,-1,0,1,-1,-1,1];

ey=[0,0,1,0,-1,1,1,-1,-1];

S=[1,1.4,1.4,1,1.2,1,1.2,1/tau,1/tau];

M=[1, 1, 1, 1, 1, 1, 1, 1, 1;-4,-1,-1,-1,-1, 2, 2, 2, 2;

4,-2,-2,-2,-2, 1, 1, 1, 1; 0, 1, 0,-1, 0, 1,-1,-1, 1;

0,-2, 0, 2, 0, 1,-1,-1, 1; 0, 0, 1, 0,-1, 1, 1,-1,-1;

0, 0,-2, 0, 2, 1, 1,-1,-1; 0, 1,-1, 1,-1, 0, 0, 0, 0;

0, 0, 0, 0, 0, 1,-1, 1,-1];

InvMS=M\diag(S);

for

t=1:200000

%%

%%%%%%% collision %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for

i=1:lx

for j=1:ly

jx=rho(i,j) * u(i,j);jy=rho(i,j) * v(i,j);

mEq(1,i,j)= rho(i,j);

mEq(2,i,j)= -2*rho(i,j) + 3 * (jx^2 + jy^2);

mEq(3,i,j)= rho(i,j) - 3 * (jx^2 + jy^2);

mEq(4,i,j)= jx;

mEq(5,i,j)= -jx;

mEq(6,i,j)= jy;

mEq(7,i,j)= -jy;

mEq(8,i,j)= jx^2 - jy^2;

mEq(9,i,j)= jx * jy;

end

m = reshape(M * reshape(f,9,lx*ly),9,lx,ly);

f = f - reshape(InvMS * reshape(m-mEq,9,lx*ly),9,lx,ly);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%%% streaming %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for

i=1:9

f(i,:, = circshift(f(i,:,:), [0,ex(i),ey(i)]);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%% Boundary Condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% Left Wall (bounce back)

rhoInput=f(1,1,lyInlet)+f(3,1,lyInlet)+f(5,1,lyInlet)+2*(f(4,1,lyInlet)+

…

f(7,1,lyInlet)+f(8,1,lyInlet))/(1-Uo);

f(2,1,lyInlet)=f(4,1,lyInlet) + 2/3 * rhoInput * Uo;

f(6,1,lyInlet)=f(8,1,lyInlet) + rhoInput * Uo / 6 ;

f(9,1,lyInlet)=f(7,1,lyInlet) + rhoInput * Uo / 6 ;

% %%%% Right Wall (bounce back)

% f(2,1,lyInlet)= f(2,lx,lyOutlet) ;

% f(6,1,lyInlet)= f(6,lx,lyOutlet) ;

% f(9,1,lyInlet)= f(9,lx,lyOutlet);

%

% f(2,lx+1,lyOutlet)= f(2,1,lyInlet) ;

% f(6,lx+1,lyOutlet)= f(6,1,lyInlet) ;

% f(9,lx+1,lyOutlet)= f(9,1,lyInlet);

f(2,lx,lyOutlet)= f(2,1,lyInlet) ;

f(6,lx,lyOutlet)= f(6,1,lyInlet) ;

f(9,lx,lyOutlet)= f(9,1,lyInlet);

% -------------------------------------------------------------------------

% for i=1:9

% rho(lx,lyOutlet)= rho(lx-1,lyOutlet);

% v(lx,lyOutlet)=0;

% u(lx,lyOutlet)=u(lx-1,lyOutlet);

% % Right boundary

% % f(i,lx,lyOutlet) = mEq(i,lx,lyOutlet) + …

% % 18*w(i)*ex(i)ey(i) ( f(6,lx,lyOutlet) - …

% % f(9,lx,lyOutlet)-mEq(6,lx,lyOutlet)+mEq(9,lx,lyOutlet) );

% end

%%%% Bottom Wall (bounce back)

f(3,:,1)=f(5,:,1);

f(6,:,1)=f(8,:,1);

f(7,:,1)=f(9,:,1);

%%%% Top Wall

f(5,:,ly)=f(3,:,ly);

f(9,:,ly)=f(7,:,ly);

f(8,:,ly)=f(6,:,ly);

%%%% Bottom First Obstacle wall

f(3,lxInlet,hight1)=f(5,lxInlet,hight1);

f(6,lxInlet,hight1)=f(8,lxInlet,hight1);

f(7,lxInlet,hight1)=f(9,lxInlet,hight1);

%%%% Bottom Second Obstacle wall

f(3,lxOutlet,hight2)=f(5,lxOutlet,hight2);

f(6,lxOutlet,hight2)=f(8,lxOutlet,hight2);

f(7,lxOutlet,hight2)=f(9,lxOutlet,hight2);

%%%% Left Cavity wall

f(2,lenght1,lyLeftCavity)=f(4,lenght1,lyLeftCavity);

f(6,lenght1,lyLeftCavity)=f(8,lenght1,lyLeftCavity);

f(9,lenght1,lyLeftCavity)=f(7,lenght1,lyLeftCavity);

%%%% Right Cavity Wall

f(4,lenght1+lenght2,lyRightCavity)=f(2,lenght1+lenght2,lyRightCavity);

f(8,lenght1+lenght2,lyRightCavity)=f(6,lenght1+lenght2,lyRightCavity);

f(7,lenght1+lenght2,lyRightCavity)=f(9,lenght1+lenght2,lyRightCavity);

%%%% Bottom Cavity Wall

f(3,lxBottomCavity,1)=f(5,lxBottomCavity,1);

f(6,lxBottomCavity,1)=f(8,lxBottomCavity,1);

f(7,lxBottomCavity,1)=f(9,lxBottomCavity,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

%%%%%%%%% rho and Velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%

rho = reshape(sum(f),lx,ly);

u(:, = reshape(ex * reshape(f,9,lx*ly),lx,ly)./rho;

v(:, = reshape(ey * reshape(f,9,lx*ly),lx,ly)./rho;

u(1,lyInlet) = Uo;u(lxInlet,lyLeftCavity)= 0;u(lxOutlet,lyRightCavity)= 0;

v(1,lyInlet) = 0 ;v(lxInlet,lyLeftCavity)= 0;v(lxOutlet,lyRightCavity)= 0;

%

% v(lx,lyOutlet)=0;

% u(lx,lyOutlet)=u(lx-1,lyOutlet);

figure(2);

Uaverage=sqrt(u(:,:).^2+v(:,:).^2);

imagesc(Uaverage(:,ly:-1:1)’./Uo);

drawnow

int2str(t)

end

pleb01 · October 16, 2011, 9:44am

Hi Narender,

Zou and He’s method does not explicitly conserve mass. It is possible to show that there is a small mass flux at the wall at each time step (rho_in is not exactly equal to rho_out at the wall). However, I think this “error” is proportional to dv/dy, so it shouldn’t have a big effect for Poiseuliie flow because this term should be zero. But if your code is predicting a (non-constant) vertical component of velocity to some order in grid-spacing then you could well see a mass change that decreases with mesh refinement (although I would have thought it a mass loss rather than a mass gain).

Perhaps there are some other things to consider. You say you are simulating incompressible flow in a channel, but the LBE is weakly compressible and only approximates the incompressible Navier-Stokes equations in the limit of small Mach number. If you use the He and Luo “incompressible” model then then error is an order of Ma better (for time independent flows) at approximating the incompressible Navier Stokes equations. However, the Zou and He method may need to be modified before being used with this model. To implement the wall velocity boundary condition you need to know the density, which is found from the known distributions at the wall. For example, at the south wall rho=f0+f1+f3+2(f4+f7+f7). So rho is not constant. Moreover, if we say rho=rho_0+deltarho, where rho_0 is constant and deltarho is the departure from the incompressible density then delta*rho is of the same order as in the “normal” LBE and not the incompressible LBE. These are just thoughts and I don’t know for sure, but you may want to check that your code is consistent.

I’d be very careful with the corners too. I’d argue that there are 5 unknowns here: if we take the south west corner as an example then, using standard notation, we don’t know f1, f2, f5, f6, f8 (f8 and f6 never enter the flow, but they are still needed to compute the moments). If you use Zou and He velocity boundary conditions at the inlet and the south wall, for example, then I’m not sure you have enough information to impose the conditions consistently (because the vertical boundary treats f1 as unknown but the horizontal boundary treats it as known - and similarly for other directions, and at the outlet…I’m also assuming that by uniform inlet you mean you’ve imposed the analytical solution here, and not u=constant everywhere, because this would cause and issue with the velocity at the boundary).

An alternative would be to impose a pressure gradient, i.e. impose an inflow density and a slightly smaller outflow density. I say this thinking about the “normal” LBE, but perhaps it could be generalised to the incompressible LBE if you play with the inlet/outlet pressure. Note also that since the flow is basically one dimensional then you could specialise your code to 1D. If your boundary conditions are consistent then your results for this flow should converge with dx^2.

Two other very quick points related to this thread. Firstly, half-way bounce back doesn’t give you zero velocity exactly half way between the nodes. We can show that the velocity at this place is proportional to dx^2. In other words, it is accurate to second order. Moreover, the exact position of the wall (where the velocity vanishes) is a function of the viscosity, which is a bit strange - it’s a numerical error. However, for applications where you need a fine mesh and high Re numbers (small relaxation time), this error is small.

Secondly, I haven’t looked at your code, Mohammad, but Zou and He’s method is less stable than bounce back so your code may crash for high Re flows.

Good luck!