n-s.c

四阶龙格库塔法求解流体力学-- 关于N-S方程的串行求解源程序

！！！为保护知识产权，由于本文的源程序也有我导师的辛劳，所以有些地方做了修改！！！-----因此此程序是运行不起来的，要你自己看看然后做些修改才可以使用

 

/*----------header file----------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*----------define paramant----------*/
#define NIB 1
#define NJB 1

#define NL  4
#define NI 45
#define NJ 50

#define ETA 0.001
#define XMA  3.5f
#define GAMA 1.4f

#define FJETRR 7.01217
#define FJETRV 11.88855
#define FJETPP 20.16412
#define EPU 0.0001f
#define MAXKI 2250

#define CL 0.4469
#define CMIU 0.09
#define PRL 0.72
#define PRT 0.9
#define REN 1823.4


/*----------general paramment ----------*/
double XX[NI][NJ];
double YY[NI][NJ];

double U[NL][NI][NJ],R[NL][NI][NJ];
double TE[NI][NJ];
double Q[NL][NI][NJ];

int flag = 1;
int KI = 0;
double RES,DDQ;
int IS,JS;

double CX[NI][NJ],CY[NI][NJ],AX[NI][NJ],AY[NI][NJ],AA[NI][NJ],CC[NI][NJ];
double FJ[NI][NJ];

double TRK[NI][NJ],PP[NI][NJ],TRE[NI][NJ];

double TVL[NI][NJ],TV[NI][NJ];

FILE *Res,*DatFs2,*DatGrid,*PltFs2;
/*----------sub function----------*/
void grid();
void jacob();
void init_fs2();
void rksol();
void bond();
void turbmod();
void dissp();
void spacial();
void disk();
void fs2();

/*---------- main ----------*/
int main()
{
 
 grid();
 fs2();
 printf("\n thank you !");
 return (0);
}

/*---------- grid ----------*/
void grid()
{
 int i,j;
 double x,y;
 if((DatGrid = fopen("grid.dat","r+")) == NULL)
 {
  printf("\n can not open this file !");
 }
 for(j=0;j<NJ;j++)
  for(i=0;i<NI;i++)
  {
   fscanf(DatGrid,"%lf%lf",&x,&y);
   XX[i][j] = x;YY[i][j] = y;
  }
 fclose(DatGrid);
}

/*---------- fs2 ----------*/
void fs2()
{ 
 jacob();
 
 init_fs2();
 
 if((Res = fopen("res.dat","w+")) == NULL)
 {
  printf("\n can not create res.dat !");
  exit(1);
  
 }
 printf("\n\tKI\tMax Residue\tPosition\n");
 fprintf(Res,"\tKI\tMax Residue\tPosition\n");

 do
 {  
  rksol();
 
  KI = KI + 1;
  
  printf("\t%d\t%16.9f\t%d\t%d\n",KI,RES,IS,JS);
  fprintf(Res,"\t%d\t%16.9f\t%d\t%d\n",KI,RES,IS,JS);
  if(RES <= EPU)
  {
   printf("\tKI=%d\tReidue<Epsilon=%f\n\t!!STOP!!\n",KI,RES);
   fprintf(Res,"\tKI=%d\tReidue<Epsilon=%f\n\t!!STOP!!\n",KI,RES);
   flag = 0;
   break;
  }
  if(KI >= MAXKI)
  {
   printf("\tKI>%d\tReidue=%f\n\t!!STOP!!\n",MAXKI,RES);
   fprintf(Res,"\tKI>%d\tReidue=%f\n\t!!STOP!!\n",MAXKI,RES);  
   flag = 0;
   break;
  }
  if(flag == 0)break; 
 }while(flag != 0);
 
 disk();
}
/*---------- jacob ----------*/
/*----------------------------------------------------------------
 求坐标变换的雅可比矩阵
----------------------------------------------------------------*/

void jacob()
{
 double XI,YI,XJ,YJ,DI,DJ,COB;
 int I1,I2,J1,J2,i,j;
 
 for(i=0;i<NI;i++)
 {
  I1 = i - 1;
  I2 = i + 1;
  DI = 2.0;
  if(i == 0)
  {
   I1 = i;
   DI = 1.0;
  }
  if(i == (NI-1))
  {
   I2 = i;
   DI = 1.0;   
  }
  for( j=0;j<NJ;j++)
  {
   J1 = j - 1;
   J2 = j + 1;
   DJ = 2.0;
   if(j == 0)
   {
    J1 = j;
    DJ = 1.0;
   }
   if(j == (NJ-1))
   {
    J2 = j;
    DJ = 1.0;   
   }
   XI = (XX[I2][j] - XX[I1][j]) / DI;
   YI = (YY[I2][j] - YY[I1][j]) / DI;
   XJ = (XX[i][J2] - XX[i][J1]) / DJ;
   YJ = (YY[i][J2] - YY[i][J1]) / DJ;
   CX[i][j] = YJ;
   CY[i][j] = -XJ;
   AX[i][j] = -YI;
   AY[i][j] = XI;
   COB = XI * YJ - XJ * YI;
   FJ[i][j] = (double)(1.0 / COB);
  }
 }
 for( i=0;i<NI;i++)
  for( j=0;j<NJ;j++)
  {
   CX[i][j] = CX[i][j] * FJ[i][j];
   CY[i][j] = CY[i][j] * FJ[i][j];
   AX[i][j] = AX[i][j] * FJ[i][j];
   AY[i][j] = AY[i][j] * FJ[i][j];
   CC[i][j] = sqrt(CX[i][j] * CX[i][j] + CY[i][j] * CY[i][j]);
   AA[i][j] = sqrt(AX[i][j] * AX[i][j] + AY[i][j] * AY[i][j]);
  }
 for( i=0;i<NI;i++)
  for( j=0;j<NJ;j++)
  {
   TE[i][j] = ETA;
  }
 
}
/*---------- init_fs2 ----------*/
/*---------------------------------------------------------------
 赋初场
----------------------------------------------------------------*/
void init_fs2()
{ 
 double RRF = 1.0;
 double VXF = XMA;
 double VYF = 0.0;
 double PPF = 1.0 / GAMA;
 double ENF = PPF / (GAMA - 1.0) + 0.5f * (XMA * XMA);
 double RKF = 0.01f * (VXF * VXF + VYF * VYF);
 double REF = RRF * CMIU * pow((RKF / RRF),1.5f) / 0.03f;
 int i,j;
 
 for( i=0;i<NI;i++)
  for( j=0;j<NJ;j++)
  {
   U[0][i][j] = RRF;
   U[1][i][j] = VXF;
   U[2][i][j] = VYF;    
   U[3][i][j] = ENF;
   TRK[i][j] = RKF;
   TRE[i][j] = REF;
   PP[i][j] = PPF;
  }
}


/*---------- bond ----------*/
/*---------------------------------------------------------------
 边界条件
 物面取无滑移条件；
 上、后远场外推；
 前远场取来流值；
----------------------------------------------------------------*/
void bond()
{
 int l,i,j,k1,k2;
 for( i=0;i<NI;i++)
 {
  U[0][i][0] = U[0][i][1];
  U[1][i][0] = 0.0;
  U[2][i][0] = 0.0;
  U[3][i][0] = U[3][i][1] - 0.5f * (U[1][i][1] * U[1][i][1] + U[2][i][1] * U[2][i][1]) / U[0][i][1];
  if(XX[i][0] > -0.5 && XX[i][0] < 0.5)
  {
   U[0][i][0] = FJETRR;
   U[1][i][0] = 0.0;
   U[2][i][0] = FJETRV;
   U[3][i][0] = FJETPP / GAMA / (GAMA - 1.0) + 0.5f * U[2][i][0] * U[2][i][0] / FJETRR;
  }
 }
 k1 = NJ -1;k2 = NJ - 2;
 for( l=0;l<NL;l++)
  for( i=0;i<NI;i++)
  {
   
   U[l][i][k1] = U[l][i][k2];
  }
 k1 = NI -1;k2 = NI -2;
 for( l=0;l<NL;l++)
  for( j=0;j<NJ;j++)
  {
   
   U[l][k1][j] = U[l][k2][j];
  }
}

/*---------- turbmod ----------*/
void turbmod()
{
 int i,j;
 double TT;
 for(i=0;i<NI;i++)
  for(j=0;j<NJ;j++)
  {
   TT = GAMA * PP[i][j] / U[0][i][j];if(TT == 0){printf("\t%f\n",TT);getchar();}
   TVL[i][j] = pow(TT,1.5f) * (1.0f + CL) / (TT + CL);
   //TV[i][j] = TVL[i][j];
   TV[i][j] = 1.0;
  }

}

/*---------- rksol ----------*/
/*---------------------------------------------------------------
 三阶Runge-Kutta法求解N-S方程
----------------------------------------------------------------*/
void rksol()
{
 int l,i,j;
 
 bond();
 for(l=0;l<NL;l++)
  for(i=0;i<NI;i++)
   for(j=0;j<NJ;j++)
   {
    R[l][i][j] = U[l][i][j];
   }   
 
 turbmod();
 dissp();
 spacial();
 
 for( l=0;l<NL;l++)
  for( i=NIB;i<(NI-1);i++)