C++LBM顶盖驱动流

这一段代码源自何雅玲的格子BOLTZMANN方法的理论及应用其后面源代码，经过调试可用，输出文件可以导入到TECPLOT中。使用时，在VISUAL STUDIO中新建一个C+工程，代码赋值进去调试即可。#include#include#include#include#include#include#includeusing namespace std; const int Q=9; /D2Q9模型const int NX=256; /X方向const int NY=256; /Y方向const double U=0.1; /顶盖速度int eQ2=0,0,1,0,0,1,-1,0,0,-1,1,1,-1,1,-1,-1,1,-1;double wQ=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 rhoNX+1NY+1,uNX+1NY+12,u0NX+1NY+12,fNX+1NY+1Q,FNX+1NY+1Q;int i,j,k,ip,jp,n;double c,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error;void init();double feq(int k,double rho,double u2);void evolution();void output(int m);void Error();int main()using namespace std;init();for(n=0;n+)evolution();if(n%100=0)Error();coutThenth computation result:endlThe u,v of point(NX/2,NY/2) is:setprecision(6)uNX/2NY/20,uNX/2NY/21endl;coutThe max relative error of uv is:setiosflags(ios:scientific)error=1000)if(n%1000=0) output(n);if(error1.0e-6) break;return 0;void init()dx=1.0;dy=1.0;Lx=dx*double(NY);Ly=dy*double(NX);dt=dx;c=dx/dt; /1.0rho0=1.0;Re=1000;niu=U*Lx/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+)uij0=0;uij1=0;rhoij=rho0;uiNY0=U;for(k=0;kQ;k+)fijk=feq(k,rhoij,uij);double feq(int k,double rho,double u2) double eu,uv,feq;eu=(ek0*u0+ek1*u1);uv=(u0*u0+u1*u1);feq=wk*rho*(1.0+3.0*eu+4.5*eu*eu-1.5*uv);return feq;void evolution()for(i=1;iNY;i+) /演化for(j=1;jNY;j+)for(k=0;kQ;k+)ip=i-ek0;jp=j-ek1;Fijk=fipjpk+(feq(k,rhoipjp,uipjp)-fipjpk)/tau_f;for(i=1;iNX;i+) /计算宏观量for(j=1;jNY;j+)u0ij0=uij0;u0ij1=uij1;rhoij=0;uij0=0;uij1=0;for(k=0;kQ;k+)fijk=Fijk;rhoij+=fijk;uij0+=ek0*fijk;uij1+=ek1*fijk;uij0/=rhoij;uij1/=rhoij;/边界处理for(j=1;jNY;j+) /左右边界for(k=0;kQ;k+)rhoNXj=rhoNX-1j;fNXjk=feq(k,rhoNXj,uNXj)+fNX-1jk-feq(k,rhoNX-1j,uNX-1j);rho0j=rho1j;f0jk=feq(k,rho0j,u0j)+f1jk-feq(k,rho1j,u1j);for(i=0;i=NX;i+)/上下边界for(k=0;kQ;k+)rhoi0=rhoi1;fi0k=feq(k,rhoi0,ui0)+fi1k-feq(k,rhoi1,ui1); rhoiNY=rhoiNY-1; uiNY0=U; fiNYk=feq(k,rhoiNY,uiNY)+fiNY-1k-feq(k,rhoiNY-1,uiNY-1);void output(int m) /输出ostringstream name;namecavity_m.dat;ofstream out(name.str().c_str();outTitle=LBM Lid Driven FlownVARIABLES = X,Y,U,VnZONE T=BOX,I=NX+1,J=NY+1,F=POINTendl;for(j=0;j=NY;j+)for(i=0;i=NX;i+)outdouble(i)/Lx double(j)/Ly uij0 uij1endl;void Error()double temp1,temp2;temp1=0;temp2=0;for(i=1;iNX;i+)for(j=1;jNY;j+)temp1 +=(uij0-u0ij0)*(uij0-u0i