Lid Driven Cavity, Boundary Conditions

Lattice Boltzmann Method LBM: Theory
#1

Hi every one
I have written a program for solving Lid Driven Cavity for Re=1000 in MATLAB. But unfortunately I can get a correct result.
I think I have problem in North boundary condition! when I run the program, I face lots of NAN values. I don’t know why!

I will be so thankful if anybody can help me.
Here is my code:


% 2D Lattice Boltzmann model for Lid-Driven Cavity.
%  c4  c3   c2
%    \  |  /
%  c5 -c9 - c1
%    /  |  \
%  c6  c7   c8

tic; hold on;
clear, clc;
nx=100; ny=100; tstep=400; alpha=0.01; omega=1/(3*alpha+0.5);
u_ini=0.1; v_ini=0; Re=u_ini*nx/alpha;
w1=4/9; w2=1/9; w3=1/36; f=zeros(nx,ny,9); f_eq=f; density=5;


u=(sum(f(:,:,[1 2 8]),3)-sum(f(:,:,[4 5 6]),3))./density;
v=(sum(f(:,:,[2 3 4]),3)-sum(f(:,:,[6 7 8]),3))./density;



for ii=1:tstep
        
            
    f_eq(:,:,1)=w2*density.*(1+3*u+9/2*u.^2-3/2*(u.^2+v.^2));
    f_eq(:,:,3)=w2*density.*(1+3*v+9/2*v.^2-3/2*(u.^2+v.^2));
    f_eq(:,:,5)=w2*density.*(1-3*u+9/2*u.^2-3/2*(u.^2+v.^2));
    f_eq(:,:,7)=w2*density.*(1-3*v+9/2*v.^2-3/2*(u.^2+v.^2));
    
    f_eq(:,:,2)=w3*density.*(1+3*(u+v)+9/2*(u+v).^2-3/2*(u.^2+v.^2));
    f_eq(:,:,4)=w3*density.*(1+3*(-u+v)+9/2*(-u+v).^2-3/2*(u.^2+v.^2));
    f_eq(:,:,6)=w3*density.*(1-3*(u+v)+9/2*(u+v).^2-3/2*(u.^2+v.^2));
    f_eq(:,:,8)=w3*density.*(1+3*(u-v)+9/2*(u-v).^2-3/2*(u.^2+v.^2));
    
    f_eq(:,:,9)=w1*density.*(1-3/2*(u.^2+v.^2));
    
    f=omega*f_eq+(1-omega)*f;
    
    % Propegate (This part of code [propegate] is always constant for all LBM
    % problems.)
    f(:,:,4)=f([2:nx 1],[ny 1:ny-1],4); f(:,:,3)=f(:,[ny 1:ny-1],3);
    f(:,:,2)=f([nx 1:nx-1],[ny 1:ny-1],2); f(:,:,5)=f([2:nx 1],:,5);
    f(:,:,1)=f([nx 1:nx-1],:,1); f(:,:,6)=f([2:nx 1],[2:ny 1],6);
    f(:,:,7)=f(:,[2:ny 1],7); f(:,:,8)=f([nx 1:nx-1],[2:ny 1],8);
    
    % Boundary Conditions

    %At i=1, Bounceback
    f(1,:,1)=f(1,:,5); f(1,:,2)=f(1,:,6); f(1,:,8)=f(1,:,4);

    %At i=nx, Bounceback
    f(nx,:,4)=f(nx,:,8); f(nx,:,5)=f(nx,:,1); f(nx,:,6)=f(nx,:,2);

    %At j=1, Bounceback
    f(:,1,2)=f(:,1,6); f(:,1,3)=f(:,1,7); f(:,1,4)=f(:,1,8);

    %At j=ny, Know Velocity
    densityN=f(:,ny,9)+f(:,ny,1)+f(:,ny,5)+2*(f(:,ny,3)+f(:,ny,4)+f(:,ny,2));
    f(:,ny,7)=f(:,ny,3);
    f(:,ny,6)=f(:,ny,2)+0.5*(f(:,ny,1)-f(:,ny,5))-u_ini.*densityN/2;
    f(:,ny,8)=f(:,ny,4)+0.5*(f(:,ny,5)-f(:,ny,1))+u_ini.*densityN/2;
    
    
    density=sum(f,3);
    u=(sum(f(:,:,[1 2 8]),3)-sum(f(:,:,[4 5 6]),3))./density;
    v=(sum(f(:,:,[2 3 4]),3)-sum(f(:,:,[6 7 8]),3))./density;
            
end

contour(u');
toc;
#2

Hi Emad,

Quickly scanning through your code, I see you initiate f at zero everywhere. This might cause your problem, since density, as the sum of all f-values, is used as a dividing factor later on in your code. See what happens if you initiate f to f_eq. You can use zero u,v velocity and initial density to calculate this, but at least the nodes will be populated.

Also, try running your code step by step for each time step. This way you can see where the NaN values first appear and point you in the right direction.

Hope this helps.

Elrohir

-Edit

One more thing I just noticed. Even if your code would be correct, the combination of Zou-He and bounceback conditions is not very stable. From my own experience with LBM on LDC simulations, this will blow up the solution at higher values of omega. See what happens if you lower omega to, say, 1.5. Either lower your Reynolds number or use more nodes to discretise the domain.

#3

Dear Elrohir

first of all, thank you very much for your suggestions.

I have changed the omega to 1.5 and I could get much better result. but still I could not get the correct figure.
About omega, I want to know if there is any formula for calculating it.

#4

Thanks for your code.

#5

i am working on a project fdlbm(finite difference latiice boltezmam )
and i need a code for lid cavity
and i need it very much
best regards
mahmoud