Télécharger pr12f.eso

Retour à la liste

Numérotation des lignes :

  1. C PR12F SOURCE KK2000 14/04/10 21:15:27 8032
  2. SUBROUTINE PR12f(IM1, ICH1, ICH2, ICH3,
  3. & ICH4, ICH5, ICH6, ICH7,
  4. & ICH8, ICH9,
  5. & Cp, Cvm,
  6. & OCH1, OCH2, OCH3,
  7. & OCH4, OCH5, OCH6,
  8. & OCH7, OCH8, OCH9,
  9. & LOGNEG, LOGBOR, MESERR,
  10. & VALER, VAL1, VAL2)
  11. C************************************************************************
  12. C
  13. C PROJET : CASTEM 2000
  14. C
  15. C NOM : PR12f
  16. C
  17. C DESCRIPTION : Primitive variables from conserved variables
  18. C
  19. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
  20. C
  21. C AUTEUR : Jose R. Garcia-Cascales,
  22. C Universidad Politecnica de Cartagena,
  23. C jr.garcia@upct.es
  24. C
  25. C************************************************************************
  26. C
  27. C APPELES (E/S) :
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES :
  32. C
  33. C IM1 : MELEME contenant les centres des ELTs
  34. C
  35. C ICH1 : CHPOINT contenant la masse volumique (mass, gas/vapour)
  36. C
  37. C ICH2 : CHPOINT contenant la masse volumique (mass, liquid)
  38. C
  39. C ICH3 : CHPOINT contenant les dèbits (momentum, gas/vapour)
  40. C ( NDIM composantes);
  41. C
  42. C ICH4 : CHPOINT contenant les dèbits (momentum, liquid)
  43. C ( NDIM composantes);
  44. C
  45. C ICH5 : CHPOINT contenat l'énergie totale per
  46. C unité de volume (alpha rho et) (gas/vapour);
  47. C
  48. C ICH6 : CHPOINT contenat l'énergie totale per
  49. C unité de volume (alpha rho et) (liquid);
  50. C
  51. C ICH7 : CHPOINT containing void fraction of the previous
  52. C time step, in order to avoid problems when a
  53. C phase disappears
  54. C
  55. C ICH8 : CHPOINT containing vapour temperature of the previous
  56. C time step, in order to avoid problems when a
  57. C phase disappears
  58. C
  59. C ICH9 : CHPOINT containing liquid temperature of the previous
  60. C time step, in order to avoid problems when the liquid
  61. C phase disappears
  62. C
  63. C Cp : Parameter characterizing the pressure
  64. C correction term
  65. C
  66. C Cvm : Parameter characterizing virtual mass correction,
  67. C it is included in the new CATHARE pressure
  68. C conrrection term
  69. C
  70. C SORTIES : OCH1 : CHPOINT void fraction (alpha);
  71. C
  72. C OCH2 : CHPOINT gas/vapour velocity (uvx, uvy, uvz);
  73. C
  74. C OCH3 : CHPOINT liquid velocity (ulx, uly, ulz);
  75. C
  76. C OCH4 : CHPOINT pressure;
  77. C
  78. C OCH5 : CHPOINT vapour temperature (Tv);
  79. C
  80. C OCH6 : CHPOINT liquid temperature (Tl);
  81. C
  82. C OCH7 : CHPOINT vapour density (rv);
  83. C
  84. C OCH8 : CHPOINT liquid density (rl);
  85. C
  86. C OCH9 : CHPOINT interfacial pressure correction term (pi)
  87. C
  88. C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité
  89. C negative a été detectée -> le programme s'arrete
  90. C (sa valeur stockée en MESERR(1) et VALER(1))
  91. C
  92. C LOGBOR : (LOGICAL)
  93. C gamma a été detecté dehor GAMMIN et GAMMAX
  94. C (sa valeur stockée en MESERR(2) et VALER(2),VAL1,VAL2)
  95. C
  96. C MESERR,
  97. C VALER,
  98. C VAL1,
  99. C VAL2 : pour message d'erreur
  100. C
  101. C
  102. C
  103. C************************************************************************
  104. C
  105. C HISTORIQUE (Anomalies et modifications éventuelles)
  106. C
  107. C HISTORIQUE : Créée le 20.11.05.
  108. C
  109. C************************************************************************
  110. C
  111. IMPLICIT INTEGER(I-N)
  112. IMPLICIT REAL*8(A-H,O-Z)
  113.  
  114.  
  115. C**** Les variables
  116. C
  117. INTEGER NLCE,N1,IGEOMC
  118. C
  119. INTEGER IM1
  120. & ,ICH1, ICH2, ICH3
  121. & ,ICH4, ICH5, ICH6
  122. & ,ICH7, ICH8, ICH9
  123. & ,OCH1,OCH2,OCH3
  124. & ,OCH4,OCH5,OCH6
  125. & ,OCH7,OCH8,OCH9
  126. C
  127. REAL*8 VALER(2),VAL1,VAL2
  128. & ,alpha, uvx, uvy, uvz, uv2, ulx, uly, ulz, ul2
  129. & ,p, Tv, Tl, rv, rl, ev, el
  130. & ,oldav, oldTv
  131. & ,A, B, delta
  132. & ,gamv, cpv, gaml, cpl, pil
  133. & ,FFV, FFL
  134. & ,rm, rvm, rlm, Cp, cvm, pi
  135. PARAMETER(gamv = 1.4D0, cpv = 1008.D0,
  136. & gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8)
  137. C
  138. CHARACTER*(8) TYPE
  139. CHARACTER*(40) MESERR(2)
  140. LOGICAL LOGNEG, LOGBOR
  141. C
  142. INTEGER NAT, NSOUPO, N, NC
  143. C
  144. C**** Les includes
  145. C
  146. -INC CCOPTIO
  147. -INC SMCHPOI
  148. POINTEUR MCHPO5.MCHPOI, MCHPO6.MCHPOI,
  149. & MCHPO7.MCHPOI, MCHPO8.MCHPOI, MCHPO9.MCHPOI
  150. C
  151. POINTEUR MSOUP6.MSOUPO, MSOUP7.MSOUPO, MSOUP8.MSOUPO,
  152. & MSOUP9.MSOUPO
  153. C
  154. POINTEUR MPOW1.MPOVAL, MPOW2.MPOVAL, MPOW3.MPOVAL,
  155. & MPOW4.MPOVAL, MPOW5.MPOVAL, MPOW6.MPOVAL,
  156. & MPOW7.MPOVAL, MPOW8.MPOVAL, MPOW9.MPOVAL
  157. C
  158. POINTEUR MPOV1.MPOVAL, MPOV2.MPOVAL, MPOV3.MPOVAL,
  159. & MPOV4.MPOVAL, MPOV5.MPOVAL, MPOV6.MPOVAL,
  160. & MPOV7.MPOVAL, MPOV8.MPOVAL, MPOV9.MPOVAL
  161. -INC SMELEME
  162. C
  163. C
  164. C**** Initialisation des variables pour la gestion des erreurs pas ici,
  165. C mais avant, i.e.
  166. C
  167. LOGNEG = .FALSE.
  168. LOGBOR = .FALSE.
  169. MESERR(1) = ' '
  170. MESERR(2) = ' '
  171. VALER(1) = 0.0D0
  172. VALER(2) = 0.0D0
  173. VAL1 = 0.0D0
  174. VAL2 = 0.0D0
  175.  
  176. C
  177. C**** Activation du MELEME "CENTRE"
  178. C
  179. IPT1 = IM1
  180. SEGACT IPT1
  181. N1 = IPT1.NUM(/2)
  182. SEGDES IPT1
  183. C
  184. C**** Reading MPOVA of each MCHAMPOIN
  185. C MPOW1 = alpha*rv
  186. C MPOW2 = (1 - alpha)*rl
  187. C MPOW3 = alpha*rv*uvx & alpha*rv*uvy & alpha*rv*uvz
  188. C MPOW4 = (1 - alpha)*rl*ulx & (1 - alpha)*rl*uly & (1 - alpha)*rl*ulz
  189. C MPOW5 = alpha*rv*Ev
  190. C MPOW6 = (1 - alpha)*rl*El
  191. C
  192. C MPOW7: It is the old void fraction, we will use it when
  193. C alpha approaches to zero
  194. C MPOW8: It is the old vapour temperature, we will use it when
  195. C alpha approaches to zero
  196. C MPOW9: It is the old liquid temperature, we will use it when
  197. C alpha approaches to one
  198. C
  199. CALL LICHT(ICH1,MPOW1,TYPE,IGEOMC)
  200. CALL LICHT(ICH2,MPOW2,TYPE,IGEOMC)
  201. CALL LICHT(ICH3,MPOW3,TYPE,IGEOMC)
  202. CALL LICHT(ICH4,MPOW4,TYPE,IGEOMC)
  203. CALL LICHT(ICH5,MPOW5,TYPE,IGEOMC)
  204. CALL LICHT(ICH6,MPOW6,TYPE,IGEOMC)
  205. CALL LICHT(ICH7,MPOW7,TYPE,IGEOMC)
  206. CALL LICHT(ICH8,MPOW8,TYPE,IGEOMC)
  207. CALL LICHT(ICH9,MPOW9,TYPE,IGEOMC)
  208. C
  209. C SEGACT MPOW1
  210. C SEGACT MPOW2
  211. C SEGACT MPOW3
  212. C SEGACT MPOW4
  213. C SEGACT MPOW5
  214. C SEGACT MPOW6
  215. C SEGACT MPOW7
  216. C SEGACT MPOW8
  217. C SEGACT MPOW9
  218. CC
  219. C**** LICHT active les MPOVALs en *MO
  220. C
  221. C
  222. C**** Privitive variables
  223. C
  224. C void fraction: alpha
  225.  
  226. NAT = 1
  227. NSOUPO=1
  228. N = N1
  229. NC = 1
  230. SEGINI MCHPO1
  231. MCHPO1.JATTRI(1)=2
  232. MCHPO1.IFOPOI=IFOMOD
  233. SEGINI MSOUP1
  234. MCHPO1.IPCHP(1)=MSOUP1
  235. SEGINI MPOV1
  236. MSOUP1.NOCOMP(1)='SCAL'
  237. MSOUP1.IGEOC = IM1
  238. MSOUP1.IPOVAL = MPOV1
  239.  
  240. C vapour velocities: uvx, uvy, uvz
  241.  
  242.  
  243. NC = IDIM
  244. SEGINI MCHPO2
  245. MCHPO2.JATTRI(1)=2
  246. MCHPO2.IFOPOI=IFOMOD
  247. SEGINI MSOUP2
  248. MCHPO2.IPCHP(1)=MSOUP2
  249. SEGINI MPOV2
  250. MSOUP2.NOCOMP(1)='UVX'
  251. MSOUP2.NOCOMP(2)='UVY'
  252. IF(IDIM .EQ. 3)THEN
  253. MSOUP2.NOCOMP(3)='UVZ'
  254. ENDIF
  255. MSOUP2.IGEOC = IM1
  256. MSOUP2.IPOVAL = MPOV2
  257.  
  258. C liquid velocities: ulx, uly, ulz
  259.  
  260.  
  261. NC = IDIM
  262. SEGINI MCHPO3
  263. MCHPO3.JATTRI(1)=2
  264. MCHPO3.IFOPOI=IFOMOD
  265. SEGINI MSOUP3
  266. MCHPO3.IPCHP(1)=MSOUP3
  267. SEGINI MPOV3
  268. MSOUP3.NOCOMP(1)='ULX'
  269. MSOUP3.NOCOMP(2)='ULY'
  270. IF(IDIM .EQ. 3)THEN
  271. MSOUP3.NOCOMP(3)='ULZ'
  272. ENDIF
  273. MSOUP3.IGEOC = IM1
  274. MSOUP3.IPOVAL = MPOV3
  275.  
  276. C pressure: p
  277.  
  278. NC = 1
  279. SEGINI MCHPO4
  280. MCHPO4.JATTRI(1)=2
  281. MCHPO4.IFOPOI=IFOMOD
  282. SEGINI MSOUP4
  283. MCHPO4.IPCHP(1)=MSOUP4
  284. SEGINI MPOV4
  285. MSOUP4.NOCOMP(1)='SCAL'
  286. MSOUP4.IGEOC = IM1
  287. MSOUP4.IPOVAL = MPOV4
  288.  
  289. C vapour temperature: Tv
  290.  
  291. NC = 1
  292. SEGINI MCHPO5
  293. MCHPO5.JATTRI(1)=2
  294. MCHPO5.IFOPOI=IFOMOD
  295. SEGINI MSOUP5
  296. MCHPO5.IPCHP(1)=MSOUP5
  297. SEGINI MPOV5
  298. MSOUP5.NOCOMP(1)='SCAL'
  299. MSOUP5.IGEOC = IM1
  300. MSOUP5.IPOVAL = MPOV5
  301.  
  302. C liquid temperature: Tl
  303.  
  304. NC = 1
  305. SEGINI MCHPO6
  306. MCHPO6.JATTRI(1)=2
  307. MCHPO6.IFOPOI=IFOMOD
  308. SEGINI MSOUP6
  309. MCHPO6.IPCHP(1)=MSOUP6
  310. SEGINI MPOV6
  311. MSOUP6.NOCOMP(1)='SCAL'
  312. MSOUP6.IGEOC = IM1
  313. MSOUP6.IPOVAL = MPOV6
  314.  
  315. C vapour density: rv
  316.  
  317. NC = 1
  318. SEGINI MCHPO7
  319. MCHPO7.JATTRI(1)=2
  320. MCHPO7.IFOPOI=IFOMOD
  321. SEGINI MSOUP7
  322. MCHPO7.IPCHP(1)=MSOUP7
  323. SEGINI MPOV7
  324. MSOUP7.NOCOMP(1)='SCAL'
  325. MSOUP7.IGEOC = IM1
  326. MSOUP7.IPOVAL = MPOV7
  327.  
  328. C liquid density: rl
  329.  
  330. NC = 1
  331. SEGINI MCHPO8
  332. MCHPO8.JATTRI(1)=2
  333. MCHPO8.IFOPOI=IFOMOD
  334. SEGINI MSOUP8
  335. MCHPO8.IPCHP(1)=MSOUP8
  336. SEGINI MPOV8
  337. MSOUP8.NOCOMP(1)='SCAL'
  338. MSOUP8.IGEOC = IM1
  339. MSOUP8.IPOVAL = MPOV8
  340.  
  341. C pressure correction term, pi
  342.  
  343. NC = 1
  344. SEGINI MCHPO9
  345. MCHPO9.JATTRI(1)=2
  346. MCHPO9.IFOPOI=IFOMOD
  347. SEGINI MSOUP9
  348. MCHPO9.IPCHP(1)=MSOUP9
  349. SEGINI MPOV9
  350. MSOUP9.NOCOMP(1)='SCAL'
  351. MSOUP9.IGEOC = IM1
  352. MSOUP9.IPOVAL = MPOV9
  353. C
  354. C**** "Recapitulando"
  355. C
  356. C We activate MPOV1 - MPOV9
  357. C
  358. C**** "bucle" for the calculation of the primitives at
  359. C each centre
  360.  
  361.  
  362. DO NLCE = 1, N1
  363.  
  364. oldav = MPOW7.VPOCHA(NLCE,1)
  365. oldTv = MPOW8.VPOCHA(NLCE,1)
  366.  
  367. IF(IDIM .EQ. 3)THEN
  368. if(oldav .LT. 0.01d0)then
  369. A = (gamv - 1.D0)*cpv*oldTv*MPOW1.VPOCHA(NLCE,1)/gamv
  370. else
  371. A = (gamv - 1.D0)*(MPOW5.VPOCHA(NLCE,1) - (MPOW3
  372. $ .VPOCHA(NLCE,1)**2 + MPOW3.VPOCHA(NLCE,2)**2 +
  373. & MPOW3.VPOCHA(NLCE,3)**2)/
  374. & (2.d0*(MPOW1.VPOCHA(NLCE,1)+ 1.D-20)))
  375. endif
  376. if(oldav .GT. 0.99d0)then
  377. B = (gaml - 1.D0)*MPOW2.VPOCHA(NLCE,1)*
  378. & (cpl*MPOW9.VPOCHA(NLCE,1)/gaml + pil/1000.d0)
  379. else
  380. B = (gaml - 1.D0)*(MPOW6.VPOCHA(NLCE,1) - (MPOW4
  381. $ .VPOCHA(NLCE,1)**2 + MPOW4.VPOCHA(NLCE,2)**2)+
  382. & MPOW4.VPOCHA(NLCE,3)**2/(2.D0
  383. $ *(MPOW2.VPOCHA(NLCE,1)+ 1.D-20)))
  384. endif
  385. ELSE
  386. if(oldav .LT. 0.01d0)then
  387. A = (gamv - 1.D0)*cpv*oldTv*MPOW1.VPOCHA(NLCE,1)/gamv
  388. else
  389. A = (gamv - 1.D0)*(MPOW5.VPOCHA(NLCE,1) - (MPOW3
  390. $ .VPOCHA(NLCE,1)**2 + MPOW3.VPOCHA(NLCE,2)**2)/(2.d0
  391. $ *(MPOW1.VPOCHA(NLCE,1)+ 1.D-20)))
  392. endif
  393. if(oldav .GT. 0.99d0)then
  394. B = (gaml - 1.D0)*MPOW2.VPOCHA(NLCE,1)*
  395. & (cpl*MPOW9.VPOCHA(NLCE,1)/gaml + pil/1000.d0)
  396. else
  397.  
  398. B = (gaml - 1.D0)*(MPOW6.VPOCHA(NLCE,1) - (MPOW4
  399. $ .VPOCHA(NLCE,1)**2 + MPOW4.VPOCHA(NLCE,2)**2)/(2.D0
  400. $ *(MPOW2.VPOCHA(NLCE,1)+ 1.D-20)))
  401. endif
  402. ENDIF
  403. delta = (- A - B + gaml*pil)**2.D0 + 4.D0*gaml*pil*A
  404.  
  405. p = 1.D0/2.D0*(A + B - gaml*pil + sqrt(delta))
  406.  
  407. alpha = A/p
  408. call FUNCF(alpha, FFV)
  409. uvx = FFV*MPOW3.VPOCHA(NLCE,1)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
  410. uvy = FFV*MPOW3.VPOCHA(NLCE,2)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
  411. IF(IDIM .EQ. 3)THEN
  412. uvz = FFV*MPOW3.VPOCHA(NLCE,3)/(MPOW1.VPOCHA(NLCE,1) + 1.D
  413. & -20)
  414. uv2 = uvx*uvx + uvy*uvy + uvz*uvz
  415. ELSE
  416. uv2 = uvx*uvx + uvy*uvy
  417. ENDIF
  418.  
  419. call FUNCF((1.d0 - alpha),FFL)
  420. ulx = FFL*MPOW4.VPOCHA(NLCE,1)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
  421. uly = FFL*MPOW4.VPOCHA(NLCE,2)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
  422. IF(IDIM .EQ. 3)THEN
  423. ulz = FFL*MPOW4.VPOCHA(NLCE,3)/(MPOW2.VPOCHA(NLCE,1) + 1.D
  424. & -20)
  425. ul2 = ulx*ulx + uly*uly + ulz*ulz
  426. ELSE
  427. ul2 = ulx*ulx + uly*uly
  428. ENDIF
  429.  
  430. if((alpha .GT. 0.999d0) .OR. (alpha .LT. 0.01d0))then
  431. C
  432. C Water faucet test
  433. C Tv = 323.15D0
  434. C Tl = 323.15D0
  435. C
  436. C Shock tube test
  437. C Tv = 308.15d0
  438. C Tl = 308.15d0
  439. C
  440. C Phase separation test or slosing of a water column
  441. Tv = 298.15D0
  442. Tl = 298.15D0
  443. C For the oscillating manometer should be 323 but to be practical I have left it in 298.15
  444. ev = cpv*Tv/gamv
  445. el = cpl*Tl/gaml*(p + gaml*pil)/(p + pil)
  446. else
  447. ev = MPOW5.VPOCHA(NLCE,1)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
  448. & - uv2/2.D0
  449. el = MPOW6.VPOCHA(NLCE,1)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
  450. & - ul2/2.D0
  451. Tv = gamv*ev/cpv
  452. Tl = gaml*el/(cpl*(1.D0 + pil*(gaml - 1.D0)/(p + pil)))
  453. endif
  454. rv = gamv*p/((gamv - 1.D0)*cpv*Tv)
  455. rl = gaml*(p + pil)/((gaml - 1.D0)*cpl*Tl)
  456. C
  457. C Bestion´s formula
  458. C
  459. rm = alpha*rv + (1.d0 - alpha)*rl
  460. rvm = rv + Cvm*rm
  461. rlm = rl + Cvm*rm
  462.  
  463. pi = p - Cp*(alpha*rvm*(1.D0- alpha)*rlm)/
  464. & (alpha*rlm + (1.D0- alpha)*rvm)*
  465. & (sqrt(uv2) - sqrt(ul2))**2
  466.  
  467. MPOV1.VPOCHA(NLCE,1) = alpha
  468. MPOV2.VPOCHA(NLCE,1) = uvx
  469. MPOV2.VPOCHA(NLCE,2) = uvy
  470. MPOV3.VPOCHA(NLCE,1) = ulx
  471. MPOV3.VPOCHA(NLCE,2) = uly
  472. IF(IDIM .EQ. 3)THEN
  473. MPOV2.VPOCHA(NLCE,3) = uvz
  474. MPOV3.VPOCHA(NLCE,3) = ulz
  475. ENDIF
  476. MPOV4.VPOCHA(NLCE,1) = p
  477. MPOV5.VPOCHA(NLCE,1) = Tv
  478. MPOV6.VPOCHA(NLCE,1) = Tl
  479. MPOV7.VPOCHA(NLCE,1) = rv
  480. MPOV8.VPOCHA(NLCE,1) = rl
  481. MPOV9.VPOCHA(NLCE,1) = pi
  482.  
  483. ENDDO
  484. C
  485. OCH1 = MCHPO1
  486. OCH2 = MCHPO2
  487. OCH3 = MCHPO3
  488. OCH4 = MCHPO4
  489. OCH5 = MCHPO5
  490. OCH6 = MCHPO6
  491. OCH7 = MCHPO7
  492. OCH8 = MCHPO8
  493. OCH9 = MCHPO9
  494. C
  495. SEGDES MPOV1
  496. SEGDES MPOV2
  497. SEGDES MPOV3
  498. SEGDES MPOV4
  499. SEGDES MPOV5
  500. SEGDES MPOV6
  501. SEGDES MPOV7
  502. SEGDES MPOV8
  503. SEGDES MPOV9
  504.  
  505. SEGDES MPOW1
  506. SEGDES MPOW2
  507. SEGDES MPOW3
  508. SEGDES MPOW4
  509. SEGDES MPOW5
  510. SEGDES MPOW6
  511. SEGDES MPOW7
  512. SEGDES MPOW8
  513. SEGDES MPOW9
  514.  
  515. SEGDES MSOUP1
  516. SEGDES MSOUP2
  517. SEGDES MSOUP3
  518. SEGDES MSOUP4
  519. SEGDES MSOUP5
  520. SEGDES MSOUP6
  521. SEGDES MSOUP7
  522. SEGDES MSOUP8
  523. SEGDES MSOUP9
  524.  
  525. SEGDES MCHPO1
  526. SEGDES MCHPO2
  527. SEGDES MCHPO3
  528. SEGDES MCHPO4
  529. SEGDES MCHPO5
  530. SEGDES MCHPO6
  531. SEGDES MCHPO7
  532. SEGDES MCHPO8
  533. SEGDES MCHPO9
  534.  
  535. RETURN
  536. END
  537.  
  538.  
  539.  
  540.  
  541.  
  542.  
  543.  
  544.  
  545.  
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  

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