Poiseulle flow with Velocity and Outflow Boundary Conditions

Lattice Boltzmann Method LBM: Theory
#1

Hi All,

I am trying to code a poiseulle flow with Zou-He Velocity boundary conditions at the inlet, NBC (Outflow) boundary conditions at the outlet and Zou-He No-Slip on the walls. Did anybody do something similar? My rho keeps on increasing with time. Here is my code at inlet and outlet.


% Inlet Velocity Boundary Condition

if (x==1 && y==1)           % Bottom Left Corner                 

    f(x,y,1) = f(x,y,3);
    f(x,y,2) = f(x,y,4);
    f(x,y,5) = f(x,y,7);
    f(x,y,6) = 0.5*(rho(2,2) - ( f(x,y,9) + 2*(f(x,y,3)+f(x,y,7)+f(x,y,4)) ) );
    f(x,y,8) = 0.5*(rho(2,2)- ( f(x,y,9) + 2*(f(x,y,3)+f(x,y,7)+f(x,y,4)) ) );

elseif (x==1 && y==YMAX)    % Top Left Corner

    f(x,y,1) = f(x,y,3);
    f(x,y,4) = f(x,y,2);
    f(x,y,8) = f(x,y,6);
    f(x,y,5) = 0.5*(rho(2,YMAX-1) - ( f(x,y,9) + 2*(f(x,y,2)+f(x,y,6)+f(x,y,3)) ) );
    f(x,y,7) = 0.5*(rho(2,YMAX-1) - ( f(x,y,9) + 2*(f(x,y,2)+f(x,y,6)+f(x,y,3)) ) );

else

    u(x,y)=u_in;
    v(x,y)=0;
    rho(x,y) = rho0*u(x,y) + (f(x,y,2)+f(x,y,4)+f(x,y,9)+2*(f(x,y,3)+f(x,y,6)+f(x,y,7)));
    f(x,y,1)=f(x,y,3)+2/3*rho0*u(x,y);
    f(x,y,5)=f(x,y,7)+1/6*rho0*u(x,y)+1/2*(f(x,y,4)-f(x,y,2))+(1/2)*rho0*v(x,y);
    f(x,y,8)=f(x,y,6)+1/6*rho0*u(x,y)+1/2*(f(x,y,2)-f(x,y,4))-(1/2)*rho0*v(x,y);

end


% Outlet NBC Outflow Boundary Condition

  % Velocity at outlet (Outflow condition. NBC)

  uvel = u(x-2,y);

  % Imposing this velocity as Dirichlet Boundary condition


if (x==XMAX && y==1)        % Bottom right corner
    f(x,y,3) = f(x,y,1);
    f(x,y,2) = f(x,y,4);
    f(x,y,6) = f(x,y,8);
    f(x,y,5) = 0.5*(rho(x-1,y) - ( f(x,y,9) + 2*(f(x,y,4)+f(x,y,8)+f(x,y,1)) ) );
    f(x,y,7) = 0.5*(rho(x-1,y) - ( f(x,y,9) + 2*(f(x,y,4)+f(x,y,8)+f(x,y,1)) ) );
elseif (x==XMAX && y==YMAX) %Top right corner
    f(x,y,3) = f(x,y,1);
    f(x,y,4) = f(x,y,2);
    f(x,y,7) = f(x,y,5);
    f(x,y,6) = 0.5*(rho(x-1,y) - ( f(x,y,9) + 2*(f(x,y,1)+f(x,y,5)+f(x,y,2)) ) );
    f(x,y,8) = 0.5*(rho(x-1,y) - ( f(x,y,9) + 2*(f(x,y,1)+f(x,y,5)+f(x,y,2)) ) );
else

    u(x,y)=uvel;
    v(x,y)=0;
    rho(x,y) = ((f(x,y,9)+f(x,y,2)+f(x,y,4))+2*(f(x,y,5)+f(x,y,1)+f(x,y,8))) - u(x,y)*rho0;
    f(x,y,3)=f(x,y,1)-(2/3)*rho0*u(x,y);
    f(x,y,7)=f(x,y,5)-(1/6)*rho0*u(x,y)-(1/2)*rho0*v(x,y)+(1/2)*(f(x,y,2)-f(x,y,4));
    f(x,y,6)=f(x,y,8)-(1/6)*rho0*u(x,y)+(1/2)*rho0*v(x,y)+(1/2)*(f(x,y,4)-f(x,y,2));

end

I don’t know where I am going wrong or am I missing something simple? Please correct me.

Thanks,
Narender