格子玻尔兹曼方法计算顶盖驱动流源程序（C++版）

#include
#include
#include
#include
#include
#include
#include

using namespace std;        
const int Q=9;          //D2Q9模型
const int NX=256;		//x方向
const int NY=256;		//y方向
const double U=0.1;		//顶盖速度

int e[Q][2]={{00}{10}{01}{-10}{0-1}{11}{-11}{-1-1}{1-1}};
double w[Q]={4.0/91.0/91.0/91.0/91.0/91.0/361.0/361.0/361.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 ijkipjpn;
double cRedxdyLxLydtrho0P0tau_fniuerror;

void init（）;
double feq（int kdouble 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“<			cout<<“The uv of point（NX/2NY/2） is: “<				<			cout<<“The max relative error of uv is: “
				<			if（n>=1000）
			{
				if（n%1000==0） output（n）;
				if（error<1.0e-6） break;
			}
		}
	}
	return 0;
}

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=1000;
	niu=U*Lx/Re;
	tau_f=3.0*niu+0.5;
	std::cout<<“tau_f= “<
	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;
			u[i][NY][0]=U;
			for（k=0; k				f[i][j][k]=feq（krho[i][j]u[i][j]）;
		}
}

double feq（int kdouble rho double u[2]）
{
	double euuvfeq;
	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		for（j=1;j			for（k=0;k			{
				ip=i-e[k][0];
				jp=j-e[k][1];