方柱绕流的lbm代码(c++).pdfVIP

下载本文档

489
0
约5.23千字
约 7页
2017-11-21 发布于浙江
举报
版权申诉

方柱绕流的lbm代码(c++).pdf

1、原创力文档（book118）网站文档一经付费（服务费），不意味着购买了该文档的版权，仅供个人/单位学习、研究之用，不得用于商业用途，未经授权，严禁复制、发行、汇编、翻译或者网络传播等，侵权必究。。
2、本站所有内容均由合作方或网友上传，本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺！文档内容仅供研究参考，付费前请自行鉴别。如您付费，意味着您自己接受本站规则且自行承担风险，本站不退款、不进行额外附加服务；查看《如何避免下载的几个坑》。如果您已付费下载过本站文档，您可以点击这里二次下载。
3、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等，请点击“版权申诉”（推荐），也可以打举报电话：400-050-0827(电话支持时间：9:00-18:30)。
4、该文档为VIP文档，如果想要下载，成为VIP会员后，下载免费。
5、成为VIP后，下载本文档将扣除1次下载权益。下载后，不支持退款、换文档。如有疑问请联系我们。
6、成为VIP后，您将拥有八大权益，权益包括：VIP文档下载权益、阅读免打扰、文档格式转换、高级专利检索、专属身份标志、高级客服、多端互通、版权登记。
7、VIP文档为合作方或网友上传，每下载1次，网站将根据用户上传文档的质量评分、类型等，对文档贡献者给予高额补贴、流量扶持。如果你也想贡献VIP文档。上传文档

方柱绕流的lbm代码(c)

//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],uv[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++) { evolution(); statistic(); boundary(); 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 evolution() { for(i=1;iNX;i++)//边界不参与碰撞（边界条件除外） for(j=1;jNY;j++) {