Télécharger rieco2.eso

Retour à la liste

Numérotation des lignes :

rieco2
  1. C RIECO2 SOURCE CHAT 05/01/13 02:56:21 5004
  2. C RIECOM SOURCE BOB2 96/10/01 22:21:28 2307
  3. SUBROUTINE RIECO2(NITER,EPSI,ZERO,
  4. & VGRIL,
  5. & GAMG,RG,PG,UNG,UTG,UVG,
  6. & GAMD,RD,PD,UND,UTD,UVD,
  7. & RR,UNR,UTR,UVR,RETR,PR,GAMR,LOGETD,
  8. & LOGAN,LOGNC,MESERR)
  9. C
  10. C************************************************************************
  11. C
  12. C PROJET : CASTEM 2000
  13. C
  14. C NOM : RIECOM
  15. C
  16. C DESCRIPTION : Voir FLURIE
  17. C
  18. C Solution du problème de Riemann dans le
  19. C repaire (n,t) in x/t = VGRIL
  20. C
  21. C Parametrisation de Smoller
  22. C (voir J. SMOLLER, "Shock Waves and Reaction
  23. C Diffusion Equations", Springer Verlag, 1983)
  24. C
  25. C LANGAGE : FORTRAN 77
  26. C
  27. C AUTEUR : P. GALON DRN/DMT/SEMT/TTMF
  28. C
  29. C************************************************************************
  30. C
  31. C APPELES
  32. C
  33. C RIECOM ---- RACC ---- WNVXC ---- VLH1
  34. C |
  35. C |
  36. C -------- VLH1
  37. C |
  38. C |
  39. C -------- VLF1
  40. C |
  41. C |
  42. C -------- VLF3
  43. C
  44. C************************************************************************
  45. C
  46. C**** Entrées:
  47. C
  48. C NMAX = nombre maximum d'itérations pour Newton-Rapson
  49. C
  50. C EPSI = erreur tolérée pour Newton-Rapson
  51. C
  52. C ZERO = tolérance d'egalite pour REAL*8
  53. C
  54. C VGRIL = vitesse de la surface (ALE)
  55. C
  56. C GAMG, GAMD = les "gamma" du gaz (gauche et droite)
  57. C
  58. C RG, RD = les densités
  59. C
  60. C PG, PD = les pressions
  61. C
  62. C UNG, UND = les vitesses normales
  63. C
  64. C UTG, UTD = les vitesses tangentielles
  65. C
  66. C**** Sorties:
  67. C
  68. C RR = rho(x/t = VGRIL)
  69. C
  70. C UNR = un(x/t = VGRIL)
  71. C
  72. C UTR = ut(x/t = VGRIL)
  73. C
  74. C RETR = rho*et(x/t = VGRIL)
  75. C
  76. C PR = p(x/t = VGRIL) (indispensable en sortie dans
  77. C le cas du vide)
  78. C
  79. C LOGETD = .TRUE. -> on est à droite de la discontinuité
  80. C de contact.
  81. C
  82. C .FALSE. -> on est à gauche de la discontinuité
  83. C de contact.
  84. C
  85. C GAMR = gamma(x/t = VRILL)
  86. C
  87. C LOGAN = .TRUE. -> anomalie de programmation.
  88. C
  89. C LOGNC = .TRUE. -> Newton Rapson ne converge pas
  90. C
  91. C MESERR = message d'erreur
  92. C
  93. C
  94. C************************************************************************
  95. C
  96. C HISTORIQUE (Anomalies et modifications éventuelles)
  97. C
  98. C HISTORIQUE : 6.1.98 , Modifié par A. BECCANTINI, DRN/DMT/SEMT/TTMF
  99. C pour pouvoir traiter la formation du vide.
  100. C
  101. C
  102. C************************************************************************
  103. C
  104. C**** N.B.: toutes les variables sont DECLAREES!
  105. C
  106. C
  107. IMPLICIT INTEGER(I-N)
  108. INTEGER NITER
  109. REAL*8 EPSI,ZERO,VGRIL
  110. & ,GAMG,RG,PG,UNG,UTG,UVG
  111. & ,GAMD,RD,PD,UND,UTD,UVD
  112. & ,RR,UNR,UTR,UVR,RETR,PR,GAMR
  113. & ,GM1G,GP1G,USGM1G,GM1D,GP1D,USGM1D
  114. & ,AUX,RETG,RETD,CG,CD,DC,DP,U1G,U1D,DU,CEL1,CEL2
  115. & ,U2G,U2D,CEL3
  116. & ,A1,A2,A3,X1,X2,X3
  117. & ,H1,H1P,VLF1,VLF3
  118. & ,PI,UI,RA1,R1,RA2,R2,C1,C2
  119. & ,FR1G, FR1D, FR2, FR3G, FR3D, TAU, XX
  120.  
  121. LOGICAL LOGAN,LOGNC,LOGETD
  122. CHARACTER*(40) MESERR
  123. C
  124. C**** Initialisation de LOGNC, LOGAN,MESERR ne doit pas etre faite ici,
  125. C mais avant, i.e.
  126. C
  127. C LOGNC = .FALSE.
  128. C LOGAN = .FALSE.
  129. C MESERR = ' '
  130. C
  131. C
  132. C**** N.B. On suppose que la positivité de RG, RD, PG, PD a été
  133. C déjà verifiée.
  134. C Ceci est très important si on travaille en "NaNQ"
  135. C
  136. GM1G = GAMG - 1.0D0
  137. GP1G = GAMG + 1.0D0
  138. USGM1G = 1.0D0 / GM1G
  139. GM1D = GAMD - 1.0D0
  140. GP1D = GAMD + 1.0D0
  141. USGM1D = 1.0D0 / GM1D
  142.  
  143. C
  144. C**** Calcul des energies totales.
  145. C
  146. AUX = 0.5D0 * RG * (UNG*UNG + UTG*UTG + UVG*UVG)
  147. RETG = PG * USGM1G + AUX
  148. AUX = 0.5D0 * RD * (UND*UND + UTD*UTD + UVD*UVD)
  149. RETD = PD * USGM1D + AUX
  150. C
  151. C**** Trois cas possibles:
  152. C
  153. C a) etats egaux ou discontinuité de contact;
  154. C
  155. C b) formation du vide
  156. C
  157. C c) choc + discontinuité de contact + détente
  158. C
  159. CG = SQRT(GAMG*PG/RG)
  160. CD = SQRT(GAMD*PD/RD)
  161. DC = MAX(CG,CD)
  162. DP = MAX(RG,RD) * DC * DC
  163. C
  164. U1G = UNG + 2.0D0 * USGM1G * CG
  165. U1D = UND - 2.0D0 * USGM1D * CD
  166. DU = UNG - UND
  167. CEL1 = ABS(DU)/DC
  168. CEL2 = ABS(PG-PD)/DP
  169. C
  170. C**** NB : control adimensionalizé
  171. C
  172. IF( (CEL1 .LT. ZERO) .AND. (CEL2 .LT. ZERO) ) THEN
  173. C
  174. C
  175. C******* Cas a) : PG = PD et UNG = UND
  176. C
  177. C
  178. IF(VGRIL .LE. UND) THEN
  179. RR=RG
  180. UNR=UNG
  181. UTR=UTG
  182. UVR=UVG
  183. RETR=RETG
  184. PR = PG
  185. GAMR = GAMG
  186. LOGETD = .FALSE.
  187. ELSE
  188. RR=RD
  189. UNR=UND
  190. UTR=UTD
  191. UVR=UVD
  192. RETR=RETD
  193. PR = PD
  194. GAMR = GAMD
  195. LOGETD = .TRUE.
  196. ENDIF
  197. ELSEIF((U1G-U1D) .LE. 0.0D0)THEN
  198. C
  199. C
  200. C******* Cas b) : formation de vide, connecté aux etats initiaux par
  201. C deux ondes de détente.
  202. C
  203. U2G = UNG -CG
  204. U2D = UND + CD
  205. IF(VGRIL .LT. U2G)THEN
  206. RR=RG
  207. UNR=UNG
  208. UTR=UTG
  209. UVR=UVG
  210. RETR=RETG
  211. PR = PG
  212. GAMR = GAMG
  213. LOGETD = .FALSE.
  214. ELSEIF(VGRIL .LT. U1G)THEN
  215. CEL3 = (U1G-VGRIL)*GM1G/GP1G/CG
  216. C
  217. C********** Si VGRIL = U1G alors CEL3 = 0 (vide)
  218. C
  219. C Si VGRIG = U2G alors U1G-VGRIL = GP1G / GM1G * CG
  220. C
  221. C CEL3 = 1
  222. C
  223. UNR = U1G - ( 2.0D0 * USGM1G* CEL3 * CG)
  224. CEL3 = CEL3**(2.0D0*USGM1G)
  225. RR = CEL3 * RG
  226. CEL3 = CEL3**GAMG
  227. PR = CEL3* PG
  228. UTR= UTG
  229. UVR=UVG
  230. AUX = 0.5D0 * RR * (UNR*UNR + UTR*UTR + UVR*UVR)
  231. RETR = USGM1G * PR + AUX
  232. GAMR = GAMG
  233. LOGETD = .FALSE.
  234. ELSEIF(VGRIL .LT. U1D)THEN
  235. RR = 0.0D0
  236. C
  237. C********** UNR, UTR non definies: on fait une approximation,
  238. C mais de toute façon RR=0 -> RR*UNR=RR*UTR=0
  239. C
  240.  
  241. UNR = 0.5D0*(U1D+U1G)
  242. IF(VGRIL .LT. UNR)THEN
  243. UTR = UTG
  244. UVR=UVG
  245. LOGETD = .FALSE.
  246. GAMR = GAMG
  247. ELSE
  248. UTR = UTD
  249. UVR = UVD
  250. LOGETD = .TRUE.
  251. GAMR = GAMD
  252. ENDIF
  253. RETR = 0.0D0
  254. PR = 0.0D0
  255. ELSEIF(VGRIL .LT. U2D)THEN
  256. CEL3 = (VGRIL-U1D)*GM1D/GP1D/CD
  257. C
  258. C********** Si VGRIL = U1D alors CEL3 = 0 (vide)
  259. C
  260. C Si VGRIG = U2D alors VGRIL-U1D = GP1D / GM1D * CD
  261. C
  262. C CEL3 = 1
  263. C
  264. UNR = U1D + ( 2.0D0 * USGM1D* CEL3 * CD)
  265. CEL3 = CEL3**(2.0D0*USGM1D)
  266. RR = CEL3 * RD
  267. CEL3 = CEL3**GAMD
  268. PR = CEL3* PD
  269. UTR = UTD
  270. UVR = UVD
  271. AUX = 0.5D0 * RR * (UNR*UNR + UTR*UTR + UVR*UVR)
  272. RETR = USGM1D * PR + AUX
  273. LOGETD = .TRUE.
  274. GAMR = GAMD
  275. ELSE
  276. RR = RD
  277. UNR = UND
  278. UTR = UTD
  279. UVR = UVD
  280. RETR = RETD
  281. PR = PD
  282. LOGETD = .TRUE.
  283. GAMR = GAMD
  284. ENDIF
  285. ELSE
  286. C
  287. C
  288. C******* Cas c) : il n'y a pas la formation de vide.
  289. C
  290. C
  291. C
  292. C******* X1 = 1-ONDE : racine de la resolution Riemann
  293. C
  294. A1=RD/RG
  295. A2=PD/PG
  296. A3=(UND-UNG)/CG
  297. CALL RACC(EPSI,NITER,GAMG,GAMD,A1,A2,A3,X1,LOGNC,LOGAN,MESERR)
  298. IF(LOGNC .OR. LOGAN) THEN
  299. GOTO 9999
  300. ENDIF
  301. X3=X1+LOG(A2)
  302. AUX=A1 / (VLF1(X1,GAMG)*VLF3(X3,GAMD))
  303. X2=LOG(AUX)
  304. IF(ABS(X1).LT.ZERO) X1=0.D0
  305. IF(ABS(X2).LT.ZERO) X2=0.D0
  306. IF(ABS(X3).LT.ZERO) X3=0.D0
  307. C
  308. PI=PG*EXP(-X1)
  309. CALL VLH1(H1,H1P,X1,GAMG)
  310. UI=UNG+CG*H1
  311. RA1 = VLF1(X1,GAMG)
  312. RA2 = RA1 * EXP(X2)
  313. R1 = RA1 * RG
  314. R2 = RA2 * RG
  315. C1=SQRT(GAMG*PI/R1)
  316. C2=SQRT(GAMD*PI/R2)
  317. C
  318. C******* Les differentes pentes frontieres
  319. C
  320. IF( (ABS(RA1-1.0D0) .LE. ZERO) .AND.
  321. & (X1.LT.0.0D0)) X1=0.D0
  322. IF( (ABS(RA2-A1) .LE. ZERO) .AND.
  323. & (X3.LT.0.)) X3=0.D0
  324. C
  325. IF(X1.LT.0.D0) THEN
  326. C
  327. C********** 1-ONDE : choc
  328. C
  329. AUX = (R1*UI-RG*UNG) / (R1-RG)
  330. FR1G = AUX
  331. FR1D = AUX
  332. ELSE
  333. C
  334. C********* 1-ONDE : detente
  335. C
  336. FR1G=UNG-CG
  337. FR1D=UI-C1
  338. ENDIF
  339. FR2 = UI
  340. IF(X3.LT.0.D0) THEN
  341. C
  342. C******* 3-ONDE : choc
  343. C
  344. AUX = (RD*UND-R2*UI) / (RD-R2)
  345. FR3G = AUX
  346. FR3D = AUX
  347. ELSE
  348. C
  349. C******* 3-ONDE : detente
  350. C
  351. FR3G=UI+C2
  352. FR3D=UND+CD
  353. ENDIF
  354. C
  355. C******* Position de la droite x/t=VGRIL P/R aux differents cas
  356. C
  357. IF(FR1G .GT. VGRIL) THEN
  358. C
  359. C******* ETAT : GAUCHE
  360. C
  361. RR=RG
  362. UNR = UNG
  363. UTR = UTG
  364. UVR = UVG
  365. RETR=RETG
  366. PR = PG
  367. LOGETD = .FALSE.
  368. GAMR = GAMG
  369. ELSEIF(FR1D.GT.VGRIL) THEN
  370. C
  371. C********** ETAT : 1-DETENTE
  372. C
  373. TAU=0.5D0*GM1G/GAMG
  374. AUX = GM1G*(UNG-VGRIL)/(CG*GP1G) + 2.D0/GP1G
  375. XX = -LOG(AUX)/TAU
  376. IF(XX .LT. 0.0D0) THEN
  377. LOGAN = .TRUE.
  378. MESERR= 'RIEMAN, subroutine riecom.eso '
  379. GOTO 9999
  380. ENDIF
  381. PI=PG*EXP(-XX)
  382. RR=RG*EXP(-XX/GAMG)
  383. UNR = UNG + 2.D0*CG*USGM1G * (1.D0-EXP(-TAU*XX))
  384. UTR = UTG
  385. UVR = UVG
  386. RETR = PI*USGM1G + 0.5D0*RR*(UNR*UNR+UTR*UTR+UVR*UVR)
  387. PR = PI
  388. LOGETD = .FALSE.
  389. GAMR = GAMG
  390. ELSEIF(FR2.GT.VGRIL) THEN
  391. C
  392. C********** ETAT : INTERMEDIAIRE 1
  393. C
  394. RR = R1
  395. UNR = UI
  396. UTR = UTG
  397. UVR = UVG
  398. RETR = PI*USGM1G + 0.5D0*RR*(UNR*UNR+UTR*UTR
  399. & +UVR*UVR)
  400. PR = PI
  401. LOGETD = .FALSE.
  402. GAMR = GAMG
  403. ELSEIF(FR3G.GT.VGRIL) THEN
  404. C
  405. C********** ETAT : INTERMEDIAIRE 2
  406. C
  407. RR=R2
  408. UNR = UI
  409. UTR = UTD
  410. UVR = UVD
  411. RETR = PI*USGM1D + 0.5D0*RR*(UNR*UNR+UTR*UTR
  412. & +UVR*UVR)
  413. PR = PI
  414. LOGETD = .TRUE.
  415. GAMR = GAMD
  416. ELSEIF(FR3D.GT.VGRIL) THEN
  417. C
  418. C********** ETAT : 3-DETENTE
  419. C
  420. TAU=0.5D0*GM1D/GAMD
  421. AUX = 2.D0/GP1D - GM1D*(UI-VGRIL)/(C2*GP1D)
  422. XX = LOG(AUX) / TAU
  423. IF(XX.LT.0) THEN
  424. LOGAN = .TRUE.
  425. MESERR =
  426. & 'RIEMAN, subroutine riecom.eso '
  427. GOTO 9999
  428. ENDIF
  429. PI=PI*EXP(XX)
  430. RR=R2*EXP(XX/GAMD)
  431. UNR = UI + 2.D0*C2*USGM1D*(EXP(TAU*XX)-1.D0)
  432. UTR = UTD
  433. UVR = UVD
  434. RETR = PI*USGM1D + 0.5D0*RR*(UNR*UNR+UTR*UTR
  435. & +UVR*UVR)
  436. PR = PI
  437. LOGETD = .TRUE.
  438. GAMR = GAMD
  439. ELSE
  440. C
  441. C********** ETAT : DROITE
  442. C
  443. RR = RD
  444. UNR = UND
  445. UTR = UTD
  446. UVR = UVD
  447. RETR = RETD
  448. PR = PD
  449. LOGETD = .TRUE.
  450. GAMR = GAMD
  451. ENDIF
  452. ENDIF
  453. C
  454. 9999 CONTINUE
  455. C
  456. RETURN
  457. END
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  

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