C++LBM顶盖驱动流

1、精品文档1欢迎下载这一段代码源自何雅玲的格子BOLTZMANN 法的理论及应用其后面源代码，经过调试可用，输出文件可以导入到 TECPLO中。使用时，在 VISUAL STUDIO 中新建一个 C+工程，代码赋 值进去调试即可。 using 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,

2、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 namesp

3、acestd;init();for (n=0;n+)evolution();if (n%100=0)Error();cout The nth computation result: endl The u,v of point(NX/2,NY/2) is:setprecision(6)uNX/2NY/20, uNX/2NY/21endl;cout The max relative error of uvis: setiosflags(ios:scientific)error=1000)if (n%1000=0) output(n);if (error1.0e-6) break ;return 0

4、;#include#include#include#include#include#include#include精品文档2欢迎下载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:cout tau_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=

5、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)-fipj

6、pk)/tau_f;精品文档3欢迎下载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,u

7、NX-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 ;of

8、stream out(name.str().c_str();精品文档4欢迎下载out Title=LBM Lid Driven FlownVARIABLES = X,Y,U,VnT=BOX,I= NX+1,J= NY+1,F=POINT endl;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-u0ij0)+(uij1-u0ij1)*(uij1- u0ij1);temp2 +=(uij0*uij0+uij1*uij1