求助！！！新手的poiseuille流算出来老是发散，求大神指导
#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 = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
double w = {
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,u,u0,f,F,force;
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);
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<<","<<u<<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 = 0;
u = 0;
rho = rho0;
u = 4.0*U/(Ly*Ly)*(Ly*j-j*j);
for(k=0; k<Q; k++)
{
f=feq(k,rho,u);
}
}
}
double feq(int k,double rho,double u)
{
double eu,uv,feq;
eu = (e*u + e*u);
uv = (u*u + u*u);
feq = w*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=f+(feq(k,rho,u)-f)/tau_f;
}
for (i = 0; i <= NX; i++)
for (j = 1; j < NY; j++)
for (k = 0; k < Q; k++)
{
ip=int(i-e+NX)%NX ;
jp=j-e;
f = F;
}
i=0;
for(j=1; j<=NY-1; j++)
for(k=0; k<Q; k++)
{
u = 4.0*U/(Ly*Ly)*(Ly*j-j*j);
u = 0.0;
rho = 1.0/(1.0-u)*(f+f+f+2.0*(f+f+f));
f = f+2.0/3.0*rho*u;
f = f+1.0/6.0*rho*u+0.5*(f-f)+1.0/2.0*rho*u;
f = f+1.0/6.0*rho*u-0.5*(f-f)-1.0/2.0*rho*u;
}
for (i=0; i<=NX; i++)
for(k=0;k<Q;k++)
{
rho=rho;
rho=rho;
f = feq(k,rho,u)+f-feq(k,rho,u);
f = feq(k,rho,u)+f - feq(k, rho, u);
}
for (i = 0; i <= NX; i++)
for (j=1; j <NY; j++)
{
u0 = u;
u0 = u;
dense = 0;
u_x = 0;
u_y = 0;
for (k = 0; k < Q; k++)
{
dense += f;
u_x += e * f;
u_y += e * f;
}
rho = dense;
u = u_x;
u = u_y;
u /= rho;
u /= rho;
}
}
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 << " " << u << endl;
}
}
void Error()
{
double temp1,temp2;
temp1 = 0;
temp2 = 0;
for(i=0; i<=NX; i++)
for(j=0; j<=NY; j++)
{
temp1 += ((u - u0)*(u - u0) + (u - u0)*(u - u0));
temp2 += (u*u + u*u);
}
temp1 = sqrt(temp1);
temp2 = sqrt(temp2);
error = temp1/(temp2);
}
试过用zou/he处理角点还是不行。。。现在也不知道问题出在哪里
好像有两个问题
第一，松弛时间太大，二点几不靠谱吧
第二，为什么左右既有周期边界又有速度边界？
最后建议，有问题 上图先 看程序太累了
孤独的漫步 发表于 2017-1-9 09:36
好像有两个问题
第一，松弛时间太大，二点几不靠谱吧
第二，为什么左右既有周期边界又有速度边界？
谢谢大神，我修改下，现在图有点离谱，下次会注意的:lol
页:
[1]