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

``````

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);

``````

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

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

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

Min.

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.

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.

``````

``````

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)

``````

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.

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