C FLA015    SOURCE    TTMF3     12/07/05    21:15:37     7425           
      FUNCTION FLA015(T,TC,MH2,TRAV)
C---------------------------------------------------------------------
C Calcul du débit au travers du PAR pour une température de plaque de 
C Tc et une consommation d'hydrogène mh2=mh2max
C---------------------------------------------------------------------
C
C---------------------------
C Parametres Entree/Sortie :
C---------------------------
C
C E/  T       : flottant : Temps courant (s)
C E/  TC      : flottant : Température des plaques (K)
C E/  MH2     : flottant : Hydrogène consommé (kg/s)
C  /S FLA015  : flottant : Débit total traversant le PAR (kg/s) 
C
C------------------------------
C Variables de TRAV utilisées :
C------------------------------
C
C E/  PRESSION : Pression à l'entrée du PAR (Pa)
C E/  TEMPENT  : Température à l'entrée du PAR (K)
C E/  XiMOY    : Fraction molaire à l'entrée du PAR
C E/  M(7)     : Masse molaire des constituants du mélange (kg/mol) 
C E/  U        : Vitesse minimale dans le PAR (= 0.01 m/s)
C E/  S        : Surface des plaques (m2) 
C E/  SP       : Section de passage entrée/sortie fluide du PAR (m2)
C E/  L        : Hauteur des plaques (m)
C E/  LCH      : Hauteur de la cheminée (m)
C E/  DH       : Diamètre hydraulique entre les plaques (m) 
C E/  CK       : Demi somme des coefficients de perte de charge à
C                l'entrée et à la sortie du PAR
C E/  G        : Gravité (=9.81 m/s2)
C E/  AL       : Cste fonction de CP(i)
C E/  EPS_CON  : Seuil de convergence pour Newton (|f(x)|<EPS_CON)
C
C--------------------
C Algorithme utilisé
C--------------------
C 
C L'équation de Bernouilli (qdm dans le PAR) se met sous la forme
C F(mp,ros) = a mp^2 + b mp + c (ro^2 - ros^2) = 0
C ros(mp)   = ro / (1 + D/(mp-E))/ (1 - F/mp)
C On recherche la plus grande racine de ce polynome par une méthode 
C de dichotomie (c'est pourquoi on privilégie le coté Xmax).
C
C---------------------------------------------------------------------
C
C Langage : ESOPE + FORTRAN 77
C
C Mise en oeuvre : H. Paillère (1997, TTMF)
C
C---------------------------------------------------------------------
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 FLA015
      REAL*8 XE(7),OME(7)
      REAL*8 MUE,ME,MH2,MP
-INC CCOPTIO
      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
      MAXNEWT = 10000
      COFNEWT = 1.D0 
C
      PE    = PRESSION
      TE    = TEMPENT
      XE(1) = XN2MOY
      XE(2) = XO2MOY
      XE(3) = XH2MOY
      XE(4) = XH2OMOY
      XE(5) = XHEMOY
      XE(6) = XCO2MOY
      XE(7) = XCOMOY
      CALL FLA012(XE,OME,TRAV)
      ROE   = FLA001(PE,TE,XE,TRAV)
      MUE   = FLA004(TE,XE,TRAV)
      CPE   = FLA003(OME,TRAV)
      ibou = cpi(/1)
      sum  = 0.D0
      do 5 i=1,ibou
         sum = sum + xe(i)*m(i)
 5    continue
      ME = sum
C
      ROEXT = ROE
C
C Coefficient de l'équation de Bernouilli et de la densité en sortie
      A  = Ck*COFNEWT*COFNEWT
      B  = COFNEWT*48.D0*mue*sp*L/(Dh*Dh)
      C  = -0.5D0*sp*sp*g*(L/2.D0+Lch)
      DD = S*fla009(pe,Te,Tc,Te,xe,trav)*(Tc-Te)/(Te*Cpe)*COFNEWT
      EE = al*mh2/Cpe*COFNEWT
      FF = mh2*Me/(2.D0*M(3))*COFNEWT
C
C Méthode de Dichotomie : la plus grande des racines est au-delà 
C des poles du denominateur (EE-DD) et FF
      I0   = 0
      XMIN = MAX(EE-DD,FF)
      IF (XMIN .LE. 0.D0) THEN
         XMIN = 0.D0 
      ENDIF
      XMAX  = MAX(1.D0,U*ROE*SP*500.D0)/COFNEWT
      RESUP = FLA016(XMAX,roext,A,B,C,DD,EE,FF)
      RESU  = RESUP
      RESUM = FLA016(XMIN,roext,A,B,C,DD,EE,FF)
      EPS_2 = MIN(ABS(RESUM),EPS_CON)
C
C Recherche de la racine sur l'intervalle (XMIN,XMAX)
 10   CONTINUE
      X = (XMAX + XMIN) / 2.D0
      IF (I0 .GE. MAXNEWT) GOTO 20
      IF ((XMAX-XMIN) .LT. EPS_CON) GOTO 30
C      IF (ABS(RESU) .LT. (EPS_CON*EPS_2)) GOTO 40
CC      IF (ABS(RESU) .LT. (RESUP*EPS_CON*EPS_2)) GOTO 40
         I0   = I0 + 1
         RESU = FLA016(x,roext,A,B,C,DD,EE,FF)
         IF ((RESU*RESUP) .LE. 0.D0) THEN
            XMIN  = X
         ELSE
            XMAX = X
         ENDIF
         GOTO 10
 20   CONTINUE
C
C I0 .GE. MAXNEWT : non convergence du Newton
C Se produit à l'usage lorsque la racine et les poles sont quasi nuls.
C Par suite, le débit est mis à zero.
      INTERR(1) = MAXNEWT
      CALL ERREUR(151)
      X = 0.D0
      GOTO 40
 30   CONTINUE
C 
C (XMAX-XMIN) .LT. EPS_CON : les 2 bornes de dichotomie sont égales.
C Si la fonction est de même signe en XMAX et XMIN, il n'y a pas de
C racine réelle au dela des poles. Par suite, le débit est nul.
      RESUP = FLA016(XMAX,roext,A,B,C,DD,EE,FF)
      RESUM = FLA016(XMIN,roext,A,B,C,DD,EE,FF)
      IF ((RESUM*RESUP) .GT. 0.D0) THEN
         X = 0.D0
      ENDIF
 40   CONTINUE
C
C On maintient quoiqu'il arrive un débit résiduel dans le PAR 
C de vitesse débitante u
      MP = X * COFNEWT
      IF (U*ROE*SP .GT. MP) THEN
          MP = U*ROE*SP
      ENDIF
      FLA015 = MP
C
      RETURN
      END