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
```