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.  
  147. -INC PPARAM
  148. -INC CCOPTIO
  149. -INC SMCHPOI
  150. POINTEUR MCHPO5.MCHPOI, MCHPO6.MCHPOI,
  151. & MCHPO7.MCHPOI, MCHPO8.MCHPOI, MCHPO9.MCHPOI
  152. C
  153. POINTEUR MSOUP6.MSOUPO, MSOUP7.MSOUPO, MSOUP8.MSOUPO,
  154. & MSOUP9.MSOUPO
  155. C
  156. POINTEUR MPOW1.MPOVAL, MPOW2.MPOVAL, MPOW3.MPOVAL,
  157. & MPOW4.MPOVAL, MPOW5.MPOVAL, MPOW6.MPOVAL,
  158. & MPOW7.MPOVAL, MPOW8.MPOVAL, MPOW9.MPOVAL
  159. C
  160. POINTEUR MPOV1.MPOVAL, MPOV2.MPOVAL, MPOV3.MPOVAL,
  161. & MPOV4.MPOVAL, MPOV5.MPOVAL, MPOV6.MPOVAL,
  162. & MPOV7.MPOVAL, MPOV8.MPOVAL, MPOV9.MPOVAL
  163. -INC SMELEME
  164. C
  165. C
  166. C**** Initialisation des variables pour la gestion des erreurs pas ici,
  167. C mais avant, i.e.
  168. C
  169. LOGNEG = .FALSE.
  170. LOGBOR = .FALSE.
  171. MESERR(1) = ' '
  172. MESERR(2) = ' '
  173. VALER(1) = 0.0D0
  174. VALER(2) = 0.0D0
  175. VAL1 = 0.0D0
  176. VAL2 = 0.0D0
  177.  
  178. C
  179. C**** Activation du MELEME "CENTRE"
  180. C
  181. IPT1 = IM1
  182. SEGACT IPT1
  183. N1 = IPT1.NUM(/2)
  184. SEGDES IPT1
  185. C
  186. C**** Reading MPOVA of each MCHAMPOIN
  187. C MPOW1 = alpha*rv
  188. C MPOW2 = (1 - alpha)*rl
  189. C MPOW3 = alpha*rv*uvx & alpha*rv*uvy & alpha*rv*uvz
  190. C MPOW4 = (1 - alpha)*rl*ulx & (1 - alpha)*rl*uly & (1 - alpha)*rl*ulz
  191. C MPOW5 = alpha*rv*Ev
  192. C MPOW6 = (1 - alpha)*rl*El
  193. C
  194. C MPOW7: It is the old void fraction, we will use it when
  195. C alpha approaches to zero
  196. C MPOW8: It is the old vapour temperature, we will use it when
  197. C alpha approaches to zero
  198. C MPOW9: It is the old liquid temperature, we will use it when
  199. C alpha approaches to one
  200. C
  201. CALL LICHT(ICH1,MPOW1,TYPE,IGEOMC)
  202. CALL LICHT(ICH2,MPOW2,TYPE,IGEOMC)
  203. CALL LICHT(ICH3,MPOW3,TYPE,IGEOMC)
  204. CALL LICHT(ICH4,MPOW4,TYPE,IGEOMC)
  205. CALL LICHT(ICH5,MPOW5,TYPE,IGEOMC)
  206. CALL LICHT(ICH6,MPOW6,TYPE,IGEOMC)
  207. CALL LICHT(ICH7,MPOW7,TYPE,IGEOMC)
  208. CALL LICHT(ICH8,MPOW8,TYPE,IGEOMC)
  209. CALL LICHT(ICH9,MPOW9,TYPE,IGEOMC)
  210. C
  211. C SEGACT MPOW1
  212. C SEGACT MPOW2
  213. C SEGACT MPOW3
  214. C SEGACT MPOW4
  215. C SEGACT MPOW5
  216. C SEGACT MPOW6
  217. C SEGACT MPOW7
  218. C SEGACT MPOW8
  219. C SEGACT MPOW9
  220. CC
  221. C**** LICHT active les MPOVALs en *MO
  222. C
  223. C
  224. C**** Privitive variables
  225. C
  226. C void fraction: alpha
  227.  
  228. NAT = 1
  229. NSOUPO=1
  230. N = N1
  231. NC = 1
  232. SEGINI MCHPO1
  233. MCHPO1.JATTRI(1)=2
  234. MCHPO1.IFOPOI=IFOMOD
  235. SEGINI MSOUP1
  236. MCHPO1.IPCHP(1)=MSOUP1
  237. SEGINI MPOV1
  238. MSOUP1.NOCOMP(1)='SCAL'
  239. MSOUP1.IGEOC = IM1
  240. MSOUP1.IPOVAL = MPOV1
  241.  
  242. C vapour velocities: uvx, uvy, uvz
  243.  
  244.  
  245. NC = IDIM
  246. SEGINI MCHPO2
  247. MCHPO2.JATTRI(1)=2
  248. MCHPO2.IFOPOI=IFOMOD
  249. SEGINI MSOUP2
  250. MCHPO2.IPCHP(1)=MSOUP2
  251. SEGINI MPOV2
  252. MSOUP2.NOCOMP(1)='UVX'
  253. MSOUP2.NOCOMP(2)='UVY'
  254. IF(IDIM .EQ. 3)THEN
  255. MSOUP2.NOCOMP(3)='UVZ'
  256. ENDIF
  257. MSOUP2.IGEOC = IM1
  258. MSOUP2.IPOVAL = MPOV2
  259.  
  260. C liquid velocities: ulx, uly, ulz
  261.  
  262.  
  263. NC = IDIM
  264. SEGINI MCHPO3
  265. MCHPO3.JATTRI(1)=2
  266. MCHPO3.IFOPOI=IFOMOD
  267. SEGINI MSOUP3
  268. MCHPO3.IPCHP(1)=MSOUP3
  269. SEGINI MPOV3
  270. MSOUP3.NOCOMP(1)='ULX'
  271. MSOUP3.NOCOMP(2)='ULY'
  272. IF(IDIM .EQ. 3)THEN
  273. MSOUP3.NOCOMP(3)='ULZ'
  274. ENDIF
  275. MSOUP3.IGEOC = IM1
  276. MSOUP3.IPOVAL = MPOV3
  277.  
  278. C pressure: p
  279.  
  280. NC = 1
  281. SEGINI MCHPO4
  282. MCHPO4.JATTRI(1)=2
  283. MCHPO4.IFOPOI=IFOMOD
  284. SEGINI MSOUP4
  285. MCHPO4.IPCHP(1)=MSOUP4
  286. SEGINI MPOV4
  287. MSOUP4.NOCOMP(1)='SCAL'
  288. MSOUP4.IGEOC = IM1
  289. MSOUP4.IPOVAL = MPOV4
  290.  
  291. C vapour temperature: Tv
  292.  
  293. NC = 1
  294. SEGINI MCHPO5
  295. MCHPO5.JATTRI(1)=2
  296. MCHPO5.IFOPOI=IFOMOD
  297. SEGINI MSOUP5
  298. MCHPO5.IPCHP(1)=MSOUP5
  299. SEGINI MPOV5
  300. MSOUP5.NOCOMP(1)='SCAL'
  301. MSOUP5.IGEOC = IM1
  302. MSOUP5.IPOVAL = MPOV5
  303.  
  304. C liquid temperature: Tl
  305.  
  306. NC = 1
  307. SEGINI MCHPO6
  308. MCHPO6.JATTRI(1)=2
  309. MCHPO6.IFOPOI=IFOMOD
  310. SEGINI MSOUP6
  311. MCHPO6.IPCHP(1)=MSOUP6
  312. SEGINI MPOV6
  313. MSOUP6.NOCOMP(1)='SCAL'
  314. MSOUP6.IGEOC = IM1
  315. MSOUP6.IPOVAL = MPOV6
  316.  
  317. C vapour density: rv
  318.  
  319. NC = 1
  320. SEGINI MCHPO7
  321. MCHPO7.JATTRI(1)=2
  322. MCHPO7.IFOPOI=IFOMOD
  323. SEGINI MSOUP7
  324. MCHPO7.IPCHP(1)=MSOUP7
  325. SEGINI MPOV7
  326. MSOUP7.NOCOMP(1)='SCAL'
  327. MSOUP7.IGEOC = IM1
  328. MSOUP7.IPOVAL = MPOV7
  329.  
  330. C liquid density: rl
  331.  
  332. NC = 1
  333. SEGINI MCHPO8
  334. MCHPO8.JATTRI(1)=2
  335. MCHPO8.IFOPOI=IFOMOD
  336. SEGINI MSOUP8
  337. MCHPO8.IPCHP(1)=MSOUP8
  338. SEGINI MPOV8
  339. MSOUP8.NOCOMP(1)='SCAL'
  340. MSOUP8.IGEOC = IM1
  341. MSOUP8.IPOVAL = MPOV8
  342.  
  343. C pressure correction term, pi
  344.  
  345. NC = 1
  346. SEGINI MCHPO9
  347. MCHPO9.JATTRI(1)=2
  348. MCHPO9.IFOPOI=IFOMOD
  349. SEGINI MSOUP9
  350. MCHPO9.IPCHP(1)=MSOUP9
  351. SEGINI MPOV9
  352. MSOUP9.NOCOMP(1)='SCAL'
  353. MSOUP9.IGEOC = IM1
  354. MSOUP9.IPOVAL = MPOV9
  355. C
  356. C**** "Recapitulando"
  357. C
  358. C We activate MPOV1 - MPOV9
  359. C
  360. C**** "bucle" for the calculation of the primitives at
  361. C each centre
  362.  
  363.  
  364. DO NLCE = 1, N1
  365.  
  366. oldav = MPOW7.VPOCHA(NLCE,1)
  367. oldTv = MPOW8.VPOCHA(NLCE,1)
  368.  
  369. IF(IDIM .EQ. 3)THEN
  370. if(oldav .LT. 0.01d0)then
  371. A = (gamv - 1.D0)*cpv*oldTv*MPOW1.VPOCHA(NLCE,1)/gamv
  372. else
  373. A = (gamv - 1.D0)*(MPOW5.VPOCHA(NLCE,1) - (MPOW3
  374. $ .VPOCHA(NLCE,1)**2 + MPOW3.VPOCHA(NLCE,2)**2 +
  375. & MPOW3.VPOCHA(NLCE,3)**2)/
  376. & (2.d0*(MPOW1.VPOCHA(NLCE,1)+ 1.D-20)))
  377. endif
  378. if(oldav .GT. 0.99d0)then
  379. B = (gaml - 1.D0)*MPOW2.VPOCHA(NLCE,1)*
  380. & (cpl*MPOW9.VPOCHA(NLCE,1)/gaml + pil/1000.d0)
  381. else
  382. B = (gaml - 1.D0)*(MPOW6.VPOCHA(NLCE,1) - (MPOW4
  383. $ .VPOCHA(NLCE,1)**2 + MPOW4.VPOCHA(NLCE,2)**2)+
  384. & MPOW4.VPOCHA(NLCE,3)**2/(2.D0
  385. $ *(MPOW2.VPOCHA(NLCE,1)+ 1.D-20)))
  386. endif
  387. ELSE
  388. if(oldav .LT. 0.01d0)then
  389. A = (gamv - 1.D0)*cpv*oldTv*MPOW1.VPOCHA(NLCE,1)/gamv
  390. else
  391. A = (gamv - 1.D0)*(MPOW5.VPOCHA(NLCE,1) - (MPOW3
  392. $ .VPOCHA(NLCE,1)**2 + MPOW3.VPOCHA(NLCE,2)**2)/(2.d0
  393. $ *(MPOW1.VPOCHA(NLCE,1)+ 1.D-20)))
  394. endif
  395. if(oldav .GT. 0.99d0)then
  396. B = (gaml - 1.D0)*MPOW2.VPOCHA(NLCE,1)*
  397. & (cpl*MPOW9.VPOCHA(NLCE,1)/gaml + pil/1000.d0)
  398. else
  399.  
  400. B = (gaml - 1.D0)*(MPOW6.VPOCHA(NLCE,1) - (MPOW4
  401. $ .VPOCHA(NLCE,1)**2 + MPOW4.VPOCHA(NLCE,2)**2)/(2.D0
  402. $ *(MPOW2.VPOCHA(NLCE,1)+ 1.D-20)))
  403. endif
  404. ENDIF
  405. delta = (- A - B + gaml*pil)**2.D0 + 4.D0*gaml*pil*A
  406.  
  407. p = 1.D0/2.D0*(A + B - gaml*pil + sqrt(delta))
  408.  
  409. alpha = A/p
  410. call FUNCF(alpha, FFV)
  411. uvx = FFV*MPOW3.VPOCHA(NLCE,1)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
  412. uvy = FFV*MPOW3.VPOCHA(NLCE,2)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
  413. IF(IDIM .EQ. 3)THEN
  414. uvz = FFV*MPOW3.VPOCHA(NLCE,3)/(MPOW1.VPOCHA(NLCE,1) + 1.D
  415. & -20)
  416. uv2 = uvx*uvx + uvy*uvy + uvz*uvz
  417. ELSE
  418. uv2 = uvx*uvx + uvy*uvy
  419. ENDIF
  420.  
  421. call FUNCF((1.d0 - alpha),FFL)
  422. ulx = FFL*MPOW4.VPOCHA(NLCE,1)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
  423. uly = FFL*MPOW4.VPOCHA(NLCE,2)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
  424. IF(IDIM .EQ. 3)THEN
  425. ulz = FFL*MPOW4.VPOCHA(NLCE,3)/(MPOW2.VPOCHA(NLCE,1) + 1.D
  426. & -20)
  427. ul2 = ulx*ulx + uly*uly + ulz*ulz
  428. ELSE
  429. ul2 = ulx*ulx + uly*uly
  430. ENDIF
  431.  
  432. if((alpha .GT. 0.999d0) .OR. (alpha .LT. 0.01d0))then
  433. C
  434. C Water faucet test
  435. C Tv = 323.15D0
  436. C Tl = 323.15D0
  437. C
  438. C Shock tube test
  439. C Tv = 308.15d0
  440. C Tl = 308.15d0
  441. C
  442. C Phase separation test or slosing of a water column
  443. Tv = 298.15D0
  444. Tl = 298.15D0
  445. C For the oscillating manometer should be 323 but to be practical I have left it in 298.15
  446. ev = cpv*Tv/gamv
  447. el = cpl*Tl/gaml*(p + gaml*pil)/(p + pil)
  448. else
  449. ev = MPOW5.VPOCHA(NLCE,1)/(MPOW1.VPOCHA(NLCE,1) + 1.D-20)
  450. & - uv2/2.D0
  451. el = MPOW6.VPOCHA(NLCE,1)/(MPOW2.VPOCHA(NLCE,1) + 1.D-20)
  452. & - ul2/2.D0
  453. Tv = gamv*ev/cpv
  454. Tl = gaml*el/(cpl*(1.D0 + pil*(gaml - 1.D0)/(p + pil)))
  455. endif
  456. rv = gamv*p/((gamv - 1.D0)*cpv*Tv)
  457. rl = gaml*(p + pil)/((gaml - 1.D0)*cpl*Tl)
  458. C
  459. C Bestion´s formula
  460. C
  461. rm = alpha*rv + (1.d0 - alpha)*rl
  462. rvm = rv + Cvm*rm
  463. rlm = rl + Cvm*rm
  464.  
  465. pi = p - Cp*(alpha*rvm*(1.D0- alpha)*rlm)/
  466. & (alpha*rlm + (1.D0- alpha)*rvm)*
  467. & (sqrt(uv2) - sqrt(ul2))**2
  468.  
  469. MPOV1.VPOCHA(NLCE,1) = alpha
  470. MPOV2.VPOCHA(NLCE,1) = uvx
  471. MPOV2.VPOCHA(NLCE,2) = uvy
  472. MPOV3.VPOCHA(NLCE,1) = ulx
  473. MPOV3.VPOCHA(NLCE,2) = uly
  474. IF(IDIM .EQ. 3)THEN
  475. MPOV2.VPOCHA(NLCE,3) = uvz
  476. MPOV3.VPOCHA(NLCE,3) = ulz
  477. ENDIF
  478. MPOV4.VPOCHA(NLCE,1) = p
  479. MPOV5.VPOCHA(NLCE,1) = Tv
  480. MPOV6.VPOCHA(NLCE,1) = Tl
  481. MPOV7.VPOCHA(NLCE,1) = rv
  482. MPOV8.VPOCHA(NLCE,1) = rl
  483. MPOV9.VPOCHA(NLCE,1) = pi
  484.  
  485. ENDDO
  486. C
  487. OCH1 = MCHPO1
  488. OCH2 = MCHPO2
  489. OCH3 = MCHPO3
  490. OCH4 = MCHPO4
  491. OCH5 = MCHPO5
  492. OCH6 = MCHPO6
  493. OCH7 = MCHPO7
  494. OCH8 = MCHPO8
  495. OCH9 = MCHPO9
  496. C
  497. SEGDES MPOV1
  498. SEGDES MPOV2
  499. SEGDES MPOV3
  500. SEGDES MPOV4
  501. SEGDES MPOV5
  502. SEGDES MPOV6
  503. SEGDES MPOV7
  504. SEGDES MPOV8
  505. SEGDES MPOV9
  506.  
  507. SEGDES MPOW1
  508. SEGDES MPOW2
  509. SEGDES MPOW3
  510. SEGDES MPOW4
  511. SEGDES MPOW5
  512. SEGDES MPOW6
  513. SEGDES MPOW7
  514. SEGDES MPOW8
  515. SEGDES MPOW9
  516.  
  517. SEGDES MSOUP1
  518. SEGDES MSOUP2
  519. SEGDES MSOUP3
  520. SEGDES MSOUP4
  521. SEGDES MSOUP5
  522. SEGDES MSOUP6
  523. SEGDES MSOUP7
  524. SEGDES MSOUP8
  525. SEGDES MSOUP9
  526.  
  527. SEGDES MCHPO1
  528. SEGDES MCHPO2
  529. SEGDES MCHPO3
  530. SEGDES MCHPO4
  531. SEGDES MCHPO5
  532. SEGDES MCHPO6
  533. SEGDES MCHPO7
  534. SEGDES MCHPO8
  535. SEGDES MCHPO9
  536.  
  537. RETURN
  538. END
  539.  
  540.  
  541.  
  542.  
  543.  
  544.  
  545.  
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  
  553.  
  554.  

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