Lid Driven Flow.zip_圆柱绕流_平板流动_方腔流_非恒定流_非稳定流

#include<iostream> #include<cmath> #include<cstdlib> #include<iomanip> #include<fstream> #include<sstream> #include<string> using namespace std; const int Q = 9; const int NX = 256; const int NY = 256; 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]; int i, j, k, ip, jp, n; double c, Re, dx, dy, Lx, Ly, dt, rho0, p0, tau_f, niu, error; void init(); 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++) { evolution(); if (n % 100 == 0) { 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; cout << "The max relative error of uv is:" << setiosflags(ios::scientific) << error << endl; if (n >= 1000) { if (n % 1000 == 0) output(n); if (error < 1.0e-6) break; } } } return 0; } void init() { dx = 1.0; //x方向的网格步长 dy = 1.0; //y方向的网格步长 Lx = dx * double(NX); //x方向的长度 Ly = dy * double(NY); //y方向的长度 dt = dx; //时间步长 c = dx / dt; //格子速度 rho0 = 1.0; //初始密度 Re = 400; //雷诺数 niu = U * Lx / Re; //计算运动粘性系数 tau_f = 3.0*niu + 0.5; //无量纲松弛时间 std::cout << "tau_f=" << tau_f << endl; for (i = 0; i <= NX; i++) for (j = 0; j <= NY; j++) { u[i][j][0] = 0; //初始X方向速度为0 u[i][j][1] = 0; //初始Y方向速度为0 rho[i][j] = rho0; //初始密度各点均为0 u[i][NY][0] = U; //顶盖速度为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++) //左右边界 for (k = 0; k < Q; k++) { rho[NX][j] = rho[NX - 1][j]; f[NX][j][k] = feq(k, rho[NX][j], u[NX][j]) + f[NX - 1][j][k] - feq(k, rho[NX - 1][j], u[NX - 1][j]); 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++) //上下边界 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]); rho[i][NY] = rho[i][NY - 1]; u[i][NY][0] = U; 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]); } } void output(int m) //输出 { ostringstream name; name << "cavity_" << m << ".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 + 1 << ",J=" << NY + 1 << ",F=POINT" << endl; for (j = 0; j <= NY; j++) for (i = 0; i <= NX; i++) { out << double(i) / Lx << " " << double(j) / Ly << " " << 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); }