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