About lid-driven cavity flow

gn02317298 · November 28, 2013, 10:03pm

Hi everyone,
I am a beginner in LBM and study the book “Lattice Boltzmann Method” by A.A. Mohamad.
For the lid driven cavity problem in p.81 , there is a fortran code in appendix (p.133~p.139).
I tried to rewrite it using matlab but can’t get a correct result (the density increase with iterations).
Can anyone help me to debug the code ? I spent threes days on it but still can’t identity the problem.
Thanks.
Here is my matlab code:


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LBM computer code for lid-driven cavity, D2Q9
% change the index 0 to index 9
%          6   2   5
%           \  |  /
%         3 -- 9 -- 1
%           /  |  \
%          7   4   8
clc;clear;
n=10;m=10;
uo=0.10; % velocity of lid for LBM (not for reality)
dx=1.0;
dy=dx;
dt=1.0;
alpha=0.01;
Re=uo*m/alpha; % What we are concerned is the Reynolds number Re=1000
omega=1.0/(3.*alpha+0.5); % derive from eq. 3.82
%mstep=40000;
mstep=100;
w(1:4)=1/9;w(5:8)=1/36;w(9)=4/9;
% fortran index start from 0, I change it to 9 in matlab code
cx=[1 0 -1  0 1 -1 -1  1 0 ]; % cx(1)~cx(9)
cy=[0 1  0 -1 1  1 -1 -1 0 ]; % cy(1)~cy(9)
rho=ones(m,n); % set the initial value to the value of rhoo
u=zeros(m,n);  
v=zeros(m,n);
% set up the initial velocity at top boundary (lid velocity)
for i=1:n
    u(i,m)=uo; % lid moves to the right direction, here we only set up the interior region, why?
    v(i,m)=0;
end
% initail values set up
f1=zeros(n,m);
f2=zeros(n,m);
f3=zeros(n,m);
f4=zeros(n,m);
f5=zeros(n,m);
f6=zeros(n,m);
f7=zeros(n,m);
f8=zeros(n,m);
f9=zeros(n,m);
% main loop
for kk=1:mstep
    % step 1: collision 
    for j=1:m
        for i=1:n
            t1(i,j)=u(i,j)*u(i,j)+v(i,j)*v(i,j); 
            t2_1(i,j)=u(i,j)*cx(1)+v(i,j)*cy(1);
            t2_2(i,j)=u(i,j)*cx(2)+v(i,j)*cy(2);
            t2_3(i,j)=u(i,j)*cx(3)+v(i,j)*cy(3);
            t2_4(i,j)=u(i,j)*cx(4)+v(i,j)*cy(4);
            t2_5(i,j)=u(i,j)*cx(5)+v(i,j)*cy(5);
            t2_6(i,j)=u(i,j)*cx(6)+v(i,j)*cy(6);
            t2_7(i,j)=u(i,j)*cx(7)+v(i,j)*cy(7);
            t2_8(i,j)=u(i,j)*cx(8)+v(i,j)*cy(8);
            t2_9(i,j)=u(i,j)*cx(9)+v(i,j)*cy(9);
            feq_1(i,j)=rho(i,j)*w(1)*(1.0+3.0*t2_1(i,j)+4.50*t2_1(i,j)*t2_1(i,j)-1.50*t1(i,j));
            feq_2(i,j)=rho(i,j)*w(2)*(1.0+3.0*t2_2(i,j)+4.50*t2_2(i,j)*t2_2(i,j)-1.50*t1(i,j));
            feq_3(i,j)=rho(i,j)*w(3)*(1.0+3.0*t2_3(i,j)+4.50*t2_3(i,j)*t2_3(i,j)-1.50*t1(i,j));
            feq_4(i,j)=rho(i,j)*w(4)*(1.0+3.0*t2_4(i,j)+4.50*t2_4(i,j)*t2_4(i,j)-1.50*t1(i,j));
            feq_5(i,j)=rho(i,j)*w(5)*(1.0+3.0*t2_5(i,j)+4.50*t2_5(i,j)*t2_5(i,j)-1.50*t1(i,j));
            feq_6(i,j)=rho(i,j)*w(6)*(1.0+3.0*t2_6(i,j)+4.50*t2_6(i,j)*t2_6(i,j)-1.50*t1(i,j));
            feq_7(i,j)=rho(i,j)*w(7)*(1.0+3.0*t2_7(i,j)+4.50*t2_7(i,j)*t2_7(i,j)-1.50*t1(i,j));
            feq_8(i,j)=rho(i,j)*w(8)*(1.0+3.0*t2_8(i,j)+4.50*t2_8(i,j)*t2_8(i,j)-1.50*t1(i,j));
            feq_9(i,j)=rho(i,j)*w(9)*(1.0+3.0*t2_9(i,j)+4.50*t2_9(i,j)*t2_9(i,j)-1.50*t1(i,j));
            f1(i,j)=omega*feq_1(i,j)+(1.-omega)*f1(i,j);
            f2(i,j)=omega*feq_2(i,j)+(1.-omega)*f2(i,j);
            f3(i,j)=omega*feq_3(i,j)+(1.-omega)*f3(i,j);
            f4(i,j)=omega*feq_4(i,j)+(1.-omega)*f4(i,j);
            f5(i,j)=omega*feq_5(i,j)+(1.-omega)*f5(i,j);
            f6(i,j)=omega*feq_6(i,j)+(1.-omega)*f6(i,j);
            f7(i,j)=omega*feq_7(i,j)+(1.-omega)*f7(i,j);
            f8(i,j)=omega*feq_8(i,j)+(1.-omega)*f8(i,j);
            f9(i,j)=omega*feq_9(i,j)+(1.-omega)*f9(i,j);
        end
    end
    % step 2: streaming
    for j=1:m
        for i=1:n-1
            f1(n-i+1,j)=f1(n-i,j);       % f1 basic, from right to left. no f1(1,j)
            f3(i,j)=f3(i+1,j);           % f3 basic, from left to right. no f3(n,j)                 
        end
    end
    for j=1:m-1
        for i=1:n
            f2(i,m-j+1)=f2(i,m-j);       % f2 basic, from top to bottom. no f2(i,1)
            f4(i,j)=f4(i,j+1);           % f4 basic,from bottom to top, no f4(i,j=m)
        end
    end   
    for j=1:m-1
        for i=1:n-1
            f5(n-i+1,m-j+1)=f5(n-i,m-j); % f5 related to f1,f2,from top to bottom, no f5(i,j=1),f5(i=1,j)
            f6(i,m-j+1)=f6(i+1,m-j);     % f6 related to f3,f2,from top to bottom, no f6(i=n,j),f6(i,j=1)
            f7(i,j)=f7(i+1,j+1);         % f7 related to f3,f4,from bottom to top, no f7(i=n,j),f7(i,j=m)
            f8(n-i+1,j)=f8(n-i,j+1);     % f8 related to f1,f4,from bottom to top, no f8(i=1,j),f8(i,j=m)
        end
    end
    % step 3: calculate the distribution functions at surfaces (boundaries)
    % see fig. 5.6 % (1) bounce back on west boundary 
    for j=1:m
        f1(1,j)=f3(1,j);
        f5(1,j)=f7(1,j);
        f8(1,j)=f6(1,j);
    %(2) bounce back on east boundary
        f3(n,j)=f1(n,j);
        f6(n,j)=f8(n,j);
        f7(n,j)=f5(n,j);    
    end
    %(3)bounce back on south boundary
    for i=1:n
        f2(i,1)=f4(i,1);
        f5(i,1)=f7(i,1);
        f6(i,1)=f8(i,1);
    end
    %(4)moving lid , north bundary (boundary with known velocity)
    % Apply the eqns 5.29 ~ 5.32
    %for i=2:n-1 %( only the interior grids )
    for i=1:n
        rhon(i)=f9(i,m)+f1(i,m)+f3(i,m)+2*(f2(i,m)+f6(i,m)+f5(i,m));
        f4(i,m)=f2(i,m);
        f7(i,m)=f5(i,m)+0.5*(f1(i,m)-f3(i,m))-0.5*rhon(i)*uo;
        f8(i,m)=f6(i,m)+0.5*(f3(i,m)-f1(i,m))+0.5*rhon(i)*uo; 
    end
    % step 4: calculate density rho and velocity components u,v
    % 1. calculate updated density rho(i,j) for each grid
    rho=f1+f2+f3+f4+f5+f6+f7+f8+f9;
    for i=1:n
        rho(i,m)=f9(i,m)+f1(i,m)+f3(i,m)+2*(f2(i,m)+f6(i,m)+f5(i,m));
    end
    % 2. calculate updated velocity components u,v
    u=(f1*cx(1)+f2*cx(2)+f3*cx(3)+f4*cx(4)+f5*cx(5)+f6*cx(6)+f7*cx(7)+f8*cx(8)+f9*cx(9))./rho;
    v=(f1*cy(1)+f2*cy(2)+f3*cy(3)+f4*cy(4)+f5*cy(5)+f6*cy(6)+f7*cy(7)+f8*cy(8)+f9*cy(9))./rho;
end

Philippe · November 29, 2013, 9:35am

The major issue in your code is probably that you are trying to resolve cavity flow at Re=1000 on a 10x10 lattice. Try for example Re=50 at 30x30, this should work. The implementation looks ok to me but I have not thoroughly checked each line. Apart from that … u0 could be smaller (~0.02) and you should check your units.

Another remark is on matlab: Try to get rid of the for loops, they are usually much slower than matlab’s matrix operations. for example


for i=1:n
        f2(i,1)=f4(i,1);
        f5(i,1)=f7(i,1);
        f6(i,1)=f8(i,1);
end

could be written as


f2(:,1)=f4(:,1);
f5(:,1)=f7(:,1);
f6(:,1)=f8(:,1);

gn02317298 · November 29, 2013, 10:44pm

Hi Philippe,

Thank you for the response. I have tried smaller u0 (=0.01) and used 30 x 30 lattice (Re=30) but still encountered the same problem.
If I increased the number of iteration to 100, the density became 1.13x10^7. Although I tried iterations=1, I got a denstiy rho = ~1.88 ,which is
much larger than my setting (rho=1). I carefully checked my code but still failed to find the problem.

Min

Philippe · December 3, 2013, 9:17am

Where do you get these big densities? At the inlet or somewhere inside the domain?

gn02317298 · December 3, 2013, 4:14pm

These big densities distributed on the whole domain. I use surf(rho) command in matlab to see the distribution of density.

Min.

Emad · January 27, 2014, 7:32pm

Hi,
I have the same problem. First of all, the Fortran codes for lid-driven cavity and other problems in this book have lots of errors.
Second, the major source of error in this book comes from omega relationship. At this moment, try omega=1.
Although changing omega, I still have problem and can not get the correct figure.
Check this page (http://www.palabos.org/forum/read.php?3,6497) to see my codes.

aaguirr2 · February 1, 2014, 4:30pm

Hello gn02317298

I did a not very detailed review of your code, but I noticed that you’re using Zou and He velocity condition for the lid and bounce back condition for the walls. None of this boundary conditions is fixing the density on the walls, so I would think that the density inside domain will be growing ad infinitum, although results of the velocity field should be fine.

However, checking the cavity example by Jonas Latt, which is also written for matlab, I notice that the inner density is actually converging. I don’t know, maybe I have to check my own code too.

Regards,

A. Aguirre.

gn02317298 · February 22, 2014, 5:04am

Hi Emad and aaguirr2,

Thanks for your response.
Eventually I rewrote the “streaming process”. In my original matlab code, I tried to combine some streaming process in one for-loop but now
each streaming has their own for-loop (exactly the same structure in the book). Although I don’t think there’s any difference, the results turn out to be quite different. Density can converge to 1.83~2.0. I don’t kown why it didn’t stay at one. Any suggestions ?
I would attach my code at the end.
Now, I am practicing the example " flow in a channel ". In the transient process, the density would gradually decrease from inlet to outlet.
Is there any material you can suggest me so I can change the setting of boundary conditions to make the density stay at a constant value ?


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LBM computer code for lid-driven cavity, D2Q9
% change the index 0 to index 9
%          6   2   5
%           \  |  /
%         3 -- 9 -- 1
%           /  |  \
%          7   4   8
clc;clear;
n=40;m=40;
uo=0.10; % velocity of lid for LBM (not for reality)
rhoo=1;
dx=1.0;
dy=dx;
dt=1.0;
alpha=0.01;
Re=uo*m/alpha; % What we are concerned is the Reynolds number Re=1000
omega=1.0/(3.*alpha+0.5); % derive from eq. 3.82
mstep=1000;
w(1:4)=1/9;w(5:8)=1/36;w(9)=4/9;
% fortran index start from 0, I change it to 9 in matlab code
cx=[1 0 -1  0 1 -1 -1  1 0 ]; % cx(1)~cx(9)
cy=[0 1  0 -1 1  1 -1 -1 0 ]; % cy(1)~cy(9)
rho=ones(m,n); % set the initial value to the value of rhoo
u=zeros(m,n);
v=zeros(m,n);
% set up the initial velocity at top boundary (lid velocity)
for i=1:n
    u(i,m)=uo; % lid moves to the right direction, here we only set up the interior region, why?
    v(i,m)=0;
end
% initail values set up
for j=1:m
    for i=1:n
        f1(i,j)=0;
        f2(i,j)=0;
        f3(i,j)=0;
        f4(i,j)=0;
        f5(i,j)=0;
        f6(i,j)=0;
        f7(i,j)=0;
        f8(i,j)=0;
        f9(i,j)=0;
    end
end
% main loop
for kk=1:mstep
    % step 1: collision 
    for i=1:n
        for j=1:m
            t1(i,j)=u(i,j)*u(i,j)+v(i,j)*v(i,j);
            t2_1(i,j)=u(i,j)*cx(1)+v(i,j)*cy(1);
            t2_2(i,j)=u(i,j)*cx(2)+v(i,j)*cy(2);
            t2_3(i,j)=u(i,j)*cx(3)+v(i,j)*cy(3);
            t2_4(i,j)=u(i,j)*cx(4)+v(i,j)*cy(4);
            t2_5(i,j)=u(i,j)*cx(5)+v(i,j)*cy(5);
            t2_6(i,j)=u(i,j)*cx(6)+v(i,j)*cy(6);
            t2_7(i,j)=u(i,j)*cx(7)+v(i,j)*cy(7);
            t2_8(i,j)=u(i,j)*cx(8)+v(i,j)*cy(8);
            t2_9(i,j)=u(i,j)*cx(9)+v(i,j)*cy(9);
            feq_1(i,j)=rho(i,j)*w(1)*(1.0+3.0*t2_1(i,j)+4.50*t2_1(i,j)*t2_1(i,j)-1.50*t1(i,j));
            feq_2(i,j)=rho(i,j)*w(2)*(1.0+3.0*t2_2(i,j)+4.50*t2_2(i,j)*t2_2(i,j)-1.50*t1(i,j));
            feq_3(i,j)=rho(i,j)*w(3)*(1.0+3.0*t2_3(i,j)+4.50*t2_3(i,j)*t2_3(i,j)-1.50*t1(i,j));
            feq_4(i,j)=rho(i,j)*w(4)*(1.0+3.0*t2_4(i,j)+4.50*t2_4(i,j)*t2_4(i,j)-1.50*t1(i,j));
            feq_5(i,j)=rho(i,j)*w(5)*(1.0+3.0*t2_5(i,j)+4.50*t2_5(i,j)*t2_5(i,j)-1.50*t1(i,j));
            feq_6(i,j)=rho(i,j)*w(6)*(1.0+3.0*t2_6(i,j)+4.50*t2_6(i,j)*t2_6(i,j)-1.50*t1(i,j));
            feq_7(i,j)=rho(i,j)*w(7)*(1.0+3.0*t2_7(i,j)+4.50*t2_7(i,j)*t2_7(i,j)-1.50*t1(i,j));
            feq_8(i,j)=rho(i,j)*w(8)*(1.0+3.0*t2_8(i,j)+4.50*t2_8(i,j)*t2_8(i,j)-1.50*t1(i,j));
            feq_9(i,j)=rho(i,j)*w(9)*(1.0+3.0*t2_9(i,j)+4.50*t2_9(i,j)*t2_9(i,j)-1.50*t1(i,j));
            
            f1(i,j)=omega*feq_1(i,j)+(1.-omega)*f1(i,j);
            f2(i,j)=omega*feq_2(i,j)+(1.-omega)*f2(i,j);
            f3(i,j)=omega*feq_3(i,j)+(1.-omega)*f3(i,j);
            f4(i,j)=omega*feq_4(i,j)+(1.-omega)*f4(i,j);
            f5(i,j)=omega*feq_5(i,j)+(1.-omega)*f5(i,j);
            f6(i,j)=omega*feq_6(i,j)+(1.-omega)*f6(i,j);
            f7(i,j)=omega*feq_7(i,j)+(1.-omega)*f7(i,j);
            f8(i,j)=omega*feq_8(i,j)+(1.-omega)*f8(i,j);
            f9(i,j)=omega*feq_9(i,j)+(1.-omega)*f9(i,j);
       
        end
    end
    % step 2: streaming
    for j=1:m
        for i=n:-1:2
            f1(i,j)=f1(i-1,j);       % f1 basic, from right to left. no f1(1,j)                  
        end
    end
    
    for j=1:m
        for i=1:n-1
            f3(i,j)=f3(i+1,j);       % f3 basic, from left to right. no f3(n,j) 
        end
    end
    
    for j=m:-1:2
        for i=1:n
            f2(i,j)=f2(i,j-1); % from top to bottom. no f2(i,1)
        end
    end
    
    for j=m:-1:2
        for i=n:-1:2
            f5(i,j)=f5(i-1,j-1); % from top to bottom, no f5(i,j=1),f5(i=1,j)
        end
    end
    
    for j=m:-1:2
        for i=1:n-1
            f6(i,j)=f6(i+1,j-1); % from top to bottom, no f5(i,j=1),f5(i=1,j)
        end
    end
    
    for j=1:m-1
        for i=1:n
            f4(i,j)=f4(i,j+1);
        end
    end
    
    for j=1:m-1
        for i=1:n-1
            f7(i,j)=f7(i+1,j+1);
        end
    end
    
    for j=1:m-1
        for i=n:-1:2
            f8(i,j)=f8(i-1,j+1);
        end
    end
    % step 3: calculate the distribution functions at surfaces (boundaries)
    % see fig. 5.6 
    % (1) bounce back on west boundary (bounce back can be used to simulate non-slip bc)
    for j=1:m
        f1(1,j)=f3(1,j);
        f5(1,j)=f7(1,j);
        f8(1,j)=f6(1,j);
    %(2) bounce back on east boundary
        f3(n,j)=f1(n,j);
        f7(n,j)=f5(n,j);
        f6(n,j)=f8(n,j);    
    end
    %(3)bounce back on south boundary
    for i=1:n
        f2(i,1)=f4(i,1);
        f5(i,1)=f7(i,1);
        f6(i,1)=f8(i,1);
    end
    %(4)moving lid , north bundary (boundary with known velocity)
    % Apply the eqns 5.29 ~ 5.32
    for i=2:n-1 %( only the interior grids )
        rhon(i)=f9(i,m)+f1(i,m)+f3(i,m)+2*(f2(i,m)+f6(i,m)+f5(i,m));
        f4(i,m)=f2(i,m);
        f8(i,m)=f6(i,m)+rhon(i)*uo/6;
        f7(i,m)=f5(i,m)-rhon(i)*uo/6; 
    end
    % step 4: calculate density rho and velocity components u,v
    % 1. calculate updated density rho(i,j) for each grid
    for j=1:m
        for i=1:n
            ssum=zeros(m,n);
            ssum=f1+f2+f3+f4+f5+f6+f7+f8+f9;
            rho(i,j)=ssum(i,j);
        end
    end
    % 2. calculate updated velocity components u,v
    % interior points
    % eq. 5.10 ; eq. 5.11
    for i=1:n
        for j=1:m
            usum=zeros(m,n); 
            vsum=zeros(m,n);
            usum=f1*cx(1)+f2*cx(2)+f3*cx(3)+f4*cx(4)+f5*cx(5)+f6*cx(6)+f7*cx(7)+f8*cx(8)+f9*cx(9);
            vsum=f1*cy(1)+f2*cy(2)+f3*cy(3)+f4*cy(4)+f5*cy(5)+f6*cy(6)+f7*cy(7)+f8*cy(8)+f9*cy(9);
            u(i,j)=usum(i,j)/rho(i,j);
            v(i,j)=vsum(i,j)/rho(i,j);
        end
    end
end

UU=rot90(u);
VV=rot90(v);
U=flipud(UU);
V=flipud(VV);
[x,y] = meshgrid(1:m,1:n);
figure
quiver(x,y,U,V)

qasedak700 · February 24, 2014, 10:37am

Hi dear gn02317298,

I think density fluctuation is OK (I am not sure, if I am wrong or any reason, please tell me).
But I wonder why you have used such for-loops. It takes lots of CPU time. Try to use vector-operations.
Another thing, I ran your code and I got a strange figure. I am not sure, but I think there is problem with your result.
Check aaguirr2 last post for the MATLAB code.

Arashmohammadi · December 26, 2017, 5:35pm

Hi
I have a problem in write a code of lid driven cavity.i write it for three times.but every time i cant solve the problem, and take a correct results.is any body here that can help me and send me a fortran code of lid driven cavitythat i compare it with mine. Mrarashmohammadi@yahoo.com

mahiya · August 26, 2021, 6:55am

Could you please help me out with the North side boundary condition.
As per the formula it should be like.
f4(i,m) = f2(i,m);
f8(i,m) = f6(i,m) + 0.5*(f3(i,m) - f1(i,m)) + 0.5rhon(i)uo;
f7(i,m) = f5(i,m) + 0.5(f1(i,m) - f3(i,m)) - 0.5rhon(i)*uo;