方柱绕流的LBM代码(C++)

//9. 判定条件代码
void Error() { double temp1,temp2; temp1=0; temp2=0; for(i=1;i<NX;i++) for(j=1;j<NY;j++) { if(i>=DX&&i<=DX+H&&j>=DY&&j<=DY+H) continue; temp1 +=( (u[i][j][0]-u0[i][j][0])*(u[i][j][0]-u0[i][j][0]) +(u[i][j][1]-u0[i][j][1])*(u[i][j][1]-u0[i][j][1])); temp2 +=(u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]); } temp1=sqrt(temp1);
//2. 主程序文件
int main() { using namespace std;
init(); for(n=0;n<40000;n++) { evolution(); statistic(); boundary(); if(n%10==0) { Error(); cout<<"The "<<n<<"th conputation 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; output(n); if(n>=1000) { if(error<1.0e-6) break; } } } return 0; }
//7. 边界条件代码
void boundary() { for(j=0;j<=NY;j++)//入口采用 ZOU/HE 速度边界 { i=0; u[i][j][0]=U; u[i][j][1]=0; rho[i][j]=(f[i][j][0]+f[i][j][2]+f[i][j][4]+2.0*(f[i][j][3]+f[i][j][6]+f[i][j][7]))/(1-u[i][j][0]); f[i][j][1]=f[i][j][3]+2.0*rho[i][j]*u[i][j][0]/3.0; f[i][j][5]=f[i][j][7]+rho[i][j]*u[i][j][0]/6.0-0.5*(f[i][j][2]-f[i][j][4])+rho[i][j]*u[i][j][1]/2.0;
f[i][j][8]=f[i][j][6]+rho[i][j]*u[i][j][0]/6.0+0.5*(f[i][j][2]-f[i][j][4])-rho[i][j]*u[i][j][1]/2.0; } for(j=0;j<=NY;j++)//出口采用 ZOU/HE 压力边界 { i=NX; u[i][j][0]=0; u[i][j][1]=0; rho[i][j]=1.0; u[i][j][0]=(f[i][j][0]+f[i][j][2]+f[i][j][4]+2*(f[i][j][1]+f[i][j][5]+f[i][j][8]))/rho[i][j]-1.0; f[i][j][3]=f[i][j][1]-2.0*rho[i][j]*u[i][j][0]/3.0; f[i][j][7]=f[i][j][5]-rho[i][j]*u[i][j][0]/6.0+0.5*(f[i][j][2]-f[i][j][4]); f[i][j][6]=f[i][j][8]-rho[i][j]*u[i][j][0]/6.0-0.5*(f[i][j][2]-f[i][j][4]); } for(i=1;i<NX;i++) //上下边界采用标准反弹 { j=0; f[i][j][2]=f[i][j][4]; f[i][j][5]=f[i][j][7]; f[i][j][6]=f[i][j][8]; j=NY; f[i][j][4]=f[i][j][2]; f[i][j][7]=f[i][j][5]; f[i][j][8]=f[i][j][6]; } for(i=DX+1;i<DX+H;i++) //圆柱上下边界采用标准反弹 { j=DY+H; f[i][j][2]=f[i][j][4]; f[i][j][5]=f[i][j][7]; f[i][j][6]=f[i][j][8]; j=DY; f[i][j][4]=f[i][j][2]; f[i][j][7]=f[i][j][5]; f[i][j][8]=f[i][j][6]; } for(j=DY;j<=DY+H;j++) //圆柱左右边界采用标准反弹 { i=DX; f[i][j][3]=f[i][j][1]; f[i][j][6]=f[i][j][8];
//6. 宏观统计量代码
void statistic() { for(i=1;i<NX;i++) //不统计边界上的点 for(j=1;j<NY;j++) { if(i>=DX&&i<=DX+H&&j>=DY&&j<=DY+H) continue;//抛去圆柱边界格点 u0[i][j][0]=u[i][j][0]; u0[i][j][1]=u[i][j][1]; //不需要收敛判断时可去除 rho[i][j]=0; u[i][j][0]=0; u[i][j][1]=0; for(k=0;k<Q;k++) { f[i][j][k]=F[i][j][k]; rho[i][j]+=f[i][j][k]; //求密度 u[i][j][0]+=e[k][0]*f[i][j][k]; u[i][j][1]+=e[k][1]*f[i][j][k]; } u[i][j][0]/=rho[i][j]; u[i][j][1]/=rho[i][j]; //求 xy 方向的速度分量 uv[i][j]=sqrt(u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]); p[i][j]=rho[i][j]*c*c/3;//计算压力 } }
//3. 初始化文件
void init() //初始函数 { dx=1.0; dy=1.0; Lx=dx*double(NX); Ly=dy*double(NY); dt=dx; c=dx/dt; rho0=1.0; Re=100; niu=U*H/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[i][j][0]=0; u[i][j][1]=0;
//1. 头文件及声明文件
#include<iostream> #include<cmath> #include<cstdlib> #include<iomanip> #include<fstream> #include<sstream> #include<string> using namespace std; const int Q=9; const int NX=200; const int NY=45; const int DX=20; const int DY=19; const int H=5; 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],p[NX +1][NY+1],uv[NX+1][NY+1]; int i,j,k,ip,jp,n; double c,Re,dx,dy,Lx,Ly,dt,rho0,tau_f,niu,error; void init(); double feq(int k,double rho,double u[2]); void evolution(); void statistic(); void boundary(); void output(int m); void Error();
u[0][j][0]=U; rho[i][j]=rho0; for(k=0;k<Q;k++) { f[i][j][k]=feq(k,rho[i][j],u[i][j]); } } }
//4. 平衡态分布函数代码
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; }
temp2=sqrt(temp2); error=temp1/(temp2+1e-30); }
f[i][j][7]=f[i][j][5]; i=DX+H; f[i][j][1]=f[i][j][3]; f[i][j][5]=f[i][j][6]; f[i][j][8]=f[i][j][7]; } }
//8. 输出文件代码
void output(int m) { ostringstream name; name<<"obstacle_"<<m<<".dat"; ofstream out(name.str().c_str()); out<<"Title=\"YUANZHURAOLIU\"\n"<<"VARIABLES=\"X\",\"Y\",\"U\",\"V\",\"P\",\"UV\"\n"<<"Z ONE T=\"BOX\",I=" <<NY+1<<",J=来自百度文库<<NX+1<<",F=POINT"<<endl; for(i=0;i<=NX;i++) for(j=0;j<=NY;j++) { out<<i<<" "<<j<<" "<<u[i][j][0]<<" "<<u[i][j][1]<<" "<<p[i][j]<<" "<<uv[i][j]<<endl; } }
//5. 碰撞以及迁移代码
void evolution() { for(i=1;i<NX;i++) //边界不参与碰撞（边界条件除外） for(j=1;j<NY;j++) { if(i>=DX&&i<=DX+H&&j>=DY&&j<=DY+H) continue; for(k=0;k<Q;k++) { ip=i-e[k][0]; jp=j-e[k][1]; F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau_f; } } }