Télécharger froeb3.eso

Retour à la liste

Numérotation des lignes :

froeb3
  1. C FROEB3 SOURCE KLOCZKO 05/06/14 21:15:06 5111
  2. C FROEB3 SOURCE BECC 02/06/18 21:15:55 4365
  3. SUBROUTINE FROEB3(NESP,
  4. & GAMG,ROG,PG,UNG,UTG,UVG,
  5. & GAMD,ROD,PD,UND,UTD,UVD,
  6. & YG,YD,VINF,FLU,
  7. & CELLT)
  8. C
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : FROEB3
  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 UVG, UVD = vitesses tangentielles
  48. C
  49. C VINF = vitesse de cut-off
  50. C
  51. C**** Sorties:
  52. C
  53. C FLU = table du flux a l'interface dans le repère
  54. C (n,t), i.e.
  55. C (rho*un, rho*un*un + p, rho*un*ut, rho*un*uv,
  56. C rho*un*ht, rho*un*y1, ...)
  57. C
  58. C CELLT = condition de stabilité, i.e.
  59. C
  60. C dT/diamax < cellt
  61. C
  62. C************************************************************************
  63. C
  64. C HISTORIQUE (Anomalies et modifications éventuelles)
  65. C
  66. C HISTORIQUE : Créé le 02.05.05
  67. C
  68. C
  69. C
  70. C************************************************************************
  71. C
  72. C N.B.: Toutes les variables sont DECLAREES
  73. C
  74. C
  75. C IMPLICIT NONE
  76. C
  77. C*** Déclaration des indices de boucles et du nombre d'espèces
  78. INTEGER NESP,I1,I2,I3
  79. C
  80. C*** Déclaration des variables primitives
  81. REAL*8 ROG,PG,UNG,UTG,UVG,RETG,GAMG,
  82. & ROD,PD,UND,UTD,UVD,RETD,GAMD
  83. C
  84. C*** Déclaration des variables de cacul des incréments conservatifs
  85. REAL*8 XP(5)
  86. C
  87. C*** Déclaration des variables de l'état moyen de Roe
  88. REAL*8 SQROG,SQROD,ROM,ROS,
  89. & UNM,UTM,UVM,
  90. & HG,HD,HM,
  91. & CM,GAMM,GAM1
  92. C
  93. C*** Déclaration de la condition de stabilité
  94. REAL*8 CELLT
  95. C
  96. C*** Délaration des variables de calcul du Mach de référence
  97. REAL*8 CG,CD,PHIM,VINF
  98. C
  99. C*** Déclaration des variables de calcul des valeurs propres
  100. REAL*8 VP(5),XU
  101. C
  102. C*** Déclaration des variables de calcul des matrices de passage
  103. REAL*8 R,S,TU,Q2,
  104. & MPI(5,5),MP(5,5)
  105. C
  106. C*** Déclaration des variables de calcul de la dissipation
  107. REAL*8 VDIS(5),DS(5)
  108. C
  109. C*** Déclaration des variables de calcul du flux numérique
  110. REAL*8 FC(5),FLU(*)
  111. C
  112. C*** Déclaration de la partie multiespèce
  113. REAL*8 YG(*),YD(*),FLUM
  114. C
  115. C
  116. C
  117. C*** DEBUT DES CALCULS ***
  118. C
  119. C
  120. C**** We include in the cut-off UNG,UND,UTG,UTD in order to
  121. C avoid low diffusivity in stagnation regions
  122. C
  123. VINF=MAX(VINF,((UNG**2+UTG**2+UVG**2)**0.5D0))
  124. VINF=MAX(VINF,((UND**2+UTD**2+UVD**2)**0.5D0))
  125. C
  126. C**** Calcul des Deltas XP
  127. C
  128. RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG)
  129. & + (0.5D0*ROG*(UNG**2+UTG**2+UVG**2))
  130. RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD)
  131. & + (0.5D0*ROD*(UND**2+UTD**2+UVD**2))
  132. C
  133. XP(1) = ROG - ROD
  134. XP(2) = (ROG * UNG) - (ROD * UND)
  135. XP(3) = (ROG * UTG) - (ROD * UTD)
  136. XP(4) = (ROG * UVG) - (ROD * UVD)
  137. XP(5) = RETG - RETD
  138. C
  139. C*** Calcul de l'état moyen de Roe
  140. C
  141. GAMM = SQRT(GAMG) * SQRT(GAMD)
  142. GAM1 = GAMM - 1.D0
  143. C
  144. SQROG = SQRT(ROG)
  145. SQROD = SQRT(ROD)
  146. ROM = SQROG * SQROD
  147. ROS = SQROG + SQROD
  148. C
  149. UNM = ((SQROG * UNG) + (SQROD * UND)) / ROS
  150. UTM = ((SQROG * UTG) + (SQROD * UTD)) / ROS
  151. UVM = ((SQROG * UVG) + (SQROD * UVD)) / ROS
  152. C
  153. HG = (1.D0 / ROG) * (RETG + PG)
  154. HD = (1.D0 / ROD) * (RETD + PD)
  155. HM = ((SQROG * HG) + (SQROD * HD)) / ROS
  156. C
  157. CM = SQRT(GAM1 * (HM - 0.5D0 * (UNM**2 + UTM**2+UVM**2)))
  158. C
  159. C*** Calcul de la condition de stabilité
  160. C
  161. CELLT = 1.D0 / (CM + ABS(UNM))
  162. C
  163. C*** Calcul du nombre de Mach de référence PHIM
  164. C
  165. CG = SQRT((GAMG * PG) / ROG)
  166. CD = SQRT((GAMD * PD) / ROD)
  167. IF(VINF.GE.MAX(CD,CG))THEN
  168. PHIM = 1.D0
  169. ELSE
  170. PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD))
  171. ENDIF
  172. C
  173. C*** Calcul des valeurs propres
  174. C
  175. XU = ((1.D0 - PHIM**2) * UNM)**2 + 4.D0 * (PHIM * CM)**2
  176. C
  177. VP(1) = UNM
  178. VP(2) = UNM
  179. VP(3) = UNM
  180. VP(4) = 0.5D0 * ((1.D0 + PHIM**2) * UNM + SQRT(XU))
  181. VP(5) = 0.5D0 * ((1.D0 + PHIM**2) * UNM - SQRT(XU))
  182. C
  183. C*** Calcul des matrices de passages
  184. C
  185. R = VP(4) - UNM * PHIM**2
  186. S = VP(5) - UNM * PHIM**2
  187. TU = 0.5D0 * (VP(5) - VP(4))
  188. Q2 = 0.5D0 * (UNM**2 + UTM**2+UVM**2)
  189. C
  190. MPI(1,1) = 1.D0 - GAM1 * Q2 / (CM**2)
  191. MPI(1,2) = GAM1 * UNM / (CM**2)
  192. MPI(1,3) = GAM1 * UTM / (CM**2)
  193. MPI(1,4) = GAM1 * UVM / (CM**2)
  194. MPI(1,5) = -GAM1 / (CM**2)
  195. MPI(2,1) = UVM
  196. MPI(2,2) = 0.D0
  197. MPI(2,3) = 0.D0
  198. MPI(2,4) = -1.D0
  199. MPI(2,5) = 0.D0
  200. MPI(3,1) = UTM
  201. MPI(3,2) = 0.D0
  202. MPI(3,3) = -1.D0
  203. MPI(3,4) = 0.D0
  204. MPI(3,5) = 0.D0
  205. MPI(4,1) = (S * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU
  206. MPI(4,2) = -(S * UNM * GAM1 + (PHIM*CM)**2) / TU
  207. MPI(4,3) = -(S * UTM * GAM1) / TU
  208. MPI(4,4) = -(S * UVM * GAM1) / TU
  209. MPI(4,5) = S * GAM1 /TU
  210. MPI(5,1) = -(R * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU
  211. MPI(5,2) = (R * UNM * GAM1 + (PHIM*CM)**2) / TU
  212. MPI(5,3) = (R * UTM * GAM1) / TU
  213. MPI(5,4) = (R * UVM * GAM1) / TU
  214. MPI(5,5) = -R * GAM1 / TU
  215. C
  216. MP(1,1) = 1.D0
  217. MP(1,2) = 0.D0
  218. MP(1,3) = 0.D0
  219. MP(1,4) = 0.5D0 / (PHIM*CM)**2
  220. MP(1,5) = 0.5D0 / (PHIM*CM)**2
  221. MP(2,1) = UNM
  222. MP(2,2) = 0.D0
  223. MP(2,3) = 0.D0
  224. MP(2,4) = 0.5D0 * (UNM + R) / (PHIM*CM)**2
  225. MP(2,5) = 0.5D0 * (UNM + S) / (PHIM*CM)**2
  226. MP(3,1) = UTM
  227. MP(3,2) = 0.D0
  228. MP(3,3) = -1.D0
  229. MP(3,4) = 0.5D0 * UTM / (PHIM*CM)**2
  230. MP(3,5) = 0.5D0 * UTM / (PHIM*CM)**2
  231. MP(4,1) = UVM
  232. MP(4,2) = -1.D0
  233. MP(4,3) = 0.D0
  234. MP(4,4) = 0.5D0 * UVM / (PHIM*CM)**2
  235. MP(4,5) = 0.5D0 * UVM / (PHIM*CM)**2
  236. MP(5,1) = Q2
  237. MP(5,2) = -UVM
  238. MP(5,3) = -UTM
  239. MP(5,4) = 0.5D0 * (HM + R * UNM) / (PHIM*CM)**2
  240. MP(5,5) = 0.5D0 * (HM + S * UNM) / (PHIM*CM)**2
  241. C
  242. C*** Calcul de la dissipation
  243. C
  244. DO I1 = 1,5
  245. VDIS(I1) = 0.5D0 * ABS(VP(I1))
  246. ENDDO
  247. C
  248. DO I1 = 1,5
  249. DS(I1) = 0.D0
  250. DO I2 = 1,5
  251. DO I3 = 1,5
  252. DS(I1) = DS(I1) + MP(I1,I2) * VDIS(I2) * MPI(I2,I3) * XP(I3)
  253. ENDDO
  254. ENDDO
  255. ENDDO
  256. C
  257. C*** Calcul du flux convectif
  258. C
  259. FC(1) = 0.5D0 * ((ROG * UNG) + (ROD * UND))
  260. FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD))
  261. FC(3) = 0.5D0 * ((ROG * UNG * UTG) + (ROD * UND * UTD))
  262. FC(4) = 0.5D0 * ((ROG * UNG * UVG) + (ROD * UND * UVD))
  263. FC(5) = 0.5D0 * ((UNG * (RETG + PG)) + (UND * (RETD + PD)))
  264. C
  265. C*** Calcul du flux numérique
  266. C
  267. DO I1 = 1,5
  268. FLU(I1) = FC(I1) + DS(I1)
  269. ENDDO
  270. C
  271. C*** Partie multiespèce
  272. C
  273. FLUM = FLU(1)
  274. IF(FLUM .GT. 0)THEN
  275. DO I1 = 1, NESP, 1
  276. FLU(5+I1)=FLUM * YG(I1)
  277. ENDDO
  278. ELSE
  279. DO I1 = 1, NESP, 1
  280. FLU(5+I1)=FLUM * YD(I1)
  281. ENDDO
  282. ENDIF
  283. C
  284. RETURN
  285. END
  286.  
  287.  
  288.  

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