# 2D Poiseuille flow error

Hi

I have implemented the LBM (D2Q9) for the 2D Poiseuille flow with periodic boundary conditions at the in-/out-let and no-slip at the boundaries.

The problem is that I can’t reproduce a parabolic velocity profile. The principle is very simple:

1. Find velocity and rho
2. Determine eq. distribution function
3. Perform collision. For nodes at the pipe boundary, they are reflected
4. Stream: Periodicity at in- and out-let and no streaming occurs at the pipe boundary

Can anyone give me a hint to what is wrong in the following?

``````
clear all

ny=7; nx=19;

c(1, :) = [+0 +0];
c(2, :) = [+1 +0];
c(3, :) = [+0 +1];
c(4, :) = [-1 +0];
c(5, :) = [+0 -1];
c(6, :) = [+1 +1];
c(7, :) = [-1 +1];
c(8, :) = [-1 -1];
c(9, :) = [+1 -1];

rho=ones(ny, nx);
solid_nodes=zeros(ny, nx);

f(:, :, 1) = (4./9. )*rho;
f(:, :, 2) = (1./9. )*rho;
f(:, :, 3) = (1./9. )*rho;
f(:, :, 4) = (1./9. )*rho;
f(:, :, 5) = (1./9. )*rho;
f(:, :, 6) = (1./36.)*rho;
f(:, :, 7) = (1./36.)*rho;
f(:, :, 8) = (1./36.)*rho;
f(:, :, 9) = (1./36.)*rho;

for i=1:nx
solid_nodes(1, i)  = +1;
solid_nodes(ny, i) = -1;
end

for time_step=1:1000
for j=1:ny
for i=1:nx
v_x(j, i) = 0.0;
v_y(j, i) = 0.0;
rho(j, i) = 0.0;
if ~solid_nodes(j, i)
for ctr=1:9
rho(j, i) = rho(j, i) + f(j, i, ctr);
v_x(j, i) = v_x(j, i) + c(ctr, 1)*f(j, i, ctr);
v_y(j, i) = v_y(j, i) + c(ctr, 2)*f(j, i, ctr);
end
v_x(j, i) = v_x(j, i)/rho(j, i);
v_y(j, i) = v_y(j, i)/rho(j, i);
end
end
end

coef = [3.0 9.0/2.0 3.0/2.0];
for j=1:ny
for i=1:nx
if ~solid_nodes(j, i)
v_x_eq = v_x(j, i) + pressure_grad;
v_y_eq = v_y(j, i);

v_x_squared = v_x_eq*v_x_eq;
v_y_squared = v_y_eq*v_y_eq;
v_squared   = v_x_squared + v_y_squared;

vx_vy_5 = v_x_eq + v_y_eq;
vx_vy_6 = -v_x_eq +  v_y_eq;
vx_vy_7 = -v_x_eq + -v_y_eq;
vx_vy_8 =  v_x_eq + -v_y_eq;

f_eq(j, i, 1) = rho(j, i)*(4.0/9.0 )*(1.0                                              - coef(3)*v_squared);
f_eq(j, i, 2) = rho(j, i)*(1.0/9.0 )*(1.0 + coef(1)*v_x_eq   + coef(2)*v_x_squared     - coef(3)*v_squared);
f_eq(j, i, 3) = rho(j, i)*(1.0/9.0 )*(1.0 + coef(1)*v_y_eq   + coef(2)*v_y_squared     - coef(3)*v_squared);
f_eq(j, i, 4) = rho(j, i)*(1.0/9.0 )*(1.0 - coef(1)*v_x_eq   + coef(2)*v_x_squared     - coef(3)*v_squared);
f_eq(j, i, 5) = rho(j, i)*(1.0/9.0 )*(1.0 - coef(1)*v_y_eq   + coef(2)*v_y_squared     - coef(3)*v_squared);
f_eq(j, i, 6) = rho(j, i)*(1.0/36.0)*(1.0 + coef(1)*vx_vy_5  + coef(2)*vx_vy_5*vx_vy_5 - coef(3)*v_squared);
f_eq(j, i, 7) = rho(j, i)*(1.0/36.0)*(1.0 + coef(1)*vx_vy_6  + coef(2)*vx_vy_6*vx_vy_6 - coef(3)*v_squared);
f_eq(j, i, 8) = rho(j, i)*(1.0/36.0)*(1.0 + coef(1)*vx_vy_7  + coef(2)*vx_vy_7*vx_vy_7 - coef(3)*v_squared);
f_eq(j, i, 9) = rho(j, i)*(1.0/36.0)*(1.0 + coef(1)*vx_vy_8  + coef(2)*vx_vy_8*vx_vy_8 - coef(3)*v_squared);
end
end
end

for j=1:ny
for i=1:nx
if solid_nodes(j, i)==-1
f(j, i, 3) = f(j, i, 5);
f(j, i, 6) = f(j, i, 8);
f(j, i, 7) = f(j, i, 9);

f(j, i, 5) = 0;
f(j, i, 8) = 0;
f(j, i, 9) = 0;

elseif solid_nodes(j, i)==+1
f(j, i, 5) = f(j, i, 3);
f(j, i, 8) = f(j, i, 6);
f(j, i, 9) = f(j, i, 7);

f(j, i, 3) = 0;
f(j, i, 6) = 0;
f(j, i, 7) = 0;
else
for ctr=1:9
f(j, i, ctr) = f(j, i, ctr) - (f(j, i, ctr) - f_eq(j, i, ctr))/tau;
end
end
end
end

for j=2:ny-1
j_lower = j-1;
j_upper = j+1;

for i=2:nx-1
i_lower = i-1;
i_upper = i+1;

f_temp(j,             i, 1) = f(j, i, 1);
f_temp(j,       i_upper, 2) = f(j, i, 2);
f_temp(j_upper,       i, 3) = f(j, i, 3);
f_temp(j,       i_lower, 4) = f(j, i, 4);
f_temp(j_lower,       i, 5) = f(j, i, 5);
f_temp(j_upper, i_upper, 6) = f(j, i, 6);
f_temp(j_upper, i_lower, 7) = f(j, i, 7);
f_temp(j_lower, i_lower, 8) = f(j, i, 8);
f_temp(j_lower, i_upper, 9) = f(j, i, 9);
end
end
f(:,:,:)=f_temp(:,:,:);
end

figure(1)
plot(v_x(:, 1))

``````

Hi niles_m,

As far as I am concerned, the streaming step is done for all of the nodes, including the wall nodes.

hope this concept helps your concern.

Puigar