求助！！！新手的poiseuille流算出来老是发散，求大神指导

muzhiyu

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>

using namespace std;
const int Q = 9;  
const int NX = 192;
const int NY = 32;
const double U = 0.1;

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],f[NX+1][NY+1][Q],F[NX+1][NY+1][Q],force[NX+1][NY+1][Q];

int i,j,k,ip,jp,n;
double c,Re,dx,dy,Lx,Ly,dt,rho0,P0,tau_f,niu,error,dense,u_x,u_y;

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=1; ;n++)
{

evolution();
if(n%1000 == 0)
{

Error();
cout <<"The"<<ends<<n<<"th computation result:"<<endl<<"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>= 1000)
{

if(n%1000==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; //1.0
rho0 = 1.0;
tau_f=2.41;
niu=(tau_f-0.5)/3.0;
std::cout << "tau_f=" << tau_f << endl;

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] = 4.0*U/(Ly*Ly)*(Ly*j-j*j);
for(k=0; k<Q; k++)
{
f[i][j][k]=feq(k,rho[i][j],u[i][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);
return feq;
}
void evolution()
{
        for (i=0; i<=NX;i++)
                for (j=0;j<=NY;j++)
                        for (k=0;k<Q;k++)
                        {

                                F[i][j][k]=f[i][j][k]+(feq(k,rho[i][j],u[i][j])-f[i][j][k])/tau_f;
                               
                        }
for (i = 0; i <= NX; i++)
                for (j = 1; j < NY; j++)
                        for (k = 0; k < Q; k++)
                        {
                        ip=int(i-e[k][0]+NX)%NX ;                        
                        jp=j-e[k][1];
                        f[i][j][k] = F[ip][jp][k];
                        }
        i=0;
for(j=1; j<=NY-1; j++)
for(k=0; k<Q; k++)
{

u[i][j][0] = 4.0*U/(Ly*Ly)*(Ly*j-j*j);
u[i][j][1] = 0.0;


rho[i][j] = 1.0/(1.0-u[i][j][0])*(f[i][j][0]+f[i][j][2]+f[i][j][4]+2.0*(f[i][j][3]+f[i][j][6]+f[i][j][7]));
f[i][j][1] = f[i][j][3]+2.0/3.0*rho[i][j]*u[i][j][0];
f[i][j][5] = f[i][j][7]+1.0/6.0*rho[i][j]*u[i][j][0]+0.5*(f[i][j][4]-f[i][j][2])+1.0/2.0*rho[i][j]*u[i][j][1];
f[i][j][8] = f[i][j][6]+1.0/6.0*rho[i][j]*u[i][j][0]-0.5*(f[i][j][4]-f[i][j][2])-1.0/2.0*rho[i][j]*u[i][j][1];


}
        for (i=0; i<=NX; i++)
                for(k=0;k<Q;k++)
                {
            
                        rho[i][0]=rho[i][1];
                        rho[i][NY]=rho[i][NY-1];
                        f[i][0][k] = feq(k,rho[i][0],u[i][0])+f[i][1][k]-feq(k,rho[i][1],u[i][1]);
                    f[i][NY][k] = feq(k,rho[i][NY],u[i][NY])+f[i][NY-1][k] - feq(k, rho[i][NY-1], u[i][NY-1]);       
                }
                        for (i = 0; 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];
                        dense = 0;
                        u_x = 0;
                        u_y = 0;

                        for (k = 0; k < Q; k++)
                        {

                               
                                dense += f[i][j][k];
                                u_x += e[k][0] * f[i][j][k];
                                u_y += e[k][1] * f[i][j][k];

                        }

                        rho[i][j] = dense;
                        u[i][j][0] = u_x;
                        u[i][j][1] = u_y;

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

                }

       




       
}

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\"\n" << "ZONE T=\"BOX\",I=" << NX + 1 << ",J=" << NY + 1 << ",F=POINT" << endl;
                for (i = 0; i <= NX; i++)
                        for (j = 0; j <= NY; j++)
                {
                        out << double(i) / Lx << " "
                                << double(j) / Ly/6.0 << " "   << u[i][j][0] << " " << u[i][j][1] << endl;
                }
}




void Error()
{

double temp1,temp2;
temp1 = 0;
temp2 = 0;
for(i=0; i<=NX; i++)
   for(j=0; 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);

}

muzhiyu

试过用zou/he处理角点还是不行。。。现在也不知道问题出在哪里

孤独的漫步

好像有两个问题
第一，松弛时间太大，二点几不靠谱吧
第二，为什么左右既有周期边界又有速度边界？

最后建议，有问题 上图先 看程序太累了

muzhiyu

孤独的漫步 发表于 2017-1-9 09:36
好像有两个问题
第一，松弛时间太大，二点几不靠谱吧
第二，为什么左右既有周期边界又有速度边界？

谢谢大神，我修改下，现在图有点离谱，下次会注意的