n-s.c

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

   for( j=NJB;j<(NJ-1);j++)
    {
     U[l][i][j] = R[l][i][j] + TE[i][j] * Q[l][i][j] * FJ[i][j];
    }   
 bond();
 turbmod();
 dissp();
 spacial();
 for( l=0;l<NL;l++)
  for( i=NIB;i<(NI-1);i++)
   for( j=NJB;j<(NJ-1);j++)
   {
    U[l][i][j] = 0.75f * R[l][i][j] + 0.25f * (U[l][i][j] + TE[i][j] * Q[l][i][j] * FJ[i][j]);
   }   

 bond();
 turbmod();
 dissp();
 spacial();
 for( l=0;l<NL;l++)
  for( i=NIB;i<(NI-1);i++)
   for( j=NJB;j<(NJ-1);j++)
   {
    U[l][i][j] = R[l][i][j]/3.0f + (U[l][i][j] + TE[i][j] * Q[l][i][j] * FJ[i][j]) * 2.0f / 3.0f;
   } 


 RES = 0.0;
 for( i=NIB;i<(NI-1);i++)
  for( j=NJB;j<(NJ-1);j++)
  {
   DDQ=fabs(R[0][i][j]-U[0][i][j]);
   if(DDQ > RES)
   {
    RES = DDQ;
    IS = i;
    JS = j;
   }
  }
}
/*---------- disk ----------*/
void disk()
{
 int l,i,j;
 int KI = 0;
 double VELU,VELV,FMACH,TT,TK;
 if((DatFs2 = fopen("fs2.dat","ab")) == NULL)
 {
  printf("\n can not create fs2.dat !");
  exit(1);
 }
 
 for(l=0;l<NL;l++)
  for(i=0;i<NI;i++)
   for(j=0;j<NJ;j++)
   {
    fprintf(DatFs2,"\t%f\t%d\n",U[l][i][j],KI++);
   }
 fclose(DatFs2);
 if((PltFs2 = fopen("fs2.plt","w+")) == NULL)
 {
  printf("\n can not create fs2.plt !");
  exit(1);
 }
 fprintf(PltFs2,"\tVariables = X,Y,U,V,R,E,P,M,T,K\n");
 fprintf(PltFs2,"\n\tZONE T = \"GRID1\" I = %d, J = %d\n",NI,NJ);

 
 for(j=0;j<NJ;j++)
  for(i=0;i<NI;i++)
  {
   VELU = U[1][i][j] / U[0][i][j];
   VELV = U[2][i][j] / U[0][i][j];
   FMACH = sqrt(fabs((VELU * VELU + VELV * VELV) * U[0][i][j] / GAMA / PP[i][j]));
   TT = GAMA * PP[i][j] / U[0][i][j];
   TK = TRK[i][j] / U[0][i][j];
   fprintf(PltFs2,"\t%12.5f\t%12.5f\t%12.5f\t%12.5f\t%12.5f\t%12.5f\t%12.5f\t%12.5f\t%12.5f\t%12.5f\n",XX[i][j],YY[i][j],VELU,VELV,U[0][i][j],U[3][i][j],PP[i][j],FMACH,TT,TK);
  }
 fclose(PltFs2);
}
/*---------- dissp ----------*/
void dissp()
{
 int l,i,j,I1,I2,J1,J2;
 double VF[NI][NJ],F1[NL][NI],F2[NL][NJ];
 double PRPR,DI,DJ,UU,VV,UI,VI,UJ,VJ,BMA0,BMA1,BMA2,BMA3,BMA4,BMC0,BMC1,BMC2,BMC3,BMC4;
 for(i=0;i<NI;i++)
 {
  for(j=0;j<NJ;j++)
  {
   VF[i][j] = ((U[1][i][j] * U[1][i][j]) + (U[2][i][j] * U[2][i][j])) / (U[0][i][j] * U[0][i][j]);
   PP[i][j] = (GAMA - 1.0f) * U[3][i][j] - 0.5f * (VF[i][j] * U[0][i][j]);
  }
 }
 for(l=0;l<NL;l++)
  for(i=0;i<NI;i++)
   for(j=0;j<NJ;j++)
    Q[l][i][j] = 0.0;

 for(j=NJB;j<(NJ-1);j++)
 {
  for(i=0;i<NI;i++)
  {
   PRPR = (TVL[i][j] / PRL) / TV[i][j];
   F1[0][i] = 0.5f * VF[i][j] + PP[i][j] / U[0][i][j] * PRPR * GAMA / (GAMA - 1.0);
  }
  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;   
   }
   J1 = j - 1;
   J2 = j + 1;
   UU = (CX[i][j] * U[1][i][j] + CY[i][j] * U[2][i][j]) / U[0][i][j];
   UI = (U[1][I2][j] / U[0][I2][j] - U[1][I1][j] / U[0][I1][j]) / DI;
   VI = (U[2][I2][j] / U[0][I2][j] - U[2][I1][j] / U[0][I1][j]) / DI;
   UJ = (U[1][i][J2] / U[0][i][J2] - U[1][i][J1] / U[0][i][J1]) / 2.0;
   VJ = (U[2][i][J2] / U[0][i][J2] - U[2][i][J1] / U[0][i][J1]) / 2.0;
   BMC1 = CC[i][j] * CC[i][j];
   BMC2 = CX[i][j] * UI + CY[i][j] * VI;
   BMC3 = (F1[0][I2] - F1[0][I1]) / DI;
   BMA0 = (CX[i][j] * AX[i][j] + CY[i][j] * AY[i][j]) * PRPR / (GAMA - 1.0);
   BMA1 = CX[i][j] * AX[i][j] * 4.0f / 3.0f + CY[i][j] * AY[i][j];
   BMA2 = CY[i][j] * AX[i][j] - CX[i][j] * AY[i][j] * 2.0f / 3.0f;
   BMA3 = CX[i][j] * AY[i][j] - CY[i][j] * AX[i][j] * 2.0f / 3.0f;
   BMA4 = CX[i][j] * AX[i][j] + CY[i][j] * AY[i][j] * 4.0f / 3.0f;
   F1[1][i] = (BMC1 * UI + BMC2 * CX[i][j] / 3.0f + BMA1 * UJ + BMA2 * VJ) * TV[i][j] / FJ[i][j];
   F1[2][i] = (BMC1 * VI + BMC2 * CY[i][j] / 3.0f + BMA3 * UJ + BMA4 * VJ) * TV[i][j] / FJ[i][j];
   F1[3][i] = (BMC1 * BMC3 + BMC2 * UU / 3.0f + ((BMA1 * U[1][i][j] + BMA3 * U[2][i][j]) * UJ + (BMA2 * U[1][i][j] + BMA4 * U[2][i][j]) * VJ) / U[0][i][j] + 0.5f * BMA0 * GAMA * (PP[i][J2] / U[0][i][J2] - PP[i][J1] / U[0][i][J1])) * TV[i][j] / FJ[i][j];
  
  }
  
  for(l=1;l<NL;l++)
  {
   for(i=NIB;i<(NI-1);i++)
   {
    if(i == NIB || i == (NI - 1))
    {
     Q[l][i][j] = Q[l][i][j] + 0.5f * (F1[l][i+1] - F1[l][i-1]);
    }
    else{
    Q[l][i][j] = Q[l][i][j] + (8.0f * (F1[l][i+1] - F1[l][i-1]) + F1[l][i-2]-F1[l][i+2])/12.0f;
   }
   }
  }
 }

 for(i=NIB;i<(NI-1);i++)
 {
  for(j=0;j<NJ;j++)
  {
   PRPR = (TVL[i][j] / PRL) / TV[i][j];
   F2[0][j] = 0.5f * VF[i][j] + PP[i][j] / U[0][i][j] * PRPR * GAMA / (GAMA - 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;   
   }
   I1 = i - 1;
   I2 = i + 1;
   VV = (AX[i][j] * U[1][i][j] + AY[i][j] * U[2][i][j]) / U[0][i][j];
   UI = (U[1][I2][j] / U[0][I2][j] - U[1][I1][j] / U[0][I1][j]) / 2.0;
   VI = (U[2][I2][j] / U[0][I2][j] - U[2][I1][j] / U[0][I1][j]) / 2.0;
   UJ = (U[1][i][J2] / U[0][i][J2] - U[1][i][J1] / U[0][i][J1]) / DJ;
   VJ = (U[2][i][J2] / U[0][i][J2] - U[2][i][J1] / U[0][i][J1]) / DJ;
   BMA1 = AA[i][j] * AA[i][j];
   BMA2 = AX[i][j] * UI + AY[i][j] * VI;
   BMA3 = (F2[0][J2] - F2[0][J1]) / DJ;
   BMC0 = (AX[i][j] * CX[i][j] + AY[i][j] * CY[i][j]) * PRPR / (GAMA - 1.0);
   BMC1 = AX[i][j] * CX[i][j] * 4.0f / 3.0f + AY[i][j] * CY[i][j];
   BMC2 = AY[i][j] * CX[i][j] - AX[i][j] * CY[i][j] * 2.0f / 3.0f;
   BMC3 = AX[i][j] * CY[i][j] - AY[i][j] * CX[i][j] * 2.0f / 3.0f;
   BMC4 = AX[i][j] * CX[i][j] + AY[i][j] * CY[i][j] * 4.0f / 3.0f;
   F2[1][j] = (BMC1 * UI + BMC2 * VI + BMA1 * UJ + BMA2 * AX[i][j] / 3.0f) * TV[i][j] / FJ[i][j];
   F2[2][j] = (BMC3 * UI + BMC4 * VI + BMA1 * VJ + BMA2 * AY[i][j] / 3.0f) * TV[i][j] / FJ[i][j];
   F2[3][j] = (BMA1 * BMA3 + BMA2 * VV / 3.0f + ((BMC1 * U[1][i][j] + BMC3 * U[2][i][j]) * UI + (BMC2 * U[1][i][j] + BMC4 * U[2][i][j]) * VI) / U[0][i][j] + 0.5f * BMC0 * GAMA * (PP[I2][j] / U[0][I2][j] - PP[I1][j] / U[0][I1][j])) * TV[i][j] / FJ[i][j];
  }
  for(l=1;l<NL;l++)
  {
   for(j=NJB;j<(NJ-1);j++)
   {
    if(j == NJB || j == (NJ - 1))
    {
     Q[l][i][j] = Q[l][i][j] + 0.5f * (F2[l][j+1] - F2[l][j-1]);
    }
    else{
    Q[l][i][j] = Q[l][i][j] + (8.0f * (F2[l][j+1] - F2[l][j-1]) + F2[l][j-2]-F2[l][j+2])/12.0f;
   }
   }
  }
 }

 for(l=0;l<NL;l++)
  for(i=NIB;i<(NI-1);i++)
   for(j=NJB;j<(NJ-1);j++)
    Q[l][i][j] = Q[l][i][j] / REN;
}
/*---------- spacial ----------*/
void spacial()
{ 
 double UU,VV,CF,FLAMDA;
 double EPSLON = 0.01f * EPU * EPU;
 double F1[NL][NI],F2[NL][NJ];
 int l,i,j;
 
 for( j=NJB;j<(NJ-1);j++)
  {
   for( i=0;i<NI;i++)
    {
     UU = (CX[i][j] * U[1][i][j] + CY[i][j] * U[2][i][j]) / U[0][i][j];
     CF = sqrt(GAMA * PP[i][j] / U[0][i][j]);
     FLAMDA = CC[i][j] * (sqrt((UU * UU) / (CC[i][j] * CC[i][j]) + (EPSLON * CF) * (EPSLON * CF)) + CF);
     F1[0][i] = 0.5f * (UU * U[0][i][j] + FLAMDA * U[0][i][j]) / FJ[i][j];
     F1[1][i] = 0.5f * (UU * U[1][i][j] + PP[i][j] * CX[i][j] + FLAMDA * U[1][i][j]) / FJ[i][j];
     F1[2][i] = 0.5f * (UU * U[2][i][j] + PP[i][j] * CY[i][j] + FLAMDA * U[2][i][j]) / FJ[i][j];
     F1[3][i] = 0.5f * (UU * (U[3][i][j] + PP[i][j]) + FLAMDA * U[3][i][j]) / FJ[i][j];
     
    }
   for( i=NIB;i<(NI-1);i++)
    {
     for( l=0;l<NL;l++)
     {
      Q[l][i][j] = Q[l][i][j] - (F1[l][i] - F1[l][i-1]);
     } 
    }
  }

 for( j=NJB;j<(NJ-1);j++)
  {
   for( i=0;i<NI;i++)
   {
    UU = (CX[i][j] * U[1][i][j] + CY[i][j] * U[2][i][j]) / U[0][i][j];
    CF = sqrt(GAMA * PP[i][j] / U[0][i][j]);
    FLAMDA = CC[i][j] * (sqrt((UU * UU) / (CC[i][j] * CC[i][j]) + (EPSLON * CF) * (EPSLON * CF)) + CF);
    F1[0][i] = 0.5f * (UU * U[0][i][j] - FLAMDA * U[0][i][j]) / FJ[i][j];
    F1[1][i] = 0.5f * (UU * U[1][i][j] + PP[i][j] * CX[i][j] - FLAMDA * U[1][i][j]) / FJ[i][j];
    F1[2][i] = 0.5f * (UU * U[2][i][j] + PP[i][j] * CY[i][j] - FLAMDA * U[2][i][j]) / FJ[i][j];
    F1[3][i] = 0.5f * (UU * (U[3][i][j] + PP[i][j]) - FLAMDA * U[3][i][j]) / FJ[i][j];
   }
   for( i=NIB;i<(NI-1);i++)
    {
     for( l=0;l<NL;l++)
     {
      Q[l][i][j] = Q[l][i][j] - (F1[l][i+1] - F1[l][i]);
     }
    }
  }
 
 for( i=NIB;i<(NI-1);i++)
  {
   for( j=0;j<NJ;j++)
    {
     VV = (AX[i][j] * U[1][i][j] + AY[i][j] * U[2][i][j]) / U[0][i][j];
     CF = sqrt(GAMA * PP[i][j] / U[0][i][j]);
     FLAMDA = AA[i][j] * (sqrt((VV * VV) / (AA[i][j] * AA[i][j]) + (EPSLON * CF) * (EPSLON * CF)) + CF);
     F2[0][j] = 0.5f * (VV * U[0][i][j] + FLAMDA * U[0][i][j]) / FJ[i][j];
     F2[1][j] = 0.5f * (VV * U[1][i][j] + PP[i][j] * AX[i][j] + FLAMDA * U[1][i][j]) / FJ[i][j];
     F2[2][j] = 0.5f * (VV * U[2][i][j] + PP[i][j] * AY[i][j] + FLAMDA * U[2][i][j]) / FJ[i][j];
     F2[3][j] = 0.5f * (VV * (U[3][i][j] + PP[i][j]) + FLAMDA * U[3][i][j]) / FJ[i][j];
    }
   for( j=NJB;j<(NJ-1);j++)
    {
     for( l=0;l<NL;l++)
     {
      Q[l][i][j] = Q[l][i][j] - (F2[l][j] - F2[1][j-1]);
     }
    }
  }

 for( i=NIB;i<(NI-1);i++)
  { for( j=0;j<NJ;j++)
   {
    VV = (AX[i][j] * U[1][i][j] + AY[i][j] * U[2][i][j]) / U[0][i][j];
    CF = sqrt(GAMA * PP[i][j] / U[0][i][j]);
    FLAMDA = AA[i][j] * (sqrt((VV * VV) / (AA[i][j] * AA[i][j]) + (EPSLON * CF) * (EPSLON * CF)) + CF);
    F2[0][j] = 0.5f * (VV * U[0][i][j] - FLAMDA * U[0][i][j]) / FJ[i][j];
    F2[1][j] = 0.5f * (VV * U[1][i][j] + PP[i][j] * AX[i][j] - FLAMDA * U[1][i][j]) / FJ[i][j];
    F2[2][j] = 0.5f * (VV * U[2][i][j] + PP[i][j] * AY[i][j] - FLAMDA * U[2][i][j]) / FJ[i][j];
    F2[3][j] = 0.5f * (VV * (U[3][i][j] + PP[i][j]) - FLAMDA * U[3][i][j]) / FJ[i][j];
   }
   for( j=NJB;j<(NJ-1);j++)
   { 
    for( l=0;l<NL;l++)
    {
      Q[l][i][j] = Q[l][i][j] - (F2[l][j+1] - F2[l][j]);
    }
   }
  }
}
/*---------- The End ----------*/