Welcome! Log In Create A New Profile

Advanced

About lid-driven cavity flow

Posted by gn02317298gn02317298  
About lid-driven cavity flow
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
Re: About lid-driven cavity flow
November 29, 2013 09: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);

------------------------------------------------------
Philippe Seil
CD Lab on Particulate Flow Modelling
JKU Linz, Austria
Re: About lid-driven cavity flow
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
Re: About lid-driven cavity flow
December 03, 2013 09:17AM
Where do you get these big densities? At the inlet or somewhere inside the domain?

------------------------------------------------------
Philippe Seil
CD Lab on Particulate Flow Modelling
JKU Linz, Austria
Re: About lid-driven cavity flow
December 03, 2013 04:14PM
These big densities distributed on the whole domain. I use surf(rho) command in matlab to see the distribution of density.


Min.
Re: About lid-driven cavity flow
January 27, 2014 07: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.
Re: About lid-driven cavity flow
February 01, 2014 04: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.



Edited 1 time(s). Last edit at 02/03/2014 03:56AM by aaguirr2.
Re: About lid-driven cavity flow
February 22, 2014 05: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)



Edited 1 time(s). Last edit at 02/22/2014 05:07AM by gn02317298.
Re: About lid-driven cavity flow
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.
Re: About lid-driven cavity flow
December 26, 2017 05: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
Sorry, you do not have permission to post/reply in this forum.