dgqdl_三维_C++_方腔_LBM_源码

// 三维方腔.cpp : Defines the entry point for the console application. // #include<iostream> #include<cmath> #include<cstdlib> #include<iomanip> #include<fstream> #include<sstream> #include<string> using namespace std; //变量名定义 const int Q=15; const int NX=40;//尺度 const int NY=40; const int NZ=40; const int Re=400; const double RHO=1.0; const double U=0.1; double niu=U*NX/Re; double tau=3.0*niu+0.5; int e[Q][3]={ {0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},{1,1,1},{-1,-1,-1},{1,1,-1},{-1,-1,1},{1,-1,1},{-1,1,-1},{1,-1,-1},{-1,1,1} }; double w[Q]={2.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/72,1.0/72,1.0/72,1.0/72,1.0/72,1.0/72,1.0/72,1.0/72}; //速度场权函数 double ux[NX][NY][NZ],uy[NX][NY][NZ],uz[NX][NY][NZ]; //速度 double ux0[NX][NY][NZ],uy0[NX][NY][NZ],uz0[NX][NY][NZ]; double rho[NX][NY][NZ],F[NX][NY][NZ][Q],f[NX][NY][NZ][Q]; //速度场相关 void init(); double feq(int i,int j,int k,int a); void evolution(); void outputuxy(int n); void outputuyz(int n); void outputuxz(int n); int i,j,k,a,ip,jp,kp,n; void main() //主函数 { using namespace std; init(); for(n=0;n<5000;n++) { evolution(); } outputuxz(n); } void init() //初始化 { double feq(int i,int j,int k,int a); for(i=0;i<NX;i++) for(j=0;j<NY;j++) for(k=0;k<NZ;k++) { ux[i][j][k]=0.0; uy[i][j][k]=0.0; uz[i][j][k]=0.0; rho[i][j][k]=RHO; ux[i][j][NZ-1]=U; for(a=0;a<Q;a++) //速度场 { f[i][j][k][a]=feq(i,j,k,a); } } } double feq(int i,int j,int k,int a) //速度场平衡态分布函数 { double eu,u2,feq; eu=e[a][0]*ux[i][j][k]+e[a][1]*uy[i][j][k]+e[a][2]*uz[i][j][k]; u2=ux[i][j][k]*ux[i][j][k]+uy[i][j][k]*uy[i][j][k]+uz[i][j][k]*uz[i][j][k]; feq=w[a]*rho[i][j][k]*(1.0+3.0*eu+4.5*eu*eu-1.5*u2); return feq; } void evolution() //演化函数 { for(i=1;i<NX-1;i++) //速度场格子方程 for(j=1;j<NY-1;j++) for(k=1;k<NZ-1;k++) for(a=0;a<Q;a++) { ip=i-e[a][0]; jp=j-e[a][1]; kp=k-e[a][2]; F[i][j][k][a]=f[ip][jp][kp][a]+( feq(ip,jp,kp,a)-f[ip][jp][kp][a] )/tau; } for(i=1;i<NX-1;i++) for(j=1;j<NY-1;j++) for(k=1;k<NZ-1;k++) { ux0[i][j][k]=ux[i][j][k]; uy0[i][j][k]=uy[i][j][k]; uz0[i][j][k]=uz[i][j][k]; rho[i][j][k]=0.0; ux[i][j][k]=0.0; uy[i][j][k]=0.0; uz[i][j][k]=0.0; for(a=0;a<Q;a++) { f[i][j][k][a]=F[i][j][k][a]; rho[i][j][k]+=f[i][j][k][a]; ux[i][j][k]+=e[a][0]*f[i][j][k][a]; uy[i][j][k]+=e[a][1]*f[i][j][k][a]; uz[i][j][k]+=e[a][2]*f[i][j][k][a]; } ux[i][j][k]/=rho[i][j][k]; uy[i][j][k]/=rho[i][j][k]; uz[i][j][k]/=rho[i][j][k]; } for(j=0;j<NY;j++) for(k=0;k<NZ;k++) for(a=0;a<Q;a++) { rho[NX-1][j][k]=rho[NX-2][j][k]; f[NX-1][j][k][a]=feq(NX-1,j,k,a)+f[NX-2][j][k][a]-feq(NX-2,j,k,a); rho[0][j][k]=rho[1][j][k]; f[0][j][k][a]=feq(0,j,k,a)+f[1][j][k][a]-feq(1,j,k,a); } for(i=0;i<NX;i++) for(k=0;k<NZ;k++) for(a=0;a<Q;a++) { rho[i][0][k]=rho[i][1][k]; f[i][0][k][a]=feq(i,0,k,a)+f[i][1][k][a]-feq(i,1,k,a); rho[i][NY-1][k]=rho[i][NY-2][k]; f[i][NY-1][k][a]=feq(i,NY-1,k,a)+f[i][NY-2][k][a]-feq(i,NY-2,k,a); } for(i=0;i<NX;i++) for(j=0;j<NY;j++) for(a=0;a<Q;a++) { rho[i][j][0]=rho[i][j][1]; f[i][j][0][a]=feq(i,j,0,a)+f[i][j][1][a]-feq(i,j,1,a); rho[i][j][NZ-1]=rho[i][j][NZ-2]; ux[i][j][NZ-1]=U; f[i][j][NZ-1][a]=feq(i,j,NZ-1,a)+f[i][j][NZ-2][a]-feq(i,j,NZ-2,a); } } void outputuxy(int n) //二维速度输出函数 { ostringstream name; name<<"uxy_"<<n<<".dat"; ofstream out(name.str().c_str()); out<<"title=\"LBM Lid Driven Flow\"\n"<<"VARIABLES=\"X\",\"Y\",\"U\",\"V\"\n"<<"ZONE T=\"BOX\",I="<<NX<<",J="<<NY<<",F=POINT"<<endl; for(j=0;j<NY;j++) for(i=0;i<NX;i++) { out<<double(i)/NX<<" "<<double(j)/NY<<" "<<ux[i][j][NZ/2]<<" "<<uy[i][j][NZ/2]<<endl; } } void outputuxz(int n) //二维速度输出函数 { ostringstream name; name<<"uxz_"<<n<<".dat"; ofstream out(name.str().c_str()); out<<"title=\"LBM Lid Driven Flow\"\n"<<"VARIABLES=\"X\",\"Z\",\"U\",\"W\"\n"<<"ZONE T=\"BOX\",I="<<NX<<",K="<<NZ<<",F=POINT"<<endl; for(k=0;k<NZ;k++) for(i=0;i<NX;i++) { out<<double(i)/NX<<" "<<double(k)/NZ<<" "<<ux[i][NY/2][k]<<" "<<uz[i][NY/2][k]<<endl; } } void outputuyz(int n) //二维速度输出函数 { ostringstream name; name<<"uyz_"<<n<<".dat"; ofstream out(name.str().c_str()); out<<"title=\"LBM Lid Driven Flow\"\n"<<"VARIABLES=\"Y\",\"Z\",\"V\",\"W\"\n"<<"ZONE T=\"BOX\",J="<<NY<<",K="<<NZ<<",F=POINT"<<endl; for(k=0;k<NZ;k++) for(j=0;j<NY;j++) { out<<double(j)/NY<<" "<<double(k)/NZ<<" "<<uy[NX/2][j][k]<<" "<<uz[NX/2][j][k]<<endl; } } /* for(j=0;j<=NY;j++) //速度场边界处理（非平衡）(yz面) for(k=0;k<=NZ;k++) for(a=0;a<Qu;a++) { p[NX][j][k]=p[NX-1][j][k]; f[NX][j][k][a]=feq_u(NX,j,k,a)+f[NX-1][j][k][a]-feq_u(NX-1,j,k,a); p[0][j][k]=p[1][j][k]; f[0][j][k][a]=feq_u(0,j,k,a)+f[1][j][k][a]-feq_u(1,j,k,a); } for(i=1;i<NX;i++) //(xz面) for(k=1;k<NZ;k++) for(a=0;a<Qu;a++) { p[i][0][k]=p[i][1][k]; f[i][0][k][a]=feq_u(i,0,k,a)+f[i][1][k][a]-feq_u(i,1,k,a); p[i][NY][k]=p[i][NY-1][k]; f[i][NY][k][a]=feq_u(i,NY,k,a)+f[i][NY-1][k][a]-feq_u(i,NY-1,k,a); } for(i=1;i<NX;i++) //(xy面) for(j=0;j<=NY;j++) for(a=0;a<Qu;a++) { p[i][j][0]=p[i][j][1]; f[i][j][0][a]=feq_u(i,j,0,a)+f[i][j][1][a]-feq_u(i,j,1,a); p[i][j][NZ]=p[i][j][NZ-1]; f[i][j][NZ][a]=feq_u(i,j,NZ,a)+f[i][j][NZ-1][a]-feq_u(i,j,NZ-1,a); } */