Hi Dear all,
I study LBM in very detail but the computational work is new for me. I recently develop a code for flows around a circular cylinder using the extrapolation boundary conditions proposed by Guo et al. and the Lattice BGK Model for incompressible Navier-Stokes Equation also proposed by Guo et al. (2000). I would like that some one check this code that it’s correct or wrong. I try to calculate the drag and lift coefficient using this code for different Reynolds number. Thanks in advance.
code: The code is a complete code and i compile and run it no error but after running the code i want to calculate the drag and lift coefficient but the results that i want i can’t get. So especially Jonal Latt, Thomas Morton and any one like to check this that the code is correct or not.
//Shams
/* lattice direction
x y
{ 0 0 } rest
{ + 0 } dir 1
{ - 0 } dir 3
{ 0 + } dir 2
{ 0 - } dir 4
{ + + } dir 5
{ - - } dir 7
{ + - } dir 8
{ - + } dir 6
*/
#include
#include
#include
#include
#include
#include
using namespace std;
#define L 16
#define M (20L) //grid points in y-direction //rectangular area 40D20D
#define N (40*L) //grid points in x-direction
#define M1 (M+1)
#define N1 (N+1)
#define D 1.0 //diameter of the cylinder
#define Ly (MD/L) //position of cylinder in x-direction
#define Lx (ND/L) //position of cylinder in y-direction
#define xc (10.0*D) //centre of a cylinder in x-direction
#define yc (Ly/2.0) //centre of a cylinder in y-direction
#define P0 1.0
#define U0 0.1
#define Q 9
#define Mv 1000
#define PI (3.1415926)
#define Re 100.0 //Reynolds numbers
void Lbini(void);
void Evol(void);
void exbd(void);
double feq(int k,double u,double v);
void datadeal(int m);
void data_read(void);
void time_data(int m);
void drag_lift(void);
long filesize(void);
void Remesh(void);
void DF_BD(void);
void vortex(void);
double w,w1,dx,dt;
double f[M1][N1][Q],u[M1][N1],v[M1][N1],p[M1][N1],ff[M1][N1][Q],vor[M1][N1];
double eps[M1][N1][Q];
int e[Q][2],r[Q],XC0;
int flag[M1][N1],oldflag[M1][N1];
double FD,FL; //FD for drag force on cylinder and FL lift force on cylinder
double Fxy[M1][2];
double tp[Q];
double utmp[M1][N1],vtmp[M1][N1];
double filter0=0.02, filter1=0.96;
double xcat,ycat,ucat,vcat;
void main()
{
FILE *fp;
int kmax,I,j,i;
long m;
double TEND=0, rRe,t=0.0;
int readdata;
double tau;
double oldD,oldL,newD,newL,oldFD,oldFL;
//the following are lattice velocity vectors
e[0][0]=e[0][1]=0;
e[1][0]=1; e[1][1]=0;
e[2][0]=0; e[2][1]=1;
e[3][0]=-1; e[3][1]=0;
e[4][0]=0; e[4][1]=-1;
e[5][0]=1; e[5][1]=1;
e[6][0]=-1; e[6][1]=1;
e[7][0]=-1; e[7][1]=-1;
e[8][0]=1; e[8][1]=-1;
r[1]=3; r[3]=1;
r[2]=4; r[4]=2;
r[5]=7; r[7]=5;
r[6]=8; r[8]=6;
//the following are the weighting coefficients
tp[0]=4.0/9;
tp[1]=tp[2]=tp[3]=tp[4]=1.0/9;
tp[5]=tp[6]=tp[7]=tp[8]=1.0/36;
printf("Read data? (yes=1 or no=0)\n");
scanf("%d",&readdata);
dt=dx=D/L;
XC0=(int)(xc/dx+0.1);
rRe=D*U0/Re;
tau=0.5+3*rRe/dx; w=1.0/tau; w1=1.0-w;
Lbini();
if(readdata) data_read();
for(m=0;m<Mv;m++)
Fxy[m][0]=Fxy[m][1]=0.0;
ucat=0.0; vcat=0.0;
xcat=xc; ycat=yc;
Remesh();
for(j=0;j<M1;j++) for(i=0;i<N1;i++)
oldflag[j][i]=flag[j][i];
DF_BD();
Evol();
exbd();
drag_lift();
oldFD=FD; oldFL=FL; oldD=FD; oldL=FL;
m=filesize();
if(!readdata)
{
printf("This the first run\n");
m=0;
}
printf("m=%d\n",m);
AA:
printf(“input kmax:\n”);
scanf("%d",&kmax);
TEND=TEND+kmax*dt;
printf(“w=%f U0=%f\n”,w,U0);
while((t<TEND))
{
m++;
t+=dt;
if(m<500)
{
xcat=xc; ycat=yc; ucat=0; vcat=0;
}
Remesh();
DF_BD();
Evol();
exbd();
drag_lift();
newD=FD; newL=FL;
FD=filter0*(newD+oldD)+filter1*oldFD;
FL=filter0*(newL+oldL)+filter1*oldFL;
oldFD=FD; oldFL=FL;
oldD=newD; oldL=newL;
for(j=0;j<M1;j++) for(i=0;i<N1;i++) oldflag[j][i]=flag[j][i];
if(m%100==0)
{
printf("m=%d, Drag=%f, Lift=%f\n",m,FD,FL);
}
Fxy[(m-1)%Mv][0]=FD; Fxy[(m-1)%Mv][1]=FL;
if(m%Mv==0)
{
if((fp=fopen("Force.dat","a"))==NULL)
{
printf("File Open Error\n");
exit(1);
}
for(I=0;I<Mv;I++)
fprintf(fp,"%f %f %f %f %f\n",((m/Mv-1)*Mv+I)*dt,Fxy[I][0],Fxy[I][1]);
fclose(fp);
datadeal(m/Mv);
}
}
datadeal(0);
printf("Continue? (yes=1 or no=0)\n");
scanf("%d",&readdata);
if(readdata) goto AA;
}
void Lbini()
{
int i,j,k;
double EDF;
for(j=0;j<M1;j++) for(i=0;i<N1;i++)
{
u[j][i]=0.0; v[j][i]=0.0;
p[j][i]=P0;
}
for(j=0;j<M1;j++) for(i=0;i<N1;i++) for(k=0;k<Q;k++)
{
EDF=feq(k,u[j][i],v[j][i]);
f[j][i][k]=p[j][i]*EDF;
}
}
void Remesh()
{
int i,j,k,jd,id;
double xr,yr,r,d;
double a,b,c;
for(j=0;j<M1;j++) for(i=0;i<N1;i++) flag[j][i]=1;
r=D/2.0;
for(j=0;j<M1;j++) for(i=0;i<N1;i++){
if(i<XC0)
{
xr=i*dx-xcat; yr=j*dx-ycat; }
d=xr*xr+yr*yr;
if(d<=r*r) { flag[j][i]=0; }
}
for(j=1;j<M;j++) for(i=1;i<N;i++)
{
if(flag[j][i]==1)
for(k=1;k<Q;k++)
{
if(i<XC0)
{
xr=i*dx-xcat; yr=j*dx-ycat; }
jd=j+e[k][1]; id=i+e[k][0];
if(flag[jd][id]==0)
{
a=e[k][0]*e[k][0]+e[k][1]*e[k][1];
b=-(e[k][0]*xr+e[k][1]*yr);
c=(e[k][0]*yr-e[k][1]*xr);
c=r*r*a-c*c;
eps[j][i][k]=(b*b-c)/(a*(-b+sqrt(c)));
}
}
else for(k=0;k<Q;k++)
eps[j][i][k]=1.0;
}
}
void DF_BD()
{
int i,j,k,jd,id;
int num,I;
double EDF;
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
if(flag[j][i]!=0 && oldflag[j][i]==0)
{
num=0;
u[j][i]=v[j][i]=p[j][i]=0.0;
for(I=0;I<Q;I++)
f[j][i][I]=0;
for(k=1;k<Q;k++)
{
jd=j-e[k][1]; id=i-e[k][0];
if(flag[jd][id] && oldflag[jd][id])
{
u[j][i]+=u[jd][id];
v[j][i]+=v[jd][id];
p[j][i]+=p[jd][id];
for(I=0;I<Q;I++)
{ f[j][i][I]+=f[jd][id][I];
}
num++;
}
}
if(num>0)
{
u[j][i]/=num;
v[j][i]/=num;
p[j][i]/=num;
for(I=0;I<Q;I++)
{
f[j][i][I]/=num;
}
}
else
{
u[j][i]=ucat; v[j][i]=vcat; p[j][i]=P0;
for(I=0;I<Q;I++)
{
EDF=feq(I,u[j][i],v[j][i]);
f[j][i][I]=p[j][i]*EDF;
}
}
}
}
void exbd()
{
int i,j,k;
double FMb, FMf;
for(i=0;i<=N;i++) for(k=0;k<Q;k++)
{
u[0][i]=u[1][i];
v[0][i]=v[1][i];
u[M][i]=u[M-1][i];
v[M][i]=v[M-1][i];
p[0][i]=p[1][i];
p[M][i]=p[M-1][i];
}
for(j=1;j<M;j++)
{
u[j][0]=U0; v[j][0]=0.0; p[j][0]=p[j][1];
u[j][N]=u[j][N-1]; v[j][N]=v[j][N-1]; p[j][N]=p[j][N-1];
}
for(j=1;j<M;j++) for(k=0;k<Q;k++)
{
FMb=feq(k,u[j][0],v[j][0]);
FMf=feq(k,u[j][1],v[j][1]);
f[j][0][k]=(f[j][1][k]-p[j][1]*FMf)+p[j][0]*FMb;
FMb=feq(k,u[j][N],v[j][N]);
FMf=feq(k,u[j][N-1],v[j][N-1]);
f[j][N][k]=(f[j][N-1][k]-p[j][N-1]*FMf)+p[j][N]*FMb;
}
// the following is for top and down
for(i=0;i<=N;i++) for(k=0;k<Q;k++)
{
FMb=feq(k,u[0][i],v[0][i]);
FMf=feq(k,u[1][i],v[1][i]);
f[0][i][k]=(f[1][i][k]-p[1][i]*FMf)+FMb*p[0][i];
FMb=feq(k,u[M][i],v[M][i]);
FMf=feq(k,u[M-1][i],v[M-1][i]);
f[M][i][k]=(f[M-1][i][k]-FMf*p[M-1][i])+FMb*p[M][i];
}
}
void drag_lift()
{
int i,j,k,jd,id;
double rrho;
FD=FL=0.0;
for(j=1;j<M;j++) for(i=1;i<N;i++) for(k=1;k<Q;k++)
{
jd=j-e[k][1]; id=i-e[k][0];
if((flag[j][i]==1) && (flag[jd][id]==0))
{
FD-=(e[k][0]*((f[j][i][k]+ff[j][i][r[k]])));
FL-=(e[k][1]*((f[j][i][k]+ff[j][i][r[k]])));
}
}
rrho=2./(L*U0*U0);
FD *=rrho; FL *=rrho;
}
void Evol(void)
{
int i,j,k,jd,id,idd,jdd;
double FM,FMb,FMf,FMff,u_t,v_t;
double u_t1,v_t1,u_t2,v_t2;
double ub,vb;
//collision for fluid nodes
for(j=0;j<=M;j++)
{
for(i=0;i<=N;i++)
{
for(k=0;k<Q;k++)
{
FM=feq(k,u[j][i],v[j][i]);
ff[j][i][k]=f[j][i][k]+w*(FM*p[j][i]-f[j][i][k]);
}
}
}
//collision for wall nodes (outer points)
for(j=1;j<M;j++)
{
for(i=1;i<N;i++)
{
if(flag[j][i]==0) for(k=0;k<Q;k++)
{
jd=j+e[k][1]; id=i+e[k][0];
jdd=jd+e[k][1]; idd=id+e[k][0];
if(i<XC0)
{
ub=ucat; vb=vcat;
}
if(flag[jd][id]==1)
{
u_t1=ub+(eps[jd][id][k]-1)*u[jd][id];
v_t1=vb+(eps[jd][id][k]-1)*v[jd][id];
if(eps[jd][id][k]<0.75)
{
u_t2=ub+ub+(eps[jd][id][k]-1)*u[jdd][idd];
v_t2=vb+vb+(eps[jd][id][k]-1)*v[jdd][idd];
u_t=u_t1+u_t2*(1-eps[jd][id][k])/(1+eps[jd][id][k]);
v_t=v_t1+v_t2*(1-eps[jd][id][k])/(1+eps[jd][id][k]);
FMb=feq(k,u_t,v_t);
FMf=feq(k,u[jd][id],v[jd][id]);
FMff=feq(k,u[jdd][idd],v[jdd][idd]);
ff[j][i][k]=w1*((1-eps[jd][id][k])*(f[jdd][idd][k]-FMff*p[jdd][idd])+eps[jd][id][k]*(f[jd][id][k]-FMf*p[jd][id]))+FMb*p[jd][id];
}
else
{
u_t=u_t1/eps[jd][id][k];
v_t=v_t1/eps[jd][id][k];
FMb=feq(k,u_t,v_t);
FMf=feq(k,u[jd][id],v[jd][id]);
ff[j][i][k]=w1*(f[jd][id][k]-FMf*p[jd][id])+FMb*p[jd][id];
}
}
}
}
}
//for streaming
for(j=1;j<M;j++) for(i=1;i<N;i++) if(flag[j][i]) for(k=0;k<Q;k++)
{
{
jd=j-e[k][1]; id=i-e[k][0];
f[j][i][k]=ff[jd][id][k];
}
}
for(j=1;j<M;j++)
{
for(i=1;i<N;i++)
{
u[j][i]=v[j][i]=p[j][i]=0.0;
for(k=0;k<Q;k++)
{
p[j][i]+=f[j][i][k];
}
for(k=1;k<Q;k++)
{
u[j][i]+=(f[j][i][k]*e[k][0]);
v[j][i]+=(f[j][i][k]*e[k][1]);
}
u[j][i]/=p[j][i];
v[j][i]/=p[j][i];
}
}
}
double feq(int k,double U,double V)
{
double su,uv,eu,x;
eu=(e[k][0]*U+e[k][1]*V);
uv=U*U+V*V;
su=1.0+3*eu+4.5*eu*eu-1.5*uv;
if(k==0) x=su*4.0/9;
if(k>0 && k<5) x=su/9;
if(k>4) x=su/36;
return x;
}
void datadeal(int m)
{
FILE *fp;
char fname[20];
int i,j;
vortex();
for(j=0;j<=M;j++) for(i=0;i<=N;i++) if(flag[j][i]==0)
{
utmp[j][i]=0; vtmp[j][i]=0;
}
else
{
utmp[j][i]=u[j][i]; vtmp[j][i]=v[j][i];
}
sprintf(fname,"%s %d %s","vor",m,".dat");
fp=fopen(fname,"w");
fprintf(fp,"Zone I=%d, J=%d, F=POINT\n",N1,M1);
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
fprintf(fp,"%f %f %f\n",i*dx,j*dx,vor[j][i]);
fclose(fp);
if((fp=fopen("Result.dat","w"))==NULL)
{ printf("File Open Error\n"); exit(1); }
fprintf(fp,"TITLE=%s\n","Flow over circular cylinder");
fprintf(fp,"Variables=%s,%s,%s,%s,%s\n","X","Y","Ux","Uy","Pressure");
fprintf(fp,"Zone I=%d, J=%d, F=POINT\n",N1,M1);
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
fprintf(fp,"%f %f %f %f %f \n",i*dx,j*dx,utmp[j][i]/U0,vtmp[j][i]/U0, p[j][i]);
fclose(fp);
fp=fopen("ut","wb");
fwrite(u,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("vt","wb");
fwrite(v,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("pt","wb");
fwrite(p,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("ft","wb");
fwrite(f,sizeof(double),M1*N1*Q,fp);
fclose(fp);
}
void data_read()
{
FILE *fp;
fp=fopen("ut","rb");
fread(u,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("vt","rb");
fread(v,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("pt","rb");
fread(p,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("ft","rb");
fread(f,sizeof(double),M1*N1*Q,fp);
fclose(fp);
}
void time_data(int m)
{
FILE *fp;
char fname[10];
sprintf(fname,"%s %d","uu",m);
fp=fopen(fname,"wb");
fwrite(u,sizeof(double),M1*N1,fp);
fclose(fp);
sprintf(fname,"%s %d","vv",m);
fp=fopen(fname,"wb");
fwrite(v,sizeof(double),M1*N1,fp);
fclose(fp);
sprintf(fname,"%s %d","pp",m);
fp=fopen(fname,"wb");
fwrite(p,sizeof(double),M1*N1,fp);
fclose(fp);
fp=fopen("ft","wb");
fwrite(f,sizeof(double),M1*N1*Q,fp);
fclose(fp);
}
long filesize(void)
{
FILE *stream;
double a;
long n=0;
if((stream=fopen("force.dat","r"))==NULL)
{
return 0;
}
else
{
while(1)
{
fscanf(stream,"%f",&a);
n++;
if(feof(stream))
{
break;
}
}
}
fclose(stream);
return (n-1)/5;
}
void vortex(void)
{
int i,j,k,id,jd,num;
for(j=1;j<M;j++)
for(i=1;i<N;i++)
if(flag[j][i]==1)
vor[j][i]=((u[j+1][i]-u[j-1][i])-(v[j][i+1]-v[j][i-1]))/(2*dx);
for(j=1;j<M;j++)
for(i=1;i<N;i++)
if(flag[j][i]==0)
{
vor[j][i]=0; num=0;
for(k=1;k<Q;k++)
{
jd=j+e[k][1];
id=i+e[k][0];
if(flag[jd][id]==1)
{
vor[j][i]+=vor[jd][id];
num++;
}
}
vor[j][i]/=num;
}
}