24小时热门版块排行榜

返回列表

花开时节0931

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 21
帖子: 19
在线: 5.8小时
虫号: 3159261
注册: 2014-04-24

[求助] 求助格子boltzmann方法模拟平板间气体流动已有1人参与

空气以声速从平板间进入，平板厚度10微米，想采用格子boltzmann模拟气体流动，没学过C语言，边学边改。问题还是很大，毕业论文着急，求指导。不胜感激。
#include <cmath>
#include <iostream>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
using namespace std;
const int Q=9; //D2Q9模型
const int NX=120;//X方向
const int NY=40;//Y方向
//const double U=0.1; //顶盖速度
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};
double rho[NX+1][NY+1],u[NX+1][NY+1][2],u0[NX+1][NY+1][2],f[NX+1][NY+1][Q],F[NX+1][NY+1][Q];
int i,j,k,ip,jp,n;
double c,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error,Kn;
void init();
double feq(int k,double rho,double u[2]);
double rhosum[NX][NY] = {0};
double Psum[NX][NY] = {0};

void evolution();
void output(int m) ;
void Error();
int main()
{
using namespace std;
init();
for(n=0; ;n++)
{
      evolution();
      if (n%1000==0)
      {
         Error();
         cout<<"the "<<n<<" th computation result:"<<endl<<
         "The u,v of point (NX/2,NY/2)is:" <<setprecision(6)
         <<u[NX/2][NY/2][0]<<","<<u[NX/2][NY/2][1]<<endl;
         cout<<"The max relative error of uv is: "
         <<setiosflags(ios::scientific)<< error <<endl;
         if(n>=100)
         {
            if(n%1000==0) output(n);
            if(error<1.0e-12) break;
         }
      }
}
return 0;
}
void init()
{
  dx=1.0;
  dy=1.0;
  P0=0.5;
  Lx=dx*double(NY);
  Ly=dy*double(NX);
  dt=dx;
  c=dx/dt;//1.0
  rho0=1.03;
  Kn=0.5;
  tau_f=8;
//tau_f=P0*Lx*Kn*sqrt(6/3.14)+0.5;
  //Re=8;
// niu=0.1*Lx/Re;
//tau_f=3.0*niu+0.5;

  std::cout<<"tau_f="<<tau_f<<endl;

  for(i=0;i<=NX;i++)//初始化
  for(j=0;j<=NY;j++)
  {
      u[j][0]=0;
      u[j][1]=0;
      rho[j]=rho0;
      u[0][j][0]=0.57;
      for (k=0;k<Q;k++)
      {
         f[j][k]=feq(k,rho[j],u[j]);

      }
  }
}
double feq(int k,double rho,double u[2])//计算平衡态分布函数
{
   double eu,uv,feq;
   eu=(e[k][0]*u[0]+e[k][1]*u[1]);
   uv=(u[0]*u[0]+u[1]*u[1]);
   feq=w[k]*rho*(1.0+3.0*eu+4.5*eu*eu-1.5*uv);
   return feq;
}

void evolution()
{
for(i=1;i<NX;i++)//演化
for(j=1;j<NY;j++)
for(k=0;k<Q;k++)
{
   ip=i-e[k][0];
   jp=j-e[k][1];
   F[j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],
   u[ip][jp])-f[ip][jp][k])/tau_f;
}

for(i=1;i<NX;i++)//计算宏观量
  for(j=1;j<NY;j++)
  {
   u0[j][0]=u[j][0];
   u0[j][1]=u[j][1];
   rho[j]=0;
   u[j][0]=0;
   u[j][1]=0;
   for(k=0;k<Q;k++)
   {
         f[j][k]=F[j][k];
         rho[j]+=f[j][k];
         u[j][0]+=e[k][0]*f[j][k];
         u[j][1]+=e[k][1]*f[j][k];
   rhosum[j] += f[j][k];
   }
   Psum[j] = rhosum[j] * c * c / 3;
   u[j][0]/=rho[j];
   u[j][1]/=rho[j];
}

//边界处理
//左右边界


for(i=0;i<=NX;i++)//上下边界
  for(k=0;k<Q;k++)
  {
   rho[0]=rho[1];
   f[0][k]=feq(k,rho[0],u[0])+f[1][k]
   -feq(k,rho[1],u[1]);
   rho[NY]=rho[NY-1];
   u[NY][0]=0;
   u[NY][1]=0;
   u[0][0]=0;
   u[0][1]=0;

   f[NY][k]=feq(k,rho[NY],u[NY])+f[NY-1]
   [k]-feq(k,rho[NY-1],u[NY-1]);
  }
}

void output(int m)//输出
{
   ostringstream name;
   name<<"cavity_"<< m<<".dat";
   ofstream out(name.str().c_str());
   out<<"Title=\"LBM Lid Driven Flow\"\n"<<
   "VARIABLES=\"X\",\"Y\", \"U\",\"V\",\"rhosum\",\"Psum\"\n"
   <<"ZONE T=\"BOX\",I="<<NX+1<<",J="<<NY+1<<",F=POINT"<<endl;
   for(j=0;j<=NY;j++)
   for(i=0;i<=NX;i++)
   {
      out<<
      double(i)/Lx<<","<<double(j)/Ly<<","<<u[j][0]<<","<<u[j][1]<<","<<rhosum[j]<<","<<Psum[j]<<"\n"
      <<endl;
   }
}
void Error()
{
   double temp1,temp2;
   temp1=0;
   temp2=0;
   for(i=1;i<NX;i++)
      for(j=1;j<NY;j++)
      {
         temp1+=(
         (u[j][0]-u0[j][0])*(u[j][0]-u0[j][0])
         +(u[j][1]-u0[j][1])*(u[j][1]-u0[j][1]));
            temp2+=
         (u[j][0]*u[j][0]+u[j][1]*u[j][1]);

      }
temp1=sqrt(temp1);
temp2=sqrt(temp2);
error=temp1/(temp2);
}