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

//1.头文件及声明文件#includeiostream#includecmath#includecstdlib#includeiomanip#includefstream#includesstream#includestringusingnamespacestd;constintQ=9;constintNX=200;constintNY=45;constintDX=20;constintDY=19;constintH=5;constdoubleU=0.1;inte[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};doublew[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};doublerho[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];inti,j,k,ip,jp,n;doublec,Re,dx,dy,Lx,Ly,dt,rho0,tau_f,niu,error;voidinit();doublefeq(intk,doublerho,doubleu[2]);voidevolution();voidstatistic();voidboundary();voidoutput(intm);voidError();//2.主程序文件intmain(){usingnamespacestd;init();for(n=0;n40000;n++){evolution();statistic();boundary();if(n%10==0){Error();coutThenthconputationresult:endlTheu,vofpoint(NX/2,NY/2)is:setprecision(6)u[NX/2][NY/2][0],u[NX/2][NY/2][1]endl;coutThemaxrelativeerrorofuvis:setiosflags(ios::scientific)errorendl;output(n);if(n=1000){if(error1.0e-6)break;}}}return0;}//3.初始化文件voidinit()//初始函数{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::couttau_f=tau_fendl;for(i=0;i=NX;i++)for(j=0;j=NY;j++){u[i][j][0]=0;u[i][j][1]=0;u[0][j][0]=U;rho[i][j]=rho0;for(k=0;kQ;k++){f[i][j][k]=feq(k,rho[i][j],u[i][j]);}}}//4.平衡态分布函数代码doublefeq(intk,doublerho,doubleu[2]){doubleeu,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);returnfeq;}//5.碰撞以及迁移代码voidevolution(){for(i=1;iNX;i++)//边界不参与碰撞（边界条件除外）for(j=1;jNY;j++){if(i=DX&&i=DX+H&&j=DY&&j=DY+H)continue;for(k=0;kQ;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;}}}//6.宏观统计量代码voidstatistic(){for(i=1;iNX;i++)//不统计边界上的点for(j=1;jNY;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;kQ;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;//计算压力}}//7.边界条件代码voidboundary(){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;iNX;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;iDX+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];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.输出文件代码voidoutput(intm){ostringstreamname;nameobstacle_m.dat;ofstreamout(name.str().c_str());outTitle=\YUANZHURAOLIU\\nVARIABLES=\X\,\Y\,\U\,\V\,\P\,\UV\\nZONET=\BOX\,I=NY+1,J=NX+1,F=POINTendl;for(i=0;i=NX;i++)for(j=0;j=NY;j++){outiju[i][j][0]u[i][j][1]p[i][j]uv[i][j]endl;}}//9.判定条件代码voidError(){doubletemp1,temp2;temp1=0;temp2=0;for(i=1;iNX;i++)for(j=1;jNY;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);temp2=sqrt(temp2);error=temp1/(temp2+1e-30);}