[Aldou]

Hi, I am newcomer in CFD. I try to write script for channel flow using Lattice Bolzmann method in MATLAB. My problem is I cannot get the right figure of stream lines (picture enclosed). I keep getting strange behaviour at the inlet and dont know why. As BC on the walls I use Zou-He non-slip condition. As inlet BC I use velocity condition (not fully developed velocity profile...only parametr I know is inlet velocity in x direction along the inlet, y direction is supposed to be 0). As outlet BC I use zero velocity gradient. If someone knew why is it so please give me an advice.

D2Q9 scheme
6 2 5
3 9 1
7 4 8

(500=M x 20=N) lattice nodes
.____________________ outlet
.
.
.
. <- inlet + corner nodes
.
.
.____________________ outlet


just to be precise in the distribution function: f("row","column","number of dist.fcn.")

%north bounce back (Zou-He non-slip)

 
f(N,2:M,4) = f(N,2:M,2);
f(N,2:M,7) = f(N,2:M,5) + 0.5*(f(N,2:M,1) - f(N,2:M,3));
f(N,2:M,8) = f(N,2:M,6) - 0.5*(f(N,2:M,1) - f(N,2:M,3));

%south bounce back (Zou-He non-slip)

f(1,2:M,2) = f(1,2:M,4);
f(1,2:M,5) = f(1,2:M,7) - 0.5*(f(1,2:M,1) - f(1,2:M,3));
f(1,2:M,6) = f(1,2:M,8) + 0.5*(f(1,2:M,1) - f(1,2:M,3));

% inlet velocity BC with corner node treatment

Ro_in(:) = (f(2:N-1,1,9)+f(2:N-1,1,2)+f(2:N-1,1,4)+2*(f(2:N-1,1,3)+f(2:N-1,1,6)+f(2:N-1,1,7)))/(1-u_lb);
f(2:N-1,1,1) = f(2:N-1,1,3) + (2*Ro_in(:)*u_lb)/3;
f(2:N-1,1,5) = f(2:N-1,1,7) - 0.5*(f(2:N-1,1,2) - f(2:N-1,1,4)) + (Ro_in(:)*u_lb)/6;
f(2:N-1,1,8) = f(2:N-1,1,6) + 0.5*(f(2:N-1,1,2) - f(2:N-1,1,4)) + (Ro_in(:)*u_lb)/6;

%left bottom corner

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

%left top corner

f(N,1,1) = f(N,1,3);
f(N,1,4) = f(N,1,2);
f(N,1,8) = f(N,1,6);
f(N,1,5) = 0.5*(Ro_in(length(Ro_in)) - (f(N,1,9)+f(N,1,1)+f(N,1,2)+f(N,1,3)+f(N,1,4)+f(N,1,6)+f(N,1,8)));
f(N,1,7) = f(N,1,5);

%outlet zero velocity gradient

f(:,M,3) = f(:,M-1,3);
f(:,M,6) = f(:,M-1,6);
f(:,M,7) = f(:,M-1,7);

...then i compute velocity and density from every lattice node, define new distribution functions and repeat cycle.

I hope my problem is cleared and thanks for helping.