C FLA014    SOURCE    TTMF3     12/07/05    21:15:36     7425
      FUNCTION FLA014(T,TC,INDIC,TRAV)
C---------------------------------------------------------------------
C Calcul du débit au travers du PAR pour une température de plaque de
C TC et un controle de mh2 par un déficit en O2/H2 à l'entrée du PAR
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/  INDIC   : entier   : Déficit en O2 (2) ou en H2 (3)
C  /S FLA014  : 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                (dans l'ordre N2, O2, H2 et H2O)
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 Comme il y a un déficit en O2/H2, mh2/mp est connu et égal au min de
C Yh2 et de (2Wh2/Wo2)*Yo2. Par suite mh2/mp est une variable muette.
C
C L'équation de Bernouilli (qdm dans le PAR) se met sous la forme
C F(mp) = a mp^2 + b mp - ro^2 (1 - ( mp /((mp +D)E) )^2 ) = 0
C On recherche la plus grande racine de ce polynome par un Newton.
C On choisit xp = mp * 10^3 pour lisser les gradients.
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 FLA014
      REAL*8 XE(7),OME(7)
      REAL*8 MUE,ME

-INC PPARAM
-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
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)
C
      ibou = cpi(/1)
      sum  = 0.D0
      do 10 i=1,ibou
         sum = sum + xe(i)*m(i)
 10   continue
      ME = sum
C
C Coefficient de la partie linéaire en x
      C = -0.5D0*sp*sp*g*(L/2.D0+Lch)
      A = -Ck/C*1.D-6
      B = -48.D0*mue*sp*L/(Dh*Dh)/C*1.D-3
C
C Coefficient de la partie non-linéaire en x
C Initialisation de DD et EE avec mh2 controlé par le déficit en H2 à
C l'entrée du PAR (indic=3) (resp. par le déficit en O2, indic=2)
      IF (INDIC .EQ. 3) THEN
         DD = 1.D3*S*fla009(pe,Te,Tc,Te,xe,trav)*
     &        (Tc-Te)/(Te*(Cpe-al*ome(3)))
         EE = 1.D0 - Me*ome(3)/(2*M(3))
      ELSE
         DD = 1.D3*S*fla009(pe,Te,Tc,Te,xe,trav)*(Tc-Te)/
     &        (Te*(Cpe-2*al*M(3)*ome(2)/M(2)))
         EE = 1.D0 - Me*ome(2)/M(2)
      ENDIF
C
C Méthode de Newton : x_(k+1) = x_k - f(x_k)/f'(x_k)
C (l'initialisation X=DD correspond d'après l'équation de bilan
C d'énergie du gaz au débit qu'on aurait pour une température
C en sortie du PAR égale à 2*TE)
      RESU = 1.D0
      X    = DD
      I0   = 0
200   CONTINUE
      IF (I0 .GE. MAXNEWT) GOTO 201
      IF ((ABS(RESU)) .LT. (EPS_CON*100.D0)) GOTO 201
         I0   = I0 + 1
         CALL FLA017(x,y,yp,roe,A,B,DD,EE)
         RESU = Y
         X    = X - Y/YP
         GOTO 200
201   CONTINUE
C
C On maintient quoiqu'il arrive un débit résiduel dans le PAR
C de vitesse débitante u
      IF (I0 .GE. MAXNEWT) THEN
         WRITE(IOIMP,*) 'fla014 : non convergence : ',(X*1.D-3),TC
         INTERR(1) = MAXNEWT
         CALL ERREUR(151)
         X = 0.D0
      ENDIF
      BOBO = X * 1.D-3
      IF (U*ROE*SP .GT. BOBO) THEN
         BOBO = U*ROE*SP
      ENDIF
      FLA014 = BOBO
C
      RETURN
      END