#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];