Lid cavity flow

Shaq · March 23, 2008, 1:58pm

Hello!

I’m a newbie in the LBM world, I’m trying to simulate a lid cavity flow (square domain), the sliding wall is the left one. To do that I modified the cylinder.m example code, in particular I applied bounce back BC on the right wall and the Zou/He BC on the left sliding wall. However I can’t get the right results. Can anybody help me?
Here I paisted the code, any suggestions will be very useful!

Thank you

P.S Mr Latt please help me!

% GENERAL FLOW CONSTANTS
lx = 100;
ly = 100;

uMax = 0.02; % vertical left-wall velocity
vMax = 0; % horizontal left-wall velocity
Re = 100; % Reynolds number
nu = uMax * 2.ly / Re; % kinematic viscosity
omega = 1. / (3nu+1./2.); % relaxation parameter
maxT = 400000; % total number of iterations
tPlot = 5; % cycles

% D2Q9 LATTICE CONSTANTS
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];
col = [2: (ly-1)];

[y,x] = meshgrid(1:ly,1:lx);
obst = zeros(lx,ly);
obst(:,[1,ly]) = 1;
obst(lx,: ) = 1;
bbRegion = find(obst);

% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)
fIn = reshape( t’ * ones(1,lx*ly), 9, lx, ly);

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT

% MACROSCOPIC VARIABLES
rho = sum(fIn);
ux  = reshape ( ...
      (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
uy  = reshape ( ...
      (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
  
% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS

uy(:,1,col) = uMax;  %left-wall y - velocity
ux(:,1,col) = 0;        %left-wall x - velocity
rho(:,1,col) = 1 ./ (1-ux(:,1,col)) .* ( ...
    sum(fIn([1,3,5],1,col)) + ...
    2*sum(fIn([4,7,8],1,col)) );

% COLLISION STEP
for i=1:9
   cu = 3*(cx(i)*ux+cy(i)*uy);
   fEq(i,:,: )  = rho .* t(i) .* ...
     ( 1 + cu + 1/2*(cu.*cu) ...
          - 3/2*(ux.^2+uy.^2) );
   fOut(i,:,: ) = fIn(i,:,: ) - ...
     omega .* (fIn(i,:,: )-fEq(i,:,: ));
end

% MICROSCOPIC BOUNDARY CONDITIONS
for i=1:9
     % Left boundary (Zou/He BC)
     fOut(2,1,col) = fIn(4,1,col) + 2/3*rho(:,1,col)*vMax;
     fOut(6,1,col) = fIn(8,1,col) + 1/2*(fIn(5,1,col)-fIn(3,1,col))+...
         1/2*rho(:,1,col)*uMax + 1/6*rho(:,1,col)*vMax;
     fOut(9,1,col) = fIn(7,1,col) + 1/2*(fIn(3,1,col)-fIn(5,1,col))-...
         1/2*rho(:,1,col)*uMax + 1/6*rho(:,1,col)*vMax;
     
     % Bounce back region
     fOut(i,bbRegion) = fIn(opp(i),bbRegion);
end

 % STREAMING STEP
for i=1:9
   fIn(i,:,: ) = ...
     circshift(fOut(i,:,: ), [0,cx(i),cy(i)]);
end

end

orestis · March 25, 2008, 11:11am

Hi,
in what sense the results are not those expected? Is the code unstable? Or does it give results uncompatbile with results obtained with other methods?

If the problem is instability maybe you should check your relaxation frequency (omega). This must not be bigger than 1.8 because the zou-he boundary condition is unstable above that number…

Orestis

jlatt · March 27, 2008, 3:13pm

Hello,

There are two points in your code which seem a little odd to me:

The loop "for i=1:9 " which comes right after the comment “% MICROSCOPIC BOUNDARY CONDITIONS” is not needed for the Zou/He boundary condition. It is only needed for bounce-back nodes and should be postponed correspondingly.
Be aware of the fact that Zou/He boundary nodes need to execute a normal BGK collision. I know that many people get this point wrong, and I therefore emphasize this once more: Zou/He boundary nodes are required to execute a normal BGK collision step! In your code, this simply means that implementation of Zou/He should come before collision, and that the values of the Zou/He algorithm should be assigned to fIn([2,6,9],1,col) instead of fOut.

Another point you should be aware of if you compare the results with the literature on lid-driven cavity is that the physical boundary is located half-way between nodes for bounce-back, but right on top of nodes for Zou/He. That means that the physical height of your cavity is ly-2, but the width is lx-3/2.

I haven’t tried it out, but this might fix your code. Btw, I like you adaptation of the Matlab code to lid-driven cavity. Once the code works, I encourage you to send it to me, and we could post it on the lbmethod.org web page along with other LB samples.

Shaq · March 29, 2008, 5:24pm

Thank you very much Latt, now the code works!
As you suggested in the previous post I can send the code to you, then you can put it with other LB samples on lbmethod.org, if you want. I will be very happy about that. Let me know…

jlatt · March 30, 2008, 4:36pm

Let’s go for it. You’ll get my e-mail address by clicking on my profile name above.

lbmagis · June 26, 2012, 10:49pm

hello every one,
i have to simulate a lid driven cavity 2D flow with LBM ,i’d like to ask u how to write a matlab"s program that draw the vorticity contour;please could u write it to me because i am new on matlab,
thank you very much

orestis · June 27, 2012, 7:26am

Do you want to know how to compute the vorticity or to aplly a contour plot in matlab?

I think it’s the latter, so you need to use the “contour” function on the vorticity field.

lbmagis · June 27, 2012, 3:47pm

hello Malaspin,
thank u so much Malaspin
i want to compute the vorticity and view the contour of it “those lines”
thank u for sending me both the matlab program of how to compute vorticity and contour’s vorticity thnx a lot
haave a great day

lbmagis · June 27, 2012, 3:53pm

hello,
i want both matlab program that compute and draw the contour of vorticicty
thank u so much and have a nice day

shikhabhuyan · November 9, 2015, 12:29pm

I am also new to this field, please if anyone can help me out. My code is not giving desired output.
here is the code
#include<stdio.h>
#include"crun4.h"

void main()
{
FILE *a,*ifp,fp3,fp4;
a=fopen(“lid.dat”,“w”);
ifp = fopen(“grid.dat”, “w”);
fp3=fopen(“uvel.dat”,“w”);
fp4 = fopen(“vvel.dat”, “w”);
fprintf(ifp,“ZONE I=200\t J=200”);
fprintf(ifp,"\n");
dx=1.0/nx;
dy=1.0/ny;
//dt=dx;
//tau=0.7;
nu=1/Rex;
tau=(3nu+0.5);
dt=(3nu)/(tau-0.5);
u1=(dt/dx);
x[0][0]=0.0;
y[0][0]=0.0;
x[nx][j]=1.0;
y[i][ny]=1.0;

/*----------------------------grid generation----------------------------*/
	for(i=1;i<nx;i++)
	{
            for(j=1;j<ny;j++)
	  {	

	
		x[i][j]=i*dx;
		y[i][j]=j*dy;
                fprintf(ifp,"%lf\t %lf\n", x[i][j],y[i][j]);
		//printf("grids are: %f\t%f\n",x[i],y[j]);
	}
}

/*----------------------------initial condition----------------------------*/

for(i=0;i<=nx;i++)
{
	for(j=0;j<ny;j++)
		{
			rho[i][j]=1.0;
			ux[i][j]=0.0;
			uy[i][j]=0.0;
			for(k=0;k<9;k++)
				{
					f[i][j][k]=w_f[k]*rho[i][j];

// //printf(“initialised PDF\n”);
}
}
}
for(i=1;i<nx;i++)
{
rho[i][ny]=1.0;
ux[i][ny]=u1;
uy[i][ny]=0.0;
// F2=0.01;
// for(k=0;k<9;k++)
// {
// F1=ex[k]0.1;
// feq=w_f[k]rho[i][ny](1.0+(3.0F1)+(4.5F1F1)-(1.5*F2));
// f[i][ny][k]=feq;
// //printf("%f\n",f[i][ny][k]);
// }
}
for(i=0;i<=nx;i++)
{

      // rhon=f[i][ny][0]+f[i][ny][1]+f[i][ny][3]+2.0*((f[i][ny][2])+(f[i][ny][6])+(f[i][ny][5]));
	F2=u1*u1;
	for(k=0;k<9;k++)
	{
		F1=ex[k]*u1;
		feq=w_f[k]*rho[i][ny]*(1.0+(3.0*F1)+(4.5*F1*F1)-(1.5*F2));			
		f[i][ny][k]=feq;

// //printf("%f\n",f[i][ny][k]);
}
/f[i][ny][4]=f[i][ny][2];
f[i][ny][8]=f[i][ny][6]+(rho[i][ny](ux[i][ny]/6.0));
f[i][ny][7]=f[i][ny][5]+(rho[i][ny](ux[i][ny]/6.0));/

}


/*-------------------------------main loop---------------------------------*/

error=10.0;
do
{
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
E2=ux[i][j]ux[i][j]+uy[i][j]uy[i][j];
for(k=0;k<9;k++)
{
E1=ex[k]ux[i][j]+ey[k]uy[i][j];
feq=w_f[k]rho[i][j](1.0+(3.0E1)+(4.5E1E1)-(1.5E2));
f[i][j][k]=f[i][j][k]+((feq-f[i][j][k])/tau);
//printf("%f\n",fp[i][j][k]);
}
}
}
/-----------------------------Streaming step-----------------------------/
for(j=0;j<=ny;j++)
{
for(i=nx;i>=1;i–)
{
f[i][j][1]=f[i-1][j][1];

	}

	for(i=0;i<nx;i++)
	{
		f[i][j][3]=f[i+1][j][3];
	}
}
for(j=ny;j>=1;j--)							/*At j=0,we have pdfs from Boundary conditions*/
{
	for(i=0;i<=nx;i++)
	{
		f[i][j][2]=f[i][j-1][2];
	}
	for(i=nx;i>0;i--)
	{
		f[i][j][5]=f[i-1][j-1][5];
	}
	for(i=0;i<nx;i++)
	{
		f[i][j][6]=f[i+1][j-1][6];
	}
}
for(j=0;j<ny;j++)							/*At j=ny,pdfs are defined*/
{
	for(i=0;i<=nx;i++)
	{
		f[i][j][4]=f[i][j+1][4];
	}
	for(i=0;i<nx;i++)
	{
		f[i][j][7]=f[i+1][j+1][7];
	}
	for(i=nx;i>=1;i--)
	{
		f[i][j][8]=f[i-1][j+1][8];
	}
}


/*---------------------------------------Boundary Condition--------------------------------------*/

/*---------------------left wall boundary condition-----------------------*/
for(j=1;j<ny;j++)
{
	f[0][j][1]=f[0][j][3];
	f[0][j][5]=f[0][j][7];
	f[0][j][8]=f[0][j][6];
}
/*---------------------right wall boundary condition-----------------------*/
for(j=1;j<ny;j++)
{
	f[nx][j][3]=f[nx][j][1];
	f[nx][j][6]=f[nx][j][8];
	f[nx][j][7]=f[nx][j][5];
}
/*---------------------bottom wall boundary condition-----------------------*/

for(i=0;i<=nx;i++)
{
	f[i][0][2]=f[i][0][4];
	f[i][0][5]=f[i][0][7];
	f[i][0][6]=f[i][0][8];
}
/*---------------------top wall boundary condition-----------------------*/




/*-----------------------------Calculation of velocity and density-----------------------------*/
for(i=1;i<nx;i++)
	{
		for(j=1;j<ny;j++)
			{
			  sum=0.0;
			   for(k=0;k<9;k++)
			   {
			     sum=sum+f[i][j][k];
			   }
			   rho[i][j]=sum;

// rho[i][j]=f[i][j][1]+f[i][j][2]+f[i][j][3]+f[i][j][4]+f[i][j][5]+f[i][j][6]+f[i][j][7]+f[i][j][8]+f[i][j][0];
// ux[i][j]=((f[i][j][1]+f[i][j][5]+f[i][j][8])-(f[i][j][3]+f[i][j][6]+f[i][j][7]))/rho[i][j];
// uy*/[i][j]=((f[i][j][2]+f[i][j][5]+f[i][j][6])-(f[i][j][4]+f[i][j][7]+f[i][j][8]))/rho[i][j];
}
}

	for(i=1;i<nx;i++)
	{
		for(j=1;j<ny;j++)
			{
			  uxsum=0.0;
			  uysum=0.0;
			   for(k=0;k<9;k++)
			   {
			     uxsum=uxsum+(f[i][j][k]*ex[k]);
			     uysum=uysum+(f[i][j][k]*ey[k]);
			   }
			   ux[i][j]=uxsum/rho[i][j];
			   uy[i][j]=uysum/rho[i][j];

			}
	}
	
/*-------------------------------------error checking condition------------------------------------------*/

e1=e2=0.0;
for(i=1;i<nx;i++)
{
	for(j=1;j<ny;j++)
	{
		e1=e1+(sqrt((ux[i][j]-ux0[i][j])*(ux[i][j]-ux0[i][j])+(uy[i][j]-uy0[i][j])*(uy[i][j]-uy0[i][j])));
		e2=e2+(sqrt((ux[i][j]*ux[i][j])+(uy[i][j]*uy[i][j])));
		ux0[i][j]=ux[i][j];
		uy0[i][j]=uy[i][j];
	}
}
error=e1/e2;
if(timestep%1000==0)
{
	printf("%d\t%0.13f\n",timestep,error);
}
	timestep++;

// if(timestep>=100000) error=0.0;
}
while(error>=1e-8);
printf("%d\t%0.13f\n",timestep,error);
/---------------------------------------------output------------------------------------------------------/
fprintf(a,“Zone I=%d\tJ=%d \n”,nx+1,ny+1 );
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
fprintf(a,"%f\t%f\t%f\t %f\n",x[i],y[j],ux[i][j],uy[i][j]);
}
}
for(j=1;j<=ny;j++)
{
fprintf(fp3,"\n%f\t%f",y[100][j],ux[100][j]);
}
for(i=1;i<=nx;i++)
{
fprintf(fp4,"\n%f\t%f",x[i][100],uy[i][100]);
}
fclose(ifp);
fclose(a);
fclose(fp3);
fclose(fp4);
getch();
}
crun4.h

#include<stdio.h>
#include<math.h>
#define nx 200
#define ny 200
#define arr 201

int i,j,k,timestep=0;
double x[arr][arr],y[arr][arr],rho[arr][arr],ux[arr][arr],uy[arr][arr],f[arr][arr][10],fp[arr][arr][10],feq,ux0[arr][arr],uy0[arr][arr],rhon,sum,uxsum,uysum;
double w_f[9]={(4.0/9.0),(1.0/9.0),(1.0/9.0),(1.0/9.0),(1.0/9.0),(1.0/36.0),(1.0/36.0),(1.0/36.0),(1.0/36.0)};
int ex[9]= {0,1,0,-1,0,1,-1,-1,1};
int ey[9]= {0,0,1,0,-1,1,1,-1,-1};
double dx,dy,E1,E2,tau,nu,dt,F1,F2,e1,e2,error,u1;
double Rex=400;

And please suggest me any links from where I can get the LBM concepts clear.