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

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