请教关于Shan-Chen模型模拟两相分离的问题 程序越界【附源代码】

stusos

小弟刚刚接触LBM，基础较差，最近在尝试用S-C模型模拟两相流分离，但是出现很奇怪的结果  好像是越界了，也找不出哪里的问题，速度和密度变为 1.#QNAN

很是郁闷 找不到问题出在哪里，恳请各位前辈赐教，不胜感激！


下面是程序源码

#include<iostream>
#include<cmath>
#include<cstdlib>
#include<iomanip>
#include<fstream>
#include<sstream>
#include<string>
#include <math.h>

using namespace std;        
const int Q=9;          //D2Q9模型
const int NX=61;                //x方向
const int NY=61;                //y方向


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],phi[NX+1][NY+1],temrho[NX+1][NY+1],temu[NX+1][NY+1][2],u[NX+1][NY+1][2],u0[NX+1][NY+1][2],ueq[NX+1][NY+1][2],f[NX+1][NY+1][Q],F[NX+1][NY+1][Q];
int i,j,k,ip,jp,n;
int bot,top,lef,rig;
double c,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error,grad_phi_x,grad_phi_y;
double  G = -6.0;  

void init();
double feq(int k,double rho, double u[2]);
void evolution();
void output(int m);
void Error();

int main()
{
        using namespace std;
        init();
        for(n=0;;n++)
        {
                evolution();
                if(n%100==0)
                {
                        Error();
                        cout<<"The"<<n<<"th computation result:"<<endl;
                        cout<<"The u,v of point(NX/2,NY/2) is: "<<setprecision(6)
                                <<u[NX/2][NY/2][0]<<", "<<u[NX/2][NY/2][1]<<endl;
                        cout<<"The max relative error of uv is: "
                                <<setiosflags(ios::scientific)<<error<<endl;
                        if(n>=100)
                        {
                                if(n%100==0) output(n);
                                if(error<1.0e-6) break;
                        }
                }
        }
        return 0;
}

void init()
{
        dx=1.0;
        dy=1.0;
        Lx=dx*double(NX);
        Ly=dy*double(NY);
        dt=dx;
        c=dx/dt;
        rho0=1.0;
       
        tau_f=1;
        std::cout<<"tau_f= "<<tau_f<<endl;

        for(i=0;i<=NX;i++)        //初始化
                for(j=0;j<=NY;j++)
                {
                        u[j][0]=0;
                        u[j][1]=0;
                        //rho[j]=rho0;
                        rho[j]=rho0+rand()%10*0.1;  //初始状态   各点随机密度 以便能够运动
                        //cout<<rho[j]<<"\n";
                       
                        for(k=0; k<Q;k++)
                                f[j][k]=feq(k,rho[j],u[j]);
                }
                //output(0);
}
//计算平衡态分布函数
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);
        return feq;
}

void evolution()       
{
        //①计算有效密度
        for(int i=0;i<NX;i++)
                for(int j=0;j<NY;j++)
                {
                        temrho[j]=0;
                        temu[j][0]=0;
                        temu[j][1]=0;
                        for (int k=0;k<Q;k++)
                        {
                                temrho[j]+=f[j][k];                //计算每一点的宏观密度
                                temu[j][0]+=e[k][0]*f[j][k];
                                temu[j][1]+=e[k][1]*f[j][k];
                        }
                        temu[j][0]/=temrho[j];//计算每一点的宏观速度
                        temu[j][1]/=temrho[j];
                        phi[j]= 1.0-exp(-temrho[j]); //计算有效密度   自由参数rho0取1
                }
               

                //②计算平衡态速度
                for(int i=0;i<NX;i++)
                        for(int j=0;j<NY;j++)
                        {
                                bot = (j+NY-1)%NY;
                                top = (j+1)%NY;
                                lef = (i+NX-1)%NX;
                                rig = (i+1)%NX;


                                //计算相邻的格点的有效密度与速度的乘积
                                grad_phi_x = (phi[rig][j]-phi[lef][j])*w[1];
                                grad_phi_y = (phi[top]-phi[bot])*w[1];
                                grad_phi_x+= (phi[rig][top]-phi[lef][top]+phi[rig][bot]-phi[lef][bot])*w[2];
                                grad_phi_y+= (phi[rig][top]+phi[lef][top]-phi[lef][bot]-phi[rig][bot])*w[2];

                                // interparticule potential in equilibrium velocity 计算平衡态速度
                                //tmp_phi为这一点处的有效密度 grad_phi_x 为邻近格子点的有效密度与速度的乘积 tmp_rho 为这一点的宏观密度
                                //ueq[j][0] 为x方向的平衡态速度,ueq[j][1]为y方向的平衡态速度
                                ueq[j][0] =temu[j][0]+ tau_f*(-G*phi[j]*grad_phi_x)/temrho[j];
                                ueq[j][1] =temu[j][1]+tau_f*(-G*phi[j]*grad_phi_y)/temrho[j];
                               

                                //output(j);
                        }


        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];
                                F[j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp],ueq[ip][jp])-f[ip][jp][k])/tau_f;  //碰撞之后的分布函数
                        }

                        for(i=1;i<NX;i++)        //计算宏观量
                                for(j=1;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];
                                        //cout<<u[j][0]<<" "<<u[j][1]<<"\n";
                                }

                                //边界处理，采用非平衡态外推格式
                                for(j=1;j<NY;j++)        //左右边界
                                        for(k=0;k<Q;k++)
                                        {
                                                rho[NX][j]=rho[NX-1][j];
                                                f[NX][j][k]=feq(k,rho[NX][j],ueq[NX][j])+f[NX-1][j][k]-feq(k,rho[NX-1][j],ueq[NX-1][j]);
                                                rho[0][j]=rho[1][j];
                                                f[0][j][k]=feq(k,rho[0][j],ueq[0][j])+f[1][j][k]-feq(k,rho[1][j],ueq[1][j]);
                                        }

                                        for(i=1;i<NX;i++)        //上下边界
                                                for(k=0;k<Q;k++)
                                                {
                                                        rho[0]=rho[1];
                                                        f[0][k]=feq(k,rho[0],ueq[0])+f[1][k]-feq(k,rho[1],ueq[1]);

                                                        rho[NY]=rho[NY-1];
                                                       
                                                        f[NY][k]=feq(k,rho[NY],ueq[NY])+f[NY-1][k]-feq(k,rho[NY-1],ueq[NY-1]);
                                                }

}

void output(int m)        //输出
{
        ostringstream name;
        name<<"cavity_"<<m<<".dat";
        ofstream out(name.str().c_str());
        out<<"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;
        for(j=0; j<=NY; j++)
                for(i=0; i<=NX; i++)
                {
                        out<<i<<" "<<j<<" "
                                <<u[j][0]<<" "<<u[j][1]<<" "<<rho[j]<<endl;
                }
}

void Error()
{
        double temp1,temp2;
        temp1=0;
        temp2=0;
        for(i=1; i<NX; i++)
                for(j=1;j<NY;j++)
                {
                        temp1 += (u[j][0]-u0[j][0])*(u[j][0]-u0[j][0])+(u[j][1]-u0[j][1])*(u[j][1]-u0[j][1]);
                        temp2 += u[j][0]*u[j][0]+u[j][1]*u[j][1];
                }
                temp1=sqrt(temp1);
                temp2=sqrt(temp2);
                error=temp1/(temp2+1e-30);
}

孤独的漫步

LZ，你的i,j循环里面都是只有j参数啊，是不是这个有问题~~~

stusos

首先非常感谢！
但是不是那个原因，你可以下载我的附件查看源代码，源代码里面是有i的，可能是发布的时候给过滤了 [i] 当成 ubb代码 不显示了

longzaijt

你算点间作用力的时候，为什么要算临近点的有效密度的差，应该直接由有效密度计算。而且算两相分离的话可以用两套分布函数

upc_mao

longzaijt 发表于 2016-3-9 19:22
你算点间作用力的时候，为什么要算临近点的有效密度的差，应该直接由有效密度计算。而且算两相分离的话可以 ...

我也有这样的疑惑

qq281362022

问题解决了吗

kdlianglin

请问 您的问题解决了吗？