关于三阶紧致格式编写方腔流问题

liuheshan1223

本人最近在编写关于三阶紧致差分计算方腔流的问题，编写的程序检查了好几遍了，感觉没什么错误，但就是不收敛，结果也不对，快崩溃了！哪位牛人有时间帮忙看下，小弟不甚感激！
   程序如下：
   
//三阶迎风紧致差分计算方腔流，人工伪压缩算法，交错网格
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define C 1.5
#define Re 100.0
#define Lx 1.0  //方腔流x方向的长度
#define Ly 1.0  //方腔流y方向的长度
#define dt 0.0004
#define Ix 100  //x方向网格数
#define Jy 100  //y方向网格数
#define xl 0.1//收敛控制参数


//全局变量

double Fux[2][Ix+1][Jy+2],Fuy[2][Ix+1][Jy+2],Fvx[2][Ix+2][Jy+1],Fvy[2][Ix+2][Jy+1];//对流项的F的三阶紧致差分矩阵，2代表正负，0为正，1为负
double u[Ix+1][Jy+2],v[Ix+2][Jy+1],p[Ix+2][Jy+2],dx,dy;//流场的各种参数




/*边界条件*/

void bounduv()
{

  int i,j;


  for(i=1;i<=Ix-1;i++)
{

        u[0]=-u[1];           //下边界
    u[Jy+1]=2.0-u[Jy];    //上边界

}

    for(j=1;j<=Jy;j++)
{
    u[0][j]=u[Ix][j]=0.0;//左边界 ,右边界
}

for(i=1;i<=Ix;i++)
{
  v[0]=v[Jy]=0.0;//下边界 ,上边界
}

for(j=1;j<=Jy-1;j++)
{
   v[0][j]=-v[1][j];          //左边界
   v[Ix][j]=-v[Ix-1][j]; //左边界
}
}


//压力边界条件
void boundp()
{
int i,j;
        for(i=1;i<=Ix;i++)
  {
    p[0]=p[1];
    p[Jy+1]=p[Jy];
  }

  for(j=1;j<=Jy;j++)
  {
         p[0][j]=p[1][j];
     p[Ix+1][j]=p[Ix][j];
  }  

}

/*初始化*/

void Init()
{
int i,j;
//速度全部初始化为0
for(i=0;i<=Ix;i++)
   for(j=0;j<=Jy+1;j++)
      u[j]=0.0;


for(i=0;i<=Ix+1;i++)
   for(j=0;j<=Jy;j++)
     v[j]=0.0;

   bounduv();

//压力的初值设为1
for(i=0;i<=Ix+1;i++)
   for(j=0;j<=Jy+1;j++)
       p[j]=1.0;

    boundp();
}



//对流项的F的三阶紧致差分格式，2代表正负，0为正，1为负
void Q_FS()
{
   int i,j;
   
   for(j=0;j<Jy+2;j++)
      Fux[0][0][j]=0.5*(3*(u[1][j]-u[0][j])-(u[2][j]-u[1][j]));

   for(i=0;i<Ix+1;i++)
      Fuy[0][0]=0.5*(3*(u[1]-u[0])-(u[2]-u[1]));

   for(j=0;j<Jy+1;j++)
      Fvx[0][0][j]=0.5*(3*(v[1][j]-v[0][j])-(v[2][j]-v[1][j]));

   for(i=0;i<Ix+2;i++)
      Fvy[0][0]=0.5*(3*(v[1]-v[0])-(v[2]-v[1]));

   for(j=0;j<Jy+2;j++)
      Fux[1][Ix][j]=0.5*(3*(u[Ix][j]-u[Ix-1][j])-(u[Ix-1][j]-u[Ix-2][j]));

   for(i=0;i<Ix+1;i++)
      Fuy[1][Jy+1]=0.5*(3*(u[Jy+1]-u[Jy])-(u[Jy]-u[Jy-1]));

   for(j=0;j<Jy+1;j++)
      Fvx[1][Ix+1][j]=0.5*(3*(v[Ix+1][j]-v[Ix][j])-(v[Ix][j]-v[Ix-1][j]));

   for(i=0;i<Ix+2;i++)
      Fvy[1][Jy]=0.5*(3*(v[Jy]-v[Jy-1])-(v[Jy-1]-v[Jy-2]));


   for(j=0;j<Jy+2;j++)
           for(i=1;i<Ix;i++)
           {
                   Fux[0][j]=0.5*((0.5*u[i+1][j]+2*u[j]-2.5*u[i-1][j])-Fux[0][i-1][j]);
           Fux[0][Ix][j]=Fux[0][Ix-1][j];
           }

   for(j=0;j<Jy+1;j++)
           for(i=1;i<Ix+1;i++)
           {
                   Fvx[0][j]=0.5*((0.5*v[i+1][j]+2*v[j]-2.5*v[i-1][j])-Fvx[0][i-1][j]);
           Fvx[0][Ix+1][j]=Fux[0][Ix][j];           
           }

   for(i=0;i<Ix+1;i++)
           for(j=1;j<Jy+1;j++)
           {
                   Fuy[0][j]=0.5*((0.5*u[j+1]+2*u[j]-2.5*u[j-1])-Fuy[0][j-1]);
           Fuy[0][Jy]=Fuy[0][Jy-1];
           }

   for(i=0;i<Ix+2;i++)
           for(j=1;j<Jy;j++)
           {
                   Fvy[0][j]=0.5*((0.5*v[j+1]+2*v[j]-2.5*v[j-1])-Fvy[0][j-1]);
           Fvy[0][Jy+1]=Fuy[0][Jy];                
           }

   for(j=0;j<Jy+2;j++)
           for(i=Ix-1;i>0;i--)
           {
                   Fux[1][j]=0.5*((2.5*u[i+1][j]-2*u[j]-0.5*u[i-1][j])-Fux[1][i+1][j]);
               Fux[1][0][j]=Fux[1][1][j];
           }

   for(j=0;j<Jy+1;j++)
           for(i=Ix;i>0;i--)
           {
                   Fvx[1][j]=0.5*((2.5*v[i+1][j]-2*v[j]-0.5*v[i-1][j])-Fvx[1][i+1][j]);
               Fvx[1][0][j]=Fvx[1][1][j];
       }

   for(i=0;i<Ix+1;i++)
           for(j=Jy;j>0;j--)
           {
                   Fuy[1][j]=0.5*((2.5*u[j+1]-2*u[j]-0.5*u[j-1])-Fuy[1][j+1]);
               Fuy[1][0][j]=Fuy[1][1][j];               
           }
   for(i=0;i<Ix+2;i++)
           for(j=Jy-1;j>0;j--)
           {
                  Fvy[1][j]=0.5*((2.5*v[j+1]-2*v[j]-0.5*v[j-1])-Fvy[1][j+1]);
              Fvy[1][0][j]=Fvy[1][1][j];
       }

}


//三阶紧致差分格式主算子，其中粘性项采用中心差分
void Compact_solver()
{
        int i,j;
    double  adu,adv;

               
                Q_FS();


        for(i=1;i<Ix;i++)
                for(j=1;j<Jy+1;j++)
                {
                        adv=0.25*(v[j+1]+v[j-1]+v[i+1][j]+v[i-1][j]);
                        u[j]=u[j]+dt*((1.0/Re)*((u[i+1][j]-2*u[j]+u[i-1][j])/(dx*dx)+(u[j+1]-2*u[j]+u[j-1])/(dy*dy))-(p[j]-p[i-1][j])/dx-
                                  (adv-fabs(adv))/2*Fuy[1][j]/dy-(adv+fabs(adv))/2*Fuy[0][j]/dy-
                                          (u[j]-fabs(u[j]))/2*Fux[1][j]/dx-(u[j]+fabs(u[j]))/2*Fux[0][j]/dx);
                }

        for(i=1;i<Ix+1;i++)
                for(j=1;j<Jy;j++)
                {
                        adu=0.25*(u[j+1]+u[j-1]+u[i+1][j]+u[i-1][j]);
                        v[j]=v[j]+dt*((1.0/Re)*((v[i+1][j]-2*v[j]+v[i-1][j])/(dx*dx)+(v[j+1]-2*v[j]+v[j-1])/(dy*dy))-(p[j]-p[j-1])/dy-
                                  (adu-fabs(adu))/2*Fvx[1][j]/dx-(adu+fabs(adu))/2*Fvx[0][j]/dx-
                                          (v[j]-fabs(v[j]))/2*Fvx[1][j]/dy-(v[j]+fabs(v[j]))/2*Fvx[0][j]/dy);
        }
     
      bounduv();

        for(i=1;i<Ix+1;i++)
                for(j=1;j<Jy+1;j++)
                        p[j]=p[j]-dt*C*C*((u[j]-u[i-1][j])/dx+(v[j]-u[j-1])/dy);

      boundp();

}

/*输出结果*/

void Output()
{
int i,j;
FILE *fp;


fp=fopen("result.txt","w");

for(i=0;i<Ix+1;i++)
        for(j=0;j<Jy+1;j++)
{
  fprintf(fp,"%20f%20f%20.10e%20.10e%20.10en",i*dx,j*dy,0.5*(u[j]+u[j+1]),0.5*(v[j]+v[i+1][j]),0.25*(p[j]+p[i+1][j+1]+p[i+1][j]+p[j+1]));
}
   fclose(fp);
}


/*主函数*/

void main()
{
    int i,j;
        double T,N,m=0,d;
        double ut[Ix+2][Jy+2],vt[Ix+2][Jy+2],pt[Ix+2][Jy+2];

        dx=Lx/Ix;
        dy=Ly/Jy;

        Init();
       
        T=0;
        N=0;

        do
        {
                N++;
                T+=dt;

                printf("N=%10g      T=%10fn",N,T);
   

        for(i=0;i<Ix+2;i++)
                for(j=0;j<Jy+2;j++)
                {
                   ut[j]=u[j];
                   vt[j]=v[j];
                   pt[j]=p[j];
                }               


        Compact_solver();
       printf("%20.10e%20.10e%20.10en",u[2][2],v[2][2],p[2][2]);      

        for(i=1;i<Ix+1;i++)
                for(j=1;j<Jy+1;j++)
                {
               
                d=max(max(fabs(u[j]-ut[j])/dt,fabs(v[j]-vt[j])/dt),fabs(pt[j]-p[j])/(C*C*dt));
        
        if(d>m)
                     m=d;               
               
                }
    printf("%20.10en",m);               
        }while(m>=xl&&N<100);
   
   Output();
}

shirazbj

c的程序在xp下用什么编译器呀？

for(i=1;i<=Ix-1;i++)
{

        u[0]=-u[1];           //下边界
    u[Jy+1]=2.0-u[Jy];    //上边界

}

这u到底是一维数组还是两维数组？还有v?

[ 本帖最后由 shirazbj 于 2010-6-2 22:31 编辑 ]

liuheshan1223

不好意思不太明白您说的是什么意思！如果是压力的算法，我采用的是人工伪压缩！

liuheshan1223

晕，怎么贴上去都成了这样了，都是二维的数组！

liuheshan1223

终于做出来了，哈哈！幸福啊！有个地方搞错了！