C++LBM頂蓋驅(qū)動流

1、這一段代碼源自何雅玲的格子BOLTZMANN法的理論及應(yīng)用其后面源代碼，經(jīng)過調(diào)試可用，輸出文件可以導(dǎo)入到 TECPLO中。使用時，在 VISUAL STUDIO中新建一個 C+工程，代碼賦 值進去調(diào)試即可。#include#include#include #include #include #include #include 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,

2、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();vo

3、id output( int m);void Error();int main()using namespacestd;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

4、 (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: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=0;kQ;k+)fi

5、jk=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

6、;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=rh

7、o1j;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(na

8、me.str().c_str();ZONEout 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*uij