Télécharger fla014.eso

Retour à la liste

Numérotation des lignes :

  1. C FLA014 SOURCE TTMF3 12/07/05 21:15:36 7425
  2. FUNCTION FLA014(T,TC,INDIC,TRAV)
  3. C---------------------------------------------------------------------
  4. C Calcul du débit au travers du PAR pour une température de plaque de
  5. C TC et un controle de mh2 par un déficit en O2/H2 à l'entrée du PAR
  6. C---------------------------------------------------------------------
  7. C
  8. C---------------------------
  9. C Parametres Entree/Sortie :
  10. C---------------------------
  11. C
  12. C E/ T : flottant : Temps courant (s)
  13. C E/ TC : flottant : Température des plaques (K)
  14. C E/ INDIC : entier : Déficit en O2 (2) ou en H2 (3)
  15. C /S FLA014 : flottant : Débit total traversant le PAR (kg/s)
  16. C
  17. C------------------------------
  18. C Variables de TRAV utilisées :
  19. C------------------------------
  20. C
  21. C E/ PRESSION : Pression à l'entrée du PAR (Pa)
  22. C E/ TEMPENT : Température à l'entrée du PAR (K)
  23. C E/ XiMOY : Fraction molaire à l'entrée du PAR
  24. C E/ M(7) : Masse molaire des constituants du mélange (kg/mol)
  25. C (dans l'ordre N2, O2, H2 et H2O)
  26. C E/ U : Vitesse minimale dans le PAR (= 0.01 m/s)
  27. C E/ S : Surface des plaques (m2)
  28. C E/ SP : Section de passage entrée/sortie fluide du PAR (m2)
  29. C E/ L : Hauteur des plaques (m)
  30. C E/ LCH : Hauteur de la cheminée (m)
  31. C E/ DH : Diamètre hydraulique entre les plaques (m)
  32. C E/ CK : Demi somme des coefficients de perte de charge à
  33. C l'entrée et à la sortie du PAR
  34. C E/ G : Gravité (=9.81 m/s2)
  35. C E/ AL : Cste fonction de CP(i)
  36. C E/ EPS_CON : Seuil de convergence pour Newton (|f(x)|<EPS_CON)
  37. C
  38. C--------------------
  39. C Algorithme utilisé
  40. C--------------------
  41. C
  42. C Comme il y a un déficit en O2/H2, mh2/mp est connu et égal au min de
  43. C Yh2 et de (2Wh2/Wo2)*Yo2. Par suite mh2/mp est une variable muette.
  44. C
  45. C L'équation de Bernouilli (qdm dans le PAR) se met sous la forme
  46. C F(mp) = a mp^2 + b mp - ro^2 (1 - ( mp /((mp +D)E) )^2 ) = 0
  47. C On recherche la plus grande racine de ce polynome par un Newton.
  48. C On choisit xp = mp * 10^3 pour lisser les gradients.
  49. C
  50. C---------------------------------------------------------------------
  51. C
  52. C Langage : ESOPE + FORTRAN 77
  53. C
  54. C Mise en oeuvre : H. Paillère (1997, TTMF)
  55. C
  56. C---------------------------------------------------------------------
  57. IMPLICIT INTEGER(I-N)
  58. IMPLICIT REAL*8 (A-H,O-Z)
  59. REAL*8 FLA014
  60. REAL*8 XE(7),OME(7)
  61. REAL*8 MUE,ME
  62. -INC CCOPTIO
  63. segment trav
  64. integer iKALP,iMODEL
  65. real*8 e,L,Lch,Dh,S,sp,Ck
  66. real*8 mc,Cpc
  67. real*8 g,R,deltah
  68. real*8 M(nbesp),cstmod(ncst)
  69. real*8 Cpi(nbesp),al
  70. real*8 eps_mh2,eps_dt,eps_con,u
  71. real*8 XH2MOY,XO2MOY,XN2MOY,XH2OMOY,PRESSION,TEMPENT
  72. real*8 XHEMOY,XCO2MOY,XCOMOY
  73. endsegment
  74. C
  75. MAXNEWT = 10000
  76. C
  77. PE = PRESSION
  78. TE = TEMPENT
  79. XE(1) = XN2MOY
  80. XE(2) = XO2MOY
  81. XE(3) = XH2MOY
  82. XE(4) = XH2OMOY
  83. XE(5) = XHEMOY
  84. XE(6) = XCO2MOY
  85. XE(7) = XCOMOY
  86. CALL FLA012(XE,OME,TRAV)
  87. ROE = FLA001(PE,TE,XE,TRAV)
  88. MUE = FLA004(TE,XE,TRAV)
  89. CPE = FLA003(OME,TRAV)
  90. C
  91. ibou = cpi(/1)
  92. sum = 0.D0
  93. do 10 i=1,ibou
  94. sum = sum + xe(i)*m(i)
  95. 10 continue
  96. ME = sum
  97. C
  98. C Coefficient de la partie linéaire en x
  99. C = -0.5D0*sp*sp*g*(L/2.D0+Lch)
  100. A = -Ck/C*1.D-6
  101. B = -48.D0*mue*sp*L/(Dh*Dh)/C*1.D-3
  102. C
  103. C Coefficient de la partie non-linéaire en x
  104. C Initialisation de DD et EE avec mh2 controlé par le déficit en H2 à
  105. C l'entrée du PAR (indic=3) (resp. par le déficit en O2, indic=2)
  106. IF (INDIC .EQ. 3) THEN
  107. DD = 1.D3*S*fla009(pe,Te,Tc,Te,xe,trav)*
  108. & (Tc-Te)/(Te*(Cpe-al*ome(3)))
  109. EE = 1.D0 - Me*ome(3)/(2*M(3))
  110. ELSE
  111. DD = 1.D3*S*fla009(pe,Te,Tc,Te,xe,trav)*(Tc-Te)/
  112. & (Te*(Cpe-2*al*M(3)*ome(2)/M(2)))
  113. EE = 1.D0 - Me*ome(2)/M(2)
  114. ENDIF
  115. C
  116. C Méthode de Newton : x_(k+1) = x_k - f(x_k)/f'(x_k)
  117. C (l'initialisation X=DD correspond d'après l'équation de bilan
  118. C d'énergie du gaz au débit qu'on aurait pour une température
  119. C en sortie du PAR égale à 2*TE)
  120. RESU = 1.D0
  121. X = DD
  122. I0 = 0
  123. 200 CONTINUE
  124. IF (I0 .GE. MAXNEWT) GOTO 201
  125. IF ((ABS(RESU)) .LT. (EPS_CON*100.D0)) GOTO 201
  126. I0 = I0 + 1
  127. CALL FLA017(x,y,yp,roe,A,B,DD,EE)
  128. RESU = Y
  129. X = X - Y/YP
  130. GOTO 200
  131. 201 CONTINUE
  132. C
  133. C On maintient quoiqu'il arrive un débit résiduel dans le PAR
  134. C de vitesse débitante u
  135. IF (I0 .GE. MAXNEWT) THEN
  136. WRITE(IOIMP,*) 'fla014 : non convergence : ',(X*1.D-3),TC
  137. INTERR(1) = MAXNEWT
  138. CALL ERREUR(151)
  139. X = 0.D0
  140. ENDIF
  141. BOBO = X * 1.D-3
  142. IF (U*ROE*SP .GT. BOBO) THEN
  143. BOBO = U*ROE*SP
  144. ENDIF
  145. FLA014 = BOBO
  146. C
  147. RETURN
  148. END
  149.  
  150.  
  151.  
  152.  
  153.  
  154.  
  155.  
  156.  
  157.  
  158.  
  159.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales