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:
- Find velocity and rho
- Determine eq. distribution function
- Perform collision. For nodes at the pipe boundary, they are reflected
- 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;
pressure_grad=0.001; tau=1;
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))