方柱绕流的lbm模块.pdfVIP

//1. 头文件及声明文件 #includeiostream #includecmath #includecstdlib #includeiomanip #includefstream #includesstream #includestring using namespace std; const int Q=9; const int NX=200; const int NY=45; const int DX=20; const int DY=19; const int H=5; const double U=0.1; int e[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}}; double w[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}; double rho[NX+1][NY+1],u[NX+1][NY+1][2],u0[NX+1][NY+1][2],f[NX+1][NY+1][Q],F[NX+1][NY+1][Q],p[NX +1][NY+1]; int i,j,k,ip,jp,n; double c,Re,dx,dy,Lx,Ly,dt,rho0,tau_f,niu,error; void init(); double feq(int k,double rho,double u[2]); void evolution(); void statistic(); void boundary(); void output(int m); void Error(); //2. 主程序文件 int main() { using namespace std; init(); for(n=0;n40000;n++) { boundary(); evolution(); statistic(); if(n%10==0) { Error(); coutThe nth conputation result:endlThe u,v of point (NX/2,NY/2)is: setprecision(6)u[NX/2][NY/2][0],u[NX/2][NY/2][1]endl; coutThe max relative error of uv is:setiosflags(ios::scientific)errorendl; output(n); if(n=1000) { if(error1.0e-6) break; } } } return 0; } //3. 初始化文件 void init() //初始函数 { dx=1.0; dy=1.0; Lx=dx*double(NX); Ly=dy*double(NY); dt=dx; c=dx/dt; rho0=1.0; Re=100; niu=U*H/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++) { u[i][j][0]=0; u[i][j][1]=0; u[0][j][0]=U; rho[i][j]=rho0; for(k=0;kQ;k++) { f[i][j][k]=feq(k,rho[i][j],u[i][j]); } } } //4. 平衡态分布函数代码 double feq(int k,double rho,double u[2]) { double eu,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); return feq; } //5. 边界条件代码 void boundary() { f[0][0][5]=f[0][0][7]; f[0][NY][8]=f[0][NY][6]; f[NX][0][6]=f[NX