Posted by gn02317298
 About lid-driven cavity flow November 28, 2013 10:03PM Registered: 5 years ago Posts: 5
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 Registered: 10 years ago Posts: 60
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 Registered: 5 years ago Posts: 5
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 Registered: 10 years ago Posts: 60
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 Registered: 5 years ago Posts: 5
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 Registered: 5 years ago Posts: 14
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.
 Re: About lid-driven cavity flow February 01, 2014 04:30PM Registered: 8 years ago Posts: 3
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 Registered: 5 years ago Posts: 5

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 Registered: 6 years ago Posts: 4
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 Registered: 1 year ago Posts: 1
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.