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:

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

------------------------------------------------------

Philippe Seil

CD Lab on Particulate Flow Modelling

JKU Linz, Austria

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

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 |

Re: About lid-driven cavity flow December 03, 2013 04:14PM |
Registered: 5 years ago Posts: 5 |

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.

Check this page (http://www.palabos.org/forum/read.php?3,6497) to see my codes.

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

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

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 |

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 ?

Edited 1 time(s). Last edit at 02/22/2014 05:07AM by gn02317298.

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

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

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.