LBM 长立柱.rar_LBM 方柱绕流_LBM方柱绕流_方 绕流_方柱绕流lbm_绕流

//1. 头文件及声明文件

#include<iostream>

#include<cmath>

#include<cstdlib>

#include<iomanip>

#include<fstream>

#include<sstream>

#include<string>

using namespace std;

int M=0;

const int m=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 Error();

//2. 主程序文件

int main()

evolution();

statistic();

if(n%100==0)

{

Error();

cout<<"The "<<n<<"th conputation 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;

output(n);

if(n>=1000)

{

//3. 初始化文件

void init() //初始函数

{

dx=1.0;

M=m*H;

Re=200;

niu=U*(m+2)*H/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;

u[i][j][1]=0;