Sources Esope | Cast3M

flamr
C FLAMR     SOURCE    TTMF3     12/07/05    21:15:46     7425
      SUBROUTINE FLAMR(YH2e,YO2e,YN2e,YH2Oe,YHEe,YCO2e,YCOe,
     &              PRESe,TEMPe,
     &              TEMPSOR,TEMPPLA,Q,
     &              YH2SOR,YO2SOR,YN2SOR,YH2OSOR,YHESOR,YCO2SOR,YCOSOR,
     &              XH2ONF,EFF,DDT,iKA,iMOD,iCST,IHI)
C------------------------------------------------------------------------
C Modélisation d'un recombineur catalytique (PAR) par une approche 0D.
*------------------------------------------------------------------------
C
C---------------------------
C Parametres Entree/Sortie :
C---------------------------
C
C E/  YH2e    : flottant : Fraction massique moyenne de H2 (PAR inlet)
C E/  YO2e    : flottant : Fraction massique moyenne de O2 (PAR inlet)
C E/  YN2e    : flottant : Fraction massique moyenne de N2 (PAR inlet)
C E/  YH2Oe   : flottant : Fraction massique moyenne de H2O (PAR inlet)
C E/  YHEe    : flottant : Fraction massique moyenne de HE (PAR inlet)
C E/  YCO2e   : flottant : Fraction massique moyenne de CO2 (PAR inlet)
C E/  YCOe    : flottant : Fraction massique moyenne de CO (PAR inlet)
C E/  PRESe   : flottant : Pression moyenne autour du PAR (Pa)
C E/  TEMPe   : flottant : Température moyenne du gaz (PAR inlet, K)
C  /S TEMPSOR : flottant : Température moyenne du gaz (PAR outlet, K)
C E/S TEMPPLA : flottant : Température moyenne des plaques du PAR (K)
C  /S Q       : flottant : Débit dans le PAR (kg/s)
C  /S YH2SOR  : flottant : Fraction massique moyenne de H2 (PAR outlet)
C  /S YO2SOR  : flottant : Fraction massique moyenne de O2 (PAR outlet)
C  /S YN2SOR  : flottant : Fraction massique moyenne de N2 (PAR outlet)
C  /S YH2OSOR : flottant : Fraction massique moyenne de H2O (PAR outlet)
C  /S YHESOR  : flottant : Fraction massique moyenne de HE (PAR outlet)
C  /S YCO2SOR : flottant : Fraction massique moyenne de CO2 (PAR outlet)
C  /S YCOSOR  : flottant : Fraction massique moyenne de CO (PAR outlet)
C  /S XH2ONF  : flottant : Flag pour l'allumage (négatif si le PAR est ON)
C  /S EFF     : flottant : Efficacité du recombineur
C E/  DDT     : flottant : Pas de temps (s)
C E/  IKA     : entier   : Flag pour correction du coefficient d'échange
C                          comme dans RALOC (voir fla009.eso)
C E/  IMOD    : entier   : Modélisation du PAR (1=SIEMENS, 2=HEATER)
C E/  ICST    : entier   : Pointeur de LISTREEL, paramètres du modèle
C  /S IHI     : entier   : Flag de convergence et gestion des erreurs
C
C------------------------------------------------------------------------
C
C Langage : ESOPE + FORTRAN 77
C
C Modèle Initial : P. Berne SERAC/LESI/96-11
C
C----------------
C Modifications :
C----------------
C
C 1997 : H. Paillère - mise en oeuvre dans CASTEM
C 2003 : L. Dada     - correction RALOC sur le coefficient d'échange
C 2005 : F. Dabbene  - modification de la loi SIEMENS
C 2010 : F. Dabbene  - differenciation allumage/extinction
C                    - arret en cas de non-convergence des Newtons
C 2012 : F. Dabbene  - ajout du modèle HEATER
C
C------------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC SMLREEL
C Ordre des espèces : 1 N2, 2 O2, 3 H2, 4 H2O, 5 HE, 6 CO2, 7 CO
      real*8 xi(7),om(7)
C
C Description du segment de travail TRAV (initialisé ici, il est
C transmis à tous les .eso du modèle (de fla001.eso à fla019.eso).
C Les valeurs numériques sont celles du recombineur SIEMENS FR/1-150
C   Hauteur totale du recombineur = 1 m
C   Longueur du recombineur (Lgr) = 0.15 m
C   Profondeur du recombineur (W) = 0.15 m
C   Nombre de plaques (n)         = 14
C
C IKALP   : Flag pour correction RALOC du coefficient d'échange (=IKA)
C IMODEL  : Modélisation du PAR (1=SIEMENS, 2=AUTRE) (=IMOD)
C E       : Distance entre deux plaques (= 0.01 m)
C           (Lgr/(n+1), W=largeur et n=nombre de plaques)
C L       : Hauteur des plaques (= 0.15 m)
C LCH     : Hauteur de la cheminée (= 0.85 m)
C DH      : Diamètre hydraulique entre les plaques (= 0.01875 m)
C           (2EW/(E+W))
C S       : Surface des plaques (= 0.675 m2)
C           (2n*LW, n=nombre de plaques)
C SP      : Section de passage entrée/sortie fluide du PAR (= 0.0225 m2)
C           (Lgr*W)
C CK      : Demi somme des coefficients de perte de charge à l'entrée et à
C           la sortie du PAR (= 3/2)
C MC      : Masse des plaques (= 0.133 kg)
C CPC     : Chaleur spécifique des plaques (= 132.6 J/kg/K)
C G       : Gravité (=9.81 m/s2)
C R       : Constante des gaz (= 8.314 J/mol/K)
C DELTAH  : Energie libérée par la combustion d'un kg d'H2 (= 120 MJ/kg)
C M(7)    : Masse molaire des constituants du mélange (kg/mol)
C ncst    : Nombre de constantes du modèle (= longueur du listreel ICST)
C CSTMOD  : Constantes du modèle  (= contenu de ICST)
C CPI(7)  : Chaleurs massiques des constituants (J/kg/K)
C AL      : Cste fonction de CP(i)
C EPS_MH2 : Seuil en deça duquel on suppose que la consommation d'H2 max
C           est mise à 0. (kg/s)
C EPS_DT  : Ecart entre la température des plaques et la température du
C           gaz (Tc-T) en deça duquel on suppose qu'il n'y a pas
C           d'échange convectif (K)
C EPS_CON : Critère de convergence
C U       : Vitesse minimale dans le PAR (= 0.01 m/s)
C
      segment trav
           integer iKALP,iMODEL
           real*8 e,L,Lch,Dh,S,sp,Ck
           real*8 mc,Cpc
           real*8 g,R,deltah
           real*8 M(nbesp),cstmod(ncst)
           real*8 Cpi(nbesp),al
           real*8 eps_mh2,eps_dt,eps_con,u
           real*8 XH2MOY,XO2MOY,XN2MOY,XH2OMOY,PRESSION,TEMPENT
           real*8 XHEMOY,XCO2MOY,XCOMOY
      endsegment
C
C      write(6,*) '************************ dans flamr : debut'
      IHI = 0
      MLREEL = ICST
      SEGACT MLREEL*MOD
      NCST = PROG(/1)
      NBESP = 7
      SEGINI TRAV
C
C Paramètres dépendant de la géométrie
C (à changer si le PAR n'est pas un FR/1-150)
      IMODEL = IMOD
      e   = 0.01D0
      L   = 0.15D0
      Lch = 0.85D0
      Dh  = 0.01875D0
      S   = 0.675D0
      sp  = 0.0225D0
      Ck  = 1.5D0
      mc  = 0.133D0
C
C Paramètres physico-chimiques
      Cpc    = 132.6D0
      g      = 9.81D0
      R      = 8.3144621D0
      deltah = 1.2D8
      M(1) = 2.D0 * 0.0140067D0
      M(2) = 2.D0 * 0.0159994D0
      M(3) = 2.D0 * 0.00100794D0
      M(4) = 2.D0 * 0.00100794D0 + 0.0159994D0
      M(5) = 0.004002602D0
      M(6) = 2.D0 * 0.0159994D0 + 0.0120107D0
      M(7) = 0.0120107D0 + 0.0159994D0
C
      Cpi(1) = 1046.D0
      Cpi(2) = 942.D0
      Cpi(3) = 14490.D0
      Cpi(4) = 2014.D0
      Cpi(5) = 5192.D0
      Cpi(6) = 1117.D0
      Cpi(7) = 1042.D0
      al=0.5D0*(Cpi(3)+M(2)*Cpi(2)/(2*M(3))-M(4)*Cpi(4)/M(3))
C
C Transformation des fractions massiques en fractions molaires
      OM(1) = YN2e
      OM(2) = YO2e
      OM(3) = YH2e
      OM(4) = YH2Oe
      OM(5) = YHEe
      OM(6) = YCO2e
      OM(7) = YCOe
C      write(6,*) '************************ dans flamr : fla013'
c      write(6,*) 'om(1) ',om(1)
c      write(6,*) 'om(2) ',om(2)
c      write(6,*) 'om(3) ',om(3)
c      write(6,*) 'om(4) ',om(4)
c      write(6,*) 'om(5) ',om(5)
c      write(6,*) 'om(6) ',om(6)
c      write(6,*) 'om(7) ',om(7)
      CALL FLA013(xi,om,trav)
c      write(6,*) 'xi(1) ',xi(1)
c      write(6,*) 'xi(2) ',xi(2)
c      write(6,*) 'xi(3) ',xi(3)
c      write(6,*) 'xi(4) ',xi(4)
c      write(6,*) 'xi(5) ',xi(5)
c      write(6,*) 'xi(6) ',xi(6)
c      write(6,*) 'xi(7) ',xi(7)
C
C Données d'entrée dans TRAV
      XN2MOY   = xi(1)
      XO2MOY   = xi(2)
      XH2MOY   = xi(3)
      XH2OMOY  = xi(4)
      XHEMOY   = xi(5)
      XCO2MOY  = xi(6)
      XCOMOY   = xi(7)
      PRESSION = PRESe
      TEMPENT  = TEMPe
      iKALP    = iKA
C
C Gestion du FR/1-150
C
C
C Efficacite du recombineur
      IF (XH2MOY .LE. XO2MOY) THEN
         EFF = 1.D0
      ELSE
         EFF = 0.6D0
      ENDIF
C
C Allumage du recombineur
      XH2ONF = PROG(3)
cc      write(6,*) 'XH2ON XH2MOY',xh2onf,xh2moy
      IF (XH2ONF .LE. 0.D0) THEN
         CHI = 1.D0
      ELSE
         IF (XH2MOY .GE. XH2ONF) THEN
            CHI    = 1.D0
            XH2ONF = -XH2ONF
         ELSE
            CHI = 0.D0
         ENDIF
      ENDIF
C
C Extinction du recombineur
      XH2OFF = PROG(4)
      IF (XH2ONF .LE. 0.D0) THEN
      IF (XH2MOY .LE. XH2OFF) THEN
         CHI    = 0.D0
         XH2ONF = -XH2ONF
      ENDIF
      ENDIF
cc      write(6,*) 'XH2ON chi   ',xh2onf,chi
      PROG(1) = PROG(1) * EFF * CHI
      PROG(2) = PROG(2) * EFF * CHI
C
      DO 10 I1=1,NCST
         CSTMOD(I1) = PROG(I1)
 10   CONTINUE
      SEGDES MLREEL
C
C Critère de convergence et vitesse minimale dans le PAR
      eps_mh2 = 1.D-10
      eps_dt  = 1.D-3
      eps_con = 1.D-8
      u       = 0.01D0
C
C Limitation du pas de temps à 0.5s
C (Découpage du pas de temps en sous pas de temps si trop grand)
      t = 0.0D0
      I = 1
      TMAX = DDT
      Tc   = TEMPPLA
      dt = 0.5D0
      IF (ddt .LT. dt) THEN
         dt = ddt
      ENDIF
 20   CONTINUE
      IF (.NOT.(t.lt.tmax)) GOTO 30
C
C Calcul de la temperature des plaques au temps t+dt
         IF ((t+dt) .GT. tmax) THEN
            dt = tmax - t
         ENDIF
         call fla011(Tc1,t,dt,Tc,trav)
         IF (TC1 .LT. 0.D0) THEN
            IHI = -100
            GOTO 30
         ENDIF
C
C Calcul du débit et des conditions de sortie
C au temps t+dt pour la nouvelle temperature de plaques Tc
C (si pb de débit négatif on sort en erreur)
         t  = t + dt
         Tc = Tc1
         call fla019(t,Tc,deb,Tgaz,xi,ihi,trav)
         IF ((DEB .LT. 0.D0) .OR. (Tgaz .LT. 0.D0)) THEN
            IHI = -1 * IHI
            GOTO 30
         ENDIF
         GOTO 20
 30   CONTINUE
C
C Sortie lorsque le temps final est atteint
c      write(6,*) '************************ dans flamr : fla012'
c      write(6,*) 'xi(1) ',xi(1)
c      write(6,*) 'xi(2) ',xi(2)
c      write(6,*) 'xi(3) ',xi(3)
c      write(6,*) 'xi(4) ',xi(4)
c      write(6,*) 'xi(5) ',xi(5)
c      write(6,*) 'xi(6) ',xi(6)
c      write(6,*) 'xi(7) ',xi(7)
      CALL FLA012(xi,om,trav)
c      write(6,*) 'om(1) ',om(1)
c      write(6,*) 'om(2) ',om(2)
c      write(6,*) 'om(3) ',om(3)
c      write(6,*) 'om(4) ',om(4)
c      write(6,*) 'om(5) ',om(5)
c      write(6,*) 'om(6) ',om(6)
c      write(6,*) 'om(7) ',om(7)
      YN2SOR  = om(1)
      YO2SOR  = om(2)
      YH2SOR  = om(3)
      YH2OSOR = om(4)
      YHESOR  = om(5)
      YCO2SOR = om(6)
      YCOSOR  = om(7)
      TEMPPLA = Tc
      Q       = deb
      TEMPSOR = Tgaz
C
      SEGSUP trav
      RETURN
      END