# include <stdio.h>
# include <iostream>
# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <fstream>
# include <sstream>
# include <string>
using namespace std; //不解!!!!!!!!!!!!!!!!!!!
const int Q=9; //D2Q9之意
const int NX=100; //X方向格子数
const int NY=100; //Y方向格子
const double T2=354,T1=273; //两端温度,这里温度定义为内能密度,即内能除以物质量
double Gg[2]={0,-9.8}; //重力
const double Lx=100; //总长度
const double Ly=100;
double e[Q][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}}; //直角方向 9个移动方向
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],vel[NX+1][NY+1]; //密度,速度,分布函数定义
double T[NX+1][NY+1],g[NX+1][NY+1][Q],G[NX+1][NY+1][Q]; //温度场和内能分布函数的定义
double wf,wg,ut[2]; //两个系数和加速度
int i,j,k,ip,jp,n;
double c,Ra,Prandtl,R=287,dx,dy,dt,rho0,P0,tau_f,tau_g,niu,error;
//#define Fa (ut[0]*(e[k][0]-u[ip][jp][0])+ut[1]*(e[k][1]-u[ip][jp][1]))*feq(k,rho[ip][jp],u[ip][jp])/R/T[ip][jp]; //第一种外力
#define Fa dt*Gg[1]*(e[k][1]-u[ip][jp][1])*feq(k,rho[ip][jp],u[ip][jp])/R/T[ip][jp]; //第二种外力
//#define Fa (T[i][j]-T1)*w[k]*(e[k][1]/3-u[ip][jp][1]/3+(e[k][0]*u[ip][jp][0]+e[k][1]*u[ip][jp][1])*e[k][1]/9)*Gg[1]
void init(); //函数声明
double feq(int k, double rho, double u[2]);
// double feq_t(int k, double T, double u[2]);
void evolution();
void output(int m);
void Error();
int main()
{
using namespace std; //变量名声明?
init(); //初始化
for(n=0; ;n++) //循环
{
evolution(); //1.演化计算
if(n%200==0) //2.每200步输出一次
{
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) //每1000步写出一次
output(n); //收敛判断
if(error<1.0e-6)
break;
}
}
}
return 0;
}
void init() //初始化函数
{
dx=Lx/double(NX);; //物理参数设置,每格长度为1.0
dy=Ly/double(NY);
dt=dx; //时间步长,默认1.0
c=dx/dt; //格子速度
rho0=1.0; //密度
Ra=10000; //Ra数
Prandtl=0.71; //Prandtl数
niu=0.00001; //运动粘性系数
tau_f=0.1*dy*pow(3*Prandtl,0.5)/pow(Ra,0.5)+0.5; //无量纲松弛时间
tau_g=tau_f/Prandtl;
wf=dt/(tau_f+0.5*dt);
wg=dt/(tau_g+0.5*dt);
std::cout<<"tau_f= "<<tau_f<<endl; //输出无量纲松弛时间
for(i=0;i<=NX;i++) //每个格子中速度和密度,分布函数f的初始化
{
for(j=0;j<=NY;j++)
{
u[i][j][0]=0;
u[i][j][1]=0;
rho[i][j]=rho0;
// T[i][j]=T2-i*(T2-T1)/NX;
T[i][j]=T1;
T[NX][j]=T1;
T[0][j]=T2;
for(k=0;k<Q;k++)
{
f[i][j][k]=feq(k,rho[i][j],u[i][j]);
g[i][j][k]=feq(k,T[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;
}
/*
double feq_t(int k,double T, 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]);
/* 完整He模型
if (k<1)
feq=-T*uv*2/3;
else if (0<k<5)
feq=T*(1.5+1.5*eu+4.5*eu*eu-1.5*uv)/9;
else if (4<k)
feq=T*(3+6*eu+4.5*eu*eu-1.5*uv)/36;
return feq;
feq=w[k]*T*(1.0+3.0*eu+4.5*eu*eu-1.5*uv);
}
*/
void evolution() //演化计算
{
for(i=1;i<NX;i++) //分布函数演化
{
for(j=1;j<NY;j++)
{
for(k=0;k<Q;k++)
{
ip=i-int(e[k][0]);
jp=j-int(e[k][1]);
ut[0]=u[ip][jp][0]-u[i][j][0];
ut[1]=u[ip][jp][1]-u[i][j][1];
// F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau_f+;
F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau_f+dt*Fa;
G[i][j][k]=g[ip][jp][k]+(feq(k,T[ip][jp],u[ip][jp])-g[ip][jp][k])/tau_g;
}
}
}
for(i=1;i<NX;i++) //计算宏观量
{
for(j=1;j<NX;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;
T[i][j]=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];
}
for(k=0;k<Q;k++)
{
g[i][j][k]=G[i][j][k];
T[i][j]+=g[i][j][k];
// ug[i][j][0]+=e[k][0]*g[i][j][k];
// ug[i][j][1]+=e[k][1]*g[i][j][k];
}
// u[i][j][0]=(u[i][j][0]+ug[i][j][0])/(rho[i][j]+T[i][j]);
// u[i][j][1]=(u[i][j][1]+ug[i][j][1])/(rho[i][j]+T[i][j]);
u[i][j][0]/=rho[i][j];
u[i][j][1]/=rho[i][j];
vel[i][j]=sqrt(u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]);
}
}
/* 边界处理 */
for(j=1;j<NY;j++) //左右边界
{
rho[NX][j]=rho[NX-1][j];rho[0][j]=rho[1][j]; //边界上的密度直接继承旁边网格
// rho[NX][j]=6*(f[NX][j][5]+f[NX][j][1]+f[NX][j][8]); //壁面质量守恒
// rho[0][j]=6*(f[0][j][6]+f[0][j][3]+f[0][j][7]);
for(k=0;k<Q;k++)
{
// u[NX][j][0]=0;
// u[NX][j][1]=0;
// u[NX][j][0]=u[NX-1][j][0];
// u[NX][j][1]=u[NX-1][j][1];
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]); //非平衡外推格式
// u[0][j][0]=0;
// u[0][j][1]=0;
// u[0][j][0]=u[1][j][0];
// u[0][j][1]=u[1][j][0];
f[0][j][k]=feq(k,rho[0][j],u[0][j])+f[1][j][k]-feq(k,rho[1][j],u[1][j]); //非平衡外推格式
T[NX][j]=T1; //温度边界
g[NX][j][k]=feq(k,T[NX][j],u[NX][j])+g[NX-1][j][k]-feq(k,T[NX-1][j],u[NX-1][j]);
T[0][j]=T2;
g[0][j][k]=feq(k,T[0][j],u[0][j])+g[1][j][k]-feq(k,T[1][j],u[1][j]);
}
}
for(i=1;i<=NY;i++) //上下边界
{
rho[i][0]=rho[i][1];rho[i][NY]=rho[i][NY-1]; //边界上的密度直接继承旁边网格
// rho[i][0]=6*(f[i][0][7]+f[i][0][4]+f[i][0][8]); //壁面质量守恒
// rho[i][NX]=6*(f[i][NX][6]+f[i][NX][2]+f[i][NX][5]);
for(k=0;k<Q;k++)
{
// u[i][0][0]=0;
// u[i][0][1]=0;
// u[i][0][0]=u[i][1][0];
// u[i][0][1]=u[i][1][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]);
// u[i][NY][0]=0;
// u[i][NY][1]=0;
// u[i][NY][0]=u[i][NY-1][0];
// u[i][NY][1]=u[i][NY-1][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]);
T[i][0]=T[i][1]; //温度边界
g[i][0][k]=feq(k,T[i][0],u[i][0])+g[i][1][k]-feq(k,T[i][1],u[i][1]);
T[i][NY]=T[i][NY-1];
g[i