zhouhe定压边界模拟泊肃叶流动

城市阳光

          进出口 采用zhouhe定压边界模拟，上下边界是反弹处理，现在的问题是流体不流动，如图所示希望高手能指点一下，这个问题困惑很久了，如果边界修改成非平衡外推，则可以流动，zhouhe边界的代码也是对的（和matlab代码对照过），现在感觉的问题是在求取进入口边界上未知的f时，用的从固壁上反弹回来的f, 感觉这个反弹回来的f有问题（反弹应该在边界处理之后还是之前？？）程序放在附件里面了，C语言代码编写的，可在VC里面直接运行，希望各位能指点一下，万分感激，下面是代码的主体部分

void init()
{

dx = 1.0;
dy = 1.0;
Lx = dx*double(NX);
Ly = dy*double(NY);
dt = dx;
c = dx/dt; //1.0
rho0 = 1.0;////初始密度
Re = 10;
Kn=0.194;
niu=0.4;///////粘度系数
tau_f = (3*niu+0.5);//////松弛时间为1.7

std::cout << "tau_f= " << tau_f << endl;
printf("omega = %8.6f\n",1/tau_f);
printf("松弛时间 = %8.6f\n",tau_f);
printf("粘度系数 = %8.6f\n",niu);

for(i=0; i<=NX ;i++) //初始化
    for(j=0; j<=NY ; j++)
        {

         u[j][0] = 0;
         u[j][1] = 0;
         rho[j] = rho0;
          
       
         for(k=0; k<Q; k++)
                 {

            f[j][k]=feq(k,rho[j],u[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);
//feq = w[k]*(rho + 3.0*eu +4.5*eu*eu - 1.5*uv);
return feq;

}
void evolution() //////////
{


           // 计算宏观变量
   
      
                // 计算宏观变量
   
       for(i=0; i<=NX; i++)
          for(j=0; j<=NY; j++)
                  {

              u0[j][0] = u[j][0];
              u0[j][1] = u[j][1];
                       rho[j]=0;
               u[j][0]=0;
                        u[j][1]=0;
             for(k=0; k<Q; k++)
                         {     
                 //  f[j][k]=F[j][k];
                   rho[j] += f[j][k];
                   u[j][0]+= e[k][0]*f[j][k];
                   u[j][1]+= e[k][1]*f[j][k];

                         }
                 u[j][0] /=rho[j];
                 u[j][1] /=rho[j];

                  }

  

       
        for (i = 0; i<=NX; i++)//碰撞方程，包括进出口和上下固壁边界
            for (j = 0; j<=NY; j++)       
                 for (k = 0; k<Q; k++)
                        {
             
                  F[j][k] = f[j][k] + (feq(k, rho[j], u[j]) - f[j][k]) / tau_f;

                        }





for (int i=0; i<=NX; i++)///////////////迁移过程，整体迁移，包含四个边界
                {
                        for (int j=0; j<=NY; j++)
                        {                       
                                for (int k=0; k<Q; k++)
                                {
                               ip=int((i+e[k][0]+NX+1)%(NX+1));                          // STREAMING
                   jp=int((j+e[k][1]+NY+1)%(NY+1));
                   f[ip][jp][k]=F[j][k];

                                       
                       
                                }
                        }
                }

//// 左右定压边界

                    for(j=1;j<NY;j++)
                         {  
                       
                                rho[0][j]=1.01;
                                u[0][j][0]=1-(f[0][j][0]+f[0][j][2]+f[0][j][4]+2*(f[0][j][3]+f[0][j][7]+f[0][j][6]))/rho[0][j];
                                f[0][j][1]=f[0][j][3]+2/3*rho[0][j]*u[0][j][0];
                                f[0][j][5]=f[0][j][7]+1/6*rho[0][j]*u[0][j][0]+1/2*(f[0][j][4]-f[0][j][2]);
                                f[0][j][8]=f[0][j][6]+1/6*rho[0][j]*u[0][j][0]+1/2*(f[0][j][2]-f[0][j][4]);

                                rho[NX][j]=1.0;
                                u[NX][j][0]=-1+(f[NX][j][0]+f[NX][j][2]+f[NX][j][4]+2*(f[NX][j][1]+f[NX][j][5]+f[NX][j][8]))/rho[0][j];
                                f[NX][j][3]=f[NX][j][1]-2/3*rho[NX][j]*u[NX][j][0];
                                f[NX][j][7]=f[NX][j][5]-1/6*rho[NX][j]*u[NX][j][0]+1/2*(f[NX][j][2]-f[NX][j][4]);
                                f[NX][j][6]=f[NX][j][8]-1/6*rho[NX][j]*u[NX][j][0]+1/2*(f[NX][j][4]-f[NX][j][2]);
                        }
         
       






                        for (i = 0; i <= NX; i++) //上下反弹边界
        {

                    f[0][2]= f[0][4];
                     f[0][5]= f[0][7];
                          f[0][6]= f[0][8];


                         f[NY][4]= f[NY][2];
                           f[NY][7]=f[NY][5];
                        f[NY][8]= f[NY][6];


         
        }
       

                
   


}

void output(int m)
{

FILE*fp=NULL;
char filename[30];
sprintf(filename,"cavity_%d.dat",m);
fp=fopen(filename,"w");
fprintf(fp,"Title = \"LBM LID DRIVEN FLOW\"\n");
fprintf(fp,"Variables = \"X\", \"Y\", \"P\", \"U\", \"V\"\n");
fprintf(fp,"Zone T = \"BOX\",I = %d, J = %d, F = POINT\n", NX+1, NY+1);

for(j=0; j<=NY; j++)
   for(i=0; i<=NX; i++)
   {

      fprintf(fp,"%d %d %12.10f %12.10f %12.10f\n", i, j, rho[j], u[j][0], u[j][1]);

   }


}

速度图

poiseuilleb不流动.rar (2.28 KB, 下载次数: 19)

2016-7-11 22:31 上传

点击文件名下载附件

Ge_Zhang

看起来是只算了一步，没继续演化。

muzhiyu

兄弟，你的问题解决了吗

城市阳光

muzhiyu 发表于 2017-1-11 09:28
兄弟，你的问题解决了吗

没有，改成非平衡外推了