Rayleigh-Benard对流格子Boltzmann代码

//naturalconvectionforclosedsqurecavity.cpp:Definestheentrypointfortheconsoleapplication.//coupleddoubledistributionfountion#includestdafx.h#includeiostream#includecmath#includecstdlib#includeiomanip#includefstream#includesstream#includestringusingnamespacestd;//===============================================Constantdeclarations==================================================constintQ=9;//D2Q9modelconstintNX=100;//latticeofthexdirectionconstintNY=50;//latticeoftheydirectionconstdoubleT_hot=2.;//ThebottomwalltemperatureconstdoubleT_cold=1.;//ThetopwalltemperaturedoubleGg[2]={0,-0.005};//accelerationofgravity//================================================Variabledeclaration====================================================//ThediscretevelocitydirectionofD2Q9modeldoublee[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};//weightcoefficientofD2Q9modeldoublew[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],//densityu[NX+1][NY+1][2],//velocityofn+1timeu0[NX+1][NY+1][2],//velocityofntimef[NX+1][NY+1][Q],//velocitydistributionfuntionbeforeevolutionF[NX+1][NY+1][Q],//velocitydietributionfuntionafterevolutionvel[NX+1][NY+1];//velocityof(i.j)doubleT[NX+1][NY+1],//tempertureofn+1timeg[NX+1][NY+1][Q],//temperturedistributionfuntionbeforeevolutionG[NX+1][NY+1][Q];//temperturedistributionfuntionafterevolutiondoublewf,wg;inti,j,k,ip,jp,n;doubleLX,LY,c,Fa,Ra,Prandtl,Rg=8.3145,dx,dy,dt,rho0,P0,tau_f,tau_g,niu,error,beita,Nu,Tav;//==================================================funtionsofthiscode===================================================voidinit();doublefeq(intk,doublerho,doubleu[2]);doublegeq(intk,doubleT,doubleu[2]);voidevolution();voidoutput(intm);voidError();//=======================================================mainloop============================================================intmain(){usingnamespacestd;init();for(n=0;;n++){evolution();if(n%5000==0){Error();coutThenerrorsetiosflags(ios::scientific)errorNU=NuTav=Tavendl;if(n=20000){if(n%20000==0)output(n);if(error1.0e-10)break;}}}return0;}//======================================================Theinitialization====================================================voidinit(){dy=1;//Latticestepofxdirectiondt=dy;//timestepLY=double(NY)*dy;LX=2*LY;c=dy/dt;//Thelatticevelocityrho0=1.0;//initialdensityRa=5.0e3;//RayleighnumberPrandtl=0.71;//Prandtlnumberniu=1.e-5;//Rg=1/(3*(T_hot+T_cold)/2);//kinematicviscositytau_f=0.1*NY*pow(3*Prandtl,0.5)/pow(Ra,0.5)+0.5;//relaxationtimeofvelocityfieldtau_g=(tau_f-0.5)/Prandtl+0.5;//relaxationtimeoftemperturefieldwf=dt/(tau_f+0.5*dt);wg=dt/(tau_g+0.5*dt);beita=Ra*niu/(abs(Gg[1])*(T_hot-T_cold)*NY*NY*NY*Prandtl);for(i=0;i=NX;i++){for(j=0;j=NY;j++){u[i][j][0]=0;u[i][j][1]=0;rho[i][j]=rho0;T[i][j]=(T_hot+T_cold)/2.;for(k=0;kQ;k++){f[i][j][k]=feq(k,rho[i][j],u[i][j]);//initializationdistributionfunctiong[i][j][k]=geq(k,T[i][j],u[i][j]);//velocityfieldandtemperture}}}}//Thevelocityequilibriumdistributionfunctiondoublefeq(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;}//Thetempretureequilibriumdistributionfunctiondoublegeq(intk,doubleT,doubleu[2]){doubleeu,uv,geq;eu=(e[k][0]*u[0]+e[k][1]*u[1]);uv=(u[0]*u[0]+u[1]*u[1]);geq=w[k]*T*(1.0+3.0*eu+4.5*eu*eu-1.5*uv);returngeq;}//evolution:collision,streaming,macro,boundary;voidevolution(){for(i=1;iNX;i++){for(j=1;jNY;j++){for(k=0;kQ;k++){ip=i-int(e[k][0]);//collisionandstreamingjp=j-int(e[k][1]);////evolutionequationofvelocityfield//Fa=(1-1/(2*tau_f))*w[k]*(3*(e[k][1]-u[ip][jp][1])+9*e[k][1]*u[ip][jp][1]*e[k][1])*(beita*(T[ip][jp]-(T_hot+T_cold)/2)*Gg[1]);Fa=dt*Gg[1]*(e[k][1]-u[ip][jp][1])*feq(k,rho[ip][jp],u[ip][jp])/Rg/T[ip][jp];F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau_f+dt*Fa;//evolutionequationoftemperturefieldG[i][j][k]=g[ip][jp][k]+(geq(k,T[ip][jp],u[ip][jp])-g[ip][jp][k])/tau_g;}}}//Themicroscopicdistributionfunctiontothemacroscopicphysicalquantityfor(i=1;iNX;i++){for(j=1;jNY;j++){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;T[i][j]=0;for(k=0;kQ;k++){f[i][j][k]=F[i][j][k];rho[i][j]+=f[i][j][k];g[i][j][k]=G[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];T[i][j]+=g[i][j][k];}u[i][j][0]/=rho[i][j];u[i][j][1]/=rho[i][j];vel[i][j]=sqrt(u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]);}}Nu=0.0;//computationalNufor(i=0;i=NX;i++){Nu=Nu+(T[i][0]-T[i][1]);}Tav=0.0;for(i=0;i=NX;i++){for(j=0;j=NY;j++){Tav=Tav+(T[i][j]-1);}}Tav=Tav/((NX+1)*(NY+1));//---------------------------------------------------------boundarytreatment---------------------------------------------------//Thenon-equilibriumextrapolationschemeforboundarytreatmentandperiodicboundary//leftandrightperiodicfor(j=0;j=NY;j++){for(k=0;kQ;k++){rho[NX][j]=rho[NX-1][j];T[NX][j]=T[NX-1][j];if(k==1||k==5||k==8){f[0][j][k]=f[NX][j][k];g[0][j