⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 unit1.cpp

📁 气体反吹模拟程序
💻 CPP
字号:
🔍 
//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "Unit1.h"
#include "math.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
#define pi 3.1415926
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------

void __fastcall TForm1::SimulationClick(TObject *Sender)
{
float PP[20],P[20],MP[20],M[20],OP[20],OM[20],D[2],F[2];
float T1[3],C1[2],C2[2];
float AB[2],AA[3][4],XA[3];
for(int i=0;i<20;i++)
{
PP[i]=0.0;
P[i]=0.0;
MP[i]=0.0;
M[i]=0.0;
OP[i]=0.0;
OM[i]=0.0;
}
for(int i=0;i<2;i++)
{
D[i]=0.0;
F[i]=0.0;
C1[i]=0.0;
C2[i]=0.0;
AB[i]=0.0;
}
for(int i=0;i<3;i++)
{
T1[i]=0.0;
XA[i]=0.0;
}
for(int i=0;i<3;i++)
{
 for(int j=0;j<4;j++)
 {
  AA[i][j]=0.0;
 }
}

float PAT;
      PAT=100000.0;//大气压
int N1,N2,N2S;
      N1=4;
      N2=8;
      N2S=N1+N2;
float RAIR;
     RAIR=287.0;//气体常数
  //rem air vessel parameter
float V0,P0,T0,EK[8],ZAIR;
 T0=300.0;//温度
 ZAIR=1.4;//压缩系数
 EK[0]=1.4;
 P0=700000.0;//进口气压
 V0=0.0302; //容积
  //rem nozzel parameter
float XT1,DD,BOUE;
 XT1=1.6;
 DD=0.006;
 BOUE=1.0;
 //rem two pipe prarmeter
 D[0]=0.025;
 D[1]=0.025;
 F[0]=0.013;
 F[1]=0.014;
 //rem pulse valve prarmeters
float TOU,TUP,TDOWN,DF,BOUV,JPR,TMAX,GR,G,TOLM,TOLP,IDG,ALM;
   TOU=0.046;
   TUP=0.15215;
   TDOWN=0.3;
   TOU=TUP+TDOWN+TOU;
   T1[0]=TUP;
   T1[1]=TOU-TDOWN;
   T1[2]=TOU;
  BOUV=1.0;
 JPR=3.0;
 TMAX=2.0;
 GR=0.6;
 G=9.81;
 TOLM=0.002;
 TOLP=40.0;
 IDG=0.0;
 ALM=50.0;
 //rem initial condition
float M0,B;
       M0=P0*V0/(RAIR*T0);//PV=MRT,求管内气体质量
       StringGrid_TMP->Cells[0][0]="M0="+FloatToStr(M0);
       B=sqrt(ZAIR*RAIR*T0)/1.2964; //等温过程,B=sqrt(ZRT)
       StringGrid_TMP->Cells[1][0]="B="+FloatToStr(B);
float  DDM;
     DDM=0.0;
float P0K1,SUS,ECR,AAE,AAV;
   EK[1]=1.0/EK[0];
   P0K1=pow(P0,EK[1]);
   EK[2]=EK[0]/(EK[0]-1.0);
   EK[3]=2.0/(EK[0]+1.0);
   EK[4]=sqrt(2.0*EK[0]/(EK[0]-1.0));
   SUS=sqrt(EK[3]*EK[0]*pow(EK[3],2*EK[2]/EK[0]));
   EK[5]=(EK[0]+1.0)/EK[0];
   EK[6]=2.0/EK[0];
   ECR=pow(2.0/(EK[0]+1.0),EK[2]);
   EK[7]=(EK[0]-1.0)/EK[0];
   AAE=0.25*XT1*pi*DD*DD;
   AAV=0.25*pi*DF*DF;
float DT=0.001;
  for(int i=0;i<N1;i++)
  {
  P[i]=P0;
  OP[i]=P0;
  M[i]=0.0;
  OM[i]=0.0;
  }
float XX=1.0;
   C1[0]=0.25*F[0]*pow(B,2)*DT*XX/(pow(0.25*pi,2)*pow(D[0],5));
   C1[1]=0.25*F[1]*pow(B,2)*DT*XX/(pow(0.25*pi,2)*pow(D[1],5));
   //print c1[0] c2[0]
    StringGrid_TMP->Cells[2][0]="C1[0]="+FloatToStr(C1[0]);
    StringGrid_TMP->Cells[3][0]="C1[1]="+FloatToStr(C1[1]);
    C2[0]=0.0;
    C2[1]=0.0;
for(int i=N1;i<N2S;i++)
    {
      M[i]=0.0;
      OM[i]=0.0;
      P[i]=PAT;
      OP[i]=PAT;
    }
//rem caculate alpha
  AB[0]=B/(0.25*pi*D[0]*D[0]);
  AB[1]=B/(0.25*pi*D[1]*D[1]);
  StringGrid_TMP->Cells[4][0]="AB[0]="+FloatToStr(AB[0]);
  StringGrid_TMP->Cells[5][0]="AB[1]="+FloatToStr(AB[1]);
float T=0.0;
int J=0;

//计算
float PLP=0.0,PKP=0.0,F1=0.0,F1P=0.0,F1M=0.0,PRAI=0.0,PED=0.0,POPP=0.0,F2=0.0;
float F2P=0.0,F2M=0.0,DP=0.0,DM=0.0,TEM=0.0,CONZ2=0.0,CONZ1=0.0;
float BOU=0.0,CVAL1=0.0,CVAL2=0.0,PJP=0.0,PUP=0.0,F1P2=0.0,F2P2=0.0,F3=0.0,F3P=0.0,F3M=0.0,F3P2=0.0;
float DP2=0.0,DEM=0.0;

Step1:
       T=T+DT;
       J=J+1;
       //上游边界条件
       if(T>TMAX) return;


       //rem upstream boundary
       MP[0]=2.0*M[0]-OM[0];
       PP[0]=2.0*P[0]-OP[0];

         for(int ki=0;ki<5;ki++)
           {
           PLP=C1[0]*(MP[0]+M[1])*fabs(MP[0]+M[1]);
           F1=AB[0]*(MP[0]-M[1])*(PP[0]+P[1])-PP[0]*PP[0]+(1.0+C2[0])*P[1]*P[1]+PLP ;
           F1P=AB[0]*(MP[0]-M[1])-2.0*PP[0];
           F1M=AB[0]*(PP[0]+P[1])+2.0*C1[0]*fabs(MP[0]+M[1]);
           F2=M0*P0K1-P0K1*(DDM+MP[0]*DT)-M0*pow(PP[0],EK[1]);
           F2P=-EK[1]*M0*pow(PP[0],(EK[1]-1.0));
           F2M=-DT*P0K1;
           DP=(F1/F1M-F2/F2M)/(F2P/F2M-F1P/F1M);
           DM=-(F1+F1P*DP)/F1M;
           PP[0]=PP[0]+DP;
           MP[0]=MP[0]+DM;
           if(fabs(DP)<TOLP&&fabs(DM)<TOLM)
                 break ;
           }

        DDM=DDM+MP[0]*DT;
        TEM=T0*pow(PP[0]/P0,EK[7]);
        CONZ1=BOUE*AAE*EK[4]/sqrt(RAIR*TEM);
        CONZ2=BOUE*AAE*SUS/sqrt(RAIR*TEM);
       //rem downstream boundary condition
        MP[N2S-1]=2*M[N2S-1]-OM[N2S-1];
        PP[N2S-1]=2*P[N2S-1]-OP[N2S-1];

    if(P[N2S-2]>PAT)
      {
         if(PP[N2S-1]<=PAT)  PP[N2S-1]=P[N2S-2];
         for(int ki=0;ki<15;ki++)
         {
          PKP=-P[N2S-2]*P[N2S-2]+C1[1]*(MP[N2S-1]+M[N2S-2]*fabs(MP[N2S-1]+M[N2S-2]));
          F1=AB[1]*(MP[N2S-1]-M[N2S-2])*(PP[N2S-1]+P[N2S-2])+(1.0+C2[1])*PP[N2S-1]*PP[N2S-1]+PKP;
          F1P=AB[1]*(MP[N2S-1]-M[N2S-2])+2.0*(1.0+C2[1])*PP[N2S-1];
          F1M=AB[1]*(PP[N2S-1]+P[N2S-2])+2.0*C1[1]*fabs(MP[N2S-1]+M[N2S-2]);
          PRAI=PAT/PP[N2S-1];

          //PP=input(pp);
          PED=pow(PRAI,EK[6])-pow(PRAI,EK[5]);
          if(PRAI>1.0&&PRAI<1.002)   PED=fabs(PED);
           POPP=sqrt(PED);
           F2=MP[N2S-1]-CONZ1*PP[N2S-1]*POPP;
           //input I
           F2P=CONZ1*(-POPP-0.5/POPP*(-EK[6]*pow(PRAI,EK[6])+EK[5]*pow(PRAI,EK[5])));
           F2M=1.0;
           if(PRAI<=ECR)
            {
             F2=MP[N2S-1]-CONZ2*PP[N2S-1];
             F2P=-CONZ2;
             F2M=1.0;
            }
           DP=(F1/F1M-F2/F2M)/(F2P/F2M-F1P/F1M);
           DM=-(F1+F1P*DP)/F1M;
           PP[N2S-1]=PP[N2S-1]+DP;
           MP[N2S-1]=MP[N2S-1]+DM;
           if(fabs(DP)<40.0&&fabs(DM)<0.002)
           break;
         }
      }
     else
     {
        PP[N2S-1]=P[N2S-2];
        MP[N2S-1]=0.0;
     }
     //1200 rem pulse valve
      BOU=BOUV;
       if(T<=T1[0]) BOU=T/T1[0]*BOUV;
       if(T>T1[1]&&T<T1[2]) BOU=BOUV*(1.0-(T-T1[1])/TDOWN);
       if(T>=T1[2]) BOU=0.0;
       CVAL1=BOU*AAV*EK[4]/sqrt(RAIR*TEM);
       CVAL2=BOU*AAV*SUS/sqrt(RAIR*TEM);
       MP[N1-1]=2*M[N1-1]-OM[N1-1];
       MP[N1]=MP[N1-1];
       PP[N1-1]=2*P[N1-1]-OP[N1-1];
       PP[N1]=2*P[N1]-OP[N1];
     for(int ki=0;ki<15;ki++)
       {
       PJP=-P[N1-2]*P[N1-2]+C1[0]*(MP[N1-1]+M[N1-2])*fabs(MP[N1-1]+M[N1-2]);
       F1=AB[0]*(MP[N1-1]-M[N1-2])*(PP[N1-1]+P[N1-2])+(1+C2[0])*PP[N1-1]*PP[N1-1]+PJP;
       F1P=AB[0]*(MP[N1-1]-M[N1-2])+2.0*(1+C2[0])*PP[N1-1];
       F1M=AB[0]*(PP[N1-1]+P[N1-2])+2.0*C1[0]*fabs(MP[N1-1]+M[N1-2]);
       F1P2=0.0;
       PUP=AB[1]*(MP[N1]-M[N1+1])*(PP[N1]+P[N1+1])-PP[N1]*PP[N1];
       F2=PUP+(1.0+C2[1])*P[N1+1]*P[N1+1]+C1[1]*(MP[N1]+M[N1+1])*fabs(MP[N1]+M[N1+1]);
       F2P2=AB[1]*(MP[N1]-M[N1+1])-2.0*PP[N1];
       F2M=AB[1]*(PP[N1]+P[N1+1])+2.0*C1[1]*fabs(MP[N1]+M[N1+1]);
       F2P=0.0;
       PRAI=PP[N1]/PP[N1-1];
       POPP=sqrt(pow(PRAI,EK[6])-pow(PRAI,EK[5]));
       F3=MP[N1-1]-CVAL1*PP[N1-1]*POPP;
       F3P=CVAL1*(-POPP-0.5/POPP*(-EK[6]*pow(PRAI,EK[6])+EK[5]*pow(PRAI,EK[5])));
       F3M=1.0;
       F3P2=-CVAL1*PP[N1-1]/POPP*(EK[6]*pow(PRAI,EK[6])-EK[5]*pow(PRAI,EK[5]))/PP[N1];
       if(PRAI<=ECR)
        {
         F3=MP[N1-1]-CVAL2*PP[N1-1];
         F3P=-CVAL2;
         F3M=1.0;
         F3P2=0.0;
        }
       AA[0][0]=F1P;AA[0][1]=F1M;AA[0][2]=F1P2;AA[0][3]=-F1;
       AA[1][0]=F2P;AA[1][1]=F2M;AA[1][2]=F2P2;AA[1][3]=-F2;
       AA[2][0]=F3P;AA[2][1]=F3M;AA[2][2]=F3P2;AA[2][3]=-F3;
      //1530goto sub3000
       for(int i=0;i<3;i++)
       {
        XA[i]=0.0;
       }
       for(int k=0;k<3;k++)
       {
         //GOSUB 5000
         DD=AA[k][k];
         int L;
         L=k;
         for(int i=k+1;i<3;i++)
          {
          if(fabs(AA[i][k])<=fabs(DD))
            continue;
          DD=AA[i][k];
          L=i;
          }
         if(DD<=0.001)
         {
          //no solution
         //ShowMessage("No Solution");
         goto Step2;
         }
        if(L!=k)
         {
          for(int j=k;j<4;j++)
          {
            float TT;
            TT=AA[L][j];
            AA[L][j]=AA[k][j] ;
            AA[k][j]=TT;
          }
         }
         for(int j=k+1;j<4;j++)
         {
         AA[k][j]=AA[k][j]/AA[k][k];
         }
         for(int i=k+1;i<3;i++)
         {
           for(int j=k+1;j<4;j++)
           {
            AA[i][j]=AA[i][j]-AA[i][k]*AA[k][j];
           }
         }
      }//k=0:3
       XA[2]=AA[2][3];
       for(int i=1;i>=0;i--)
       {
         for(int j=i+1;j<3;j++)
         {
           XA[i]=XA[i]+AA[i][j]*XA[j];
         }
         XA[i]=AA[i][3]-XA[i];
       }
     

     //ShowMessage(FloatToStr(P[N2S-1]));
     //ShowMessage(FloatToStr(MP[N2S-1]));
    //计算脉冲阀
     //1550
Step2:
     DP=XA[0];DM=XA[1];DP2=XA[2];
     PP[N1-1]=PP[N1-1]+DP;
     MP[N1-1]=MP[N1-1]+DM;
     PP[N1]=PP[N1]+DP2;
     MP[N1]=MP[N1-1];
     if(fabs(DP)<TOLP&&fabs(DM)<TOLM) break;
   }//ki=1:15

    //rem interior sections
    //rem before valve
   for(int i=1;i<N1-1;i++)
    {
    PP[i]=2*P[i]-OP[i];
    MP[i]=2*M[i]-OM[i];
     for(int ki=0;ki<15;ki++)
     {
      PJP=-P[i-1]*P[i-1]+C1[0]*(MP[i]+M[i-1])*fabs(MP[i]+M[i-1]);
      F1=AB[0]*(MP[i]-M[i-1])*(PP[i]+P[i-1])+(1+C2[0])*PP[i]*PP[i]+PJP;
      F1P=AB[0]*(MP[i]-M[i-1])+2.0*(1+C2[0])*PP[i];
      F1M= AB[0]*(PP[i]+P[i-1])+2.0*C1[0]*fabs(MP[i]+M[i-1]);
      PUP=AB[0]*(MP[i]-M[i+1])*(PP[i]+P[i+1])-PP[i]*PP[i];
      F2=PUP+(1.0+C2[0])*P[i+1]*P[i+1]+C1[0]*(MP[i]+M[i+1])*fabs(MP[i]+M[i+1]);
      F2P=AB[0]*(MP[i]-M[i+1])-2*PP[i];
      F2M=AB[0]*(PP[i]+P[i+1])+2.0*C1[0]*fabs(MP[i]+M[i+1]);
      DP=(F1/F1M-F2/F2M)/(F2P/F2M-F1P/F1M);
      DM=-(F1+F1P*DP)/F1M;
      PP[i]=PP[i]+DP;
      MP[i]=MP[i]+DM;
      if(fabs(DP)<TOLP&&fabs(DM)<TOLM) break;
     }
    }
    //rem after valve
Step3:

     for(int i=N1+1;i<N2S-2;i++)
     {
      PP[i]=2*P[i]-OP[i];
      MP[i]=2*M[i]-OM[i];
      for(int ki=0;ki<15;ki++)
       {
         PJP=-P[i-1]*P[i-1]+C1[1]*(MP[i]+M[i-1])*fabs(MP[i]+M[i-1]);
         F1=AB[1]*(MP[i]-M[i-1])*(PP[i]+P[i-1])+(1.0+C2[1])*PP[i]*PP[i]+PJP;
         F1P=AB[1]*(MP[i]-M[i-1])+2.0*(1.0+C2[1])*PP[i];
         F1M= AB[1]*(PP[i]+P[i-1])+2.0*C1[1]*fabs(MP[i]+M[i-1]);
         PUP=AB[1]*(MP[i]-M[i+1])*(PP[i]+P[i+1])-PP[i]*PP[i];
         F2=PUP+(1.0+C2[1])*P[i+1]*P[i+1]+C1[1]*(MP[i]+M[i+1])*fabs(MP[i]+M[i+1]);
         F2P=AB[1]*(MP[i]-M[i+1])-2.0*PP[i];
         F2M=AB[1]*(PP[i]+P[i+1])+2.0*C1[1]*fabs(MP[i]+M[i+1]);
         DP=(F1/F1M-F2/F2M)/(F2P/F2M-F1P/F1M);
         DM=-(F1+F1P*DP)/F1M;
         PP[i]=PP[i]+DP;
         MP[i]=MP[i]+DM;
         if(fabs(DP)<TOLP&&fabs(DM)<TOLM) break;
       }
      }
   //input "after valve"
     for(int i=0;i<N2S;i++)
     {
      OM[i]=M[i];
      M[i]=MP[i];
      OP[i]=P[i];
      P[i]=PP[i];
     }
     StringGrid_TMP->Cells[0][J]="T*1000="+FloatToStr(T*1000.0);
     StringGrid_TMP->Cells[1][J]="P[N1-1]="+FloatToStr(P[N1-1]);
     StringGrid_TMP->Cells[2][J]="P[N1]="+FloatToStr(P[N1]);
     StringGrid_TMP->Cells[2][J]="M[N1]*100000.0="+FloatToStr(M[N1]*100000.0);
     StringGrid_TMP->Cells[3][J]="P[N2S-1]="+FloatToStr(P[N2S-1]);
     StringGrid_TMP->Cells[4][J]="M[N2S-1]*100000.0="+FloatToStr(M[N2S-1]*100000.0);

   //print T*1000,P[N1-1],P[N1],P[N2S-1],M[N2S-1]*1000;
    if(T<TOU+0.2) goto Step1;
    //gosub 7670
   T=T+DT;
   J=J+1;
   //rem upstream boundary conditions
   MP[N1]=0.0;
   PP[N1]=2.0*P[N1]-OP[N1];
   for(int ki=0;ki<5;ki++)
   {
    PLP=C1[1]*(MP[N1]+M[N1+1])*fabs(MP[N1]+M[N1+1]);
    F1=AB[1]*(MP[N1]-M[N1+1])*(PP[N1]+P[N1+1])-PP[N1]*PP[N1]+(1.0+C2[1])*P[N1+1]*P[N1+1]+PLP;
    DF=AB[0] *(MP[N1]-M[N1+1])-2.0*PP[N1];
    DP=-F1/DF;
    PP[N1]=PP[N1]+DP;
    if(fabs(DP)<TOLP) break;
   }
   CONZ1=BOUE*AAE*EK[4]/sqrt(RAIR*TEM);
   CONZ2=BOUE*AAE*SUS/sqrt(RAIR*TEM);
   //rem downstream boundary conditions
   MP[N2S-1]=2.0*M[N2S-1]-OM[N2S-1];
   PP[N2S-1]=2.0*P[N2S-1]-OP[N2S-1];
   if(P[N2S-1]>PAT)
   {
    for(int ki=0;ki<15;ki++)
    {
    PKP=-P[N2S-2]*P[N2S-2]+C1[1]*(MP[N2S-1]+M[N2S-2])*fabs(MP[N2S-1]+M[N2S-2]);
    F1=AB[1]*(MP[N2S-1]-M[N2S-2])*(PP[N2S-1]+P[N2S-2])+(1.0+C2[1])*PP[N2S-1]*PP[N2S-1]+PKP;
    F1P=AB[1]*(MP[N2S-1]-M[N2S-2])+2.0*(1.0+C2[1])*PP[N2S-1];
    F1M=AB[1]*(PP[N2S-1]+P[N2S-2])+2.0*C1[1]*fabs(MP[N2S-1]+M[N2S-2]);
    PRAI=PAT/PP[N2S-1];
    if(PRAI>=1.0) break;
    POPP=sqrt(pow(PRAI,EK[6])-pow(PRAI,EK[5]));
    F2=MP[N2S-1]-CONZ1*PP[N2S-1]*POPP;
    F2P=CONZ1*(-POPP-0.5/POPP*(-EK[6]*pow(PRAI,EK[6])+EK[5]*pow(PRAI,EK[5])));
    F2M=1.0;
    if(PRAI<=ECR)
    {
     F2=MP[N2S-1]-CONZ2*PP[N2S-1];
     F2P=-CONZ2;
     F2M=1.0;
    }
    DP=(F1/F1M-F2/F2M)/(F2P/F2M-F1P/F1M);
    DM=-(F1+F1P*DP)/F1M;
    PP[N2S-1]=PP[N2S-1]+DP;
    MP[N2S-1]=MP[N2S-1]+DM;
    if(fabs(DP)<40.0&&fabs(DM<TOLM)) break;
    }
   }
   PP[N2S-1]=PAT;
   MP[N2S-1]=0.0;
  if(T>TOU&&T<(TOU+0.5))
     goto Step3;
   DEM=M0*pow((1.0-(P[0])),EK[1]);
   //print P0,DDM,P[0],DEM,TOU,TUP

}
//---------------------------------------------------------------------------

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -