自己编写的两项shanchen模型模拟通道内的流动程序有点问题，想求指教

qq281362022

初场为二维通道中有一枚液滴
#include<iostream>
#include<cmath>
#include<cstdlib>
#include<iomanip>
#include<fstream>
#include<sstream>
#include<string>
using namespace std;
const int Q=9;
const int xDim=250;
const int yDim=50;
const int tMax=1;
const double Re=10.0;
const double g11=-0.15,g22=-0.15;
const double g12=0.28,g21=0.28;
const double gw1=0.0,gw2=0.8;
const double tau1=1.0,tau2=1.0;
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};
int opposite[Q]={0,3,4,1,2,7,8,5,6};
double f1[xDim][yDim][Q],f2[xDim][yDim][Q],g1[xDim][yDim][Q],g2[xDim][yDim][Q],feq1[xDim][yDim][Q],feq2[xDim][yDim][Q],u_1[xDim][yDim],u_2[xDim][yDim],v_1[xDim][yDim],
            v_2[xDim][yDim],uSqr1[xDim][yDim],uSqr2[xDim][yDim],rho1[xDim][yDim],rho2[xDim][yDim],uProf[yDim],phi1[xDim][yDim],phi2[xDim][yDim],
                fx11[xDim][yDim],fy11[xDim][yDim],fx22[xDim][yDim],fy22[xDim][yDim],fx12[xDim][yDim],fy12[xDim][yDim],fx21[xDim][yDim],fy21[xDim][yDim],
                fwx1[xDim][yDim],fwy1[xDim][yDim],fwx2[xDim][yDim],fwy2[xDim][yDim],fintX1[xDim][yDim],fintX2[xDim][yDim],fintY1[xDim][yDim],fintY2[xDim][yDim],u0[xDim][yDim],
                v0[xDim][yDim],ueq1[xDim][yDim],veq1[xDim][yDim],ueq2[xDim][yDim],veq2[xDim][yDim],u_m[xDim][yDim],v_m[xDim][yDim],rhom[xDim][yDim],pre[xDim][yDim],
                fTmi1[xDim][yDim][Q],fTmi2[xDim][yDim][Q],fTmii1[xDim][yDim][Q],fTmii2[xDim][yDim][Q];
int x,y,k,xd,yd,n;
double c,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,omega,error,xt,yt,zx,zy;
void init();
double f_eq(int k,double rho,double u,double v);
void collide();
void stream();
void computeMacros();
void output(int m);
int main()
{
        std::cout<<"tau_f= "<<tau_f<<endl;
        for(x=1;x<=xDim;x++)
        {
                for(y=1;y<=yDim;y++)
                {
                        if(y==1||y==yDim)
                        {
                                rho1[x][y]=0;
                            rho2[x][y]=0;
                         }
                        else
                        {
                        if (x>20&&x<30&&y>20&&y<30)
                        {
                                rho1[x][y]=0.01;
                            rho2[x][y]=1.0;
                         }
                        else
                        {
                                rho1[x][y]=1.0;
                            rho2[x][y]=0.01;
                        }
                        }
                }
        }
   for(x=1;x<=xDim;x++)
        {
                for(y=1;y<=yDim;y++)
                {
                 u_1[x][y]=0;
             v_1[x][y]=0;
                 u_2[x][y]=0;
             v_2[x][y]=0;
                        for(k=0;k<Q;k++)
                        {
            f1[x][y][k]=f_eq(k,rho1[x][y],u_1[x][y],v_1[x][y]);
            f2[x][y][k]=f_eq(k,rho2[x][y],u_2[x][y],v_2[x][y]);
                        }
                }
        }
   

using namespace std;
        int iter;
        double p[xDim][yDim][Q],h[xDim][yDim][Q],m[xDim][yDim][Q];
        for(iter=0;iter<tMax;iter++)
        {                 

         //流动方程
//1
        for(x=1;x<=xDim-1;x++)
        {
            for(y=2;y<=yDim-1;y++)
            {
                    p[xDim][y][1]=f1[xDim][y][1];
                        f1[x+1][y][1]=f1[x][y][1];
                        f1[1][y][1]=p[xDim][y][1];
                        p[xDim][y][1]=f2[xDim][y][1];
                        f2[x+1][y][1]=f2[x][y][1];
                        f2[1][y][1]=p[xDim][y][1];
                }
          }
//2
        for(x=1;x<=xDim;x++)
        {            
            for(y=2;y<=yDim-2;y++)
                {
                        p[x][2][2]=f1[x][2][2];
                        f1[x][y][2]=f1[x][y+1][2];
                        f1[x][yDim-1][2]=p[x][2][2];
                        p[x][2][2]=f2[x][2][2];
                        f2[x][y][2]=f2[x][y+1][2];
                        f2[x][yDim-1][2]=p[x][2][2];
                }
        }
//3
        for(x=1;x<=xDim-1;x++)
        {
            for(y=2;y<=yDim-1;y++)
            {
                        p[1][y][3]=f1[1][y][3];
                        f1[x][y][3]=f1[x+1][y][3];
                        f1[xDim][y][3]=p[1][y][3];
                        p[1][y][3]=f2[1][y][3];
                        f2[x][y][3]=f2[x+1][y][3];
                        f2[xDim][y][3]=p[1][y][3];
                }
        }
//4
        for(x=1;x<=xDim;x++)
        {            
            for(y=2;y<=yDim-2;y++)
                {
                        p[x][yDim-1][4]=f1[x][yDim-1][4];
                        f1[x][y+1][4]=f1[x][y][4];
                        f1[x][2][4]=p[x][yDim-1][4];
                        p[x][yDim-1][4]=f2[x][yDim-1][4];
                        f2[x][y+1][4]=f2[x][y][4];
                        f2[x][2][4]=p[x][yDim-1][4];
                }
        }
//5       
        for(x=1;x<=xDim-1;x++)
        {
                 for(y=2;y<=yDim-2;y++)
                {
                        p[xDim][y+1][5]=f1[xDim][y+1][5];
                        h[x][yDim-1][5]=f1[x][yDim-1][5];
                        m[xDim][yDim-1][5]=f1[xDim][yDim-1][5];
                        f1[x+1][y+1][5]=f1[x][y][5];
                        f1[1][y][5]=p[xDim][y+1][5];
                        f1[x][yDim-1][5]=h[x][yDim-1][5];
                        f1[1][2][5]=m[xDim][yDim-1][5];
                        p[xDim][y+1][5]=f2[xDim][y+1][5];
                        h[x][yDim-1][5]=f2[x][yDim-1][5];
                        m[xDim][yDim-1][5]=f2[xDim][yDim-1][5];
                        f2[x+1][y+1][5]=f2[x][y][5];
                        f2[1][y][5]=p[xDim][y+1][5];
                        f2[x][yDim-1][5]=h[x][yDim-1][5];
                        f2[1][2][5]=m[xDim][yDim-1][5];
                 }
        }
//6
        for(x=1;x<=xDim-1;x++)
        {
                 for(y=2;y<=yDim-2;y++)
                {
                        p[1][y][6]=f1[1][y][6];
                        h[x+1][yDim-1][5]=f1[x+1][yDim-1][6];
                        m[1][yDim-1][6]=f1[1][yDim-1][6];
                        f1[x][y+1][6]=f1[x+1][y][6];
                        f1[xDim][y+1][6]=p[1][y][6];
                        f1[x][2][6]=h[x+1][yDim-1][5];
                        f1[xDim][2][6]=m[1][yDim-1][6];
                        p[1][y][6]=f2[1][y][6];
                        h[x+1][yDim-1][5]=f2[x+1][yDim-1][6];
                        m[1][yDim-1][6]=f2[1][yDim-1][6];
                        f2[x][y+1][6]=f2[x+1][y][6];
                        f2[xDim][y+1][6]=p[1][y][6];
                        f2[x][2][6]=h[x+1][yDim-1][5];
                        f2[xDim][2][6]=m[1][yDim-1][6];
                 }
        }
//7
        for(x=1;x<=xDim-1;x++)
        {
                 for(y=2;y<=yDim-2;y++)
                {
                        p[1][y+1][7]=f1[1][y+1][7];
                        h[x+1][2][7]=f1[x+1][2][7];
                        m[1][2][7]=f1[1][2][7];
                        f1[x][y][7]=f1[x+1][y+1][7];
                        f1[xDim][y][7]=p[1][y+1][7];
                        f1[x][yDim-1][7]=h[x+1][2][7];
                        f1[xDim][yDim-1][7]=m[1][2][7];
                        p[1][y+1][7]=f2[1][y+1][7];
                        h[x+1][2][7]=f2[x+1][2][7];
                        m[1][2][7]=f2[1][2][7];
                        f2[x][y][7]=f2[x+1][y+1][7];
                        f2[xDim][y][7]=p[1][y+1][7];
                        f2[x][yDim-1][7]=h[x+1][2][7];
                        f2[xDim][yDim-1][7]=m[1][2][7];
                 }
        }
//8
        for(x=1;x<=xDim-1;x++)
        {
                 for(y=2;y<=yDim-2;y++)
                {
                        p[xDim][y+1][8]=f1[xDim][y+1][8];
                        h[x][2][8]=f1[x][2][8];
                        m[xDim][2][8]=f1[xDim][2][8];
                        f1[x+1][y][8]=f1[x][y+1][8];
                        f1[1][y][8]=p[xDim][y+1][8];
                        f1[x+1][yDim-1][8]=h[x][2][8];
                        f1[1][yDim-1][8]=m[xDim][2][8];
                        p[xDim][y+1][8]=f2[xDim][y+1][8];
                        h[x][2][8]=f2[x][2][8];
                        m[xDim][2][8]=f2[xDim][2][8];
                        f2[x+1][y][8]=f2[x][y+1][8];
                        f2[1][y][8]=p[xDim][y+1][8];
                        f2[x+1][yDim-1][8]=h[x][2][8];
                        f2[1][yDim-1][8]=m[xDim][2][8];
                 }
        }
//边界条件        zouhe入口       
         for(y=2;y<=yDim-1;y++)
         {
                uProf[y]=0.001;
                u_1[1][y]=uProf[y];
                v_1[1][y]=0;
                rho1[1][y]=f1[1][y][0]+f1[1][y][2]+f1[1][y][4]+2*(f1[1][y][3]+f1[1][y][6]+f1[1][y][7])/(1-u_1[1][y]);
                fTmii1[1][y][1]=f1[1][y][3]+(2*rho1[1][y]*u_1[1][y])/3;
                fTmii1[1][y][5]=(rho1[1][y]*u_1[1][y])/6+f1[1][y][7]+(f1[1][y][4]-f1[1][y][2])/2+(rho1[1][y]*v_1[1][y])/2;
        fTmii1[1][y][8]=(rho1[1][y]*u_1[1][y])/6+f1[1][y][6]+(f1[1][y][2]-f1[1][y][4])/2-(rho1[1][y]*v_1[1][y])/2;
                f1[1][y][1]=fTmii1[1][y][1];
                f1[1][y][5]=fTmii1[1][y][5];
                f1[1][y][8]=fTmii1[1][y][8];
     }
//自由流出口
         for(y=2;y<=yDim-1;y++)
         {
                 for(k=0;k<9;k++)
                 {
                f1[xDim][y][k]=f1[xDim-1][y][k];
                f2[xDim][y][k]=f2[xDim-1][y][k];
         }       
         }
//标准反弹         
     for(x=2;x<=xDim-1;x++)
         {               
                 for(k=0;k<Q;k++)
                 {
                        fTmi1[x][2][2]=f1[x][2][4];
                        fTmi1[x][2][5]=f1[x][2][7];
                        fTmi1[x][2][6]=f1[x][2][8];
                    fTmi1[x][2][2]=f1[x][2][4];
                        fTmi1[x][2][5]=f1[x][2][7];
                        fTmi1[x][2][6]=f1[x][2][8];
              }
                 for(k=0;k<Q;k++)
                 {
                         f1[x][2][2]=fTmi1[x][2][2];
                         f1[x][2][5]=fTmi1[x][2][5];
                         f1[x][2][6]=fTmi1[x][2][6];
                         f2[x][2][2]=fTmi2[x][2][2];
                         f2[x][2][5]=fTmi2[x][2][5];
                         f2[x][2][6]=fTmi2[x][2][6];
                 }
      }

     for(x=2;x<=xDim-1;x++)
         {               
                 for(k=0;k<Q;k++)
                 {
                        fTmi1[x][yDim][4]=f1[x][yDim-1][2];
                        fTmi1[x][yDim][7]=f1[x][yDim-1][5];
                        fTmi1[x][yDim][8]=f1[x][yDim-1][6];
                    fTmi1[x][yDim][4]=f1[x][yDim-1][2];
                        fTmi1[x][yDim][7]=f1[x][yDim-1][5];
                        fTmi1[x][yDim][8]=f1[x][yDim-1][6];
              }
                 for(k=0;k<Q;k++)
                 {
                         f1[x][yDim][4]=fTmi1[x][yDim][4];
                         f1[x][yDim][7]=fTmi1[x][yDim][7];
                         f1[x][yDim][8]=fTmi1[x][yDim][8];
                         f2[x][yDim][4]=fTmi2[x][yDim][4];
                         f2[x][yDim][7]=fTmi2[x][yDim][7];
                         f2[x][yDim][8]=fTmi2[x][yDim][8];
                 }
      }
//四个角点处理
         f1[1][2][2]=f1[1][2][4];
         f2[1][2][2]=f2[1][2][4];
         f1[1][2][5]=f1[1][2][7];
         f2[1][2][5]=f2[1][2][7];
         f1[xDim][2][2]=f1[1][2][4];
         f2[xDim][2][2]=f2[1][2][4];
         f1[xDim][2][6]=f1[1][2][8];
         f2[xDim][2][6]=f2[1][2][8];
         f1[xDim][yDim-1][4]=f1[1][2][2];
         f2[xDim][yDim-1][4]=f2[1][2][2];
         f1[xDim][yDim-1][7]=f1[1][2][5];
         f2[xDim][yDim-1][7]=f2[1][2][5];
         f1[1][yDim-1][4]=f1[1][2][2];
         f2[1][yDim-1][4]=f2[1][2][2];
         f1[1][yDim-1][8]=f1[1][2][6];
         f2[1][yDim-1][8]=f2[1][2][6];
//求解宏观变量
         for(x=1;x<=xDim;x++)
        {
                for(y=2;y<=yDim-1;y++)
                {
                        rho1[x][y]=f1[x][y][0]+f1[x][y][1]+f1[x][y][2]+f1[x][y][3]+f1[x][y][4]+f1[x][y][5]+f1[x][y][6]+f1[x][y][7]+f1[x][y][8];
                        rho2[x][y]=f2[x][y][0]+f2[x][y][1]+f2[x][y][2]+f2[x][y][3]+f2[x][y][4]+f2[x][y][5]+f2[x][y][6]+f2[x][y][7]+f2[x][y][8];
            u_1[x][y]=(f1[x][y][1]+f1[x][y][5]+f1[x][y][8]-f1[x][y][3]-f1[x][y][6]-f1[x][y][7])/rho1[x][y];
                        u_2[x][y]=(f2[x][y][1]+f2[x][y][5]+f2[x][y][8]-f2[x][y][3]-f2[x][y][6]-f2[x][y][7])/rho2[x][y];
                        v_1[x][y]=(f1[x][y][5]+f1[x][y][6]+f1[x][y][2]-f1[x][y][7]-f1[x][y][8]-f1[x][y][4])/rho1[x][y];
                    v_2[x][y]=(f2[x][y][5]+f2[x][y][6]+f2[x][y][2]-f2[x][y][7]-f2[x][y][8]-f2[x][y][4])/rho2[x][y];
                    phi1[x][y]=1-exp(-rho1[x][y]);
                        phi2[x][y]=1-exp(-rho2[x][y]);
                 }
      }
         
//流体间的作用力                       
        for(x=1;x<=xDim;x++)
        {
                for(y=2;y<=yDim-1;y++)
                {
                xt=0;
                        yt=0;
                    for(k=1;k<Q;k++)
                        {
                                xt=xt+phi1[x+e[k][0]][y+e[k][1]]*e[k][0]*w[k];
                                yt=yt+phi1[x+e[k][0]][y+e[k][1]]*e[k][1]*w[k];
                        }
                        fwx1[x][y]=-phi1[x][y]*gw1*9*xt;
            fwy1[x][y]=-phi1[x][y]*gw1*9*xt;
                        xt=0;
                        yt=0;
                    for(k=1;k<Q;k++)
                        {
                                xt=xt+phi2[x+e[k][0]][y+e[k][1]]*e[k][0]*w[k];
                                yt=yt+phi2[x+e[k][0]][y+e[k][1]]*e[k][1]*w[k];
                        }
                        fwx2[x][y]=-phi1[x][y]*gw2*9*xt;
            fwy2[x][y]=-phi1[x][y]*gw2*9*xt;
                        fintX1[x][y]=fx11[x][y]+fx12[x][y]+fwx1[x][y];
                        fintY1[x][y]=fy11[x][y]+fy12[x][y]+fwy1[x][y];
                        fintX2[x][y]=fx22[x][y]+fx21[x][y]+fwx2[x][y];
                        fintY2[x][y]=fy22[x][y]+fy21[x][y]+fwy2[x][y];

                        u0[x][y]=((rho1[x][y]*u_1[x][y])/tau1+(rho2[x][y]*u_2[x][y])/tau2)/(rho1[x][y]/tau1+rho2[x][y]/tau2);
                        v0[x][y]=((rho1[x][y]*v_1[x][y])/tau1+(rho2[x][y]*v_2[x][y])/tau2)/(rho1[x][y]/tau1+rho2[x][y]/tau2);
                       
                        u_1[x][y]=u0[x][y]+tau1*fintX1[x][y]/rho1[x][y];
                        v_1[x][y]=v0[x][y]+tau1*fintY1[x][y]/rho1[x][y];
                        u_2[x][y]=u0[x][y]+tau2*fintX2[x][y]/rho2[x][y];
                        v_2[x][y]=v0[x][y]+tau2*fintY2[x][y]/rho2[x][y];
                        rhom[x][y]=rho1[x][y]+rho2[x][y];
                }
        }
        for(x=1;x<=xDim;x++)
        {
                for(y=2;y<=yDim-1;y++)
                {               
                for(k=0;k<Q;k++)
                        {
                           zx=0;
                           zy=0;
                           zx=f1[x][y][1]*e[1][0]+f1[x][y][2]*e[2][0]+f1[x][y][3]*e[3][0]+f1[x][y][4]*e[4][0]+f1[x][y][5]*e[5][0]+f1[x][y][6]*e[6][0]+f1[x][y][7]*e[7][0]+f1[x][y][8]*e[8][0]+f1[x][y][0]*e[0][0]+
                                   f2[x][y][1]*e[1][0]+f2[x][y][2]*e[2][0]+f2[x][y][3]*e[3][0]+f2[x][y][4]*e[4][0]+f2[x][y][5]*e[5][0]+f2[x][y][6]*e[6][0]+f2[x][y][7]*e[7][0]+f2[x][y][8]*e[8][0]+f2[x][y][0]*e[0][0];
                           zy=f1[x][y][1]*e[1][1]+f1[x][y][2]*e[2][1]+f1[x][y][3]*e[3][1]+f1[x][y][4]*e[4][1]+f1[x][y][5]*e[5][1]+f1[x][y][6]*e[6][1]+f1[x][y][7]*e[7][1]+f1[x][y][8]*e[8][1]+f1[x][y][0]*e[0][1]+
                                   f2[x][y][1]*e[1][1]+f2[x][y][2]*e[2][1]+f2[x][y][3]*e[3][1]+f2[x][y][4]*e[4][1]+f2[x][y][5]*e[5][1]+f2[x][y][6]*e[6][1]+f2[x][y][7]*e[7][1]+f2[x][y][8]*e[8][1]+f2[x][y][0]*e[0][1];
                    }
                           u_m[x][y]=(zx+(fintX1[x][y]+fintX2[x][y])/2)/rhom[x][y];
                           v_m[x][y]=(zy+(fintY1[x][y]+fintY2[x][y])/2)/rhom[x][y];
                  
                        pre[x][y]=(rho1[x][y]+rho2[x][y])/2+2*(g11*phi1[x][y]*phi1[x][y]+g22*phi2[x][y]*phi2[x][y]+g12*phi1[x][y]*phi2[x][y]+g21*phi1[x][y]*phi2[x][y]);
                }
        }



//碰撞方程

                for(x=1;x<=xDim;x++)
        {
                for(y=2;y<=yDim-1;y++)
                {
                        for(k=0;k<Q;k++)
                        {
                                feq1[x][y][k]=f_eq(k,rho1[x][y],u_1[x][y],v_1[x][y]);
                    feq2[x][y][k]=f_eq(k,rho2[x][y],u_2[x][y],v_2[x][y]);
                                f1[x][y][k]=f1[x][y][k]-tau1*(f1[x][y][k]-feq1[x][y][k]);
                                f2[x][y][k]=f2[x][y][k]-tau2*(f2[x][y][k]-feq2[x][y][k]);       
                        }
                 }
         }
       
         
          

  ostringstream name;
        name<<"cavity_100_"<<iter<<".dat";
        ofstream out(name.str().c_str());
        out<<"Title= \"LBM \"\n"<<"VARIABLES=\"X\",\"Y\",\"u1\",\"v1\",\"um\",\"vm\",\"Rho1\",\"Rho2\",\"f0\",\"f1\",\"f5\",\"f8\",\"4\"\n"<<"ZONE T=\"BOX\",I="<<yDim<<",J="<<xDim<<",F=POINT"<<endl;
        for(x=1;x<=xDim;x++)
        {
                for(y=1;y<=yDim;y++)
                {
                        out<<x<<" "<<y<<" "<<u_1[x][y]<<" "<<v_1[x][y]<<" "<<u_m[x][y]<<" "<<v_m[x][y]<<" "<<rho1[x][y]<<" "<<rho2[x][y]<<" "<<f1[x][y][0]<<" "<<f1[x][y][1]<<" "<<f1[x][y][5]<<" "<<f1[x][y][8]<<" "<<zx<<endl;
                }
    }
}
}
double f_eq(int k,double rho,double u,double v)  /*Equiblibrium distribution*/
{
        double eu,uv,feq;
        eu=(e[k][0]*u+e[k][1]*v);
        uv=(u*u+v*v);
        feq=w[k]*rho*(1.0+3.0*(e[k][0]*u+e[k][1]*v)+4.5*(e[k][0]*u+e[k][1]*v)*(e[k][0]*u+e[k][1]*v)-1.5*(u*u+v*v));
        return feq;
}

城市阳光

你遇到的问题是什么？