Télécharger fhllcl.eso

Retour à la liste

Numérotation des lignes :

  1. C FHLLCL SOURCE BECC 10/05/05 21:15:09 6674
  2. SUBROUTINE FHLLCL(NESP,
  3. & GAMG,ROG,PG,UNG,UTG,
  4. & GAMD,ROD,PD,UND,UTD,
  5. & YG,YD,VINF0,FLU1,
  6. & CELLT)
  7. C
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : FHLLCL
  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 +Low Mach preconditioning
  21. C
  22. C Voir: * Batten et al. SIAM J. Sci. Comp.
  23. C Vol. 18 pp 1553-6, 1997
  24. C * Batten et al. JCP, 1997
  25. C Vol. 137 pp 38--
  26. C * Luo et al., AIIA Journal, 2005
  27. C Vol. 43, 6, pp 1160--
  28. C
  29. C N.B.: SPLIT DE FRACTIONS MASSIQUES ET UNT A LA
  30. C LARROUTUROU
  31. C
  32. C LANGUAGE : FORTRAN 77
  33. C
  34. C AUTEUR : A. BECCANTINI DRN/DMT/SEMT/LTMF
  35. C
  36. C************************************************************************
  37. C
  38. C APPELES
  39. C
  40. C************************************************************************
  41. C
  42. C**** Entrées:
  43. C
  44. C NESP = nombre d'especes considérées dans les Equations
  45. C d'Euler
  46. C
  47. C GAMG, GAMD = les "gamma" du gaz (gauche et droite)
  48. C
  49. C ROG, ROD = les densités
  50. C
  51. C PG, PD = les pressions
  52. C
  53. C UNG, UND = vitesses normales
  54. C
  55. C UTG, UTD = vitesses tangentielles
  56. C
  57. C YG, YD = tables des fractiones massiques
  58. C
  59. C**** Sorties:
  60. C
  61. C FLU1 = table du flux a l'interface dans le repaire
  62. C (n,t), i.e.
  63. C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  64. C rho*un*y1, ...)
  65. C
  66. C CELLT = condition de stabilité, i.e.
  67. C
  68. C dT/diamax < cellt (temps/longueur)
  69. C
  70. C************************************************************************
  71. C
  72. C HISTORIQUE (Anomalies et modifications éventuelles)
  73. C
  74. C HISTORIQUE : Créé le 11.09.2000
  75. C
  76. C************************************************************************
  77. C
  78. C N.B.: Toutes les variables sont DECLAREES
  79. C
  80. C It exactly captures:
  81. C
  82. C Stationary shocks (test for instance
  83. C
  84. C GAMG=1.4D0
  85. C GAMD=1.4D0
  86. C ROG=1.0D0
  87. C ROD=4.5584471D0
  88. C PG=1.0D5
  89. C PD=0.18279373D7
  90. C UNG=0.14877919D4
  91. C UND=0.32638130D3
  92. C UTG=0.0D0
  93. C UTD=0.0D0
  94. C ) and stationary contacts
  95. C
  96. IMPLICIT INTEGER(I-N)
  97. INTEGER NESP
  98. & , IESP
  99. REAL*8 GAMG,ROG,PG,UNG,UTG
  100. & ,GAMD,ROD,PD,UND,UTD
  101. & ,YG(*),YD(*),FLU1(*),CELLT,VINF0
  102. & ,HALF
  103. PARAMETER (HALF=0.5D0)
  104. REAL*8 USGGM1,USGDM1,HTG,HTD,REING,REIND,SQRG,SQRD,USOSQR
  105. & ,HTAVE,UNAVE,REIAVE,USGAVE,GAMAVE,ROAVE,PAVE,ASONAV
  106. & ,ASONG,ASOND
  107. & ,VINF,PHIG,PHID,PHIAV,XU,SG,SGS,SD,SDS
  108. & ,SM,PSTAR,RETG,RETD,OMEGA
  109. & ,ROGS,RUNGS,RETGS,RODS,RUNDS,RETDS
  110. C
  111. USGGM1=1.0D0/(GAMG-1.0D0)
  112. USGDM1=1.0D0/(GAMD-1.0D0)
  113. C
  114. C Evaluation of the Roe average state to compute the non-linear wave
  115. C speeds (According to Shuye, JCP 142, pag. 217--, 1998
  116. C
  117. REING=PG*USGGM1
  118. REIND=PD*USGDM1
  119. HTG=(HALF*UNG*UNG)+(GAMG*REING/ROG)
  120. HTD=(HALF*UND*UND)+(GAMD*REIND/ROD)
  121. SQRG=ROG**0.5D0
  122. SQRD=ROD**0.5D0
  123. USOSQR=1.0D0/(SQRG+SQRD)
  124. HTAVE=((SQRG*HTG)+(SQRD*HTD))*USOSQR
  125. UNAVE=((SQRG*UNG)+(SQRD*UND))*USOSQR
  126. USGAVE=((SQRG*USGGM1)+(SQRD*USGDM1))*USOSQR
  127. GAMAVE=(1.0D0/USGAVE)+1.0D0
  128. REIAVE=((SQRG*REING)+(SQRD*REIND))*USOSQR
  129. C
  130. ROAVE=REIAVE*GAMAVE/(HTAVE-(HALF*UNAVE*UNAVE))
  131. PAVE=REIAVE*(GAMAVE-1.0D0)
  132. C
  133. ASONAV=(GAMAVE*PAVE/ROAVE)**0.5D0
  134. C
  135. C**** Evaluation of the speed of the Genuinely-Non-Linear Waves
  136. C for the Preconditioned System
  137. C
  138. ASONG=(GAMG*PG/ROG)**0.5D0
  139. ASOND=(GAMD*PD/ROD)**0.5D0
  140. C
  141. C**** We include in the cut-off UNG,UND,UTG,UTD in order to
  142. C avoid low diffusivity in stagnation regions
  143. C
  144. VINF=MAX(VINF0,((UNG*UNG+UTG*UTG)**0.5D0))
  145. VINF=MAX(VINF0,((UND*UND+UTD*UTD)**0.5D0))
  146. C
  147. C*** Calcul des nombres de Mach de référence
  148. C PHIG, PHID, PHIAV
  149. C
  150. IF(VINF.GE.ASONG)THEN
  151. PHIG = 1.D0
  152. ELSE
  153. PHIG = VINF/ASONG
  154. ENDIF
  155. IF(VINF.GE.ASOND)THEN
  156. PHID = 1.D0
  157. ELSE
  158. PHID = VINF/ASOND
  159. ENDIF
  160. IF(VINF.GE.ASONAV)THEN
  161. PHIAV = 1.D0
  162. ELSE
  163. PHIAV = VINF/ASONAV
  164. ENDIF
  165. C
  166. C*** Calcul des valeurs propres
  167. C
  168. XU = (((1.D0 - PHIG**2) * UNG)**2) + (4.D0 * ((PHIG * ASONG)**2))
  169. SG= ((1.D0 + (PHIG**2)) * UNG) - (XU**0.5D0)
  170. SG= SG*HALF
  171. C write(*,*) SG,(UNG-ASONG)
  172. C
  173. XU = (((1.D0 - PHIAV**2) * UNAVE)**2) +
  174. & (4.D0 * ((PHIAV * ASONAV)**2))
  175. SGS= ((1.D0 + (PHIAV**2)) * UNAVE) - (XU**0.5D0)
  176. SGS= SGS*HALF
  177. C write(*,*) SGS,(UNAVE-ASONAV)
  178. C
  179. SG=MIN(SG,SGS)
  180. C
  181. SDS= ((1.D0 + (PHIAV**2)) * UNAVE) + (XU**0.5D0)
  182. SDS= SDS*HALF
  183. C write(*,*) SDS,(UNAVE+ASONAV)
  184. C
  185. XU = (((1.D0 - PHID**2) * UND)**2) + (4.D0 * ((PHID * ASOND)**2))
  186. SD= ((1.D0 + (PHID**2)) * UND) + (XU**0.5D0)
  187. SD= SD*HALF
  188. C write(*,*) SD,(UND+ASOND)
  189. C
  190. SD=MAX(SD,SDS)
  191. C
  192. C**** Speed of the contact discontinuity
  193. C
  194. SM=ROD*UND*(SD-UND)-ROG*UNG*(SG-UNG)+PG-PD
  195. SM=SM/((ROD*(SD-UND))-(ROG*(SG-UNG)))
  196. C
  197. C**** Interfacial pressure
  198. C
  199. PSTAR=ROG*(UNG-SG)*(UNG-SM)+PG
  200. RETG=(USGGM1*PG)+(HALF*ROG*((UNG*UNG)+(UTG*UTG)))
  201. RETD=(USGDM1*PD)+(HALF*ROD*((UND*UND)+(UTD*UTD)))
  202. C
  203. IF(SG .GT. 0)THEN
  204. FLU1(1)=ROG*UNG
  205. FLU1(2)=ROG*UNG*UNG+PG
  206. FLU1(4)=UNG*(RETG+PG)
  207. ELSEIF(SM .GT. 0)THEN
  208. OMEGA=1.0D0/(SG-SM)
  209. ROGS=ROG*(SG-UNG)
  210. ROGS=ROGS*OMEGA
  211. RUNGS=((SG-UNG)*(ROG*UNG))+(PSTAR - PG)
  212. RUNGS=RUNGS*OMEGA
  213. RETGS=((SG-UNG)*RETG)+(PSTAR*SM)-(PG*UNG)
  214. RETGS=RETGS*OMEGA
  215. C
  216. FLU1(1)=ROGS*SM
  217. FLU1(2)=RUNGS*SM+PSTAR
  218. FLU1(4)=(RETGS+PSTAR)*SM
  219. C
  220. ELSEIF(SD .GT. 0)THEN
  221. OMEGA=1.0D0/(SD-SM)
  222. RODS=ROD*(SD-UND)
  223. RODS=RODS*OMEGA
  224. RUNDS=((SD-UND)*(ROD*UND))+(PSTAR - PD)
  225. RUNDS=RUNDS*OMEGA
  226. RETDS=((SD-UND)*RETD)+(PSTAR*SM)-(PD*UND)
  227. RETDS=RETDS*OMEGA
  228. C
  229. FLU1(1)=RODS*SM
  230. FLU1(2)=RUNDS*SM+PSTAR
  231. FLU1(4)=(RETDS+PSTAR)*SM
  232. ELSE
  233. FLU1(1)=ROD*UND
  234. FLU1(2)=ROD*UND*UND+PD
  235. FLU1(4)=UND*(RETD+PD)
  236. ENDIF
  237. C
  238. FLU1(3)=(FLU1(1)*(UTG+UTD))+(ABS(FLU1(1))*(UTG-UTD))
  239. FLU1(3)=HALF*FLU1(3)
  240. C
  241. DO IESP=1,NESP,1
  242. FLU1(4+IESP)=(FLU1(1)*(YG(IESP)+YD(IESP)))+
  243. & (ABS(FLU1(1))*(YG(IESP)-YD(IESP)))
  244. FLU1(4+IESP)=HALF*FLU1(4+IESP)
  245. ENDDO
  246. C
  247. CELLT=1.0D0/(MAX(SG,SD))
  248. C
  249. RETURN
  250. END
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  

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