Télécharger froebm.eso

Retour à la liste

Numérotation des lignes :

  1. C FROEBM SOURCE KLOCZKO 05/06/14 21:15:08 5111
  2. C FROEBM SOURCE BECC 02/06/18 21:15:55 4365
  3. SUBROUTINE FROEBM(NESP,
  4. & GAMG,ROG,PG,UNG,UTG,
  5. & GAMD,ROD,PD,UND,UTD,
  6. & YG,YD,VINF,FLU,
  7. & CELLT)
  8. C
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : FROEBM
  13. C
  14. C DESCRIPTION : Formulation Volumes Finis pour les Equations
  15. C d'Euler Multi-Espèces relatives à un mélange
  16. C de gaz parfaits.
  17. C
  18. C Calcul du flux aux interfaces avec la méthode
  19. C de Roe-Turkel (cf. Thèse Cécile Viozat) adaptée
  20. C aux écoulements à faible nombre de Mach.
  21. C
  22. C LANGUAGE : FORTRAN 77
  23. C
  24. C AUTEUR : T. KLOCZKO DM2S/SFME/LTMF
  25. C
  26. C************************************************************************
  27. C
  28. C APPELES :
  29. C
  30. C************************************************************************
  31. C
  32. C**** Entrées:
  33. C
  34. C NESP = nombre d'especes considérées dans les Equations
  35. C d'Euler
  36. C
  37. C GAMG, GAMD = les "gamma" du gaz (gauche et droite)
  38. C
  39. C ROG, ROD = les densités
  40. C
  41. C PG, PD = les pressions
  42. C
  43. C UNG, UND = vitesses normales
  44. C
  45. C UTG, UTD = vitesses tangentielles
  46. C
  47. C VINF = vitesse de cut-off
  48. C
  49. C**** Sorties:
  50. C
  51. C FLU = table du flux a l'interface dans le repaire
  52. C (n,t), i.e.
  53. C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  54. C rho*un*y1, ...)
  55. C
  56. C CELLT = condition de stabilité, i.e.
  57. C
  58. C dT/diamax < cellt
  59. C
  60. C************************************************************************
  61. C
  62. C HISTORIQUE (Anomalies et modifications éventuelles)
  63. C
  64. C HISTORIQUE : 29.04.05
  65. C
  66. C
  67. C
  68. C************************************************************************
  69. C
  70. C N.B.: Toutes les variables sont DECLAREES
  71. C
  72. C
  73. C IMPLICIT NONE
  74. C
  75. C*** Déclaration des indices de boucles et du nombre d'espèces
  76. INTEGER NESP,I1,I2,I3
  77. C
  78. C*** Déclaration des variables primitives
  79. REAL*8 ROG,PG,UNG,UTG,RETG,GAMG,
  80. & ROD,PD,UND,UTD,RETD,GAMD
  81. C
  82. C*** Déclaration des variables de cacul des incréments conservatifs
  83. REAL*8 XP(4)
  84. C
  85. C*** Déclaration des variables de l'état moyen de Roe
  86. REAL*8 SQROG,SQROD,ROM,ROS,
  87. & UNM,UTM,
  88. & HG,HD,HM,
  89. & CM,GAMM,GAM1
  90. C
  91. C*** Déclaration de la condition de stabilité
  92. REAL*8 CELLT
  93. C
  94. C*** Délaration des variables de calcul du Mach de référence
  95. REAL*8 CG,CD,PHIM,VINF
  96. C
  97. C*** Déclaration des variables de calcul des valeurs propres
  98. REAL*8 VP(4),XU
  99. C
  100. C*** Déclaration des variables de calcul des matrices de passage
  101. REAL*8 R,S,TU,Q2,
  102. & MPI(4,4),MP(4,4)
  103. C
  104. C*** Déclaration des variables de calcul de la dissipation
  105. REAL*8 VDIS(4),DS(4)
  106. C
  107. C*** Déclaration des variables de calcul du flux numérique
  108. REAL*8 FC(4),FLU(*)
  109. C
  110. C*** Déclaration de la partie multiespèce
  111. REAL*8 YG(*),YD(*),FLUM
  112. C
  113. C
  114. C
  115. C*** DEBUT DES CALCULS ***
  116. C
  117. C
  118. C**** We include in the cut-off UNG,UND,UTG,UTD in order to
  119. C avoid low diffusivity in stagnation regions
  120. C
  121. VINF=MAX(VINF,((UNG*UNG+UTG*UTG)**0.5D0))
  122. VINF=MAX(VINF,((UND*UND+UTD*UTD)**0.5D0))
  123. C
  124. C**** Calcul des Deltas XP
  125. C
  126. RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG)
  127. & + (0.5D0*ROG*((UNG*UNG)+(UTG*UTG)))
  128. RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD)
  129. & + (0.5D0*ROD*((UND*UND)+(UTD*UTD)))
  130. C
  131. XP(1) = ROG - ROD
  132. XP(2) = (ROG * UNG) - (ROD * UND)
  133. XP(3) = (ROG * UTG) - (ROD * UTD)
  134. XP(4) = RETG - RETD
  135. C
  136. C*** Calcul de l'état moyen de Roe
  137. C
  138. GAMM = SQRT(GAMG) * SQRT(GAMD)
  139. GAM1 = GAMM - 1.D0
  140. C
  141. SQROG = SQRT(ROG)
  142. SQROD = SQRT(ROD)
  143. ROM = SQROG * SQROD
  144. ROS = SQROG + SQROD
  145. C
  146. UNM = ((SQROG * UNG) + (SQROD * UND)) / ROS
  147. UTM = ((SQROG * UTG) + (SQROD * UTD)) / ROS
  148. C
  149. HG = (1.D0 / ROG) * (RETG + PG)
  150. HD = (1.D0 / ROD) * (RETD + PD)
  151. HM = ((SQROG * HG) + (SQROD * HD)) / ROS
  152. C
  153. CM = SQRT(GAM1 * (HM - 0.5D0 * (UNM**2 + UTM**2)))
  154. C
  155. C*** Calcul de la condition de stabilité
  156. C
  157. CELLT = 1.D0 / (CM + ABS(UNM))
  158. C
  159. C*** Calcul du nombre de Mach de référence PHIM
  160. C
  161. CG = SQRT((GAMG * PG) / ROG)
  162. CD = SQRT((GAMD * PD) / ROD)
  163. IF(VINF.GE.MAX(CD,CG))THEN
  164. PHIM = 1.D0
  165. ELSE
  166. PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD))
  167. ENDIF
  168. C
  169. C*** Calcul des valeurs propres
  170. C
  171. XU = ((1.D0 - PHIM**2) * UNM)**2 + 4.D0 * (PHIM * CM)**2
  172. C
  173. VP(1) = UNM
  174. VP(2) = UNM
  175. VP(3) = 0.5D0 * ((1.D0 + PHIM**2) * UNM + SQRT(XU))
  176. VP(4) = 0.5D0 * ((1.D0 + PHIM**2) * UNM - SQRT(XU))
  177. C
  178. C*** Calcul des matrices de passages
  179. C
  180. R = VP(3) - UNM * PHIM**2
  181. S = VP(4) - UNM * PHIM**2
  182. TU = 0.5D0 * (VP(4) - VP(3))
  183. Q2 = 0.5D0 * (UNM**2 + UTM**2)
  184. C
  185. MPI(1,1) = 1.D0 - GAM1 * Q2 / (CM**2)
  186. MPI(1,2) = GAM1 * UNM / (CM**2)
  187. MPI(1,3) = GAM1 * UTM / (CM**2)
  188. MPI(1,4) = -GAM1 / (CM**2)
  189. MPI(2,1) = UTM
  190. MPI(2,2) = 0.D0
  191. MPI(2,3) = -1.D0
  192. MPI(2,4) = 0.D0
  193. MPI(3,1) = (S * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU
  194. MPI(3,2) = -(S * UNM * GAM1 + (PHIM*CM)**2) / TU
  195. MPI(3,3) = -(S * UTM * GAM1) / TU
  196. MPI(3,4) = S * GAM1 /TU
  197. MPI(4,1) = -(R * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU
  198. MPI(4,2) = (R * UNM * GAM1 + (PHIM*CM)**2) / TU
  199. MPI(4,3) = (R * UTM * GAM1) / TU
  200. MPI(4,4) = -R * GAM1 / TU
  201. C
  202. MP(1,1) = 1.D0
  203. MP(1,2) = 0.D0
  204. MP(1,3) = 0.5D0 / (PHIM*CM)**2
  205. MP(1,4) = 0.5D0 / (PHIM*CM)**2
  206. MP(2,1) = UNM
  207. MP(2,2) = 0.D0
  208. MP(2,3) = 0.5D0 * (UNM + R) / (PHIM*CM)**2
  209. MP(2,4) = 0.5D0 * (UNM + S) / (PHIM*CM)**2
  210. MP(3,1) = UTM
  211. MP(3,2) = -1.D0
  212. MP(3,3) = 0.5D0 * UTM / (PHIM*CM)**2
  213. MP(3,4) = 0.5D0 * UTM / (PHIM*CM)**2
  214. MP(4,1) = Q2
  215. MP(4,2) = -UTM
  216. MP(4,3) = 0.5D0 * (HM + R * UNM) / (PHIM*CM)**2
  217. MP(4,4) = 0.5D0 * (HM + S * UNM) / (PHIM*CM)**2
  218. C
  219. C*** Calcul de la dissipation
  220. C
  221. DO I1 = 1,4
  222. VDIS(I1) = 0.5D0 * ABS(VP(I1))
  223. ENDDO
  224. C
  225. DO I1 = 1,4
  226. DS(I1) = 0.D0
  227. DO I2 = 1,4
  228. DO I3 = 1,4
  229. DS(I1) = DS(I1) + MP(I1,I2) * VDIS(I2) * MPI(I2,I3) * XP(I3)
  230. ENDDO
  231. ENDDO
  232. ENDDO
  233. C
  234. C*** Calcul du flux convectif
  235. C
  236. FC(1) = 0.5D0 * ((ROG * UNG) + (ROD * UND))
  237. FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD))
  238. FC(3) = 0.5D0 * ((ROG * UNG * UTG) + (ROD * UND * UTD))
  239. FC(4) = 0.5D0 * ((UNG * (RETG + PG)) + (UND * (RETD + PD)))
  240. C
  241. C*** Calcul du flux numérique
  242. C
  243. DO I1 = 1,4
  244. FLU(I1) = FC(I1) + DS(I1)
  245. ENDDO
  246. C
  247. C*** Partie multiespèce
  248. C
  249. FLUM = FLU(1)
  250. IF(FLUM .GT. 0)THEN
  251. DO I1 = 1, NESP, 1
  252. FLU(4+I1)=FLUM * YG(I1)
  253. ENDDO
  254. ELSE
  255. DO I1 = 1, NESP, 1
  256. FLU(4+I1)=FLUM * YD(I1)
  257. ENDDO
  258. ENDIF
  259. C
  260. RETURN
  261. END
  262.  
  263.  
  264.  

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