C++LBM顶盖驱动流

1、这一段代码源自何雅玲的格子BOLTZMANN法的理论及应用其后面源代码，经过调试可用，输出文件可以导入到TECPLO冲。使用时，在VISUALSTUDIO中新建一个C+C程，代码赋值进去调试即可。#include<iostream>#include<cmath>#include<cstdlib>#include<iomanip>#include<fstream>#include<sstream>#include<string>usingnamespacestd;constintQ=9;/D2Q9模型consti

2、ntNX=256;/X方向constintNY=256;/Y方向constdoubleU=0.1;/顶盖速度inteQ2=0,0,1,0,0,1,-1,0,0,-1,1,1,-1,1,-1,-1,1,-1;doublewQ=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;doublerhoNX+1NY+1,uNX+1NY+12,u0NX+1NY+12,fNX+1NY+1Q,FNX+1NY+1Q;inti,j,k,ip,jp,n;doublec,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error;

3、voidinit();doublefeq(intk,doublerho,doubleu2);voidevolution();voidoutput(intm);voidError();intmain()usingnamespacestd;init();for(n=0;n+)evolution();if(n%100=0)Error();cout<<"The"<<n<<"thcomputationresult:"<<endl<<"Theu,vofpoint(NX/2,NY/2)is:&quo

4、t;<<setprecision(6)<<uNX/2NY/20<<","<<uNX/2NY/21<<endl;cout<<"Themaxrelativeerrorofuvis:"<<setiosflags(ios:scientific)<<error<<endl;if(n>=1000)(if(n%1000=0)output(n);if(error<1.0e-6)break;return0;voidinit()(dx=1.0;dy=1.0

5、;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_f<<endl;for(i=0;i<=NX;i+)/初始化for(j=0;j<=NY;j+)(uij0=0;uij1=0;rhoij=rho0;uiNY0=U;for(k=0;k<Q;k+)(fijk=feq(k,rhoij,uij);doublefeq(intk,double

6、rho,doubleu2)(doubleeu,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);returnfeq;voidevolution()for(i=1;i<NY;i+)演化for(j=1;j<NY;j+)for(k=0;k<Q;k+)ip=i-ek0;jp=j-ek1;Fijk=fipjpk+(feq(k,rhoipjp,uipjp)-fipjpk)/tau_f;for(i=1;i<NX;i+)/计算宏观量for(j=1;j<NY;j+)

7、u0ij0=uij0;u0ij1=uij1;rhoij=0;uij0=0;uij1=0;for(k=0;k<Q;k+)fijk=Fijk;rhoij+=fijk;uij0+=ek0*fijk;uij1+=ek1*fijk;uij0/=rhoij;uij1/=rhoij;/边界处理for(j=1;j<NY;j+)/左右边界for(k=0;k<Q;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(

8、k,rho1j,u1j);)for(i=0;i<=NX;i+)/上下边界for(k=0;k<Q;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);)voidoutput(intm)/输出ostringstreamname;namev&avity'<<m<<.dat”;ofstreamout(name.str().c_s

9、tr();out<<"Title="LBMLidDrivenFlow"n"'VARIABLES="X","Y","U","V"n"<VZONET="BOX",I="<<NX+1<<,J="<<NY+1<<,F=POINT"<<endl;for(j=0;j<=NY;j+)for(i=0;i<=NX;i+)out<<double(i)/Lx<<""<<double(j)/Ly<<""<<uij0<<""<<uij1<<endl;)voidError()doubletemp1,temp2;temp1=0;temp2=0;for(i=1;i<NX;i+)for(j=1;j<NY;j+)tempi+=(uij0-u0ij0)*(uij0-u0ij0)+(uij1-u0ij1)*(u