Télécharger fcolcp.eso

Retour à la liste

Numérotation des lignes :

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

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