backsteplbm.cpp

格子Boltzmann方法 格子Boltzmann方法是为了保留格子气自动机方法的优点

#include <stdio.h>
#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <iomanip.h>//C head file
#include <time.h>
#define Q 9
#define M 600//the number of x axis direction nodes 
#define N 60 //the number of y axis direction nodes
#define M1 20
#define N1 20// backforward facing height and width
//discrte direction e
int e[Q][2] = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
double w[Q]={4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36};//setting weights w
double C, p[N+1][M+1], f[N+1][M+1][Q], F[N+1][M+1][Q], u[N+1][M+1][2];
double Computefeq(int k, double p, double u[2])//equbrium distribution funciton
{
  double eu,uv,feq;
  eu = (e[k][0]*u[0]+e[k][1]*u[1])/C;
  uv = (u[0]*u[0]+u[1]*u[1])/(C*C);
  feq = w[k]*p*(1+3*eu+4.5*eu*eu-1.5*uv);
  return feq;
}
void datadeal() //write into files
{
int i,j;
FILE *fp;

	if((fp=fopen("data.dat","w"))==NULL)
	{
		printf("cannot open file\n");
		return;	
	}

  
   for(i=0;i<=M;i++)
		   for(j=0;j<=N;j++)
		   {fprintf(fp," %d %d %le %le %le\n", i, j, u[j][i][0], u[j][i][1], p[j][i]);
   
	   }
   fclose(fp);
  }  

void main()
{
double Re, w, error_p, error_u, error_v, duration;
double temp_p[1000], temp_u[1000], temp_v[1000];
int m, m0, k, i, j, id, jd;//m is the number of cylces , m0 is temp time stpes, k is direction, i is axis direction                            //node j is y axis direction node
   FILE *fp;
   time_t start, finish;
  //Computer particle velocity
   cout <<"Please input the value of Re, w:";//input Re and Omega
   cin >>Re >> w;
   c = 3*N/(Re*(1.0/w-0.5));

   cout <<"c = " <<C;
   start = clock();
//initial u,v,rho
   for(j = 0; j <= N; j++) 
       for(i = 0; i <= M; i++) {
           u[j][i][0] = u[j][i][1] = 0.0;
           p[j][i] = 1.0;
}
   for(j = N1; j <N; j++) //setting inlet velocity
           //u[j][0][0] =1;
   u[j][0][0] = 1*6.0*(((N-N1)*(j-N1))-(j-N1)*(j-N1))/((N-N1)*(N-N1));
          
   for(j = 0; j <= N; j++)
      for(i = 0; i <= M; i++)
         for(k = 0; k < Q; k++)
           f[j][i][k] = Computefeq(k, p[j][i], u[j][i]);
   for(m = 0; ;m++) {
//interpoint
   for(j = N1+1; j < N; j++)
      for(i = 1; i < M; i++)
         for(k = 0; k < Q; k++) {
          id = i-e[k][0];
          jd = j-e[k][1];
         F[j][i][k] = f[jd][id][k]*(1-w)+w*Computefeq(k, p[jd][id], u[jd][id]);
		 }
for(j = 1; j < N1+1; j++)
      for(i = M1+1; i < M; i++)
         for(k = 0; k < Q; k++) {
          id = i-e[k][0];
          jd = j-e[k][1];
         F[j][i][k] = f[jd][id][k]*(1-w)+w*Computefeq(k, p[jd][id], u[jd][id]);
		 }

	 
//calculate marcoscopic variables
   for(j = N1+1; j < N; j++)
       for(i = 1; i < M; i++) {
           p[j][i] = 0; u[j][i][0] = 0; u[j][i][1] = 0;
           for(k = 0; k < Q; k++) {
               f[j][i][k] = F[j][i][k];
               p[j][i] += f[j][i][k];
               u[j][i][0] += e[k][0]*f[j][i][k];
              u[j][i][1] += e[k][1]*f[j][i][k];
}
           u[j][i][0] *= (C/p[j][i]);
           u[j][i][1] *= (C/p[j][i]);
		 
}
	     
	   for(j = 1; j < N1+1; j++)
       for(i = M1+1; i < M; i++) {
           p[j][i] = 0; u[j][i][0] = 0; u[j][i][1] = 0;
           for(k = 0; k < Q; k++) {
               f[j][i][k] = F[j][i][k];
               p[j][i] += f[j][i][k];
              u[j][i][0] += e[k][0]*f[j][i][k];
               u[j][i][1] += e[k][1]*f[j][i][k];
}
           u[j][i][0] *= (C/p[j][i]);
           u[j][i][1] *= (C/p[j][i]);
		 
		
}

//the nodes on the boundary
    for(j = N1+1; j < N; j++) //left 
    for(k = 0; k < Q; k++) {
p[j][0]=p[j][1];
f[j][0][k] = Computefeq(k, p[j][0], u[j][0])+f[j][1][k]-Computefeq(k, p[j][1],
																   u[j][1]);}

for(j = 1; j < N; j++) //right
    for(k = 0; k < Q; k++) {
u[j][M][0]=u[j][M-1][0];
u[j][M][1]=u[j][M-1][1];
f[j][M][k] = Computefeq(k, p[j][M], u[j][M])+f[j][M-1][k]-Computefeq(k, p[j][M-1],
u[j][M-1]);}



    for(j = 1; j <N1; j++) //left real wall
		for(k = 0; k < Q; k++) {
      p[j][M1]=p[j][M1+1];
	  f[j][M1][k] = Computefeq(k, p[j][M1], u[j][M1])+f[j][M1+1][k]-Computefeq(k, p[j][M1+1],
u[j][M1+1]);}

        
	for(i = 1; i < M; i++) //top 
           for( k = 0; k < Q; k++ ){
               p[N][i] = p[N-1][i];
          f[N][i][k] = Computefeq(k, p[N][i], u[N][i])+f[N-1][i][k]-Computefeq(k, p[N-1][i],
u[N-1][i]);}


         

for(i = 1; i < M1; i++) //bottom 1stpart
for( k = 0; k < Q; k++ )
{p[N1][i] = p[N1+1][i];
f[N1][i][k] = Computefeq(k, p[N1][i], u[N1][i])+f[N1+1][i][k]-Computefeq(k, p[N1+1][i], u[N1+1][i]);
		   }

for(i = M1+1; i < M; i++) //bottom 2nd part
for( k = 0; k < Q; k++ )
{p[0][i] = p[1][i];
f[0][i][k] = Computefeq(k, p[0][i], u[0][i])+f[1][i][k]-Computefeq(k, p[1][i], u[1][i]);
		   }

//6 corner on the boundary
          for( k = 0; k < Q; k++ ) {
                p[N1][0] = p[N1+1][1];
f[N1][0][k] = Computefeq(k, p[N1][0], u[N1][0])+f[N1+1][1][k]-Computefeq(k, p[N+1][1], u[N1+1][1]);
                p[N][0] = p[N-1][1];
f[N][0][k]=Computefeq(k,p[N][0],u[N][0])+f[N-1][1][k]-Computefeq(k,p[N-1][1],u[N-1][1]);
                p[N][M] = p[N-1][M-1];
f[N][M][k] = Computefeq(k, p[N][M], u[N][M])+f[N-1][M-1][k]-Computefeq(k,
p[N-1][M-1], u[N-1][M-1]);
p[0][M] = p[1][M-1];
f[0][M][k]= Computefeq(k, p[0][M], u[0][M])+f[1][M-1][k]-Computefeq(k,
p[1][M-1], u[1][M-1]);

p[N1][M1]=p[N1+1][M1+1];
f[N1][M1][k]= Computefeq(k, p[N1][M1], u[N1][M1])+f[N1+1][M1+1][k]-Computefeq(k,
p[N1+1][M1+1], u[N1+1][M1+1]);
p[0][M1]=p[1][M1+1];
f[0][M1+1][k]= Computefeq(k, p[0][M1], u[0][M1])+f[1][M1+1][k]-Computefeq(k,
																	  p[1][M1+1], u[1][M1+1]);}
		  
		  
//computer steady state time
if (m%100 == 0)
cout <<"The "<< m<<"th computation result:"
<<endl<<"The p,u,v of point(M1+20,N1/2)is:"
<<setprecision(7)<<p[M1+20][N1/2]
<<","<<u[M1+20][N1/2][0]<<","<<u[M1+20][N1/2][1]<<endl;
if(m < 1000) {
temp_p[m] = p[M1+20][N1/2];
error_p = fabs(p[M1+20][N1/2]-temp_p[0]);
temp_u[m] = u[M1+20][N1/2][0];
error_u = fabs(u[M1+20][N1/2][0]-temp_u[0]);
temp_v[m] = u[M1+20][N1/2][1];
error_v = fabs(u[M1+20][N1/2][1]-temp_v[0]);
if (m%100 == 0)
cout <<"The error of p, u, v is:" << setiosflags(ios::scientific)
<<error_p<<","<<error_u <<","<<error_v <<"."<<endl;
}
else {
m0 = m%1000;
error_p = fabs(p[M1+20][N1/2]-temp_p[m0]);
error_u = fabs(u[M1+20][N1/2][0]-temp_u[m0]);
error_v = fabs(u[M1+20][N1/2][1]-temp_v[m0]);
if (m%100 == 0)
cout <<"The error of p, u, v is:" <<error_p <<","<<error_u
<<","<<error_v <<"."<<endl;
if(error_u < 5.0e-6 && error_v < 5.0e-6)
break;
temp_p[m0] = p[M1+20][N1/2];
temp_u[m0] = u[M1+20][N1/2][0];
temp_v[m0] = u[M1+20][N1/2][1];
}
if(!(m%2000))
datadeal();
}
finish = clock();
duration = (finish - start) / CLK_TCK;
cout <<"The computation times: "
<<m << "/nduration :" << duration;
if((fp = fopen("times_duration","w")) == NULL) { //save time duration
printf("File Open Error\n"); exit(1);
}
fprintf(fp,"times: %d \nduration: %f", m, duration);
fclose(fp);
datadeal();
}