Télécharger pre52f.eso

Retour à la liste

Numérotation des lignes :

  1. C PRE52F SOURCE CB215821 16/04/21 21:18:04 8920
  2. SUBROUTINE PRE52F(LOGTEM,
  3. & ICEN,IFACE,IFACEL,INORM,
  4. & IALPH, IGRALP, IALALP,
  5. & IUVC, IGRUVC, IALUVC,
  6. & IULC, IGRULC, IALULC,
  7. & IPC, IGRPC, IALPC,
  8. & ITVC, IGRTVC, IALTVC,
  9. & ITLC, IGRTLC, IALTLC,
  10. & IRVC, IRLC,
  11. & DELTAT,
  12. & IALPHF, IUVF, IULF, IPF, ITVF, ITLF,
  13. & IRVF, IRLF,
  14. & LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
  15. C************************************************************************
  16. C
  17. C PROJET : CASTEM 2000
  18. C
  19. C NOM : PRE121
  20. C
  21. C DESCRIPTION : Voir PRE42F
  22. C
  23. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
  24. C
  25. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
  26. C Modified for two-fluid flow by
  27. C Jose R. Garcia-Cascales
  28. C
  29. C************************************************************************
  30. C
  31. C APPELES (Outils) : KRIPAD, LICHT
  32. C
  33. C APPELES (Calcul) : AUCUN
  34. C
  35. C************************************************************************
  36. C
  37. C ENTREES
  38. C
  39. C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps
  40. C sinon 1er ordre en temps
  41. C
  42. C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE
  43. C
  44. C ICEN : MELEME de 'POI1' SPG des CENTRES
  45. C
  46. C IFACE : MELEME de 'POI1' SPG des FACES
  47. C
  48. C IFACEL : MELEME de 'SEG3' avec
  49. C CENTRE d'Elt "gauche"
  50. C CENTRE de Face
  51. C CENTRE d'Elt "droite"
  52. C
  53. C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2)
  54. C
  55. C INORM : CHPOINT des cosinus directeurs de normales aux faces
  56. C
  57. C 2) Pointeurs des CHPOINTs
  58. C
  59. C IALPH : CHPOINT "CENTRE" containing void fraction
  60. C
  61. C IGRALP : CHPOINT "CENTRE" contaninig the void fraction gradient
  62. C (2 composantes)
  63. C
  64. C IALALP : CHPOINT "CENTRE" contaning the LIMITA of the
  65. C void fraction gradient
  66. C
  67. C IUVC : CHPOINT "CENTRE" containing the vapour velocity UvX, UVY ;
  68. C
  69. C IGRUVC : CHPOINT "CENTRE" containig the vapour velocity gradient
  70. C (4 composantes)
  71. C
  72. C IALUVC : CHPOINT "CENTRE" contaning the LIMITA of the vapour
  73. C velocity (2 composantes)
  74. C
  75. C IULC : CHPOINT "CENTRE" containing the liquid velocity UvX, UVY ;
  76. C
  77. C IGRULC : CHPOINT "CENTRE" containig the liquid velocity gradient
  78. C (4 composantes)
  79. C
  80. C IALULC : CHPOINT "CENTRE" contaning the LIMITA of the liquid
  81. C velocity (2 composantes)
  82. C
  83. C IPC : CHPOINT "CENTRE" containing the pressure P ;
  84. C
  85. C IGRPC : CHPOINT "CENTRE" containing the gradient of pressure
  86. C (2 composantes)
  87. C
  88. C IALPC : CHPOINT "CENTRE" containig the LIMITA of pressure
  89. C gradient
  90. C
  91. C ITVC : CHPOINT "CENTRE" containing the vapour temperature Tv ;
  92. C
  93. C IGRTVC : CHPOINT "CENTRE" containing the gradient of vapour
  94. C temperature (2 composantes)
  95. C
  96. C IALTVC : CHPOINT "CENTRE" containig the LIMITA of vapour
  97. C temperature gradient
  98. C
  99. C ITLC : CHPOINT "CENTRE" containing the liquid temperature Tl ;
  100. C
  101. C IGRTLC : CHPOINT "CENTRE" containing the gradient of liquid
  102. C temperature (2 composantes)
  103. C
  104. C IALTLC : CHPOINT "CENTRE" containig the LIMITA of liquid
  105. C temperature gradient
  106. C
  107. C IRVC : CHPOINT "CENTRE" containing the vapour density
  108. C
  109. C IRLC : CHPOINT "CENTRE" containing the liquid density
  110. C
  111. C
  112. C DELTAT : REAL*8, encrement en temps pour calculer la prediction
  113. C
  114. C SORTIES
  115. C
  116. C C SORTIES
  117. C
  118. C IALPHF : MCHAML defined en la MELEME of pointers
  119. C IFACEL, it contains the void fraction
  120. C (a gauche et a droite de chaque face).
  121. C Only one component ('SCAL')
  122. C
  123. C IUVF : MCHAML "FACEL" vapour velocity and the
  124. C director cosines (n,t) in the corresponding face;
  125. C in the 2D case 6 composantes:
  126. C 'UVN' = normal velocity
  127. C 'UVT' = tangent velocity,
  128. C these two velocities are defined in a local
  129. C reference framework defined over the MELEME
  130. C of pointers IFACE
  131. C 'NX' = n.x
  132. C 'NY' = n.y
  133. C 'TX' = t.x
  134. C 'TY' = t.y
  135. C
  136. C IULF : MCHAML "FACEL" liquid velocity
  137. C in the 2D case 2 composantes:
  138. C 'ULN' = normal velocity
  139. C 'ULT' = tangent velocity
  140. C
  141. C IPF : MCHAML "FACEL" pressure
  142. C Only one component ('SCAL')
  143. C
  144. C ITVF : MCHAML "FACEL" vapour temperature
  145. C Only one component ('SCAL')
  146. C
  147. C ITVL : MCHAML "FACEL" liquid temperature
  148. C Only one component ('SCAL')
  149. C
  150. C IRVF : MCHAML "FACEL" vapour density
  151. C Only one component ('SCAL')
  152. C
  153. C IRLF : MCHAML "FACEL" liquid density
  154. C Only one component ('SCAL')
  155. C
  156. C LOGAN : anomalie detectee (changement de la convention dans
  157. C la table domaine)
  158. C
  159. C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité
  160. C negative a été detectée -> en interactif le
  161. C programme s'arrete en GIBIANE
  162. C (erreur stocké en MESERR et VALER)
  163. C
  164. C LOGBOR : (LOGICAL): si .TRUE. un gamma a ete detecte
  165. C dehor 1 et 3 (sa valeur stockée en MESERR et VALER;
  166. C en VAL1 et en VAL2 on stocke 1.0 et 3.0)
  167. C
  168. C MESERR
  169. C VALER
  170. C VAL1,
  171. C VAL2 : pour les messages d'erreur
  172. C
  173. C************************************************************************
  174. C
  175. C HISTORIQUE (Anomalies et modifications éventuelles)
  176. C
  177. C HISTORIQUE : Créée le 11.6.98.
  178. C
  179. C************************************************************************
  180. C
  181. C
  182. C ATTENTION: Cet programme marche si le MAILLAGE est convex;
  183. C si non il faut changer l'argoritme de calcul de
  184. C l'orientation des normales aux faces.
  185. C
  186. C La positivité n'est pas controlle parce que c'est déjà fait
  187. C dans l'operateur PRIM
  188. C
  189. C
  190. C************************************************************************
  191. C
  192. IMPLICIT INTEGER(I-N)
  193. IMPLICIT REAL*8(A-H,O-Z)
  194.  
  195. C
  196. C**** Les variables
  197. C
  198. INTEGER ICEN, IFACE, IFACEL, INORM,
  199. & IALPH, IGRALP, IALALP,
  200. & IUVC, IGRUVC, IALUVC,
  201. & IULC, IGRULC, IALULC,
  202. & IPC, IGRPC, IALPC,
  203. & ITVC, IGRTVC, IALTVC,
  204. & ITLC, IGRTLC, IALTLC,
  205. & IRVC, IRLC,
  206. & IALPHF, IUVF, IULF, IPF, ITVF, ITLF,
  207. & IRVF, IRLF,
  208. & IGEOM, NFAC, NCEN,
  209. & N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE,
  210. & NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1,
  211. & IDIMP1, INDCEL,
  212. & l, m, n
  213. REAL*8 VALER, VAL1, VAL2, XG, YG, XC, YC, XD, YD, DELTAT
  214. & , DXG, DYG, DXD, DYD
  215. & , CNX, CNY, CTX, CTY, ORIENT
  216. & , alphag, uvxg, uvyg, uvng, uvtg, ulxg, ulyg, ulng, ultg
  217. & , pg, Tvg, Tlg, rvg, rlg
  218. & , alphad, uvxd, uvyd, uvnd, uvtd, ulxd, ulyd, ulnd, ultd
  219. & , pd, Tvd, Tld, rvd, rld
  220. & , VALCEL, DCEL, ALCEL
  221. & , alpha, uvx, uvy, ulx, uly, p, Tv, Tl
  222. & , rv, rl, drvdp, drvdhv, drldp, drldhl
  223. & , hv, hl, cv, cl
  224. & , A(8,8), B(8,8), denom
  225. & , dVdx(8), dVdy(8), dV(8), LIMITA(8), dVdVp(8,8)
  226. & , dVpdV(8,8), P1(8,8), P2(8,8)
  227. & , dhvdp, dhvdTv, dhldp, dhldTl
  228. CHARACTER*(40) MESERR
  229. CHARACTER*(8) TYPE
  230. LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGALP
  231. C
  232. C**** perfect and stiffened gas eos
  233. C
  234. real*8 gamv, cpv, gaml, cpl, pil
  235. PARAMETER(gamv = 1.4D0, cpv = 1008.D0,
  236. & gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8)
  237. C
  238. C**** LOGALP = .TRUE. -> Prediction avec limiteur
  239. C
  240. PARAMETER(LOGALP = .TRUE.)
  241. C
  242. C**** Les Includes
  243. C
  244. -INC SMCOORD
  245. -INC CCOPTIO
  246. -INC SMCHPOI
  247. POINTEUR MPALPC.MPOVAL, MPGRAL.MPOVAL,
  248. & MPUVC.MPOVAL, MPGRUV.MPOVAL,
  249. & MPULC.MPOVAL, MPGRUL.MPOVAL,
  250. & MPPC.MPOVAL, MPGRP.MPOVAL,
  251. & MPTVC.MPOVAL, MPGRTV.MPOVAL,
  252. & MPTLC.MPOVAL, MPGRTL.MPOVAL,
  253. & MPRVC.MPOVAL, MPRLC.MPOVAL,
  254. & MPNORM.MPOVAL,
  255. & MPALPP.MPOVAL, MPUVP.MPOVAL, MPULP.MPOVAL,
  256. & MPPP.MPOVAL, MPTVP.MPOVAL, MPTLP.MPOVAL,
  257. & MPRVP.MPOVAL, MPRLP.MPOVAL
  258. -INC SMCHAML
  259. POINTEUR MLALP.MELVAL, MLP.MELVAL,
  260. & MLTV.MELVAL, MLTL.MELVAL,
  261. & MLRV.MELVAL, MLRL.MELVAL
  262. POINTEUR MLUVN.MELVAL, MLUVT.MELVAL,
  263. & MLULN.MELVAL, MLULT.MELVAL,
  264. & MLVNX.MELVAL, MLVNY.MELVAL, MLVTX.MELVAL, MLVTY.MELVAL,
  265. & MLLNX.MELVAL, MLLNY.MELVAL, MLLTX.MELVAL, MLLTY.MELVAL
  266. C**** Some melemes for the LIMITAs
  267. POINTEUR MELLAL.MPOVAL,
  268. & MELLUV.MPOVAL,
  269. & MELLUL.MPOVAL,
  270. & MELLP.MPOVAL,
  271. & MELLTV.MPOVAL,
  272. & MELLTL.MPOVAL
  273. -INC SMLENTI
  274. -INC SMELEME
  275. C
  276. C
  277. C**** Initialisation des parametres d'erreur déjà faite, i.e.
  278. C
  279. C LOGNEG = .FALSE.
  280. C LOGBOR = .FALSE.
  281. C MESERR = ' '
  282. C MOTERR(1:40) = MESERR(1:40)
  283. C VALER = 0.0D0
  284. C VAL1 = 0.0D0
  285. C VAL2 = 0.0D0
  286. C
  287. C
  288. C**** KRIPAD pour la correspondance global/local de centre
  289. C
  290. CALL KRIPAD(ICEN,MLENT1)
  291. C
  292. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  293. C
  294. C Si i est le numero global d'un noeud de ICEN,
  295. C MLENT1.LECT(i) contient sa position, i.e.
  296. C
  297. C I = numero global du noeud centre
  298. C MLENT1.LECT(i) = numero local du noeud centre
  299. C
  300. C MLENT1 déjà activé, i.e.
  301. C
  302. C SEGACT MLENT1
  303. C
  304. C**** Activation de CHPOINTs
  305. C
  306. C void fraction + grad + limiter
  307. C vapour velocity + grad + limiter
  308. C liquid velocity + grad + limiter
  309. C pressure + grad + limiter
  310. C vapour temperature + grad + limiter
  311. C liquid temperature + grad + limiter
  312. C vapour density
  313. C liquid density
  314. C
  315. C cosinus directeurs des normales aux surface
  316. C
  317. CALL LICHT(IALPH, MPALPC, TYPE, IGEOM)
  318. CALL LICHT(IGRALP, MPGRAL, TYPE, IGEOM)
  319. CALL LICHT(IUVC, MPUVC, TYPE, IGEOM)
  320. CALL LICHT(IGRUVC, MPGRUV, TYPE, IGEOM)
  321. CALL LICHT(IULC, MPULC, TYPE, IGEOM)
  322. CALL LICHT(IGRULC, MPGRUL, TYPE, IGEOM)
  323. CALL LICHT(IPC , MPPC , TYPE, IGEOM)
  324. CALL LICHT(IGRPC, MPGRP, TYPE, IGEOM)
  325. CALL LICHT(ITVC, MPTVC, TYPE, IGEOM)
  326. CALL LICHT(IGRTVC, MPGRTV, TYPE, IGEOM)
  327. CALL LICHT(ITLC, MPTLC, TYPE, IGEOM)
  328. CALL LICHT(IGRTLC, MPGRTL, TYPE, IGEOM)
  329. CALL LICHT(IRVC, MPRVC, TYPE, IGEOM)
  330. CALL LICHT(IRLC, MPRLC, TYPE, IGEOM)
  331. CALL LICHT(INORM, MPNORM, TYPE, IGEOM)
  332. C
  333. C**** Les MPOVALs 'Prediction'
  334. C
  335. IF(LOGTEM)THEN
  336. SEGINI, MPALPP = MPALPC
  337. SEGINI, MPUVP = MPUVC
  338. SEGINI, MPULP = MPULC
  339. SEGINI, MPPP = MPPC
  340. SEGINI, MPTVP = MPTVC
  341. SEGINI, MPTLP = MPTLC
  342. SEGINI, MPRVP = MPRVC
  343. SEGINI, MPRLP = MPRLC
  344. ELSE
  345. MPALPP = MPALPC
  346. MPUVP = MPUVC
  347. MPULP = MPULC
  348. MPPP = MPPC
  349. MPTVP = MPTVC
  350. MPTLP = MPTLC
  351. MPRVP = MPRVC
  352. MPRLP = MPRLC
  353. ENDIF
  354. C
  355. C**** Les Limiteurs
  356. C
  357. CALL LICHT(IALALP, MELLAL, TYPE, IGEOM)
  358. CALL LICHT(IALUVC, MELLUV, TYPE, IGEOM)
  359. CALL LICHT(IALULC, MELLUL, TYPE, IGEOM)
  360. CALL LICHT(IALPC , MELLP, TYPE, IGEOM)
  361. CALL LICHT(IALTVC, MELLTV, TYPE, IGEOM)
  362. CALL LICHT(IALTLC, MELLTL, TYPE, IGEOM)
  363. C
  364. C
  365. C**** MPOVAL sont déjà activés i.e.:
  366. C
  367. C SEGACT MPALPP
  368. C SEGACT MPGRAL
  369. C SEGACT IALALP
  370. C SEGACT MPUVC
  371. C SEGACT IALUVC
  372. C SEGACT MPGRUV
  373. C SEGACT MPULC
  374. C SEGACT MPGRUL
  375. C SEGACT IALULC
  376. C SEGACT MPPC
  377. C SEGACT MPGRP
  378. C SEGACT IALPC
  379. C SEGACT MPTVC
  380. C SEGACT MPGRTV
  381. C SEGACT IALTVC
  382. C SEGACT MPTLC
  383. C SEGACT MPGRTL
  384. C SEGACT IALTLC
  385. C SEGACT MPRVC
  386. C SEGACT MPRLC
  387. C SEGACT MPNORM
  388. C
  389. C**** Le MELEME FACEL
  390. C
  391. IPT1 = IFACEL
  392. IPT2 = IFACE
  393. SEGACT IPT1
  394. SEGACT IPT2
  395. NFAC = IPT1.NUM(/2)
  396. C
  397. C**** Creation de MCHAMLs contenant les etats gauche et droite,
  398. C
  399. C i.e.:
  400. C
  401. C vitesse + cosinus directors du repere local
  402. C densité
  403. C pression
  404. C gamma
  405. C
  406. C**** Cosinus directors du repere local et vitesse
  407. C
  408. C Les cosinus directeurs
  409. C
  410. C VAPOUR PHASE
  411.  
  412. N1 = 2
  413. N3 = 6
  414. L1 = 28
  415. SEGINI MCHEL1
  416. IUVF = MCHEL1
  417. MCHEL1.TITCHE = 'UV '
  418. MCHEL1.IMACHE(1) = IFACE
  419. MCHEL1.IMACHE(2) = IFACEL
  420. MCHEL1.CONCHE(1) = ' (n,t) in (x,y) '
  421. MCHEL1.CONCHE(2) = ' UV in (n,t) '
  422. C
  423. C**** Valeurs des cosinus definies par respect au repair global, i.e.
  424. C
  425. MCHEL1.INFCHE(1,1) = 2
  426. MCHEL1.INFCHE(1,3) = NIFOUR
  427. MCHEL1.INFCHE(1,4) = 0
  428. MCHEL1.INFCHE(1,5) = 0
  429. MCHEL1.INFCHE(1,6) = 0
  430. MCHEL1.IFOCHE = IFOUR
  431. C
  432. C**** Valeurs de vitesse definies par respect au repair local, i.e.
  433. C
  434. MCHEL1.INFCHE(2,1) = 1
  435. MCHEL1.INFCHE(2,3) = NIFOUR
  436. MCHEL1.INFCHE(2,4) = 0
  437. MCHEL1.INFCHE(2,5) = 0
  438. MCHEL1.INFCHE(2,6) = 0
  439. C
  440. C**** Le cosinus directeurs
  441. C
  442. N1PTEL = 1
  443. N1EL = NFAC
  444. N2PTEL = 0
  445. N2EL = 0
  446. C
  447. C**** MCHAML a N2 composantes:
  448. C
  449. C cosinus directeurs du repere local (n,t1)
  450. C
  451. C IDIM = 2 -> 4 composantes
  452. C
  453. N2 = 4
  454. SEGINI MCHAM1
  455. MCHEL1.ICHAML(1) = MCHAM1
  456. MCHAM1.NOMCHE(1) = 'NVX '
  457. MCHAM1.NOMCHE(2) = 'NVY '
  458. MCHAM1.NOMCHE(3) = 'TVX '
  459. MCHAM1.NOMCHE(4) = 'TVY '
  460. MCHAM1.TYPCHE(1) = 'REAL*8 '
  461. MCHAM1.TYPCHE(2) = 'REAL*8 '
  462. MCHAM1.TYPCHE(3) = 'REAL*8 '
  463. MCHAM1.TYPCHE(4) = 'REAL*8 '
  464. SEGINI MLVNX
  465. SEGINI MLVNY
  466. SEGINI MLVTX
  467. SEGINI MLVTY
  468. MCHAM1.IELVAL(1) = MLVNX
  469. MCHAM1.IELVAL(2) = MLVNY
  470. MCHAM1.IELVAL(3) = MLVTX
  471. MCHAM1.IELVAL(4) = MLVTY
  472. SEGDES MCHAM1
  473. C
  474. C**** Vitesse
  475. C
  476. N1EL = NFAC
  477. N1PTEL = 3
  478. N2EL = 0
  479. N2PTEL = 0
  480. C
  481. C**** MCHAML a N2 composantes:
  482. C
  483. C IDIM = 2 -> 2 composantes
  484. C
  485. N2 = 2
  486. SEGINI MCHAM1
  487. MCHEL1.ICHAML(2) = MCHAM1
  488. SEGDES MCHEL1
  489. MCHAM1.NOMCHE(1) = 'UVN '
  490. MCHAM1.NOMCHE(2) = 'UVT '
  491. MCHAM1.TYPCHE(1) = 'REAL*8 '
  492. MCHAM1.TYPCHE(2) = 'REAL*8 '
  493. SEGINI MLUVN
  494. SEGINI MLUVT
  495. MCHAM1.IELVAL(1) = MLUVN
  496. MCHAM1.IELVAL(2) = MLUVT
  497. SEGDES MCHAM1
  498. C
  499. C**** Cosinus directors du repere local et vitesse
  500. C
  501. C Les cosinus directeurs
  502. C
  503. C LIQUID PHASE
  504.  
  505. N1 = 2
  506. N3 = 6
  507. L1 = 28
  508. SEGINI MCHEL1
  509. IULF = MCHEL1
  510. MCHEL1.TITCHE = 'UL '
  511. MCHEL1.IMACHE(1) = IFACE
  512. MCHEL1.IMACHE(2) = IFACEL
  513. MCHEL1.CONCHE(1) = ' (n,t) in (x,y) '
  514. MCHEL1.CONCHE(2) = ' UL in (n,t) '
  515. C
  516. C**** Valeurs des cosinus definies par respect au repair global, i.e.
  517. C
  518. MCHEL1.INFCHE(1,1) = 2
  519. MCHEL1.INFCHE(1,3) = NIFOUR
  520. MCHEL1.INFCHE(1,4) = 0
  521. MCHEL1.INFCHE(1,5) = 0
  522. MCHEL1.INFCHE(1,6) = 0
  523. MCHEL1.IFOCHE = IFOUR
  524. C
  525. C**** Valeurs de vitesse definies par respect au repair local, i.e.
  526. C
  527. MCHEL1.INFCHE(2,1) = 1
  528. MCHEL1.INFCHE(2,3) = NIFOUR
  529. MCHEL1.INFCHE(2,4) = 0
  530. MCHEL1.INFCHE(2,5) = 0
  531. MCHEL1.INFCHE(2,6) = 0
  532. C
  533. C**** Le cosinus directeurs
  534. C
  535. N1PTEL = 1
  536. N1EL = NFAC
  537. N2PTEL = 0
  538. N2EL = 0
  539. C
  540. C**** MCHAML a N2 composantes:
  541. C
  542. C cosinus directeurs du repere local (n,t1)
  543. C
  544. C IDIM = 2 -> 4 composantes
  545. C
  546. N2 = 4
  547. SEGINI MCHAM1
  548. MCHEL1.ICHAML(1) = MCHAM1
  549. MCHAM1.NOMCHE(1) = 'NLX '
  550. MCHAM1.NOMCHE(2) = 'NLY '
  551. MCHAM1.NOMCHE(3) = 'TLX '
  552. MCHAM1.NOMCHE(4) = 'TLY '
  553. MCHAM1.TYPCHE(1) = 'REAL*8 '
  554. MCHAM1.TYPCHE(2) = 'REAL*8 '
  555. MCHAM1.TYPCHE(3) = 'REAL*8 '
  556. MCHAM1.TYPCHE(4) = 'REAL*8 '
  557. SEGINI MLLNX
  558. SEGINI MLLNY
  559. SEGINI MLLTX
  560. SEGINI MLLTY
  561. MCHAM1.IELVAL(1) = MLLNX
  562. MCHAM1.IELVAL(2) = MLLNY
  563. MCHAM1.IELVAL(3) = MLLTX
  564. MCHAM1.IELVAL(4) = MLLTY
  565. SEGDES MCHAM1
  566. C
  567. C**** Vitesse
  568. C
  569. N1EL = NFAC
  570. N1PTEL = 3
  571. N2EL = 0
  572. N2PTEL = 0
  573. C
  574. C**** MCHAML a N2 composantes:
  575. C
  576. C IDIM = 2 -> 2 composantes
  577. C
  578. N2 = 2
  579. SEGINI MCHAM1
  580. MCHEL1.ICHAML(2) = MCHAM1
  581. SEGDES MCHEL1
  582. MCHAM1.NOMCHE(1) = 'ULN '
  583. MCHAM1.NOMCHE(2) = 'ULT '
  584. MCHAM1.TYPCHE(1) = 'REAL*8 '
  585. MCHAM1.TYPCHE(2) = 'REAL*8 '
  586. SEGINI MLULN
  587. SEGINI MLULT
  588. MCHAM1.IELVAL(1) = MLULN
  589. MCHAM1.IELVAL(2) = MLULT
  590. SEGDES MCHAM1
  591. C
  592. C**** Void fraction
  593. C
  594. N1 = 1
  595. N3 = 6
  596. L1 = 15
  597. SEGINI MCHEL2
  598. IALPHF = MCHEL2
  599. MCHEL2.IMACHE(1) = IFACEL
  600. MCHEL2.TITCHE = 'ALPHA '
  601. MCHEL2.CONCHE(1) = ' '
  602. C
  603. C**** Valeurs independente du repére, i.e.
  604. C
  605. MCHEL2.INFCHE(1,1) = 0
  606. MCHEL2.INFCHE(1,3) = NIFOUR
  607. MCHEL2.INFCHE(1,4) = 0
  608. MCHEL2.INFCHE(1,5) = 0
  609. MCHEL2.INFCHE(1,6) = 0
  610. MCHEL2.IFOCHE = IFOUR
  611. N2 = 1
  612. SEGINI MCHAM1
  613. MCHEL2.ICHAML(1) = MCHAM1
  614. SEGDES MCHEL2
  615. MCHAM1.NOMCHE(1) = 'SCAL '
  616. MCHAM1.TYPCHE(1) = 'REAL*8 '
  617. SEGINI MLALP
  618. MCHAM1.IELVAL(1) = MLALP
  619. SEGDES MCHAM1
  620. C
  621. C**** Pressure
  622. C
  623. MCHEL1 = IALPHF
  624. SEGINI, MCHEL2 = MCHEL1
  625. IPF = MCHEL2
  626. MCHEL2.TITCHE = 'P '
  627. C
  628. C**** MCHAM1 = MCHAML de la alpha
  629. C
  630. SEGINI, MCHAM2 = MCHAM1
  631. MCHEL2.ICHAML(1) = MCHAM2
  632. SEGDES MCHEL2
  633. SEGINI MLP
  634. MCHAM2.IELVAL(1) = MLP
  635. SEGDES MCHAM2
  636. C
  637. C**** Vapour temperature
  638. C
  639. MCHEL1 = IALPHF
  640. SEGINI, MCHEL2 = MCHEL1
  641. ITVF = MCHEL2
  642. MCHEL2.TITCHE = 'TV '
  643. C
  644. C**** MCHAM1 = MCHAML de la alpha
  645. C
  646. SEGINI, MCHAM2 = MCHAM1
  647. MCHEL2.ICHAML(1) = MCHAM2
  648. SEGDES MCHEL2
  649. SEGINI MLTV
  650. MCHAM2.IELVAL(1) = MLTV
  651. SEGDES MCHAM2
  652. C
  653. C**** Liquid temperature
  654. C
  655. MCHEL1 = IALPHF
  656. SEGINI, MCHEL2 = MCHEL1
  657. ITLF = MCHEL2
  658. MCHEL2.TITCHE = 'TL '
  659. C
  660. C**** MCHAM1 = MCHAML de la alpha
  661. C
  662. SEGINI, MCHAM2 = MCHAM1
  663. MCHEL2.ICHAML(1) = MCHAM2
  664. SEGDES MCHEL2
  665. SEGINI MLTL
  666. MCHAM2.IELVAL(1) = MLTL
  667. SEGDES MCHAM2
  668. C
  669. C**** Vapour density
  670. C
  671. MCHEL1 = IALPHF
  672. SEGINI, MCHEL2 = MCHEL1
  673. IRVF = MCHEL2
  674. MCHEL2.TITCHE = 'RV '
  675. C
  676. C**** MCHAM1 = MCHAML de la alpha
  677. C
  678. SEGINI, MCHAM2 = MCHAM1
  679. MCHEL2.ICHAML(1) = MCHAM2
  680. SEGDES MCHEL2
  681. SEGINI MLRV
  682. MCHAM2.IELVAL(1) = MLRV
  683. SEGDES MCHAM2
  684. C
  685. C**** Liquid density
  686. C
  687. MCHEL1 = IALPHF
  688. SEGINI, MCHEL2 = MCHEL1
  689. IRLF = MCHEL2
  690. MCHEL2.TITCHE = 'RL '
  691. C
  692. C**** MCHAM1 = MCHAML de la alpha
  693. C
  694. SEGINI, MCHAM2 = MCHAM1
  695. MCHEL2.ICHAML(1) = MCHAM2
  696. SEGDES MCHEL2
  697. SEGINI MLRL
  698. MCHAM2.IELVAL(1) = MLRL
  699. SEGDES MCHAM2
  700.  
  701. C
  702. C**** Recapitulatif: les MELVALs et les MPOVALs actives
  703. C
  704. C MLVNX, MLVNY,
  705. C MLVTX, MLVTY,
  706. C
  707. C MLLNX, MLLNY,
  708. C MLLTX, MLLTY
  709. C
  710. C MLUVN, MLUVT -> vapour velocities
  711. C
  712. C MLULN, MLULT -> liquid velocities
  713. C
  714. C MLALP -> void fraction
  715. C
  716. C MLP -> pressure
  717. C
  718. C MLTV -> vapour temperature
  719. C
  720. C MLTL -> liquid temperature
  721. C
  722. C MLRV -> vapour density
  723. C
  724. C MLRL -> liquid density
  725. C****
  726. C MPALPC -> void fraction
  727. C
  728. C MPUVC -> vapour velocity
  729. C
  730. C MPULC -> liquid velocity
  731. C
  732. C MPPC -> pressure
  733. C
  734. C MPTVC -> vapour temperature
  735. C
  736. C MPTLC -> liquid temperature
  737. C
  738. C MPRVC -> vapour density
  739. C
  740. C MPRLC -> liquid density
  741. C
  742. C MPNORM -> normales aux faces
  743. C
  744. C***********************************************************************
  745. C********* PREDICTION **************************************************
  746. C***********************************************************************
  747. C
  748. Cµµµµµµµµµµµµµµµµµµµµµ AQUI ME QUEDO
  749. C
  750. IF(LOGTEM)THEN
  751. C
  752. IPT3 = ICEN
  753. SEGACT IPT3
  754. NCEN = IPT3.NUM(/2)
  755. DO NLCE = 1, NCEN
  756. IF(LOGALP)THEN
  757. C
  758. C************* Gradients limités
  759. C
  760. LIMITA(1) = MELLAL.VPOCHA(NLCE,1)
  761. LIMITA(2) = MELLUV.VPOCHA(NLCE,1)
  762. LIMITA(3) = MELLUV.VPOCHA(NLCE,2)
  763. LIMITA(4) = MELLUL.VPOCHA(NLCE,1)
  764. LIMITA(5) = MELLUL.VPOCHA(NLCE,2)
  765. LIMITA(6) = MELLP.VPOCHA(NLCE,1)
  766. LIMITA(7) = MELLTV.VPOCHA(NLCE,1)
  767. LIMITA(8) = MELLTL.VPOCHA(NLCE,1)
  768. ELSE
  769. C
  770. C************* Gradients non limités
  771. C
  772. LIMITA(1) = 1.0D0
  773. LIMITA(2) = 1.0D0
  774. LIMITA(3) = 1.0D0
  775. LIMITA(4) = 1.0D0
  776. LIMITA(5) = 1.0D0
  777. LIMITA(6) = 1.0D0
  778. LIMITA(7) = 1.0D0
  779. LIMITA(8) = 1.0D0
  780. ENDIF
  781.  
  782. alpha = MPALPP.VPOCHA(NLCE,1)
  783. uvx = MPUVP.VPOCHA(NLCE,1)
  784. uvy = MPUVP.VPOCHA(NLCE,2)
  785. ulx = MPULP.VPOCHA(NLCE,1)
  786. uly = MPULP.VPOCHA(NLCE,2)
  787. p = MPPP.VPOCHA(NLCE,1)
  788. Tv = MPTVP.VPOCHA(NLCE,1)
  789. Tl = MPTLP.VPOCHA(NLCE,1)
  790. rv = MPRVP.VPOCHA(NLCE,1)
  791. rl = MPRLP.VPOCHA(NLCE,1)
  792. C
  793. C************* Value of the enthalpies and the derivatives of
  794. C the densities respect pressure and enthalpies,
  795. C for the moment only for perfect and stiffened gases.
  796. hv = Cpv*Tv
  797. hl = Cpl*Tl
  798. drvdp = rv/p
  799. drvdhv = -rv/hv
  800. cv = 1.d0/SQRT(drvdhv/rv + drvdp)
  801. drldp = rl/(p + pil)
  802. drldhl = -rl/hl
  803. cl = 1.d0/SQRT(drldhl/rl + drldp)
  804. dhvdp = 0.d0
  805. dhvdTv = Cpv
  806. dhldp = 0.d0
  807. dhldTl = Cpl
  808.  
  809. denom = drldhl*(drvdhv*p + rv*(alpha*drvdp*p + rv -
  810. & alpha*rv)) + rl*(alpha*rl*rv*cv**(-2) -
  811. & (-1 + alpha)*drldp*(drvdhv*p + rv**2))
  812.  
  813. A(1,1) = (alpha*(drldhl*p + rl**2)*rv*cv**(-2)*ulx -
  814. & (-1 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*uvx)
  815. & /denom
  816.  
  817. A(1,2) = -((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv**2)
  818. $ /denom
  819. A(1,3) = 0.d0
  820. A(1,4) = ((-1.d0 + alpha)*alpha*rl**2*rv*cv**(-2))
  821. $ /denom
  822. A(1,5) = 0.d0
  823. A(1,6) = ((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv*cv**(-2)*(ulx
  824. $ -uvx))/denom
  825. A(1,7) = 0.d0
  826. A(1,8) = 0.d0
  827.  
  828. A(2,1) = p/(alpha*rv)
  829. A(2,2) = uvx
  830. A(2,3) = 0.d0
  831. A(2,4) = 0.d0
  832. A(2,5) = 0.d0
  833. A(2,6) = 1.d0/rv
  834. A(2,7) = 0.d0
  835. A(2,8) = 0.d0
  836.  
  837. A(3,1) = 0.d0
  838. A(3,2) = 0.d0
  839. A(3,3) = uvx
  840. A(3,4) = 0.d0
  841. A(3,5) = 0.d0
  842. A(3,6) = 0.d0
  843. A(3,7) = 0.d0
  844. A(3,8) = 0.d0
  845.  
  846. A(4,1) = p/(-rl + alpha*rl)
  847. A(4,2) = 0.d0
  848. A(4,3) = 0.d0
  849. A(4,4) = ulx
  850. A(4,5) = 0.d0
  851. A(4,6) = 1.d0/rl
  852. A(4,7) = 0.d0
  853. A(4,8) = 0.d0
  854.  
  855. A(5,1) = 0.d0
  856. A(5,2) = 0.d0
  857. A(5,3) = 0.d0
  858. A(5,4) = 0.d0
  859. A(5,5) = ulx
  860. A(5,6) = 0.d0
  861. A(5,7) = 0.d0
  862. A(5,8) = 0
  863.  
  864. A(6,1) = -((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulx - uvx)
  865. $ )/denom
  866. A(6,2) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
  867. A(6,3) = 0.d0
  868. A(6,4) = -((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))
  869. $ /denom
  870. A(6,5) = 0.d0
  871. A(6,6) = (-((-1 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*ulx)
  872. $ +alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uvx)/denom
  873. A(6,7) = 0.d0
  874. A(6,8) = 0.d0
  875.  
  876. A(7,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulx - uvx))
  877. & /denom
  878. A(7,2) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))
  879. $ *rv)/denom
  880. A(7,3) = 0.d0
  881. A(7,4) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
  882. A(7,5) = 0.d0
  883. A(7,6) = ((-1.d0 + alpha)*rl*cl**(-2)*(drvdp*p - rv)*(ulx -
  884. $ uvx))/denom
  885. A(7,7) = uvx
  886. A(7,8) = 0.d0
  887.  
  888. A(8,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulx - uvx))
  889. & /denom
  890. A(8,2) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
  891. A(8,3) = 0.d0
  892. A(8,4) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))
  893. $ )/denom
  894. A(8,5) = 0.d0
  895. A(8,6) = (alpha*(drldp*p - rl)*rv*cv**(-2)*(ulx - uvx))
  896. & /denom
  897. A(8,7) = 0.d0
  898. A(8,8) = ulx
  899.  
  900. B(1,1) = (alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uly -
  901. & (-1.d0 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*uvy)
  902. & /denom
  903. B(1,2) = 0.d0
  904. B(1,3) = -((-1 + alpha)*alpha*rl*cl**(-2)*rv**2)/denom
  905. B(1,4) = 0.d0
  906. B(1,5) = ((-1.d0 + alpha)*alpha*rl**2*rv*cv**(-2))
  907. $ /denom
  908. B(1,6) = ((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv*cv**(-2)*(uly
  909. $ -uvy))/denom
  910. B(1,7) = 0.d0
  911. B(1,8) = 0.d0
  912.  
  913. B(2,1) = 0.d0
  914. B(2,2) = uvy
  915. B(2,3) = 0.d0
  916. B(2,4) = 0.d0
  917. B(2,5) = 0.d0
  918. B(2,6) = 0.d0
  919. B(2,7) = 0.d0
  920. B(2,8) = 0.d0
  921.  
  922. B(3,1) = p/(alpha*rv)
  923. B(3,2) = 0.d0
  924. B(3,3) = uvy
  925. B(3,4) = 0.d0
  926. B(3,5) = 0.d0
  927. B(3,6) = 1.d0/rv
  928. B(3,7) = 0.d0
  929. B(3,8) = 0.d0
  930.  
  931. B(4,1) = 0.d0
  932. B(4,2) = 0.d0
  933. B(4,3) = 0.d0
  934. B(4,4) = uly
  935. B(4,5) = 0.d0
  936. B(4,6) = 0.d0
  937. B(4,7) = 0.d0
  938. B(4,8) = 0.d0
  939.  
  940. B(5,1) = p/(-rl + alpha*rl)
  941. B(5,2) = 0.d0
  942. B(5,3) = 0.d0
  943. B(5,4) = 0.d0
  944. B(5,5) = uly
  945. B(5,6) = 1.d0/rl
  946. B(5,7) = 0.d0
  947. B(5,8) = 0.d0
  948.  
  949. B(6,1) = -((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(uly - uvy)
  950. $ )/denom
  951. B(6,2) = 0.d0
  952. B(6,3) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
  953. B(6,4) = 0.d0
  954. B(6,5) = -((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))
  955. $ /denom
  956. B(6,6) = (-((-1.d0 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)
  957. $ *uly) +alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uvy)
  958. $ /denom
  959. B(6,7) = 0.d0
  960. B(6,8) = 0.d0
  961.  
  962. B(7,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(uly - uvy))
  963. & /denom
  964. B(7,2) = 0.d0
  965. B(7,3) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))
  966. $ *rv)/denom
  967. B(7,4) = 0.d0
  968. B(7,5) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
  969. B(7,6) = ((-1.d0 + alpha)*rl*cl**(-2)*(drvdp*p - rv)*(uly -
  970. $ uvy))/denom
  971. B(7,7) = uvy
  972. B(7,8) = 0.d0
  973.  
  974. B(8,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(uly - uvy))
  975. & /denom
  976. B(8,2) = 0.d0
  977. B(8,3) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
  978. B(8,4) = 0.d0
  979. B(8,5) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))
  980. $ )/denom
  981. B(8,6) = (alpha*(drldp*p - rl)*rv*cv**(-2)*(uly - uvy))
  982. & /denom
  983. B(8,7) = 0.d0
  984. B(8,8) = uly
  985. C
  986. C dVpdV
  987. C
  988. dVpdV(1,1) = 1.d0
  989. dVpdV(1,2) = 0.d0
  990. dVpdV(1,3) = 0.d0
  991. dVpdV(1,4) = 0.d0
  992. dVpdV(1,5) = 0.d0
  993. dVpdV(1,6) = 0.d0
  994. dVpdV(1,7) = 0.d0
  995. dVpdV(1,8) = 0.d0
  996.  
  997. dVpdV(2,1) = 0.d0
  998. dVpdV(2,2) = 1.d0
  999. dVpdV(2,3) = 0.d0
  1000. dVpdV(2,4) = 0.d0
  1001. dVpdV(2,5) = 0.d0
  1002. dVpdV(2,6) = 0.d0
  1003. dVpdV(2,7) = 0.d0
  1004. dVpdV(2,8) = 0.d0
  1005.  
  1006. dVpdV(3,1) = 0.d0
  1007. dVpdV(3,2) = 0.d0
  1008. dVpdV(3,3) = 1.d0
  1009. dVpdV(3,4) = 0.d0
  1010. dVpdV(3,5) = 0.d0
  1011. dVpdV(3,6) = 0.d0
  1012. dVpdV(3,7) = 0.d0
  1013. dVpdV(3,8) = 0.d0
  1014.  
  1015. dVpdV(4,1) = 0.d0
  1016. dVpdV(4,2) = 0.d0
  1017. dVpdV(4,3) = 0.d0
  1018. dVpdV(4,4) = 1.d0
  1019. dVpdV(4,5) = 0.d0
  1020. dVpdV(4,6) = 0.d0
  1021. dVpdV(4,7) = 0.d0
  1022. dVpdV(4,8) = 0.d0
  1023.  
  1024. dVpdV(5,1) = 0.d0
  1025. dVpdV(5,2) = 0.d0
  1026. dVpdV(5,3) = 0.d0
  1027. dVpdV(5,4) = 0.d0
  1028. dVpdV(5,5) = 1.d0
  1029. dVpdV(5,6) = 0.d0
  1030. dVpdV(5,7) = 0.d0
  1031. dVpdV(5,8) = 0.d0
  1032.  
  1033. dVpdV(6,1) = 0.d0
  1034. dVpdV(6,2) = 0.d0
  1035. dVpdV(6,3) = 0.d0
  1036. dVpdV(6,4) = 0.d0
  1037. dVpdV(6,5) = 0.d0
  1038. dVpdV(6,6) = 1.d0
  1039. dVpdV(6,7) = 0.d0
  1040. dVpdV(6,8) = 0.d0
  1041.  
  1042. dVpdV(7,1) = 0.d0
  1043. dVpdV(7,2) = 0.d0
  1044. dVpdV(7,3) = 0.d0
  1045. dVpdV(7,4) = 0.d0
  1046. dVpdV(7,5) = 0.d0
  1047. dVpdV(7,6) = dhvdp
  1048. dVpdV(7,7) = dhvdTv
  1049. dVpdV(7,8) = 0.d0
  1050.  
  1051. dVpdV(8,1) = 0.d0
  1052. dVpdV(8,2) = 0.d0
  1053. dVpdV(8,3) = 0.d0
  1054. dVpdV(8,4) = 0.d0
  1055. dVpdV(8,5) = 0.d0
  1056. dVpdV(8,6) = dhldp
  1057. dVpdV(8,7) = 0.d0
  1058. dVpdV(8,8) = dhldTl
  1059. C
  1060. C dVdVp
  1061. C
  1062. dVdVp(1,1) = 1.d0
  1063. dVdVp(1,2) = 0.d0
  1064. dVdVp(1,3) = 0.d0
  1065. dVdVp(1,4) = 0.d0
  1066. dVdVp(1,5) = 0.d0
  1067. dVdVp(1,6) = 0.d0
  1068. dVdVp(1,7) = 0.d0
  1069. dVdVp(1,8) = 0.d0
  1070.  
  1071. dVdVp(2,1) = 0.d0
  1072. dVdVp(2,2) = 1.d0
  1073. dVdVp(2,3) = 0.d0
  1074. dVdVp(2,4) = 0.d0
  1075. dVdVp(2,5) = 0.d0
  1076. dVdVp(2,6) = 0.d0
  1077. dVdVp(2,7) = 0.d0
  1078. dVdVp(2,8) = 0.d0
  1079.  
  1080. dVdVp(3,1) = 0.d0
  1081. dVdVp(3,2) = 0.d0
  1082. dVdVp(3,3) = 1.d0
  1083. dVdVp(3,4) = 0.d0
  1084. dVdVp(3,5) = 0.d0
  1085. dVdVp(3,6) = 0.d0
  1086. dVdVp(3,7) = 0.d0
  1087. dVdVp(3,8) = 0.d0
  1088.  
  1089. dVdVp(4,1) = 0.d0
  1090. dVdVp(4,2) = 0.d0
  1091. dVdVp(4,3) = 0.d0
  1092. dVdVp(4,4) = 1.d0
  1093. dVdVp(4,5) = 0.d0
  1094. dVdVp(4,6) = 0.d0
  1095. dVdVp(4,7) = 0.d0
  1096. dVdVp(4,8) = 0.d0
  1097.  
  1098. dVdVp(5,1) = 0.d0
  1099. dVdVp(5,2) = 0.d0
  1100. dVdVp(5,3) = 0.d0
  1101. dVdVp(5,4) = 0.d0
  1102. dVdVp(5,5) = 1.d0
  1103. dVdVp(5,6) = 0.d0
  1104. dVdVp(5,7) = 0.d0
  1105. dVdVp(5,8) = 0.d0
  1106.  
  1107. dVdVp(6,1) = 0.d0
  1108. dVdVp(6,2) = 0.d0
  1109. dVdVp(6,3) = 0.d0
  1110. dVdVp(6,4) = 0.d0
  1111. dVdVp(6,5) = 0.d0
  1112. dVdVp(6,6) = 1.d0
  1113. dVdVp(6,7) = 0.d0
  1114. dVdVp(6,8) = 0.d0
  1115.  
  1116. dVdVp(7,1) = 0.d0
  1117. dVdVp(7,2) = 0.d0
  1118. dVdVp(7,3) = 0.d0
  1119. dVdVp(7,4) = 0.d0
  1120. dVdVp(7,5) = 0.d0
  1121. dVdVp(7,6) = -dhvdp/dhvdTv
  1122. dVdVp(7,7) = 1.d0/dhvdTv
  1123. dVdVp(7,8) = 0.d0
  1124.  
  1125. dVdVp(8,1) = 0.d0
  1126. dVdVp(8,2) = 0.d0
  1127. dVdVp(8,3) = 0.d0
  1128. dVdVp(8,4) = 0.d0
  1129. dVdVp(8,5) = 0.d0
  1130. dVdVp(8,6) = dhldp*dhvdp/(dhldTl*dhvdTv)
  1131. dVdVp(8,7) = -dhldp/(dhldTl*dhvdTv)
  1132. dVdVp(8,8) = 1.d0/dhldTl
  1133. C
  1134. C dVdVp*A*dVpdV
  1135. C
  1136. do l = 1, 8
  1137. do m = 1, 8
  1138. P1(l,m) = 0.d0
  1139. end do
  1140. end do
  1141. do l = 1, 8
  1142. do m = 1, 8
  1143. do n = 1, 8
  1144. P1(l,m) = P1(l,m) + dVdVp(l,n)*A(n,m)
  1145. end do
  1146. end do
  1147. end do
  1148. do l = 1, 8
  1149. do m = 1, 8
  1150. P2(l,m) = 0.d0
  1151. end do
  1152. end do
  1153. do l = 1, 8
  1154. do m = 1, 8
  1155. do n = 1, 8
  1156. P2(l,m) = P2(l,m) + P1(l,n)*dVpdV(n,m)
  1157. end do
  1158. end do
  1159. end do
  1160. do l = 1, 8
  1161. do m = 1, 8
  1162. A(l,m) = P2(l,m)
  1163. end do
  1164. end do
  1165. C
  1166. C dVdVp*B*dVpdV
  1167. C
  1168. do l = 1, 8
  1169. do m = 1, 8
  1170. P1(l,m) = 0.d0
  1171. end do
  1172. end do
  1173. do l = 1, 8
  1174. do m = 1, 8
  1175. do n = 1, 8
  1176. P1(l,m) = P1(l,m) + dVdVp(l,n)*B(n,m)
  1177. end do
  1178. end do
  1179. end do
  1180. do l = 1, 8
  1181. do m = 1, 8
  1182. P2(l,m) = 0.d0
  1183. end do
  1184. end do
  1185. do l = 1, 8
  1186. do m = 1, 8
  1187. do n = 1, 8
  1188. P2(l,m) = P2(l,m) + P1(l,n)*dVpdV(n,m)
  1189. end do
  1190. end do
  1191. end do
  1192. do l = 1, 8
  1193. do m = 1, 8
  1194. B(l,m) = P2(l,m)
  1195. end do
  1196. end do
  1197. C
  1198. C MEJORAR LO DE LA MULTIPLICACION DE LAS MATRICES DE CAMBIO
  1199. C
  1200. dVdx(1) = MPGRAL.VPOCHA(NLCE,1)*LIMITA(1)
  1201. dVdx(2) = MPGRUV.VPOCHA(NLCE,1)*LIMITA(2)
  1202. dVdx(3) = MPGRUV.VPOCHA(NLCE,3)*LIMITA(3)
  1203. dVdx(4) = MPGRUL.VPOCHA(NLCE,1)*LIMITA(4)
  1204. dVdx(5) = MPGRUL.VPOCHA(NLCE,3)*LIMITA(5)
  1205. dVdx(6) = MPGRP.VPOCHA(NLCE,1)*LIMITA(6)
  1206. dVdx(7) = MPGRTV.VPOCHA(NLCE,1)*LIMITA(7)
  1207. dVdx(8) = MPGRTL.VPOCHA(NLCE,1)*LIMITA(8)
  1208.  
  1209. dVdy(1) = MPGRAL.VPOCHA(NLCE,2)*LIMITA(1)
  1210. dVdy(2) = MPGRUV.VPOCHA(NLCE,2)*LIMITA(2)
  1211. dVdy(3) = MPGRUV.VPOCHA(NLCE,4)*LIMITA(3)
  1212. dVdy(4) = MPGRUL.VPOCHA(NLCE,2)*LIMITA(4)
  1213. dVdy(5) = MPGRUL.VPOCHA(NLCE,4)*LIMITA(5)
  1214. dVdy(6) = MPGRP.VPOCHA(NLCE,2)*LIMITA(6)
  1215. dVdy(7) = MPGRTV.VPOCHA(NLCE,2)*LIMITA(7)
  1216. dVdy(8) = MPGRTL.VPOCHA(NLCE,2)*LIMITA(8)
  1217.  
  1218.  
  1219. C No te olvides de afectar de la correspondiente cp las
  1220. C dVdx y dVdy correspondientes a las entalpias, no vayas
  1221. C a meter la pata con eso
  1222. do m = 1, 8
  1223. dV(m) = 0.d0
  1224. end do
  1225. do m = 1, 8
  1226. do n = 1, 8
  1227. dV(m) = dV(m) + A(m,n)*dVdx(n) + B(m,n)*dVdy(n)
  1228. end do
  1229. end do
  1230. C
  1231. C************* Prediction avec/sans gradients limités
  1232. C
  1233. MPALPP.VPOCHA(NLCE,1)= alpha - DELTAT*dV(1)
  1234. MPUVP.VPOCHA(NLCE,1) = uvx - DELTAT*dV(2)
  1235. MPUVP.VPOCHA(NLCE,2) = uvy - DELTAT*dV(3)
  1236. MPULP.VPOCHA(NLCE,1) = ulx - DELTAT*dV(4)
  1237. MPULP.VPOCHA(NLCE,2) = uly - DELTAT*dV(5)
  1238. MPPP.VPOCHA(NLCE,1) = p - DELTAT*dV(6)
  1239. MPTVP.VPOCHA(NLCE,1) = Tv - DELTAT*dV(7)
  1240. MPTLP.VPOCHA(NLCE,1) = Tl - DELTAT*dV(8)
  1241. C
  1242. C************* Updating densities(perfect and stiffened gas eos
  1243. C
  1244. MPRVP.VPOCHA(NLCE,1) = gamv*MPPP.VPOCHA(NLCE,1)/
  1245. & (Cpv*MPTVP.VPOCHA(NLCE,1)*(gamv - 1.d0))
  1246. MPRLP.VPOCHA(NLCE,1) = gaml*(MPPP.VPOCHA(NLCE,1) + pil)/
  1247. & (Cpl* MPTLP.VPOCHA(NLCE,1)*(gaml - 1.d0))
  1248.  
  1249. ENDDO
  1250. C
  1251. ENDIF
  1252. C
  1253. C
  1254. C***********************************************************************
  1255. C********* CORRECTION **************************************************
  1256. C***********************************************************************
  1257. C
  1258. C**** Boucle sur le faces
  1259. C
  1260. IDIMP1 = IDIM + 1
  1261. DO NLCF = 1, NFAC
  1262. C
  1263. C******* NLCF = numero local du centre de face
  1264. C NGCF = numero global du centre de face
  1265. C NGCEG = numero global du centre ELT "gauche"
  1266. C NLCEG = numero local du centre ELT "gauche"
  1267. C NGCED = numero global du centre ELT "droite"
  1268. C NLCED = numero local du centre ELT "droite"
  1269. C
  1270. NGCEG = IPT1.NUM(1,NLCF)
  1271. NGCF = IPT1.NUM(2,NLCF)
  1272. NGCED = IPT1.NUM(3,NLCF)
  1273. NLCEG = MLENT1.LECT(NGCEG)
  1274. NLCED = MLENT1.LECT(NGCED)
  1275. C
  1276. C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
  1277. C
  1278. NGCF1 = IPT2.NUM(1,NLCF)
  1279. IF( NGCF1 .NE. NGCF) THEN
  1280. LOGAN = .TRUE.
  1281. MESERR(1:40) = 'PRET, subroutine pre52f.eso '
  1282. GOTO 9999
  1283. ENDIF
  1284. C
  1285. C******* Cosinus directeurs des NORMALES aux faces
  1286. C
  1287. C On impose que les normales sont direct "Gauche" -> "Centre"
  1288. C
  1289. INDCEL = (NGCEG-1)*IDIMP1
  1290. XG = XCOOR(INDCEL+1)
  1291. YG = XCOOR(INDCEL+2)
  1292. INDCEL = (NGCF-1)*IDIMP1
  1293. XC = XCOOR(INDCEL + 1)
  1294. YC = XCOOR(INDCEL + 2)
  1295. INDCEL = (NGCED-1)*IDIMP1
  1296. XD = XCOOR(INDCEL+1)
  1297. YD = XCOOR(INDCEL+2)
  1298. DXG = XC - XG
  1299. DYG = YC - YG
  1300. DXD = XC - XD
  1301. DYD = YC - YD
  1302. C
  1303. C******* On calcule le sign du pruduit scalare
  1304. C (Normales de Castem) * (vecteur "gauche" -> "centre")
  1305. C
  1306. CNX = MPNORM.VPOCHA(NLCF,1)
  1307. CNY = MPNORM.VPOCHA(NLCF,2)
  1308. ORIENT = CNX * DXG + CNY * DYG
  1309. ORIENT = SIGN(1.0D0,ORIENT)
  1310. IF(ORIENT .NE. 1.0D0)THEN
  1311. LOGAN = .TRUE.
  1312. MESERR(1:30)=
  1313. & 'PRET , subroutine pre52f.eso. '
  1314. GOTO 9999
  1315. ENDIF
  1316. CNX = CNX * ORIENT
  1317. CNY = CNY * ORIENT
  1318. C
  1319. C********** Cosinus directeurs de tangent 2D
  1320. C
  1321. CTX = -1.0D0 * CNY
  1322. CTY = CNX
  1323. C
  1324. C
  1325. C******* Les autres MELVALs
  1326. C
  1327. C
  1328. C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
  1329. C GAMMA \in (1,3)
  1330. C Si non il faut le faire, en utilisant LOGBOR,
  1331. C LOGNEG, VALER, VAL1, VAL2
  1332. C
  1333. C
  1334. C
  1335. C******* NGCEG = NGCED -> Mur
  1336. C
  1337. IF( NGCEG .EQ. NGCED)THEN
  1338. C
  1339. C********** Sur le mur on fait de reconstruction sur l'etat gauche
  1340. C
  1341. C
  1342. C********** Etat gauche
  1343. C
  1344.  
  1345. VALCEL = MPALPP.VPOCHA(NLCEG,1)
  1346. ALCEL = MELLAL.VPOCHA(NLCEG,1)
  1347. DCEL = MPGRAL.VPOCHA(NLCEG,1)*DXG +
  1348. & MPGRAL.VPOCHA(NLCEG,2)*DYG
  1349. alphag = VALCEL + ALCEL*DCEL
  1350.  
  1351. VALCEL = MPUVP.VPOCHA(NLCEG,1)
  1352. ALCEL = MELLUV.VPOCHA(NLCEG,1)
  1353. DCEL = MPGRUV.VPOCHA(NLCEG,1)*DXG +
  1354. & MPGRUV.VPOCHA(NLCEG,2)*DYG
  1355. uvxg = VALCEL + ALCEL*DCEL
  1356.  
  1357. VALCEL = MPUVP.VPOCHA(NLCEG,2)
  1358. ALCEL = MELLUV.VPOCHA(NLCEG,2)
  1359. DCEL = MPGRUV.VPOCHA(NLCEG,3)*DXG +
  1360. & MPGRUV.VPOCHA(NLCEG,4)*DYG
  1361. uvyg = VALCEL + ALCEL*DCEL
  1362.  
  1363. uvng = uvxg*CNX + uvyg*CNY
  1364. uvtg = uvxg*CTX + uvyg*CTY
  1365.  
  1366. VALCEL = MPULP.VPOCHA(NLCEG,1)
  1367. ALCEL = MELLUL.VPOCHA(NLCEG,1)
  1368. DCEL = MPGRUL.VPOCHA(NLCEG,1)*DXG +
  1369. & MPGRUL.VPOCHA(NLCEG,2)*DYG
  1370. ulxg = VALCEL + ALCEL*DCEL
  1371.  
  1372. VALCEL = MPULP.VPOCHA(NLCEG,2)
  1373. ALCEL = MELLUL.VPOCHA(NLCEG,2)
  1374. DCEL = MPGRUL.VPOCHA(NLCEG,3)*DXG +
  1375. & MPGRUL.VPOCHA(NLCEG,4)*DYG
  1376. ulyg = VALCEL + ALCEL*DCEL
  1377.  
  1378. ulng = ulxg*CNX + ulyg*CNY
  1379. ultg = ulxg*CTX + ulyg*CTY
  1380.  
  1381. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  1382. ALCEL = MELLP.VPOCHA(NLCEG, 1)
  1383. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  1384. & MPGRP.VPOCHA(NLCEG, 2)*DYG
  1385. pg = VALCEL + ALCEL*DCEL
  1386.  
  1387. VALCEL = MPTVP.VPOCHA(NLCEG, 1)
  1388. ALCEL = MELLTV.VPOCHA(NLCEG, 1)
  1389. DCEL = MPGRTV.VPOCHA(NLCEG, 1)*DXG +
  1390. & MPGRTV.VPOCHA(NLCEG, 2)*DYG
  1391. Tvg = VALCEL + ALCEL*DCEL
  1392.  
  1393. VALCEL = MPTLP.VPOCHA(NLCEG, 1)
  1394. ALCEL = MELLTV.VPOCHA(NLCEG, 1)
  1395. DCEL = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
  1396. & MPGRTL.VPOCHA(NLCEG, 2)*DYG
  1397. Tlg = VALCEL + ALCEL*DCEL
  1398.  
  1399. C********** We calculate densities again, we use for
  1400. C the moment, perfect and stiffened gas eos
  1401. rvg = gamv*pg/(cpv*Tvg*(gamv - 1))
  1402. rlg = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1))
  1403. C
  1404. C********** Si l'on fait pas de prediction, ce n'est pas necessaire de
  1405. C controller la positivite' de la pression et de la densité; elle
  1406. C est déjà garantie par la proprieté LED de limiteur; mais, vu
  1407. C que le limiteur n'est pas calculé ici, mais dans un autre
  1408. C operateur, on le fait
  1409. C
  1410. LOGI1 = (pg .LT. 0.0D0) .OR. (rvg .LT. 0.0D0)
  1411. & .OR. (rlg .LT. 0.0D0)
  1412. C
  1413. IF(LOGI1)THEN
  1414. C
  1415. C************* Premier ordre en espace local
  1416. C
  1417. alphag = MPALPC.VPOCHA(NLCEG,1)
  1418. uvng = MPUVC.VPOCHA(NLCEG,1)*CNX +
  1419. & MPUVC.VPOCHA(NLCEG,2)*CNY
  1420. uvtg = MPULC.VPOCHA(NLCEG,1)*CTX +
  1421. & MPUVC.VPOCHA(NLCEG,2)*CTY
  1422. ulng = MPULC.VPOCHA(NLCEG,1)*CNX +
  1423. & MPULC.VPOCHA(NLCEG,2)*CNY
  1424. ultg = MPULC.VPOCHA(NLCEG,1)*CTX +
  1425. & MPULC.VPOCHA(NLCEG,2)*CTY
  1426. pg = MPPC.VPOCHA(NLCEG,1)
  1427. Tvg = MPTVC.VPOCHA(NLCEG,1)
  1428. Tlg = MPTLC.VPOCHA(NLCEG,1)
  1429. rvg = MPRVC.VPOCHA(NLCEG,1)
  1430. rlg = MPRLC.VPOCHA(NLCEG,1)
  1431.  
  1432. ENDIF
  1433. C
  1434. C********** Son etat droite
  1435. C
  1436. alphad = alphag
  1437. uvnd = -1.d0*uvng
  1438. uvtd = uvtg
  1439. ulnd = -1.d0*ulng
  1440. ultd = ultg
  1441. pd = pg
  1442. Tvd = Tvg
  1443. Tld = Tlg
  1444. rvd = rvg
  1445. rld = rlg
  1446. C
  1447. C************* Fin cas mur
  1448. C
  1449. ELSE
  1450. C
  1451. C************* Etat gauche
  1452. C
  1453. VALCEL = MPALPP.VPOCHA(NLCEG,1)
  1454. ALCEL = MELLAL.VPOCHA(NLCEG,1)
  1455. DCEL = MPGRAL.VPOCHA(NLCEG,1)*DXG +
  1456. & MPGRAL.VPOCHA(NLCEG,2)*DYG
  1457. alphag = VALCEL + ALCEL*DCEL
  1458.  
  1459. VALCEL = MPUVP.VPOCHA(NLCEG,1)
  1460. ALCEL = MELLUV.VPOCHA(NLCEG,1)
  1461. DCEL = MPGRUV.VPOCHA(NLCEG,1)*DXG +
  1462. & MPGRUV.VPOCHA(NLCEG,2)*DYG
  1463. uvxg = VALCEL + ALCEL*DCEL
  1464.  
  1465. VALCEL = MPUVP.VPOCHA(NLCEG,2)
  1466. ALCEL = MELLUV.VPOCHA(NLCEG,2)
  1467. DCEL = MPGRUV.VPOCHA(NLCEG,3)*DXG +
  1468. & MPGRUV.VPOCHA(NLCEG,4)*DYG
  1469. uvyg = VALCEL + ALCEL*DCEL
  1470.  
  1471. uvng = uvxg*CNX + uvyg*CNY
  1472. uvtg = uvxg*CTX + uvyg*CTY
  1473.  
  1474. VALCEL = MPULP.VPOCHA(NLCEG,1)
  1475. ALCEL = MELLUL.VPOCHA(NLCEG,1)
  1476. DCEL = MPGRUL.VPOCHA(NLCEG,1)*DXG +
  1477. & MPGRUL.VPOCHA(NLCEG,2)*DYG
  1478. ulxg = VALCEL + ALCEL*DCEL
  1479.  
  1480. VALCEL = MPULP.VPOCHA(NLCEG,2)
  1481. ALCEL = MELLUL.VPOCHA(NLCEG,2)
  1482. DCEL = MPGRUL.VPOCHA(NLCEG,3)*DXG +
  1483. & MPGRUL.VPOCHA(NLCEG,4)*DYG
  1484. ulyg = VALCEL + ALCEL*DCEL
  1485.  
  1486. ulng = ulxg*CNX + ulyg*CNY
  1487. ultg = ulxg*CTX + ulyg*CTY
  1488.  
  1489. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  1490. ALCEL = MELLP.VPOCHA(NLCEG, 1)
  1491. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  1492. & MPGRP.VPOCHA(NLCEG, 2)*DYG
  1493. pg = VALCEL + ALCEL*DCEL
  1494.  
  1495. VALCEL = MPTVP.VPOCHA(NLCEG, 1)
  1496. ALCEL = MELLTV.VPOCHA(NLCEG, 1)
  1497. DCEL = MPGRTV.VPOCHA(NLCEG, 1)*DXG +
  1498. & MPGRTV.VPOCHA(NLCEG, 2)*DYG
  1499. Tvg = VALCEL + ALCEL*DCEL
  1500.  
  1501. VALCEL = MPTLP.VPOCHA(NLCEG, 1)
  1502. ALCEL = MELLTL.VPOCHA(NLCEG, 1)
  1503. DCEL = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
  1504. & MPGRTL.VPOCHA(NLCEG, 2)*DYG
  1505. Tlg = VALCEL + ALCEL*DCEL
  1506.  
  1507. C********** We calculate densities again, we use for
  1508. C the moment, perfect and stiffened gas eos
  1509. rvg = gamv*pg/(cpv*Tvg*(gamv - 1))
  1510. rlg = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1))
  1511. C
  1512. C********** Positivite
  1513. C
  1514. LOGI1 = (pg .LT. 0.D0) .OR. (rvg .LT. 0.D0)
  1515. & .OR. (rlg .LT. 0.D0)
  1516. C
  1517. IF(LOGI1)THEN
  1518. C
  1519. C************* Premier ordre en espace local
  1520. C
  1521. alphag = MPALPC.VPOCHA(NLCEG,1)
  1522. uvng = MPUVC.VPOCHA(NLCEG,1)*CNX +
  1523. & MPUVC.VPOCHA(NLCEG,2)*CNY
  1524. uvtg = MPULC.VPOCHA(NLCEG,1)*CTX +
  1525. & MPUVC.VPOCHA(NLCEG,2)*CTY
  1526. ulng = MPULC.VPOCHA(NLCEG,1)*CNX +
  1527. & MPULC.VPOCHA(NLCEG,2)*CNY
  1528. ultg = MPULC.VPOCHA(NLCEG,1)*CTX +
  1529. & MPULC.VPOCHA(NLCEG,2)*CTY
  1530. pg = MPPC.VPOCHA(NLCEG,1)
  1531. Tvg = MPTVC.VPOCHA(NLCEG,1)
  1532. Tlg = MPTLC.VPOCHA(NLCEG,1)
  1533. rvg = MPRVC.VPOCHA(NLCEG,1)
  1534. rlg = MPRLC.VPOCHA(NLCEG,1)
  1535.  
  1536. ENDIF
  1537. C
  1538. C********** Etat droite
  1539. C
  1540. VALCEL = MPALPP.VPOCHA(NLCED, 1)
  1541. ALCEL = MELLAL.VPOCHA(NLCED, 1)
  1542. DCEL = MPGRAL.VPOCHA(NLCED, 1)*DXD +
  1543. & MPGRAL.VPOCHA(NLCED, 2)*DYD
  1544. alphad = VALCEL + ALCEL * DCEL
  1545.  
  1546. VALCEL = MPUVP.VPOCHA(NLCED, 1)
  1547. ALCEL = MELLUV.VPOCHA(NLCED, 1)
  1548. DCEL = MPGRUV.VPOCHA(NLCED, 1)*DXD +
  1549. & MPGRUV.VPOCHA(NLCED, 2)*DYD
  1550. uvxd = VALCEL + ALCEL * DCEL
  1551.  
  1552. VALCEL = MPUVP.VPOCHA(NLCED, 2)
  1553. ALCEL = MELLUV.VPOCHA(NLCED, 2)
  1554. DCEL = MPGRUV.VPOCHA(NLCED, 3)*DXD +
  1555. & MPGRUV.VPOCHA(NLCED, 4)*DYD
  1556. uvyd = VALCEL + ALCEL * DCEL
  1557.  
  1558. uvnd = uvxd*CNX + uvyd*CNY
  1559. uvtd = uvxd*CTX + uvyd*CTY
  1560.  
  1561. VALCEL = MPULP.VPOCHA(NLCED, 1)
  1562. ALCEL = MELLUL.VPOCHA(NLCED, 1)
  1563. DCEL = MPGRUL.VPOCHA(NLCED, 1)*DXD +
  1564. & MPGRUL.VPOCHA(NLCED, 2)*DYD
  1565. ulxd = VALCEL + ALCEL * DCEL
  1566.  
  1567. VALCEL = MPULP.VPOCHA(NLCED,2)
  1568. ALCEL = MELLUL.VPOCHA(NLCED, 2)
  1569. DCEL = MPGRUL.VPOCHA(NLCED, 3)*DXD +
  1570. & MPGRUL.VPOCHA(NLCED, 4)*DYD
  1571. ulyd = VALCEL + ALCEL * DCEL
  1572.  
  1573. ulnd = ulxd*CNX + ulyd*CNY
  1574. ultd = ulxd*CTX + ulyd*CTY
  1575.  
  1576. VALCEL = MPPP.VPOCHA(NLCED, 1)
  1577. ALCEL = MELLP.VPOCHA(NLCED, 1)
  1578. DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD +
  1579. & MPGRP.VPOCHA(NLCED, 2)*DYD
  1580. pd = VALCEL + ALCEL * DCEL
  1581.  
  1582. VALCEL = MPTVP.VPOCHA(NLCED, 1)
  1583. ALCEL = MELLTV.VPOCHA(NLCED, 1)
  1584. DCEL = MPGRTV.VPOCHA(NLCED, 1)*DXD +
  1585. & MPGRTV.VPOCHA(NLCED, 2)*DYD
  1586. Tvd = VALCEL + ALCEL * DCEL
  1587.  
  1588. VALCEL = MPTLP.VPOCHA(NLCED, 1)
  1589. ALCEL = MELLTL.VPOCHA(NLCED, 1)
  1590. DCEL = MPGRTL.VPOCHA(NLCED,1)*DXD +
  1591. & MPGRTL.VPOCHA(NLCED,2)*DYD
  1592. Tld = VALCEL + ALCEL * DCEL
  1593. C
  1594. C********** We calculate densities again, we use for
  1595. C the moment, perfect and stiffened gas eos
  1596. rvd = gamv*pd/(cpv*Tvd*(gamv - 1))
  1597. rld = gaml*(pd + pil)/(cpl*Tld*(gaml - 1))
  1598. C
  1599. C********** Positivite
  1600. C
  1601. LOGI1 = (pd .LT. 0.0D0) .OR. (rvd .LT. 0.0D0)
  1602. & .OR. (rld .LT. 0.0D0)
  1603. C
  1604. IF(LOGI1)THEN
  1605. C
  1606. C************* Premier ordre en espace local
  1607. C
  1608. alphad = MPALPC.VPOCHA(NLCED,1)
  1609. uvnd = MPUVC.VPOCHA(NLCED,1)*CNX +
  1610. & MPUVC.VPOCHA(NLCED,2)*CNY
  1611. uvtd = MPUVC.VPOCHA(NLCED,1)*CTX +
  1612. & MPUVC.VPOCHA(NLCED,2)*CTY
  1613. pd = MPPC.VPOCHA(NLCED,1)
  1614. Tvd = MPTVC.VPOCHA(NLCED,1)
  1615. Tld = MPTLC.VPOCHA(NLCED,1)
  1616. rvd = MPRVC.VPOCHA(NLCED,1)
  1617. rld = MPRLC.VPOCHA(NLCED,1)
  1618.  
  1619. ENDIF
  1620. ENDIF
  1621. C
  1622. C******** Les MELVALs
  1623. C
  1624. MLALP.VELCHE(1,NLCF) = alphag
  1625. MLALP.VELCHE(3,NLCF) = alphad
  1626. MLP.VELCHE(1,NLCF) = pg
  1627. MLP.VELCHE(3,NLCF) = pd
  1628. MLTV.VELCHE(1,NLCF) = Tvg
  1629. MLTV.VELCHE(3,NLCF) = Tvd
  1630. MLTL.VELCHE(1,NLCF) = Tlg
  1631. MLTL.VELCHE(3,NLCF) = Tld
  1632. MLRV.VELCHE(1,NLCF) = rvg
  1633. MLRV.VELCHE(3,NLCF) = rvd
  1634. MLRL.VELCHE(1,NLCF) = rlg
  1635. MLRL.VELCHE(3,NLCF) = rld
  1636.  
  1637. MLUVN.VELCHE(1,NLCF) = uvng
  1638. MLUVN.VELCHE(3,NLCF) = uvnd
  1639. MLUVT.VELCHE(1,NLCF) = uvtg
  1640. MLUVT.VELCHE(3,NLCF) = uvtd
  1641. MLVNX.VELCHE(1,NLCF) = CNX
  1642. MLVNY.VELCHE(1,NLCF) = CNY
  1643. MLVTX.VELCHE(1,NLCF) = CTX
  1644. MLVTY.VELCHE(1,NLCF) = CTY
  1645.  
  1646. MLULN.VELCHE(1,NLCF) = ulng
  1647. MLULN.VELCHE(3,NLCF) = ulnd
  1648. MLULT.VELCHE(1,NLCF) = ultg
  1649. MLULT.VELCHE(3,NLCF) = ultd
  1650. MLLNX.VELCHE(1,NLCF) = CNX
  1651. MLLNY.VELCHE(1,NLCF) = CNY
  1652. MLLTX.VELCHE(1,NLCF) = CTX
  1653. MLLTY.VELCHE(1,NLCF) = CTY
  1654. C
  1655. ENDDO
  1656. C
  1657. C**** Desactivation des SEGMENTs
  1658. C
  1659. SEGDES IPT1
  1660. SEGDES IPT2
  1661. C
  1662. C**** Le MPOVALs 'Prediction' sont detruits (si existentes)
  1663. C
  1664. IF(LOGTEM)THEN
  1665. SEGSUP MPALPP
  1666. SEGSUP MPUVP
  1667. SEGSUP MPULP
  1668. SEGSUP MPPP
  1669. SEGSUP MPTVP
  1670. SEGSUP MPTLP
  1671. SEGSUP MPRVP
  1672. SEGSUP MPRLP
  1673. ENDIF
  1674. C
  1675. SEGDES MPALPC
  1676. SEGDES MPGRAL
  1677. SEGDES MPUVC
  1678. SEGDES MPGRUV
  1679. SEGDES MPULC
  1680. SEGDES MPGRUL
  1681. SEGDES MPPC
  1682. SEGDES MPGRP
  1683. SEGDES MPTVC
  1684. SEGDES MPGRTV
  1685. SEGDES MPTLC
  1686. SEGDES MPGRTL
  1687. SEGDES MPRVC
  1688. SEGDES MPRLC
  1689. SEGDES MPNORM
  1690. SEGDES MELLAL
  1691. SEGDES MELLUV
  1692. SEGDES MELLUL
  1693. SEGDES MELLP
  1694. SEGDES MELLTV
  1695. SEGDES MELLTL
  1696. C
  1697. SEGDES MLALP
  1698. SEGDES MLP
  1699. SEGDES MLTV
  1700. SEGDES MLTL
  1701. SEGDES MLRV
  1702. SEGDES MLRL
  1703. SEGDES MLUVN
  1704. SEGDES MLUVT
  1705. SEGDES MLULN
  1706. SEGDES MLULT
  1707. SEGDES MLVNX
  1708. SEGDES MLVNY
  1709. SEGDES MLVTX
  1710. SEGDES MLVTY
  1711. SEGDES MLLNX
  1712. SEGDES MLLNY
  1713. SEGDES MLLTX
  1714. SEGDES MLLTY
  1715. C
  1716. C**** Destruction du MELNTI correspondance local/global
  1717. C
  1718. SEGSUP MLENT1
  1719. C
  1720. 9999 CONTINUE
  1721. C
  1722. RETURN
  1723. END
  1724.  
  1725.  
  1726.  
  1727.  
  1728.  
  1729.  
  1730.  
  1731.  
  1732.  
  1733.  

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