LBM圆柱绕流

lbm方法模拟圆柱绕流，边界条件采用平衡外推格式，有明显的涡街

#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 < Q;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 < NX;i++）
   for（j = 1;j < NY;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 < 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];
   }
   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