Sources Esope | Cast3M

deto
C DETO      SOURCE    CB215821  20/11/25    13:24:54     10792          
      SUBROUTINE DETO
C-----------------------------------------------------------------------
C Calcul les conditions Chapman-Jouguet (CJ), Adiabatic Isochore
C Complete Combustion (AICC) et Zeldovitch-Neuman-Doringts (ZND)
C d'un mélange O2/N2/H2/H2O.
C-----------------------------------------------------------------------
C
C---------------------------
C Phrase d'appel (GIBIANE) :
C---------------------------
C
C CHP2 CHP3 CHP4 = DETO CHP1
C
C------------------------
C Opérandes et résultat :
C------------------------
C
C CHP1 : CHAMPOINT, composition du mélange, pression et température
C        composantes : 'O2  ', 'N2  ', 'H2  ', 'H2O ', 'P   ', 'T   '
C CHP2 : CHAMPOINT, conditions Chapman-Jouguet (CJ)
C        composantes : 'RCJ ', 'TCJ ', 'PCJ ', 'TAUX', 'VCJ '
C CHP3 : CHAMPOINT, idem Zeldovitch-Neuman-Doringts (ZND)
C        composantes : 'RZND', 'TZND', 'PZND'
C CHP4 : CHAMPOINT, idem Adiabatic Isochore Complete Combustion (AICC)
C        composantes : 'RAIC', 'TAIC', 'PAIC'
C
C-----------------------------------------------------------------------
C
C Langage : ESOPE + FORTRAN77
C
C Auteurs : F.DABBENE et E.STUDER 02/99
C
C-----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*4 LCOMP(6)
      PARAMETER (NBC=4,NBT=14)
      DIMENSION AN(NBC),AH2(NBT),AH2O(NBT),AN2(NBT),AO2(NBT)
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMLMOTS
C
C- Création du LISTMOTS contenant le nom des composantes du CHAMPOINT :
C-   'O2  ', 'N2  ', 'H2  ', 'H2O ', 'P   ' et 'T   ' décrivant
C-   respectivement le nombre de moles des différentes especes
C-   ainsi que la pression et la température du mélange.
C
      IRET = 0
      JGN  = LOCOMP
      JGM  = 6
      SEGINI MLMOTS
      MOTS(1) = 'H2  '
      MOTS(2) = 'O2  '
      MOTS(3) = 'H2O '
      MOTS(4) = 'N2  '
      MOTS(5) = 'P   '
      MOTS(6) = 'T   '
C
C-----------------------------------------------------------------
C- Evolution des propriétés thermodynamiques des constituants en
C- fonction de la température. On retient les modèles de CHEMKIN :
C-----------------------------------------------------------------
C-  cp/R  = a1 + a2 T + a3 T**2  + a4 T**3 + a5 T**4
C-  H0/RT = a1 + a2 T/2  + a3 T**2 /3 + a4 T**3 /4 + a5 T**4 /5 + a6/T
C-  S0/R  = a1 lnT + a2 T + a3 T**2 /2 + a4 T**3 /3 + a5 T**4 /4 + a7
C
C- On a donc besoin de 7 coefficients pour caractériser la dépendance
C- à la température. De part et d'autre de la température seuil 1000K
C- deux approximations différentes sont proposées -> 14 coefs/espece
C
C   TINT --> Température seuil
C   IINT --> Nombre de coefficients nécessaire pour T donnée (7)
C   AH2  --> Coeffs a1 à a14 pour H2
C   AH2O --> Coeffs a1 à a14 pour H2O
C   AN2  --> Coeffs a1 à a14 pour N2
C   AO2  --> Coeffs a1 à a14 pour O2
C
C Référence : The CHEMKIN Thermodynamic Data Base, R.J.Kee, F.M.Rupley,
C             J.A.Miller, SANDIA Report, SAND87-8215B.
C
      TINT = 1000.D0
      IINT = 7
C
      AH2(1)   =  0.02991423D+02
      AH2(2)   =  0.07000644D-02
      AH2(3)   = -0.05633829D-06
      AH2(4)   = -0.09231578D-10
      AH2(5)   =  0.01582752D-13
      AH2(6)   = -0.08350340D+04
      AH2(7)   = -0.01355110D+02
      AH2(8)   =  0.03298124D+02
      AH2(9)   =  0.08249442D-02
      AH2(10)  = -0.08143015D-05
      AH2(11)  = -0.09475434D-09
      AH2(12)  =  0.04134872D-11
      AH2(13)  = -0.01012521D+05
      AH2(14)  = -0.03294094D+02
C
      AH2O(1)  =  0.02672146D+02
      AH2O(2)  =  0.03056293D-01
      AH2O(3)  = -0.08730260D-05
      AH2O(4)  =  0.01200996D-08
      AH2O(5)  = -0.06391618D-13
      AH2O(6)  = -0.02989921D+06
      AH2O(7)  =  0.06862817D+02
      AH2O(8)  =  0.03386842D+02
      AH2O(9)  =  0.03474982D-01
      AH2O(10) = -0.06354696D-04
      AH2O(11) =  0.06968581D-07
      AH2O(12) = -0.02506588D-10
      AH2O(13) = -0.03020811D+06
      AH2O(14) =  0.02590233D+02
C
      AN2(1)   =  0.02926640D+02
      AN2(2)   =  0.01487977D-01
      AN2(3)   = -0.05684761D-05
      AN2(4)   =  0.01009704D-08
      AN2(5)   = -0.06753351D-13
      AN2(6)   = -0.09227977D+04
      AN2(7)   =  0.05980528D+02
      AN2(8)   =  0.03298677D+02
      AN2(9)   =  0.01408240D-01
      AN2(10)  = -0.03963222D-04
      AN2(11)  =  0.05641515D-07
      AN2(12)  = -0.02444855D-10
      AN2(13)  = -0.01020900D+05
      AN2(14)  =  0.03950372D+02
C
      AO2(1)   =  0.03697578D+02
      AO2(2)   =  0.06135197D-02
      AO2(3)   = -0.01258842D-05
      AO2(4)   =  0.01775281D-09
      AO2(5)   = -0.01136435D-13
      AO2(6)   = -0.01233930D+05
      AO2(7)   =  0.03189166D+02
      AO2(8)   =  0.03212936D+02
      AO2(9)   =  0.01127486D-01
      AO2(10)  = -0.05756150D-05
      AO2(11)  =  0.01313877D-07
      AO2(12)  = -0.08768554D-11
      AO2(13)  = -0.01005249D+05
      AO2(14)  =  0.06034738D+02
C
C- Masse molaire en gramme des différents constituants
C
      AN(1) = 2.01594D0
      AN(2) = 31.9988D0
      AN(3) = 18.01534D0
      AN(4) = 28.01340D0
C
C- Création du CHPOINT contenant l'ensemble des composantes cherchées,
C- le CHAMPOINT étant laissé dans la pile afin que l'opérateur EXTR
C- fasse son travail.
C
      CALL ECROBJ('LISTMOTS',MLMOTS)
      CALL EXCOMP
      IF (IERR.NE.0) RETURN
      CALL LIROBJ('CHPOINT',ICHP1,1,IRETOU)
      IF (IERR.NE.0) RETURN
C
C- Test du CHAMPOIN :
C-  1) le CHPO n'est pas partitionné (i.e. on a toutes les composantes
C-     en tous les points), sinon erreur
C-     110 2 Le %m1:8 de pointeur %i1 n'est pas élémentaire (n&lt;>1)
C-  2) les composantes sont dans l'ordre du LISTMOTS, sinon erreur
C-     735 2 La composante %m1:4 est inexistante ou mal positionnee
C-           dans une connectivite
C
      MCHPOI = ICHP1
      SEGACT MCHPOI
      SEGACT MLMOTS
      N1 = IPCHP(/1)
      IF (N1.NE.1) THEN
         MOTERR(1:8) = 'CHAMPOIN'
         INTERR(1)   = ICHP1
         CALL ERREUR(110)
         RETURN
      ENDIF
      MSOUPO = IPCHP(1)
      SEGACT MSOUPO
      DO 10 I=1,NOCOMP(/2)
         IF (MOTS(I).NE.NOCOMP(I)) THEN
            MOTERR= MOTS(I)
            CALL ERREUR(735)
            RETURN
         ENDIF
 10   CONTINUE
      SEGSUP MCHPOI
      SEGSUP MLMOTS
C
C- Passage à la phase de calcul
C
      CALL DETO1(MSOUPO,MCHPO1,MCHPO2,MCHPO3,
     1           AN,AH2,AH2O,AN2,AO2,TINT,IINT,NBC,NBT)
      SEGSUP MSOUPO
      IF (IERR.NE.0) RETURN
C
C- Ecriture des résultats dans la pile
C
      CALL ECROBJ('CHPOINT',MCHPO3)
      CALL ECROBJ('CHPOINT',MCHPO2)
      CALL ECROBJ('CHPOINT',MCHPO1)
C
      RETURN
      END