i meet over a problem about the pressure boundary condition of Guozhaoli’s non-equilibrium extrapolationmethod for pressure boundary conditions in simulating the 2D Poiseuille. simulating solutions are very different from the analytical solutions. I do not know the cause. I write codes with c++. need helps.
#include
#include
#include
#include
#include
#include
#include
using namespace std;
const int Q = 9; //D2Q9模型
const int NX = 20; //x方向
const int NY = 10; //y方向
const double U = 0.1; //入口速度
double cr=25;
int e[Q][2] = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
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][NX+1][2];
//速度空间
double f[NX+1][NY+1][Q],feqq[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],temp[NX+1][NY+1][Q];
// 矩空间密度分布函数
double m[NX+1][NY+1][Q],meqq[NX+1][NY+1][Q],M[NX+1][NY+1][Q];
int opposite[Q]={0,3,4,1,2,7,6,5,6};
int image[NX][NY];
int i,j,k,ip,jp,n;
double Re,rho0,tau,nu,error,omega,cx,cy;
double s[Q]; //D2Q9松弛系数矩阵
//D2Q0转换矩阵
double MM[Q][Q] = {{1,1,1,1,1,1,1,1,1},{-1,-1,-1,-1,2,2,2,2,-4},{ -2,-2,-2,-2,1,1,1,1,4}
,{1,0,-1,0,1,-1,-1,1,0},{-2,0,2,0,1,-1,-1,1,0},{0,1,0,-1,1,1,-1,-1,0},{0,-2,0,2,1,1,-1,-1,0}
,{1,-1,1,-1,0,0,0,0,0},{ 0,0,0,0,1,-1,1,-1,0}};
double INM[Q][Q];//D2Q9转换矩阵的逆矩阵
//double jx,jy,jz;
//
double inrho = 1.1;
double outrho = 1.0;
double deltaP;
double umax;
double uc[NX+1][NY+1];
//
void init();
void imagec();
double feqc(int k,double rho,double u[2]);
double meqc(int k,double rho,double u[2]);
double force(int k,double u[2],double gx,double gy);
void evolution();
void outputdata();
void Error();
int main()
{
using namespace std;
init();
imagec();
for(n = 0;n<=60000 ;n++)
{
evolution();
if(n%100 == 0)
{
Error();
cout<<“it times”<<n<<" error="<<error<<endl;
if (error<0)
break;
if(n >= 100)
{
if(n%1000== 0)
outputdata();
}
}
}
return 0;
}
//
//初始子程序
void init()
{
rho0 = 1.0;
Re = 10;
// nu = 0.161; //U2.0cr/Re;
// tau = 3.0nu + 0.5;
omega= 0.9; //1.0/tau;
tau = 1.0/omega;
nu = 1.0/3.0(tau-0.5);
cx=0.2NX;
cy=0.5NY;
std::cout<<“tau=”<<tau<<endl;
std::cout<<“nu=”<<nu<<endl;
system(“pause”);
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] = 0;
for(k = 0;k < Q;k++)
{
f[i][j][k]=feqc(k,rho[i][j],u[i][j]);
feqq[i][j][k]=feqc(k,rho[i][j],u[i][j]);
m[i][j][k]=meqc(k,rho[i][j],u[i][j]);
}
}
u[0][0][0] = 0;
u[0][NY][0] = 0;
//
deltaP = (inrho-outrho)/NX;
std::cout<<“deltaP=”<<deltaP<<endl;
system(“pause”);
umax = deltaP/(2nu);
std::cout<<“umax=”<<umax<<endl;
system(“pause”);
for(i = 0; i <=NX ;i++)
for(j = 0;j <=NY ; j++)
{
uc[i][j] = umaxj*((NY-1)-j);
}
}
//
// 图像处理子程序
void imagec()
{
for(i = 1; i < NX; i++)
for(j = 1; j < NY; j++)
{
if (((i-cx)(i-cx)+(j-cy)(j-cy))<=cr*cr)
{
image[i][j]=1;
}
else
{
image[i][j]=0;
}
}
}
//
//计算平衡态函数-速度空间
double feqc(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.0eu +4.5eueu - 1.5uv);
return feq;
}
//
//计算平衡态函数-速度矩空间
double meqc(int k,double rho,double u[2])
{
double jx,jy,meq;
jx = u[1];
jy = u[2];
switch(k)
{
case 0: meq = rho;
break;
case 1: meq = -2.0rho + 3.0(jxjx + jyjy);
break;
case 2: meq = rho -3.0*(jxjx + jyjy); //energy square
break;
case 3: meq = jx; //momentum in x-dir
break;
case 4: meq = -jx; //energy flux in x-dir—/-jx
break;
case 5: meq = jy; //momentum in y-dir
break;
case 6: meq = -jy; //energy flux in y-dir—/-jy
break;
case 7: meq = (jxjx - jyjy); //diagonal comp of stress tensor
break;
case 8: meq = jxjy;
break;
}
return meq;
}
//
//计算外力
double force(int k,double u[2],double gx,double gy)
{
double fx,fy,eu;
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.0eue[k][0])gx;
fy=w[k](1.0-0.5/tau)(3.0(e[k][1]-u[1])+9.0eue[k][1])*gy;
return (fx+fy);
}
//
void evolution()
{
// 计算平衡态函数-速度矩空间
for(i = 0; i <=NX ;i++)
for(j = 0;j <=NY ; j++)
for(k = 0;k < Q;k++)
{
meqq[i][j][k]=meqc(k,rho[i][j],u[i][j]);
}
//
//
//边界处理
//左右边界处理
for(j = 1;j < NY;j++)
for(k = 0;k < Q;k++)
{
//右边界
rho[NX][j] =outrho; //rho[NX-1][j];
u[NX][j][0] =u[NX-1][j][0];
f[NX][j][k] = feqc(k,rho[NX][j],u[NX][j]) + f[NX -1][j][k] - feqc(k,rho[NX
- 1][j],u[NX - 1][j]);
//左边界
rho[0][j] = inrho; //rho[1][j];
u[0][j][0] = u[1][j][0]; //U;
f[0][j][k]= feqc(k,rho[0][j],u[0][j])+ f[1][j][k]-feqc(k,rho[1][j],u[1][j]);
}
//左右边界处理结束
//上下边界处理
for(i = 0;i <= NX;i++)
for(k = 0;k < Q;k++)
{
//下边界
rho[i][0] = rho[i][1];
f[i][0][k] = feqc(k,rho[i][0],u[i][0])+f[i][1][k]-feqc(k,rho[i][1],u[i][1]);
//上边界
rho[i][NY] = rho[i][NY - 1];
f[i][NY][k] = feqc(k,rho[i][NY],u[i][NY]) + f[i][NY - 1][k]
- feqc(k,rho[i][NY - 1],u[i][NY -1]);
}
//上下边界处理结束
//迁移开始
for(i=1;i<NX;i++)
for(j=1;j<NY;j++)
for(k=0;k<Q;k++)
{
ip = i - e[k][0];
jp = j - e[k][1];
// f0[i][j][k]= f[ip][jp][k];//只迁移，没有碰撞的f函数
// if(image[i][j]==0)
// {//if
// F[i][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau;
// F[i][j][k] = f[ip][jp][k] + (feqq[ip][jp][k]-f[ip][jp][k])/tau;
//+force(k,u[i][j],gx[i][j],gy[i][j]);
// M[i][j][k] = m[ip][jp][k] + s(k)(meqq[ip][jp][k]-m[ip][jp][k]);
// F[i][j][k] = inverMM[i][j][k];
// }//if
//else
// {
f0[i][j][k]= f[ip][jp][k];
// }
}
// 迁移结束
//
// 计算宏观变量
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] = f0[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.5gx[i][j]/rho[i][j];
u[i][j][1]= u[i][j][1]/rho[i][j];//+0.5gy[i][j]/rho[i][j];
}
// 计算平衡态函数-速度空间
for(i = 0; i <=NX ;i++)
for(j = 0; j <=NY ; j++)
for(k = 0;k < Q;k++)
{
feqq[i][j][k]=feqc(k,rho[i][j],u[i][j]);
}
//碰撞开始
for(i=1;i<NX;i++)
for(j=1;j<NY;j++)
for(k=0;k<Q;k++)
{
gx[i][j] = deltaP; //(inrho-outrho)/NX;
gy[i][j] = 0.0;
// if(image[i][j]==0)
// {//if
// F[i][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp])-f[ip][jp][k])/tau;
f[i][j][k] = f0[i][j][k] + omega*(feqq[i][j][k]-f0[i][j][k]) +3.0w[k](e[k][0]deltaP+e[k][1]0.0);
// F[i][j][k] = inverMM[i][j][k];
// f[i][j][k] = f0[i][j][k] + omega(feqq[i][j][k]-f0[i][j][k]);// +force(k,u[i][j],gx[i][j],gy[i][j]);
// M[i][j][k] = m[ip][jp][k] + s(k)*(meqq[ip][jp][k]-m[ip][jp][k]);
// }//if
}
//碰撞结束
//圆柱for(i = 1; i < NX; i++)
//for(j = 1; j < NY; j++)
//for(k = 0; k< 9; k++)
//{
//if (((i-cx)(i-cx)+(j-cy)(j-cy))<=cr*cr)
// if(image[i][j]==1)
// {
// u[i][j][0]=0;
// u[i][j][1]=0;
// f[i][j][k] = f0[i][j][k];
// temp[i][j][k]= f[i][j][k];
// f[i][j][k] = f[i][j][opposite[k]];
// f[i][j][opposite[k]]=temp[i][j][k];
//}
//}
//
}
//
void Error()
{
double temp1,temp2;
temp1 = 0;
temp2 = 0;
for(i = 1;i < NX;i++)
for(j = 1;j < NY;j++)
{
temp1 += ((u[i][j][0] - u0[i][j][0])(u[i][j][0] - u0[i][j][0]) +
(u[i][j][1] - u0[i][j][1])(u[i][j][1] - u0[i][j][1]));
temp2 += (u[i][j][0]*u[i][j][0] + u[i][j][1]*u[i][j][1]);
}
temp1 = sqrt(temp1);
temp2 = sqrt(temp2);
error = temp1/(temp2 + 1e-30);
}
void outputdata()
{
if(n%100==0)
{
ofstream file;
char temp[10];
string str;
itoa(n,temp,10);
str = string(temp);
str += “u.dat”;
file.open(str.c_str());
file<<“Title=“LBM Lid Driven Flow”\n”<<“VARIABLES=“X”,“Y”,“U”,“V”,“rho”\n”<<“ZONE T=“BOX”,I=”
<<NX+1<<",J="<<NY+1<<",F=POINT"<<endl;
if (file.is_open())
{
for ( int i = 0; i <= NX; i++ )
for ( int j = 0; j <= NY; j++ )
{
// file<<double(i)<<" “<<double(j)<<” “<<u[i][j][0]<<” “<<u[i][j][1]<<” “<<rho[i][j]<<endl;
// file<<setiosflags(ios::left);
file<<setiosflags(ios::right)<<double(i)<<” “<<double(j)<<” “<<setiosflags(ios::scientific)<<setprecision(6)<<u[i][j][0]<<” “<<uc[i][j]<<” "<<rho[i][j]<<resetiosflags(ios::scientific)<<endl;
}
}
}
}