Télécharger riecom.eso

Retour à la liste

Numérotation des lignes :

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

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