1.cpp

是格子boltzmann中的一段代码,希望对做多相流的朋友有所帮助.

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstdlib>
#include<iomanip.h>
#include<ctime>
#include"D2Q9.h"

int main()
{
    FILE* fp;
    
	SetPrameter();
    init_macro();
    init_micro();

    for(m=0; ; m++)     
    {
      evolution();
      output_screen();
      if(error_u<5.0e-6 && error_v<5.0e-6)
      break;
   }

    finish = clock();
    duration = (finish - start) / CLK_TCK;
    cout << "The computation times:" << m << "\nduration:" << duration;
   
    if((fp=fopen("times_duration","w"))==NULL)
    {
     printf("File Open Error\n");
     exit(1);
    }
    fprintf(fp, "times:%d\nduration:%f", m, duration);
    fclose(fp);
    datadeal();
	return 0;
}

void SetPrameter()
{
	cout << "Please input the value of Re,omega:";   
    cin >> Re >> omega;
    C = 3*M/(Re*(1.0/omega-0.5));       // 计算C值
    cout << "C=" << C;
    start = clock();
}

void init_macro()
{
   for(j=0; j<=N; j++)         //初始化u,rho
     for(i=0; i<=M; i++)
     {
       u[j][i][0] = u[j][i][1]=0;
       rho[j][i] = 1.0;
     }
   for(i=0; i<=M; i++)          //顶上一行u值置为1.0
      u[N][i][0] = 1.0; 
}

double Computefeq(int k,double rho,double u[2]) //计算平衡态分布函数
{
  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]*rho*(1+3*eu+4.5*eu*eu-1.5*uv);
  return feq;
}

void init_micro()
{
	for(j=0; j<=N; j++)          //分布函数初始化
     for(i=0; i<=M; i++)
       for(k=0; k<Q; k++)
         f[j][i][k] = Computefeq(k, rho[j][i], u[j][i]);
}

void evolution()
{
    for(j=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-omega) + omega * Computefeq(k, rho[jd][id], u[jd][id]);
        }

    for(j=1; j<N; j++)      // 计算宏观速度和宏观密度
      for(i=1; i<M; i++)
      {
        rho[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];
          rho[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/rho[j][i]);
        u[j][i][1] *= (C/rho[j][i]);
     }

//边界点
     for(j=1; j<N; j++)       //left and right
       for(k=0; k<Q; k++)
      {
        rho[j][0] = rho[j][1];
        f[j][0][k] = Computefeq(k, rho[j][0], u[j][0]) + f[j][1][k] - Computefeq(k, rho[j][1], u[j][1]);
        
        rho[j][N] = rho[j][N-1];
        f[j][N][k] = Computefeq(k, rho[j][N], u[j][N]) + f[j][N-1][k] - Computefeq(k, rho[j][N-1], u[j][N-1]);
      }

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

     for(k=0;k<Q;k++)         //角上四点
     {
       rho[0][0] = rho[1][1];
       f[0][0][k] = Computefeq(k, rho[0][0], u[0][0]) + f[1][1][k] - Computefeq(k, rho[1][1], u[1][1]);

       rho[0][N] = rho[1][N-1];
       f[0][N][k] = Computefeq(k, rho[0][N], u[0][N]) + f[1][N-1][k] - Computefeq(k, rho[1][N-1], u[1][N-1]);

       rho[N][0] = rho[N-1][1];
       f[N][0][k] = Computefeq(k, rho[N][0], u[N][0]) + f[N-1][1][k] - Computefeq(k, rho[N-1][1], u[N-1][1]);

       rho[N][N] = rho[N-1][N-1];
       f[N][N][k] = Computefeq(k, rho[N][N], u[N][N]) + f[N-1][N-1][k] - Computefeq(k, rho[N-1][N-1], u[N-1][N-1]);
     }
}

void output_screen()        //输出误差及收敛步骤
{
      if(m%100 == 0)    
      cout<<"The "<<m<<"th computation result:"
          <<endl<<"The rho,u,v of  point(M/2,N/2) is:"
          <<setprecision(7)<<rho[M/2][N/2]
          <<","<<u[M/2][N/2][0]<<","<<u[M/2][N/2][1]<<endl;
      
      if(m < 1000)
      {
       temp_rho[m] = rho[M/2][N/2];
       error_rho = fabs(rho[M/2][N/2] - temp_rho[0]);
       temp_u[m] = u[M/2][N/2][0];
       error_u = fabs(u[M/2][N/2][0] - temp_u[0]);
       temp_v[m] = u[M/2][N/2][1];
       error_v = fabs(u[M/2][N/2][1] - temp_u[0]);
       
       if(m%100 == 0)
        cout << "The error of rho,u,v is:" << setiosflags(ios::scientific)
            << error_rho << "," << error_u << "," << error_v << "." <<endl;
      }
      else
      {
        m0 = m % 1000;
        error_rho = fabs(rho[M/2][N/2] - temp_rho[0]);
        error_u = fabs(u[M/2][N/2][0] - temp_u[0]);
        error_v = fabs(u[M/2][N/2][1] - temp_u[0]);
        if(m%100==0)
          cout << "The error of rho,u,v is:" << error_rho << "," << error_u
           << "," << error_v << "." <<endl;
        temp_rho[m] = rho[M/2][N/2];
        temp_u[m] = u[M/2][N/2][0];
        temp_v[m] = u[M/2][N/2][1];
      }

      if(!(m % 2000))
        datadeal();
}

void datadeal()      //逐行输入数据
{
  FILE *fp;

  if((fp=fopen("u", "w"))==NULL)
  {
    printf("File Open error\n");
    exit(1);
  }
  for(j=0; j<=N; j++)
  {  
    for(i=0; i<=M; i++)
      fprintf(fp, "%e", u[j][i][0]);
    fprintf(fp, "\n");
  }
  fclose(fp);

  if((fp=fopen("v", "w"))==NULL)
  {
    printf("File Open error\n");
    exit(1);
  }
  for(j=0; j<=N; j++)
  {
    for(i=0; i<=M; i++)
      fprintf(fp, "%e", u[j][i][1]);
    fprintf(fp,"\n");
  }
  fclose(fp);

  if((fp=fopen("rho", "w"))==NULL)
  {
    printf("File Open error\n");
    exit(1);
  }
  for(j=0; j<=N; j++)
  {  
    for(i=0; i<=M; i++)
      fprintf(fp, "%e", rho[j][i]);
    fprintf(fp, "\n");
  }
  fclose(fp);
}