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 296 : CALL RACC(EPSI,NITER,GAMG,GAMD,A1,A2,A3,X1,LOGNC,LOGAN,MESERR) 297 : IF(LOGNC .OR. LOGAN) THEN 298 : GOTO 9999 299 : ENDIF 300 : X3=X1+LOG(A2) 301 : AUX=A1 / (VLF1(X1,GAMG)*VLF3(X3,GAMD)) 302 : X2=LOG(AUX) 303 : IF(ABS(X1).LT.ZERO) X1=0.D0 304 : IF(ABS(X2).LT.ZERO) X2=0.D0 305 : IF(ABS(X3).LT.ZERO) X3=0.D0 306 : C 307 : PI=PG*EXP(-X1) 307 : CALL VLH1(H1,H1P,X1,GAMG) 308 : UI=UNG+CG*H1 309 : RA1 = VLF1(X1,GAMG) 310 : RA2 = RA1 * EXP(X2) 311 : R1 = RA1 * RG 312 : R2 = RA2 * RG 313 : C1=SQRT(GAMG*PI/R1) 314 : C2=SQRT(GAMD*PI/R2) 315 : C 316 : C******* Les differentes pentes frontieres 317 : C 318 : IF( (ABS(RA1-1.0D0) .LE. ZERO) .AND. 319 : & (X1.LT.0.0D0)) X1=0.D0 320 : IF( (ABS(RA2-A1) .LE. ZERO) .AND. 321 : & (X3.LT.0.)) X3=0.D0 322 : C 323 : IF(X1.LT.0.D0) THEN 324 : C 325 : C********** 1-ONDE : choc 326 : C 327 : AUX = (R1*UI-RG*UNG) / (R1-RG) 328 : FR1G = AUX 329 : FR1D = AUX 330 : ELSE 331 : C 332 : C********* 1-ONDE : detente 333 : C 334 : FR1G=UNG-CG 335 : FR1D=UI-C1 336 : ENDIF 337 : FR2 = UI 338 : IF(X3.LT.0.D0) THEN 339 : C 340 : C******* 3-ONDE : choc 341 : C 342 : AUX = (RD*UND-R2*UI) / (RD-R2) 343 : FR3G = AUX 344 : FR3D = AUX 345 : ELSE 346 : C 347 : C******* 3-ONDE : detente 348 : C 349 : FR3G=UI+C2 350 : FR3D=UND+CD 351 : ENDIF 352 : C 353 : C******* Position de la droite x/t=VGRIL P/R aux differents cas 354 : C 355 : IF(FR1G .GT. VGRIL) THEN 356 : C 357 : C******* ETAT : GAUCHE 358 : C 359 : RR=RG 360 : UNR = UNG 361 : UTR = UTG 362 : UVR = UVG 363 : RETR=RETG 364 : PR = PG 365 : LOGETD = .FALSE. 366 : GAMR = GAMG 367 : ELSEIF(FR1D.GT.VGRIL) THEN 368 : C 369 : C********** ETAT : 1-DETENTE 370 : C 371 : TAU=0.5D0*GM1G/GAMG 372 : AUX = GM1G*(UNG-VGRIL)/(CG*GP1G) + 2.D0/GP1G 373 : XX = -LOG(AUX)/TAU 374 : IF(XX .LT. 0.0D0) THEN 375 : LOGAN = .TRUE. 376 : MESERR= 'RIEMAN, subroutine riecom.eso ' 377 : GOTO 9999 378 : ENDIF 379 : PI=PG*EXP(-XX) 380 : RR=RG*EXP(-XX/GAMG) 381 : UNR = UNG + 2.D0*CG*USGM1G * (1.D0-EXP(-TAU*XX)) 382 : UTR = UTG 383 : UVR = UVG 384 : RETR = PI*USGM1G + 0.5D0*RR*(UNR*UNR+UTR*UTR+UVR*UVR) 385 : PR = PI 386 : LOGETD = .FALSE. 387 : GAMR = GAMG 388 : ELSEIF(FR2.GT.VGRIL) THEN 389 : C 390 : C********** ETAT : INTERMEDIAIRE 1 391 : C 392 : RR = R1 393 : UNR = UI 394 : UTR = UTG 395 : UVR = UVG 396 : RETR = PI*USGM1G + 0.5D0*RR*(UNR*UNR+UTR*UTR 397 : & +UVR*UVR) 398 : PR = PI 399 : LOGETD = .FALSE. 400 : GAMR = GAMG 401 : ELSEIF(FR3G.GT.VGRIL) THEN 402 : C 403 : C********** ETAT : INTERMEDIAIRE 2 404 : C 405 : RR=R2 406 : UNR = UI 407 : UTR = UTD 408 : UVR = UVD 409 : RETR = PI*USGM1D + 0.5D0*RR*(UNR*UNR+UTR*UTR 410 : & +UVR*UVR) 411 : PR = PI 412 : LOGETD = .TRUE. 413 : GAMR = GAMD 414 : ELSEIF(FR3D.GT.VGRIL) THEN 415 : C 416 : C********** ETAT : 3-DETENTE 417 : C 418 : TAU=0.5D0*GM1D/GAMD 419 : AUX = 2.D0/GP1D - GM1D*(UI-VGRIL)/(C2*GP1D) 420 : XX = LOG(AUX) / TAU 421 : IF(XX.LT.0) THEN 422 : LOGAN = .TRUE. 423 : MESERR = 424 : & 'RIEMAN, subroutine riecom.eso ' 425 : GOTO 9999 426 : ENDIF 427 : PI=PI*EXP(XX) 428 : RR=R2*EXP(XX/GAMD) 429 : UNR = UI + 2.D0*C2*USGM1D*(EXP(TAU*XX)-1.D0) 430 : UTR = UTD 431 : UVR = UVD 432 : RETR = PI*USGM1D + 0.5D0*RR*(UNR*UNR+UTR*UTR 433 : & +UVR*UVR) 434 : PR = PI 435 : LOGETD = .TRUE. 436 : GAMR = GAMD 437 : ELSE 438 : C 439 : C********** ETAT : DROITE 440 : C 441 : RR = RD 442 : UNR = UND 443 : UTR = UTD 444 : UVR = UVD 445 : RETR = RETD 446 : PR = PD 447 : LOGETD = .TRUE. 448 : GAMR = GAMD 449 : ENDIF 450 : ENDIF 451 : C 452 : 9999 CONTINUE 453 : C 454 : RETURN 455 : END 456 : 457 : 458 : 459 : 460 : 461 :
© Cast3M 2003 - All rights reserved.
Disclaimer