#include?
#include?
#include?
#include?
#include?
#include?
#include?

using?namespace?std;
const?int?Q?=?9;?????//D2Q9模型
const?int?NX?=?800;??//x方向
const?int?NY?=?200;??//y方向
const?double?U?=?0.1;?//入口速度?
double?cr=25;

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][NX+1][2]f[NX+1][NY+1][Q]
???????F[NX+1][NY+1][Q]gx[NX+1][NY-1]gy[NX+1][NY-1]f0[NX+1][NY+1][Q];

int?ijkipjpn;
double?Rerho0taunuerrorcxcy;

void?init（）;
double?feq（int?kdouble?rhodouble?u[2]）;
double?force（int?kdouble?u[2]double?gxdouble?gy）;
void?evolution（）;
void?outputdata（）;
void?Error（）;

int?main（）
{
??using?namespace?std;
??init（）;
??for（n?=?0;n<=20000?;n++）
??{
????evolution（）;
??if（n%100?==?0）
??{
????Error（）;
????cout<<“it?times“<????if?（error<0）
???break;
??
??if（n?>=?100）
??{
????if（n%1000?==?0）
??????outputdata（）;
??}
??}
??}
?return?0;
}

void?init（）
{
?rho0?=?1.0;
?Re?=?100;
?nu?=?U*2.0*cr/Re;
?tau?=?3.0*nu?+?0.5;
??cx=0.2*NX;
??cy=0.5*NY;
??
?std::cout<<“tau=“<
?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[0][j][0]?=?U;
???for（k?=?0;k????{
????f[i][j][k]=feq（krho[i][j]u[i][j]）;
???}
??}
?u[0][0][0]?=?0;?
?u[0][NY][0]?=?0;??
}

double?feq（int?kdouble?rhodouble?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;
}
double?force（int?kdouble?u[2]double?gxdouble?gy）
{
?double?fxfyeu;
?eu?=?（e[k][0]*u[0]?+?e[k][1]*u[1]）;
?fx=w[k]*（1.0-0.5/tau）*（3.0*（e[k][0]-u[0]）+9.0*eu*e[k][0]）*gx;
?fy=w[k]*（1.0-0.5/tau）*（3.0*（e[k][1]-u[1]）+9.0*eu*e[k][1]）*gy;
?return?（fx+fy）;
}


void?evolution（）
{
???for（i=1;i???for（j=1;j???for（k=0;k?{
??

if?（i>50?&&?i<100?&&?j>40?&&?j<60）
????{
????gx[i][j]=0;
???gy[i][j]=0;
????}
???else{
????gx[i][j]=0;
???gy[i][j]=0;}?
????
?ip?=?i?-?e[k][0];
?jp?=?j?-?e[k][1];
?
?f0[i][j][k]=?f[ip][jp][k];//只遷移，沒有碰撞的f函數?
?
????F[i][j][k]?=?f0[i][j][k]?+?（feq（krho[ip][jp]u[ip][jp]）-f0[i][j][k]）/tau
??+force（ku[i][j]gx[i][j]gy[i][j]）;
??}

//?計算宏觀變量?
???for（i?=?1;i????for（j?=?1;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???????{
????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]=?u[i][j][0]/rho[i][j]+0.5*gx[i][j]/rho[i][j];
????u[i][j][1]=?u[i][j][1]/rho[i][j]+0.5*gy[i][j]/rho[i][j];
???}
?//邊界處理

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????5655??2010-11-05?09:23??LBM-cylinder.cpp

-----------?---------??----------?-----??----

?????????????????5655????????????????????1