Télécharger fcolte.eso

Retour à la liste

Numérotation des lignes :

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

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