back-step_基于LBM方法的后台阶驱动流模拟_格子玻尔兹曼_back_后台阶_

#include <iostream> #include <cmath> #include <cstdlib> //提供一些函数与符号常量 #include <iomanip> #include <fstream> //来自文件的输入输出 #include <sstream> //来自string对象的输入输出 #include <string> using namespace std; const int Q = 9; //D2Q9模型，离散速度的总个数 const int NX = 600; //X方向的格子数 const int NY = 60; //Y方向的格子数 const double U = 0.1; const double ER = 1.25; //通道高度与台阶高度的比值 int e[Q][2] = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}}; //离散速度eα 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}; //权系数wα 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]; //rho是密度，u0是n时层的速度，u是n+1时层的速度，f和F分别是演化前后的分布函数 int i,j,k,ip,jp,n; //k为演化次数 double c,Re,dx,dy,Lx,Ly,dt,rho0,p0,tau_f,niu,error; //分别为格子速度，雷诺数，x、y方向的网格步长，x、y方向的长度，时间步长，流场初始速度，p0，无量纲松弛时间，运动黏度系数，两个相邻时层速度的最大相对误差 void init(); //initialization初始化 double feq(int k,double rho,double u[2]); void evolution(); void output(int m); void Error(); int main() { using namespace std; init(); for(n = 0; ;n++) //n没有大小限制 { evolution(); if(n%100 == 0) //如果n能被100整除 { Error(); cout << "The" << n << "th computation result: " << endl << "The u,v of point(NX/2,NY/2) is: " << setprecision(6) << u[NX/2][NY/2][0] << "," << u[NX/2][NY/2][1] << endl; //setprecision(n)输出流显示浮点数数字个数 cout << "The max relative error of uv is: " << setiosflags(ios::scientific) << error << endl; //setprecision(n)与setiosflags(ios::scientific)合用,控制指数表示法的小数位数。setiosflags(ios::scientific)是用指数方式表示实数 //setprecision(n)与setiosflags(ios::fixed)合用,控制小数点右边的数字个数。setiosflags(ios::fixed)是用定点方式表示实数 if(n >= 1000) { if(n%1000 == 0) output(n); if(error<1.0e-6) break; } } } return 0; } void init() { using namespace std; dx = 1.0; dy = 1.0; Lx = dx*double(NY); //x、y方向的长度为其网格步长与格子数乘积 Ly = dy*double(NX); dt = dx; //时间步长与网格步长相等 c = dx/dt; //格子速度为1.0 rho0 = 1.0; Re = 400; niu = U*Lx/Re; //黏度系数niu=顶盖速度U*方向长度/雷诺数，即Re=vd/niu；水的切应力等于粘度系数与流速梯度乘积 tau_f = 3.0*niu+0.5; // 无量纲松弛时间=3倍的黏度系数+0.5 cout << "tau_f = " << tau_f << endl; //输出黏度系数 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; if(j > NY/ER) { u[0][j][0] = U; } for(k = 0;k < Q;k++) { f[i][j][k] = feq(k,rho[i][j],u[i][j]); } } } 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; } void evolution() { for(i = 1;i < NX;i++) for(j = 1;j < NY;j++) for(k = 0;k < Q;k++) { ip = i - e[k][0]; jp = j - e[k][1]; F[i][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau_f; } for(i = 1;i < NX;i++) //计算宏观量 for(j = 1;j < NY;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; for(k = 0;k < Q;k++) { f[i][j][k] = F[i][j][k]; rho[i][j]+= f[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]; } u[i][j][0]/= rho[i][j]; u[i][j][1]/= rho[i][j]; } //边界处理 for (j = 1; j < NY; j++) { u[NX][j][0] = u[NX - 1][j][0]; u[NX][j][1] = 0; rho[NX][j] = rho[NX - 1][j]; f[NX][j][3] = 2 * f[NX - 1][j][3] - f[NX - 2][j][3]; f[NX][j][6] = 2 * f[NX - 1][j][6] - f[NX - 2][j][6]; f[NX][j][7] = 2 * f[NX - 1][j][7] - f[NX - 2][j][7]; } for(j = NY/ER;j < NY;j++) //左边界 { u[0][j][0] = U; u[0][j][1] = 0; for(k = 0;k < Q;k++) { rho[0][j] = rho[1][j]; f[0][j][k] = feq(k,rho[0][j],u[0][j]) + f[1][j][k] - feq(k,rho[1][j],u[1][j]); } } for(i = 0;i <= NX;i++) //上边界 { u[i][NY][0] = 0; u[i][NY][1] = 0; for(k = 0;k < Q;k++) { rho[i][NY] = rho[i][NY-1]; f[i][NY][k] = feq(k,rho[i][NY],u[i][NY]) + f[i][NY-1][k] - feq(k,rho[i][NY-1],u[i][NY-1]); } } for(i = 0;i < NX;i++) //下边界 { if(i <= NX/60)//这里修改台阶的长度 { for(j = 0;j < NY/ER;j++) { u[i][j][0] = 0; u[i][j][1] = 0; for(k = 0;k < Q;k++) { f[i][j][k] = feq(k,rho[i][j],u[i][j])+f[i+1][j][k] - feq(k,rho[i+1][j],u[i+1][j]); } } } else { u[i][0][0] = 0; u[i][0][1] = 0; for(k = 0;k < Q;k++) { rho[i][0] = rho[i][1]; f[i][0][k] = feq(k,rho[i][0],u[i][0])+f[i][1][k] - feq(k,rho[i][1],u[i][1]); } } } } void output(int m) //输出 { ostringstream name; //ostringstream类通常用于执行C风格的串流的输出操作，格式化字符串，避免申请大量的缓冲区，替代sprintf。 name << "Backstep_Re="<<Re<< ",ER=" << ER << "," << m << ".dat"; ofstream out(name.str().c_str()); //ofstream是从内存到硬盘，ifstream是从硬盘到内存 out << "Tittle = \" LBM Backstep Flow\" \n " << "VARIABLES = \"X\",\"Y\",\"U\",\"V\"\n" << "ZONE T = \"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) << " " << double(j) << " " << u[i][j][0] << " " << u[i][j][1] << endl; } } void Error() { double temp1,temp2; temp1 = 0; temp2 = 0; for(i = 1;i < NX;i++) for(j = 1;j < NY;j++) { temp1 += ((u[i][j][0] - u0[i][j][0]) * (u[i][j][0] - u0[i][j][0]) + (u[i][j][1]-u0[i][j][1]) * (u[i][j][1]-u0[i][j][1])); temp2 += (u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1]); } temp1 = sqrt(temp1); temp2 = sqrt(temp2); error = temp1/(temp2 + 1e-30); }