|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
以下是本人所写的迎风型三阶紧致格式计算不可压缩流体的Poiseulle流动程序（交错网格），但是无奈是发散的，很着急，很崩溃。
请求在线的大牛帮我调试一下。我的联系方式：dingzijing09@mails.gucas.ac.cn
#include"stdlib.h"
#include"iostream.h"
#include"stdio.h"
#include"math.h"
#define R 100
#define Dt 0.0002
#define eps 1e-6
#define x 101 //x方向节点数目
#define y 21 //y方向节点数目
double U[x+2][y+2],V[x+2][y+2],P[x+2][y+2],tu[x+2][y+2],tv[x+2][y+2],fp1[x+2][y+2],ep1[x+2][y+2],fm1[x+2][y+2],em1[x+2][y+2],fp2[x+2][y+2],ep2[x+2][y+2],fm2[x+2][y+2],em2[x+2][y+2];
//Poiseulle流动的速度边界条件（交错网格，边界落在压强的网格上，U在P的左右侧，V在P的上下侧，在边界外面歌设置了一层虚拟网格）
void bounduv(double U[x+2][y+2],double V[x+2][y+2])
{
int i,j;
//水平边界的速度边界条件
for(i=0;i<=x+1;i++)
{
U[1]=0;
U[y]=0;
V[1]=-V[0];
V[y]=-V[y-1];
}
//入口和出口的压力边界条件
for(j=1;j<=y;j++)
{
U[0][j]=U[1][j];
U[x][j]=U[x-1][j];
V[1][j]=V[0][j];
V[x][j]=V[x-1][j];
}
}
void boundp(double P[x+2][y+2])
{
int i,j;
//入口和出口的压力边界条件
for(j=1;j<=y;j++)
{
P[1][j]=1.0;
P[x][j]=0.0;
}
//法向边界的压力条件
for(i=1;i<=x;i++)
{
P[2]=P[1];
P[y]=P[y-1];
}
}
//初始化
void init(double U[x+2][y+2],double V[x+2][y+2],double P[x+2][y+2])
{
int i,j;
for(i=0;i<=x+1;i++)
{
for(j=0;j<=y+1;j++)
{
U[j]=0;
V[j]=0;
}
}
for(i=1;i<=x;i++)
{
for(j=1;j<=y;j++)
{
P[j]=0;
}
}
bounduv(U,V);
boundp(P);
}
//三阶紧致格式求解F+，F-，E+，E-（分别对应fp1,fm1,ep1,em1;fp2,fm2,ep2,em2;前者是关于U的紧致差分矩阵，后者是关于V的紧致差分矩阵）
void solfe(double p1[x+2][y+2],double p2[x+2][y+2],double p3[x+2][y+2],double p4[x+2][y+2],double s[x+2][y+2])
{
int i,j;
for (j=1;j<=y;j++)
{
p1[1][j]=(3.0*(s[2][j]-s[1][j])-(s[3][j]-s[2][j]))/2.0;
p2[x][j]=(3.0*(s[x][j]-s[x-1][j])-(s[x-1][j]-s[x-2][j]))/2.0;
}
for(i=1;i<=x;i++)
{
p3[1]=(3.0*(s[2]-s[1])-(s[3]-s[2]))/2.0;
p4[y]=(3.0*(s[y]-s[y-1])-(s[y-1]-s[y-2]))/2.0;
}
for(i=2;i<=x;i++)
for(j=1;j<=y;j++)
{ p1[j]=-p1[i-1][j]/2.0+(s[i+1][j]+4.0*s[j]-5.0*s[i-1][j])/4.0;
p1[x+1][j]=p1[x][j];
}
for(i=x-1;i>=1;i--)
for(j=1;j<=y;j++)
{ p2[j]=-p2[i+1][j]/2.0+(5.0*s[i+1][j]-4.0*s[j]+s[i-1][j])/4.0;
p2[0][j]=p2[1][j];
}
for(j=2;j<=y;j++)
for(i=1;i<=x;i++)
{ p3[j]=-p3[j-1]/2.0+(s[j+1]+4.0*s[j]-5.0*s[j-1])/4.0;
p3[y]=p3[y];
}
for(j=y-1;j>=1;j--)
for(i=0;i<=x+1;i++)
{ p4[j]=-p4[j+1]/2.0+(5.0*s[j+1]-4.0*s[j]+5.0*s[j-1])/4.0;
p4[0]=p4[1];
}
}
//求解程序
void soluv(double lamb,double r,double dt,double dx,double dy,double U[x+2][y+2],double V[x+2][y+2],double P[x+2][y+2],double tu[x+2][y+2],double tv[x+2][y+2],double fp1[x+2][y+2],double fm1[x+2][y+2],double ep1[x+2][y+2],double em1[x+2][y+2],double fp2[x+2][y+2],double fm2[x+2][y+2],double ep2[x+2][y+2],double em2[x+2][y+2])
{
int i,j;
double uav,vav,adv,adv1,adv2,adv3,adv4,prs,vis;
solfe(fp1,fm1,ep1,em1,U);//求解U方向的紧致差分用的F,E
solfe(fp2,fm2,ep2,em2,V);//求解V方向的紧致差分用的F,E
//交错网格的差分算法
for(i=1;i<x;i++)
{
for(j=2;j<y;j++)
{
vav=(V[i+1][j-1]+V[i+1][j+1]+V[i-1][j-1]+V[i-1][j+1])/4.0;
adv1=(U[j]+fabs(U[j]))/2.0*fp1[j]/dx;
adv2=(U[j]-fabs(U[j]))/2.0*fm1[j]/dx;
adv3=(vav+fabs(vav))/2.0*ep1[j]/dy;
adv4=(vav-fabs(vav))/2.0*em1[j]/dy;
adv=-adv1-adv2-adv3-adv4;
vis=((U[i+1][j]-2.0*U[j]+U[i-1][j])/dx/dx+(U[j+1]-2.0*U[j]+U[j-1])/dy/dy)/r;
prs=-(P[i+1][j]-P[j])/dx;
tu[j]=(adv+vis+prs)*dt+U[j];
}
}
for(i=2;i<x;i++)
{
for(j=1;j<y;j++)
{
uav=(U[i-1][j+1]+U[i-1][j-1]+U[i+1][j+1]+U[i+1][j-1])/4.0;
adv1=(uav+fabs(uav))/2.0*fp2[j]/dx;
adv2=(uav-fabs(uav))/2.0*fm2[j]/dx;
adv3=(V[j]+fabs(V[j]))/2.0*ep2[j]/dy;
adv4=(V[j]-fabs(V[j]))/2.0*em2[j]/dy;
adv=-adv1-adv2-adv3-adv4;
vis=((V[i+1][j]-2.0*V[j]+V[i-1][j])/dx/dx+(V[j+1]-2.0*V[j]+V[j-1])/dy/dy)/r;
prs=-(P[j+1]-P[j])/dy;
tv[j]=(adv+vis+prs)*dt+V[j];
}
}
for(i=1;i<=x;i++)
{
for(j=1;j<=y;j++)
{
U[j]=tu[j];
V[j]=tv[j];
}
}
bounduv(U,V);
for(i=1;i<=x;i++)
{
for(j=1;j<=y;j++)
{
P[j]=P[j]-lamb*((U[j]-U[i-1][j])/dx+(V[j]-V[j-1])/dy);
}
}
boundp(P);
}
//收敛的判断程序
double conv(double U[x+2][y+2],double V[x+2][y+2],double dm,double dx,double dy)
{
int i,j;
double er;
dm=0;
for(i=1;i<=x;i++)
{
for(j=1;j<=y;j++)
{
er=fabs((U[j]-U[i-1][j])/dx+(V[j]-V[j-1])/dy);
if(dm<er)
{
dm=er;
}
}
}
return dm;
}
void output(char *sn,double U[x+2][y+2],double V[x+2][y+2],double P[x+2][y+2],double dx,double dy)
{
int i,j;
FILE *fp;
fp=fopen(sn,"w");
fprintf(fp,"Title=\"computation result\"\n");
fprintf(fp,"vars=\"x\",\"y\",\"u\",\"v\",\"p\"\n");
fprintf(fp,"zone t=\"zone1\",I=%d,J=%d,F=POINT\n",y,x);
for(i=1;i<=x;i++)
{
for(j=1;j<=y;j++)
{
fprintf(fp,"%16f%16f%20e%20e%20e\n",(i-1)*dx,(j-1)*dy,(U[i-1][j]+U[j])/2,(V[j-1]+V[j])/2,P[j]);
}
}
fclose(fp);
fp=fopen("result.txt","w");
i=(x+1)/2;
for(j=1;j<=y;j++)
{
fprintf(fp,"%16f%20e%20e\n",(j-1)*dy-1.0,U[j],P[j]);
}
fclose(fp);
}
void main()
{
int n;
double dx,dy,lamb,er;
double dm;
dx=1.0/(x-1);
dy=0.20/(y-1);
init(U,V,P);
lamb=1.5;
n=0;
do
{
soluv(lamb,R,Dt,dx,dy,U,V,P,tu,tv,fp1,fm1,ep1,em1,fp2,fm2,ep2,em2);
er=conv(U,V,dm,dx,dy);
n++;
cout<<n<<endl;
cout<<"er="<<er<<endl;
// cout<<fm1[10][(y-1)/2]<<endl;
}while(er>1e-3);
output("result.plt",U,V,P,dx,dy);
} |
|
IP卡
狗仔卡
发表于 2010-5-26 08:01:33
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
楼主
发表于 2010-5-26 10:11:13
