# Strange boundary condition problem (poiseuille flow)

Hi there
I´m trying to simulate a simple 2d channel using 2dq9 bgk

I took as guide the code available on this site (cylinder.m) and rewrote it in a way I could understand it.

If I simulate the flow using the inlet boundary condition the poiseuille flow it works as expected., but If I use something different like a constant speed
It shows an strange profile (NORTH WEST) and (SOUTH EAST), I´ve being trying to debug the problem but I just cant get it working right.

If you want to see what I´m talking about just change

``````  vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
%vxin = 0.04;
``````

for
% vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
vxin = 0.04;

on initialaizing the fluid and when setting the inlet boundary condition

% 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.);

% 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;

``````  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);

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
for i=2:ly-1
L = ly-2;
y_phys = i-1.5;
vxin = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);
%vxin = 0.04;

%macro (FLUX)

u_x(i,1) = vxin;
u_y(i,1) = uy;

%micro
rho0 = ((f0(i,1) + f2(i,1) + f4(i,1)) + 2*(f3(i,1) + f7(i,1) + f6(i,1))) / (1 - vxin);
rho(i,1) = rho0;

ru = rho0 * vxin;

f1(i,1) = f3(i,1) + (2/3)*ru;
f5(i,1) = f7(i,1) + (1/6)*ru + (1/2)*(f2(i,1) - f4(i,1));
f8(i,1) = f6(i,1) + (1/6)*ru + (1/2)*(f4(i,1) - f2(i,1));

end

%OUTLET
for i=2:ly-1
%macro (PRESSURE)
rho(i,lx) = rhoOut;
u_y(i,lx) = uy;
%micro

ux0 = -1 + ((f0(i,lx) + f2(i,lx) + f4(i,lx))  + 2*(f1(i,lx) + f5(i,lx) + f8(i,lx)))/ rhoOut;

ru = rhoOut * ux0;

f3(i,lx) = f1(i,lx) - (2/3) * ru;
f7(i,lx) = f5(i,lx) - (1/6) * ru + (1/2)*(f2(i,lx) - f4(i,lx));
f6(i,lx) = f8(i,lx) - (1/6) * ru + (1/2)*(f4(i,lx) - f2(i,lx));

end

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

%Compute equilibrium functions
for i=1:ly;
for j=1:lx;

%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

% 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,2)==0)
u = sqrt(u_x.^2+u_y.^2);

imagesc(u);
%contour(u);

axis equal off; drawnow
end
``````

end

Hi there
I have modified the cylinder.m example and set a inlet speed inSpeedx of 0.04;
the same phenomenon is observed!
it looks like a boundary layer, but I dont think it is one, but rather problem with the inlet boundary.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cylinder.m: Channel flow past a cylinderical
% obstacle, using a LB method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample in Matlab
% Address: EPFL, 1015 Lausanne, Switzerland
% E-mail: jonas@lbmethod.org
% http://www.lbmethod.org/_media/numerics:cylinder.m
%
% Original implementaion of Zou/He boundary condition by
% Adriano Sciacovelli (see example “cavity.m”)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% You should have received a copy of the GNU General Public
% License along with this program; if not, write to the Free
% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
% Boston, MA 02110-1301, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% GENERAL FLOW CONSTANTS
lx = 400; % number of cells in x-direction
ly = 100; % number of cells in y-direction
obst_x = lx/5+1; % position of the cylinder; (exact
obst_y = ly/2+3; % y-symmetry is avoided)
obst_r = ly/20+1; % radius of the cylinder
uMax = 0.1; % maximum velocity of Poiseuille inflow
Re = 100; % Reynolds number
nu = uMax * 2.obst_r / Re; % kinematic viscosity
omega = 1. / (3
nu+1./2.); % relaxation parameter
%omega = 0.51;
maxT = 400000; % total number of iterations
tPlot = 5; % cycles

% 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];
col = [2:(ly-1)];
in = 1; % position of inlet
out = lx; % position of outlet

%init Speeds
inSpeedx = 0.04;
inLetCol = zeros(1, ly-2);
inLetCol = inLetCol + inSpeedx;

[y,x] = meshgrid(1:ly,1:lx); % get coordinate of matrix indices

obst = … % Location of cylinder
(x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;
obst(:,[1,ly]) = 1; % Location of top/bottom boundary
bbRegion = find(obst); % Boolean mask for bounce-back cells

% INITIAL CONDITION: Poiseuille profile at equilibrium
L = ly-2; y_phys = y-1.5;
%ux = 4 * uMax / (L*L) * (y_phys.L-y_phys.y_phys);
ux = zeros(lx,ly);
ux = ux + inSpeedx;
uy = zeros(lx,ly);
rho = 1;
for i=1:9
cu = 3
(cx(i)ux+cy(i)uy);
fIn(i,:, = rho .
t(i) .

( 1 + cu + 1/2
(cu.cu) - 3/2(ux.^2+uy.^2) );
end

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

``````% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux  = reshape ( (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
uy  = reshape ( (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;

% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS
% Inlet: Poiseuille profile
y_phys = col-1.5;
%ux(:,in,col) = 4 * uMax / (L*L) * (y_phys.*L-y_phys.*y_phys);

ux(:,in,col) = inLetCol;
uy(:,in,col) = 0;
rho(:,in,col) = 1 ./ (1-ux(:,in,col)) .* ( ...
sum(fIn([1,3,5],in,col)) + 2*sum(fIn([4,7,8],in,col)) );
% Outlet: Constant pressure
rho(:,out,col) = 1;
ux(:,out,col) = -1 + 1 ./ (rho(:,out,col)) .* ( ...
sum(fIn([1,3,5],out,col)) + 2*sum(fIn([2,6,9],out,col)) );
uy(:,out,col)  = 0;

% MICROSCOPIC BOUNDARY CONDITIONS: INLET (Zou/He BC)
fIn(2,in,col) = fIn(4,in,col) + 2/3*rho(:,in,col).*ux(:,in,col);
fIn(6,in,col) = fIn(8,in,col) + 1/2*(fIn(5,in,col)-fIn(3,in,col)) ...
+ 1/2*rho(:,in,col).*uy(:,in,col) ...
+ 1/6*rho(:,in,col).*ux(:,in,col);
fIn(9,in,col) = fIn(7,in,col) + 1/2*(fIn(3,in,col)-fIn(5,in,col)) ...
- 1/2*rho(:,in,col).*uy(:,in,col) ...
+ 1/6*rho(:,in,col).*ux(:,in,col);

% MICROSCOPIC BOUNDARY CONDITIONS: OUTLET (Zou/He BC)
fIn(4,out,col) = fIn(2,out,col) - 2/3*rho(:,out,col).*ux(:,out,col);
fIn(8,out,col) = fIn(6,out,col) + 1/2*(fIn(3,out,col)-fIn(5,out,col)) ...
- 1/2*rho(:,out,col).*uy(:,out,col) ...
- 1/6*rho(:,out,col).*ux(:,out,col);
fIn(7,out,col) = fIn(9,out,col) + 1/2*(fIn(5,out,col)-fIn(3,out,col)) ...
+ 1/2*rho(:,out,col).*uy(:,out,col) ...
- 1/6*rho(:,out,col).*ux(:,out,col);

% COLLISION STEP
for i=1:9
cu = 3*(cx(i)*ux+cy(i)*uy);
fEq(i,:,:)  = rho .* t(i) .* ...
( 1 + cu + 1/2*(cu.*cu)  - 3/2*(ux.^2+uy.^2) );
fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:));
end

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

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

% VISUALIZATION
if (mod(cycle,tPlot)==1)
u = reshape(sqrt(ux.^2+uy.^2),lx,ly);
%p = reshape(rho/3,lx,ly);
u(bbRegion) = nan;
imagesc(u');
%p(bbRegion) = nan;
%imagesc(p');
axis equal off; drawnow
end
``````

end

Hi
can inlet boundary conditon be changed from velocity to pressure?, at the same way, can initial condition be changed the pressure ?.

ux = 4 * uMax / (L*L) * (y_phys.*L-y_phys.y_phys) is how to be understand. I canot understand means "ux = 4 * uMax / (LL) * (y_phys.*L-y_phys.*y_phys)". can initial condition and boundary condition that is pressure and pressure different be for example. thank a lot in advance.

ydt

that formula generates de poiseuille speed profile at the inlet. not an linear speed lets say 0,04. it generates a parabolic profile. that depends on the width of the channel.

I just want to mention that if you implement the corners as stated by GoncaolSilva, it works
it matches the analytical solution.