Télécharger flhllc.eso

Retour à la liste

Numérotation des lignes :

  1. C FLHLLC SOURCE BECC 10/05/05 21:15:10 6674
  2. SUBROUTINE FLHLLC(NESP,
  3. & GAMG,ROG,PG,UNG,UTG,
  4. & GAMD,ROD,PD,UND,UTD,
  5. & YG,YD,FLU1,
  6. & CELLT)
  7. C
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : FLHLLC
  13. C
  14. C DESCRIPTION : Formulation Volumes Finis pour les Equations
  15. C d'Euler Multi-Especes relatives à un melange
  16. C de gaz ideals.
  17. C
  18. C Calcul du flux aux interfaces avec la methode
  19. C de HLLC (avec gamma = const)
  20. C Voir: * Batten et al. SIAM J. Sci. Comp.
  21. C Vol. 18 pp 1553-6, 1997
  22. C * Batten et al. JCP, 1997
  23. C Vol. 137 pp 38--
  24. C
  25. C N.B.: SPLIT DE FRACTIONS MASSIQUES ET UNT A LA
  26. C LARROUTUROU
  27. C
  28. C LANGUAGE : FORTRAN 77
  29. C
  30. C AUTEUR : A. BECCANTINI DRN/DMT/SEMT/LTMF
  31. C
  32. C************************************************************************
  33. C
  34. C APPELES
  35. C
  36. C************************************************************************
  37. C
  38. C**** Entrées:
  39. C
  40. C NESP = nombre d'especes considérées dans les Equations
  41. C d'Euler
  42. C
  43. C GAMG, GAMD = les "gamma" du gaz (gauche et droite)
  44. C
  45. C ROG, ROD = les densités
  46. C
  47. C PG, PD = les pressions
  48. C
  49. C UNG, UND = vitesses normales
  50. C
  51. C UTG, UTD = vitesses tangentielles
  52. C
  53. C YG, YD = tables des fractiones massiques
  54. C
  55. C**** Sorties:
  56. C
  57. C FLU1 = table du flux a l'interface dans le repaire
  58. C (n,t), i.e.
  59. C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  60. C rho*un*y1, ...)
  61. C
  62. C CELLT = condition de stabilité, i.e.
  63. C
  64. C dT/diamax < cellt (temps/longueur)
  65. C
  66. C************************************************************************
  67. C
  68. C HISTORIQUE (Anomalies et modifications éventuelles)
  69. C
  70. C HISTORIQUE : Créé le 11.09.2000
  71. C
  72. C************************************************************************
  73. C
  74. C N.B.: Toutes les variables sont DECLAREES
  75. C
  76. C It exactly captures:
  77. C
  78. C Stationary shocks (test for instance
  79. C
  80. C GAMG=1.4D0
  81. C GAMD=1.4D0
  82. C ROG=1.0D0
  83. C ROD=4.5584471D0
  84. C PG=1.0D5
  85. C PD=0.18279373D7
  86. C UNG=0.14877919D4
  87. C UND=0.32638130D3
  88. C UTG=0.0D0
  89. C UTD=0.0D0
  90. C ) and stationary contacts
  91. C
  92. IMPLICIT INTEGER(I-N)
  93. INTEGER NESP
  94. & , IESP
  95. REAL*8 GAMG,ROG,PG,UNG,UTG
  96. & ,GAMD,ROD,PD,UND,UTD
  97. & ,YG(*),YD(*),FLU1(*),CELLT
  98. & ,HALF
  99. PARAMETER (HALF=0.5D0)
  100. REAL*8 USGGM1,USGDM1,HTG,HTD,REING,REIND,SQRG,SQRD,USOSQR
  101. & ,HTAVE,UNAVE,REIAVE,USGAVE,GAMAVE,ROAVE,PAVE,ASONAV
  102. & ,ASONG,ASOND,SG,SD,SM,PSTAR,RETG,RETD,OMEGA
  103. & ,ROGS,RUNGS,RETGS,RODS,RUNDS,RETDS
  104. C
  105. USGGM1=1.0D0/(GAMG-1.0D0)
  106. USGDM1=1.0D0/(GAMD-1.0D0)
  107. C
  108. C Evaluation of the Roe average state to compute the non-linear wave
  109. C speeds (According to Shuye, JCP 142, pag. 217--, 1998
  110. C
  111. REING=PG*USGGM1
  112. REIND=PD*USGDM1
  113. HTG=(HALF*UNG*UNG)+(GAMG*REING/ROG)
  114. HTD=(HALF*UND*UND)+(GAMD*REIND/ROD)
  115. SQRG=ROG**0.5D0
  116. SQRD=ROD**0.5D0
  117. USOSQR=1.0D0/(SQRG+SQRD)
  118. HTAVE=((SQRG*HTG)+(SQRD*HTD))*USOSQR
  119. UNAVE=((SQRG*UNG)+(SQRD*UND))*USOSQR
  120. USGAVE=((SQRG*USGGM1)+(SQRD*USGDM1))*USOSQR
  121. GAMAVE=(1.0D0/USGAVE)+1.0D0
  122. REIAVE=((SQRG*REING)+(SQRD*REIND))*USOSQR
  123. C
  124. ROAVE=REIAVE*GAMAVE/(HTAVE-(HALF*UNAVE*UNAVE))
  125. PAVE=REIAVE*(GAMAVE-1.0D0)
  126. C
  127. ASONAV=(GAMAVE*PAVE/ROAVE)**0.5D0
  128. C
  129. C**** Evaluation of the speed of the Genuinely-Non-Linear Waves
  130. C
  131. ASONG=(GAMG*PG/ROG)**0.5D0
  132. ASOND=(GAMD*PD/ROD)**0.5D0
  133. C
  134. SG=UNG-ASONG
  135. SG=MIN(SG,(UNAVE-ASONAV))
  136. SD=UND+ASOND
  137. SD=MAX(SD,(UNAVE+ASONAV))
  138. C
  139. C**** Speed of the contact discontinuity
  140. C
  141. SM=ROD*UND*(SD-UND)-ROG*UNG*(SG-UNG)+PG-PD
  142. SM=SM/((ROD*(SD-UND))-(ROG*(SG-UNG)))
  143. C
  144. C**** Interfacial pressure
  145. C
  146. PSTAR=ROG*(UNG-SG)*(UNG-SM)+PG
  147. RETG=(USGGM1*PG)+(HALF*ROG*((UNG*UNG)+(UTG*UTG)))
  148. RETD=(USGDM1*PD)+(HALF*ROD*((UND*UND)+(UTD*UTD)))
  149. C
  150. IF(SG .GT. 0)THEN
  151. FLU1(1)=ROG*UNG
  152. FLU1(2)=ROG*UNG*UNG+PG
  153. FLU1(4)=UNG*(RETG+PG)
  154. ELSEIF(SM .GT. 0)THEN
  155. OMEGA=1.0D0/(SG-SM)
  156. ROGS=ROG*(SG-UNG)
  157. ROGS=ROGS*OMEGA
  158. RUNGS=((SG-UNG)*(ROG*UNG))+(PSTAR - PG)
  159. RUNGS=RUNGS*OMEGA
  160. RETGS=((SG-UNG)*RETG)+(PSTAR*SM)-(PG*UNG)
  161. RETGS=RETGS*OMEGA
  162. C
  163. FLU1(1)=ROGS*SM
  164. FLU1(2)=RUNGS*SM+PSTAR
  165. FLU1(4)=(RETGS+PSTAR)*SM
  166. C
  167. ELSEIF(SD .GT. 0)THEN
  168. OMEGA=1.0D0/(SD-SM)
  169. RODS=ROD*(SD-UND)
  170. RODS=RODS*OMEGA
  171. RUNDS=((SD-UND)*(ROD*UND))+(PSTAR - PD)
  172. RUNDS=RUNDS*OMEGA
  173. RETDS=((SD-UND)*RETD)+(PSTAR*SM)-(PD*UND)
  174. RETDS=RETDS*OMEGA
  175. C
  176. FLU1(1)=RODS*SM
  177. FLU1(2)=RUNDS*SM+PSTAR
  178. FLU1(4)=(RETDS+PSTAR)*SM
  179. ELSE
  180. FLU1(1)=ROD*UND
  181. FLU1(2)=ROD*UND*UND+PD
  182. FLU1(4)=UND*(RETD+PD)
  183. ENDIF
  184. C
  185. FLU1(3)=(FLU1(1)*(UTG+UTD))+(ABS(FLU1(1))*(UTG-UTD))
  186. FLU1(3)=HALF*FLU1(3)
  187. C
  188. DO IESP=1,NESP,1
  189. FLU1(4+IESP)=(FLU1(1)*(YG(IESP)+YD(IESP)))+
  190. & (ABS(FLU1(1))*(YG(IESP)-YD(IESP)))
  191. FLU1(4+IESP)=HALF*FLU1(4+IESP)
  192. ENDDO
  193. C
  194. CELLT=1.0D0/(MAX(SG,SD))
  195. C
  196. RETURN
  197. END
  198.  
  199.  
  200.  
  201.  
  202.  
  203.  
  204.  
  205.  

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