Télécharger fcolt2.eso

Retour à la liste

Numérotation des lignes :

  1. C FCOLT2 SOURCE CHAT 05/01/12 23:57:22 5004
  2. SUBROUTINE FCOLT2(NESP,NSCA,NORDP1,XST,
  3. & GAMG,RTOTG,ACVTOG,ROG,PG,TG,UNG,UTG,UVG,ETHERG,
  4. & GAMD,RTOTD,ACVTOD,ROD,PD,TD,UND,UTD,UVD,ETHERD,
  5. & YG,YD,SCAG,SCAD,
  6. & FLU,CELLT,
  7. & LOGNC,MESERR,LOGAN)
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : FCOLT2
  13. C
  14. C DESCRIPTION : Formulation Volumes Finis pour les Equations
  15. C d'Euler Multi-Especes relatives à un melange
  16. C de gaz parfaits.
  17. C
  18. C Calcul du flux aux interfaces avec la methode
  19. C de Colella-Glaz, avec les chocs qui remplacent
  20. C les ondes de rarefaction + entropy fix.
  21. C
  22. C Cas 3D
  23. C
  24. C (voir:
  25. C 1) BECCANTINI, Thèse de doctorat)
  26. C
  27. C LANGUAGE : FORTRAN 77
  28. C
  29. C AUTEUR : A. BECCANTINI DRN/DMT/SEMT/LTMF
  30. C
  31. C************************************************************************
  32. C
  33. C APPELES : RACCOL
  34. C
  35. C************************************************************************
  36. C
  37. C**** Entrées:
  38. C
  39. C NESP = nombre d'especes considérées dans les Equations
  40. C d'Euler
  41. C
  42. C NSCA = nombre de scalaires passifs a transporter
  43. C
  44. C NORDP1 = (ordre des polynoms cv(T)) + 1
  45. C
  46. C XST = droite characteristique ou je cherche la solution
  47. C
  48. C GAMG, GAMD = les "gamma" du gaz (gauche et droite)
  49. C
  50. C RTOTG, RTOTD = les constantes des gaz
  51. C
  52. C ACVTOG, ACVTOD = les sommes des coefficients des cv,i
  53. C ACVTOT(j) = \sum_{i=1,nesp} Y_i*ACV_{i,j}
  54. C j = 1 , NORDP1
  55. C
  56. C ROG, ROD = les densités
  57. C
  58. C PG, PD = les pressions
  59. C
  60. C TG, TD = les temperatures
  61. C
  62. C UNG, UND = vitesses normales
  63. C
  64. C UTG,UVG,UTD,UVD = vitesses tangentielles
  65. C
  66. C ETHERG, ETHERD = les energies
  67. C
  68. C YG, YD = tables des fractiones massiques
  69. C
  70. C SCAG, SCAD = tables des scalaires passifs
  71. C
  72. C**** Sorties:
  73. C
  74. C FLU = table du flux a l'interface dans le repaire
  75. C (n,t), i.e.
  76. C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
  77. C rho*un*y1, ...)
  78. C
  79. C CELLT = condition de stabilité, i.e.
  80. C
  81. C dT/diamax < cellt
  82. C
  83. C LOGNC = si TRUE, la methode de newton pour le calcul
  84. C des etats intermediares ne converges pas
  85. C
  86. C MESERR = message d'erreuer
  87. C
  88. C LOGAN = si TRUE, une anomalie a été detectée
  89. C
  90. C************************************************************************
  91. C
  92. C HISTORIQUE (Anomalies et modifications éventuelles)
  93. C
  94. C HISTORIQUE : Créé le 08.02.00
  95. C
  96. C 21.02.00 transport de scalaires passifs
  97. C
  98. C************************************************************************
  99. C
  100. C N.B.: Toutes les variables sont DECLAREES
  101. C
  102. C
  103. IMPLICIT INTEGER(I-N)
  104. INTEGER NESP, NORDP1, I1, NSCA
  105. REAL*8 GAMG,RTOTG,ACVTOG(*),ROG,PG,TG,UNG,UTG,UVG,ETHERG
  106. & ,GAMD,RTOTD,ACVTOD(*),ROD,PD,TD,UND,UTD,UVD,ETHERD
  107. & ,YG(*),YD(*),SCAG(*),SCAD(*)
  108. & ,FLU(*),CELLT
  109. & ,PSROG,PSROD,ASONG,ASOND,PVAL,UVAL
  110. & ,HTHEG,HTHED
  111. & ,ROGS,PGS,TGS,ETHEGS,UGS,GAMGS
  112. & ,RODS,PDS,TDS,ETHEDS,UDS,GAMDS
  113. & ,XST,UNGXST,UNDXST,UGSXST,UDSXST,DGSXST,DDSXST,DGS,DDS
  114. & ,ROU,ECFLUX,XLIMG,XLIMD,ULIMG,ULIMD
  115. & ,UNGMCG,UNDPCD,PFLUX,ROFLUX,UNFLUX,UTFLUX,UVFLUX,HTFLUX
  116. & ,FC1,FC2,FPOIS1,FPOIS2,FPOIS3,FPOIS4,SIG12,SIG45,SIG3
  117. & ,FPOISG,FPOISD,EPSU
  118. & ,UGMCGS,UDPCDS,VITG,VITGS,VITD,VITDS,SIGCHO,SIG123,SIG456
  119. & ,SIG56
  120. & ,FUNT, TFLUX, RTOTFL, ETHFLU
  121. CHARACTER*(40) MESERR
  122. LOGICAL LOGAN, LOGNC
  123. C
  124. C**** YG, YD, FLU
  125. C
  126. C Dans le cas Euler monoespece, on doit
  127. C avoir :
  128. C YG(1) = YD(1) = 0.0D0
  129. C
  130. PSROG = PG / ROG
  131. PSROD = PD / ROD
  132. ASONG = SQRT(GAMG * PSROG)
  133. ASOND = SQRT(GAMD * PSROD)
  134. C
  135. C**** Valeurs de référence
  136. C
  137. PVAL = MAX(PG,PD)
  138. UVAL = MAX(ASONG,ASOND)
  139. EPSU = UVAL*1.0D-6
  140. C
  141. C**** États intermédiaires
  142. C
  143. HTHEG = ETHERG + PSROG
  144. HTHED = ETHERD + PSROD
  145. XLIMG = PSROG / (2.0D0 * HTHEG - PSROG)
  146. XLIMD = PSROD / (2.0D0 * HTHED - PSROD)
  147. ULIMG = UNG + SQRT( 2.0D0 * HTHEG *
  148. & (1.0D0 - XLIMG) / (1.0D0 + XLIMG))
  149. ULIMD = UND - SQRT( 2.0D0 * HTHED *
  150. & (1.0D0 - XLIMD) / (1.0D0 + XLIMD))
  151. UNGMCG = UNG - ASONG
  152. UNDPCD = UND + ASOND
  153. IF(ULIMD .GE. ULIMG)THEN
  154. C
  155. C******* Vide
  156. C
  157. C N.B.: la surface de contact n'existe pas: on mit le vide
  158. C
  159. IF((ULIMG .LE. 0.0D0) .AND. (ULIMD .GE. 0.0D0))THEN
  160. FLU(1) = 0.0D0
  161. FLU(2) = 0.0D0
  162. FLU(3) = 0.0D0
  163. FLU(4) = 0.0D0
  164. FLU(5) = 0.0D0
  165. DO I1 = 1, NESP, 1
  166. FLU(I1+5) = 0.0D0
  167. ENDDO
  168. DO I1 = 1, NSCA, 1
  169. FLU(I1+5+NESP) = 0.0D0
  170. ENDDO
  171. C
  172. ELSE
  173. ROGS = XLIMG * ROG
  174. RODS = XLIMD * ROD
  175. C
  176. C********** Les functions pois
  177. C
  178. FC1 = (ULIMG - XST)/(ULIMG - UNGMCG)
  179. SIG12 = 0.5D0 * (SIGN(1.0D0,FC1)+1.0D0)
  180. FC2 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
  181. FPOIS1 = SIG12*FC2
  182. C
  183. FC1 = (XST - ULIMD)/(UNDPCD - ULIMD)
  184. SIG45 = 0.5D0 * (SIGN(1.0D0,FC1)+1.0D0)
  185. FC2 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
  186. FPOIS4 = SIG45*FC2
  187. C
  188. SIG3 = 1.0D0 - (SIG12+SIG45)
  189. FC1 = (ULIMD - XST)/(ULIMD - ULIMG + EPSU)
  190. FC1 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
  191. FC2 = 1.0D0 - FC1
  192. FPOIS2 = SIG12 * (1.0D0 - FPOIS1) + SIG3 * FC1
  193. FPOIS3 = SIG45 * (1.0D0 - FPOIS4) + SIG3 * FC2
  194. C
  195. C********** La solution pour le calcul du flux
  196. C
  197. ROFLUX = ROG * FPOIS1 + ROGS * FPOIS2
  198. & + RODS * FPOIS3 + ROD * FPOIS4
  199. TFLUX = TG * FPOIS1 + TD * FPOIS4
  200. UNFLUX = UNG * FPOIS1 + ULIMG * FPOIS2
  201. & + ULIMD * FPOIS3 + UND * FPOIS4
  202. C
  203. FPOISG = FPOIS1 + FPOIS2
  204. FPOISD = FPOIS3 + FPOIS4
  205. RTOTFL = RTOTG * FPOISG + RTOTD * FPOISD
  206. PFLUX = ROFLUX * RTOTFL * TFLUX
  207. UTFLUX = UTG * FPOISG + UTD * FPOISD
  208. UVFLUX = UVG * FPOISG + UVD * FPOISD
  209. ECFLUX = 0.5D0 * (UNFLUX * UNFLUX + UTFLUX * UTFLUX
  210. & + UVFLUX * UVFLUX)
  211. C
  212. FUNT = 1.0D0
  213. ETHFLU = 0.0D0
  214. DO I1 = 1, NORDP1, 1
  215. FUNT = FUNT * TFLUX
  216. ETHFLU = ETHFLU + (FPOISG * ACVTOG(I1) * FUNT / I1) +
  217. & (FPOISD * ACVTOD(I1) * FUNT / I1)
  218. ENDDO
  219. HTFLUX = ETHFLU + PFLUX / ROFLUX
  220. C
  221. ROU = ROFLUX * UNFLUX
  222. FLU(1) = ROU
  223. FLU(2) = ROU * UNFLUX + PFLUX
  224. FLU(3) = ROU * UTFLUX
  225. FLU(4) = ROU * UVFLUX
  226. FLU(5) = ROU * (HTFLUX + ECFLUX)
  227. DO I1 = 1, NESP, 1
  228. FLU(I1+5) = ROU *
  229. & (FPOISG * YG(I1) + FPOISD * YD(I1))
  230. ENDDO
  231. DO I1 = 1, NSCA, 1
  232. FLU(I1+5+NESP) = ROU *
  233. & (FPOISG * SCAG(I1) + FPOISD * SCAD(I1))
  234. ENDDO
  235. ENDIF
  236. CELLT = MAX(ABS(UNGMCG),ABS(UNDPCD))
  237. C
  238. C
  239. C******* Test EF
  240. C
  241. C write(4,'(5D14.6)') XST, FPOIS1, FPOIS2, FPOIS3, FPOIS4
  242. C
  243. ELSE
  244. UNGXST=UNG-XST
  245. UNDXST=UND-XST
  246. CALL RACCOL(PVAL,UVAL,NORDP1,
  247. & RTOTG,ACVTOG,RTOTD,ACVTOD,
  248. & ROG,TG,ETHERG,UNGXST,PSROG,GAMG,ASONG,
  249. & ROD,TD,ETHERD,UNDXST,PSROD,GAMD,ASOND,
  250. & ROGS,PGS,TGS,ETHEGS,UGSXST,DGSXST,GAMGS,
  251. & RODS,PDS,TDS,ETHEDS,UDSXST,DDSXST,GAMDS,
  252. & LOGNC)
  253. C
  254. IF(LOGNC)THEN
  255. MESERR = 'Subroutine RACCOL.ESO, Newton-Raphson.'
  256. C
  257. C******* Message de warning; mais on s'arrete pas
  258. C
  259. ELSEIF((DGSXST .GT. UGSXST) .OR. (UDSXST .GT. DDSXST))THEN
  260. LOGAN = .TRUE.
  261. MESERR = 'Subroutine RACCOL.ESO'
  262. GOTO 9999
  263. ENDIF
  264. UGS=UGSXST+XST
  265. UDS=UDSXST+XST
  266. DGS=DGSXST+XST
  267. DDS=DDSXST+XST
  268. UGMCGS=UGS-SQRT(GAMGS*PGS/ROGS)
  269. UDPCDS=UDS+SQRT(GAMDS*PDS/RODS)
  270.  
  271. C
  272. C******* Tests
  273. C
  274. C write(*,*) ROGS,PGS,TGS,ETHEGS,UGS,GAMGS
  275. C write(*,*) RODS,PDS,TDS,ETHEDS,UDS,GAMDS
  276. C
  277. C
  278. C******* Les functions pois
  279. C
  280. C Les vitesses characteristiques des ondes
  281. C
  282. VITG = 0.5D0*(DGS+UNGMCG-ABS(DGS-UNGMCG))
  283. VITGS = 0.5D0*(DGS+UGMCGS+ABS(DGS-UGMCGS))
  284. VITD = 0.5D0*(DDS+UNDPCD+ABS(DDS-UNDPCD))
  285. VITDS = 0.5D0*(DDS+UDPCDS-ABS(DDS-UDPCDS))
  286. C
  287. C******* ZONES 1,2,3
  288. C
  289. FC1 = (VITGS-XST)/(ABS(VITGS - VITG)+EPSU)
  290. SIG12 = 0.5D0*(SIGN(1.0D0,FC1)+1.0D0)
  291. FC2 = 0.5D0*(FC1 + 1.0D0 - ABS(FC1 - 1.0D0))
  292. SIGCHO = 0.5D0 *(SIGN(1.0D0,(VITG -VITGS+EPSU))+1.0D0)
  293. SIG123 = 0.5D0 *(SIGN(1.0D0,((0.5D0*(UGS+UDS))-XST))+1.0D0)
  294. C
  295. FPOIS2 = (1.0D0 - SIGCHO)*(1.0D0 - FC2*SIG12)*SIG123 +
  296. & SIGCHO * (SIG123 - SIG12)
  297. FPOIS1 = (1.0D0 -FPOIS2)*SIG123
  298. C
  299. C******** ZONES 4,5,6
  300. C
  301. SIG456 = 1.0D0 - SIG123
  302. FC1 = (XST - VITDS)/(ABS(VITD-VITDS)+EPSU)
  303. SIG56 = 0.5D0*(SIGN(1.0D0,FC1)+1.0D0)
  304. FC2 = 0.5D0 * (FC1 + 1.0D0 - ABS(FC1 - 1.0D0))
  305. SIGCHO=0.5D0*(SIGN(1.0D0,VITDS-VITD+EPSU)+1.0D0)
  306. C
  307. FPOIS3= (1.0D0 - SIGCHO)*(1.0D0 - FC2*SIG56)*SIG456 +
  308. & SIGCHO * (SIG456 - SIG56)
  309. FPOIS4 = (1.0D0 - FPOIS3) * SIG456
  310. C
  311. C******* La solution pour le calcul du flux
  312. C (ro, T, w lineaire sur les detantes
  313. C p = ro R T
  314. C h = sum_i ACVTO(i)/i T^i
  315. C
  316. C
  317. ROFLUX = ROG * FPOIS1 + ROGS * FPOIS2
  318. & + RODS * FPOIS3 + ROD * FPOIS4
  319. TFLUX = TG * FPOIS1 + TGS * FPOIS2
  320. & + TDS * FPOIS3 + TD * FPOIS4
  321. UNFLUX = UNG * FPOIS1 + UGS * FPOIS2
  322. & + UDS * FPOIS3 + UND * FPOIS4
  323. C
  324. FPOISG = FPOIS1 + FPOIS2
  325. FPOISD = FPOIS3 + FPOIS4
  326. RTOTFL = RTOTG * FPOISG + RTOTD * FPOISD
  327. PFLUX = ROFLUX * RTOTFL * TFLUX
  328. UTFLUX = UTG * FPOISG + UTD * FPOISD
  329. UVFLUX = UVG * FPOISG + UVD * FPOISD
  330. ECFLUX = 0.5D0 * (UNFLUX * UNFLUX + UTFLUX * UTFLUX
  331. & + UVFLUX * UVFLUX)
  332. C
  333. FUNT = 1.0D0
  334. ETHFLU = 0.0D0
  335. DO I1 = 1, NORDP1, 1
  336. FUNT = FUNT * TFLUX
  337. ETHFLU = ETHFLU + (FPOISG * ACVTOG(I1) * FUNT / I1) +
  338. & (FPOISD * ACVTOD(I1) * FUNT / I1)
  339. ENDDO
  340. HTFLUX = ETHFLU + PFLUX / ROFLUX
  341. C
  342. ROU = ROFLUX * UNFLUX
  343. FLU(1) = ROU
  344. FLU(2) = ROU * UNFLUX + PFLUX
  345. FLU(3) = ROU * UTFLUX
  346. FLU(4) = ROU * UVFLUX
  347. FLU(5) = ROU * (HTFLUX + ECFLUX)
  348. DO I1 = 1, NESP, 1
  349. FLU(I1+5) = ROU *
  350. & (FPOISG * YG(I1) + FPOISD * YD(I1))
  351. ENDDO
  352. DO I1 = 1, NSCA, 1
  353. FLU(I1+5+NESP) = ROU *
  354. & (FPOISG * SCAG(I1) + FPOISD * SCAD(I1))
  355. ENDDO
  356. CC
  357. CELLT = MAX(ABS(VITG),ABS(VITD))
  358. C
  359. C******* Tests: les vitesses
  360. C
  361. C write(3,'(A8,D16.8)') 'VITG =', VITG
  362. C write(3,'(A8,D16.8)') 'VITGS =', VITGS
  363. C write(3,'(A8,2D16.8)') 'VITS =', UGS, UDS
  364. C write(3,'(A8,D16.8)') 'VITDS =', VITDS
  365. C write(3,'(A8,D16.8)') 'VITD =', VITD
  366. C
  367. C
  368. C******* Test EF
  369. C
  370. C write(4,'(5D14.6)') XST, FPOIS1, FPOIS2, FPOIS3, FPOIS4
  371. C
  372. ENDIF
  373. C
  374. 9999 CONTINUE
  375. RETURN
  376. END
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  

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