rieco2
C RIECO2 SOURCE CHAT 05/01/13 02:56:21 5004 C RIECOM SOURCE BOB2 96/10/01 22:21:28 2307 & VGRIL, & GAMG,RG,PG,UNG,UTG,UVG, & GAMD,RD,PD,UND,UTD,UVD, & RR,UNR,UTR,UVR,RETR,PR,GAMR,LOGETD, & LOGAN,LOGNC,MESERR) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : RIECOM C C DESCRIPTION : Voir FLURIE C C Solution du problème de Riemann dans le C repaire (n,t) in x/t = VGRIL C C Parametrisation de Smoller C (voir J. SMOLLER, "Shock Waves and Reaction C Diffusion Equations", Springer Verlag, 1983) C C LANGAGE : FORTRAN 77 C C AUTEUR : P. GALON DRN/DMT/SEMT/TTMF C C************************************************************************ C C APPELES C C RIECOM ---- RACC ---- WNVXC ---- VLH1 C | C | C -------- VLH1 C | C | C -------- VLF1 C | C | C -------- VLF3 C C************************************************************************ C C**** Entrées: C C NMAX = nombre maximum d'itérations pour Newton-Rapson C C EPSI = erreur tolérée pour Newton-Rapson C C ZERO = tolérance d'egalite pour REAL*8 C C VGRIL = vitesse de la surface (ALE) C C GAMG, GAMD = les "gamma" du gaz (gauche et droite) C C RG, RD = les densités C C PG, PD = les pressions C C UNG, UND = les vitesses normales C C UTG, UTD = les vitesses tangentielles C C**** Sorties: C C RR = rho(x/t = VGRIL) C C UNR = un(x/t = VGRIL) C C UTR = ut(x/t = VGRIL) C C RETR = rho*et(x/t = VGRIL) C C PR = p(x/t = VGRIL) (indispensable en sortie dans C le cas du vide) C C LOGETD = .TRUE. -> on est à droite de la discontinuité C de contact. C C .FALSE. -> on est à gauche de la discontinuité C de contact. C C GAMR = gamma(x/t = VRILL) C C LOGAN = .TRUE. -> anomalie de programmation. C C LOGNC = .TRUE. -> Newton Rapson ne converge pas C C MESERR = message d'erreur C C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : 6.1.98 , Modifié par A. BECCANTINI, DRN/DMT/SEMT/TTMF C pour pouvoir traiter la formation du vide. C C C************************************************************************ C C**** N.B.: toutes les variables sont DECLAREES! C C IMPLICIT INTEGER(I-N) INTEGER NITER & ,GAMG,RG,PG,UNG,UTG,UVG & ,GAMD,RD,PD,UND,UTD,UVD & ,RR,UNR,UTR,UVR,RETR,PR,GAMR & ,GM1G,GP1G,USGM1G,GM1D,GP1D,USGM1D & ,AUX,RETG,RETD,CG,CD,DC,DP,U1G,U1D,DU,CEL1,CEL2 & ,U2G,U2D,CEL3 & ,A1,A2,A3,X1,X2,X3 & ,FR1G, FR1D, FR2, FR3G, FR3D, TAU, XX LOGICAL LOGAN,LOGNC,LOGETD CHARACTER*(40) MESERR C C**** Initialisation de LOGNC, LOGAN,MESERR ne doit pas etre faite ici, C mais avant, i.e. C C LOGNC = .FALSE. C LOGAN = .FALSE. C MESERR = ' ' C C C**** N.B. On suppose que la positivité de RG, RD, PG, PD a été C déjà verifiée. C Ceci est très important si on travaille en "NaNQ" C GM1G = GAMG - 1.0D0 GP1G = GAMG + 1.0D0 USGM1G = 1.0D0 / GM1G GM1D = GAMD - 1.0D0 GP1D = GAMD + 1.0D0 USGM1D = 1.0D0 / GM1D C C**** Calcul des energies totales. C AUX = 0.5D0 * RG * (UNG*UNG + UTG*UTG + UVG*UVG) RETG = PG * USGM1G + AUX AUX = 0.5D0 * RD * (UND*UND + UTD*UTD + UVD*UVD) RETD = PD * USGM1D + AUX C C**** Trois cas possibles: C C a) etats egaux ou discontinuité de contact; C C b) formation du vide C C c) choc + discontinuité de contact + détente C CG = SQRT(GAMG*PG/RG) CD = SQRT(GAMD*PD/RD) DC = MAX(CG,CD) C U1G = UNG + 2.0D0 * USGM1G * CG U1D = UND - 2.0D0 * USGM1D * CD DU = UNG - UND CEL1 = ABS(DU)/DC C C**** NB : control adimensionalizé C C C C******* Cas a) : PG = PD et UNG = UND C C IF(VGRIL .LE. UND) THEN RR=RG UNR=UNG UTR=UTG UVR=UVG RETR=RETG PR = PG GAMR = GAMG LOGETD = .FALSE. ELSE RR=RD UNR=UND UTR=UTD UVR=UVD RETR=RETD PR = PD GAMR = GAMD LOGETD = .TRUE. ENDIF ELSEIF((U1G-U1D) .LE. 0.0D0)THEN C C C******* Cas b) : formation de vide, connecté aux etats initiaux par C deux ondes de détente. C U2G = UNG -CG U2D = UND + CD IF(VGRIL .LT. U2G)THEN RR=RG UNR=UNG UTR=UTG UVR=UVG RETR=RETG PR = PG GAMR = GAMG LOGETD = .FALSE. ELSEIF(VGRIL .LT. U1G)THEN CEL3 = (U1G-VGRIL)*GM1G/GP1G/CG C C********** Si VGRIL = U1G alors CEL3 = 0 (vide) C C Si VGRIG = U2G alors U1G-VGRIL = GP1G / GM1G * CG C C CEL3 = 1 C UNR = U1G - ( 2.0D0 * USGM1G* CEL3 * CG) CEL3 = CEL3**(2.0D0*USGM1G) RR = CEL3 * RG CEL3 = CEL3**GAMG PR = CEL3* PG UTR= UTG UVR=UVG AUX = 0.5D0 * RR * (UNR*UNR + UTR*UTR + UVR*UVR) RETR = USGM1G * PR + AUX GAMR = GAMG LOGETD = .FALSE. ELSEIF(VGRIL .LT. U1D)THEN RR = 0.0D0 C C********** UNR, UTR non definies: on fait une approximation, C mais de toute façon RR=0 -> RR*UNR=RR*UTR=0 C UNR = 0.5D0*(U1D+U1G) IF(VGRIL .LT. UNR)THEN UTR = UTG UVR=UVG LOGETD = .FALSE. GAMR = GAMG ELSE UTR = UTD UVR = UVD LOGETD = .TRUE. GAMR = GAMD ENDIF RETR = 0.0D0 PR = 0.0D0 ELSEIF(VGRIL .LT. U2D)THEN CEL3 = (VGRIL-U1D)*GM1D/GP1D/CD C C********** Si VGRIL = U1D alors CEL3 = 0 (vide) C C Si VGRIG = U2D alors VGRIL-U1D = GP1D / GM1D * CD C C CEL3 = 1 C UNR = U1D + ( 2.0D0 * USGM1D* CEL3 * CD) CEL3 = CEL3**(2.0D0*USGM1D) RR = CEL3 * RD CEL3 = CEL3**GAMD PR = CEL3* PD UTR = UTD UVR = UVD AUX = 0.5D0 * RR * (UNR*UNR + UTR*UTR + UVR*UVR) RETR = USGM1D * PR + AUX LOGETD = .TRUE. GAMR = GAMD ELSE RR = RD UNR = UND UTR = UTD UVR = UVD RETR = RETD PR = PD LOGETD = .TRUE. GAMR = GAMD ENDIF ELSE C C C******* Cas c) : il n'y a pas la formation de vide. C C C C******* X1 = 1-ONDE : racine de la resolution Riemann C A1=RD/RG A2=PD/PG A3=(UND-UNG)/CG IF(LOGNC .OR. LOGAN) THEN GOTO 9999 ENDIF X3=X1+LOG(A2) X2=LOG(AUX) C PI=PG*EXP(-X1) UI=UNG+CG*H1 RA2 = RA1 * EXP(X2) R1 = RA1 * RG R2 = RA2 * RG C C******* Les differentes pentes frontieres C & (X1.LT.0.0D0)) X1=0.D0 & (X3.LT.0.)) X3=0.D0 C IF(X1.LT.0.D0) THEN C C********** 1-ONDE : choc C AUX = (R1*UI-RG*UNG) / (R1-RG) FR1G = AUX FR1D = AUX ELSE C C********* 1-ONDE : detente C FR1G=UNG-CG ENDIF FR2 = UI IF(X3.LT.0.D0) THEN C C******* 3-ONDE : choc C AUX = (RD*UND-R2*UI) / (RD-R2) FR3G = AUX FR3D = AUX ELSE C C******* 3-ONDE : detente C FR3D=UND+CD ENDIF C C******* Position de la droite x/t=VGRIL P/R aux differents cas C IF(FR1G .GT. VGRIL) THEN C C******* ETAT : GAUCHE C RR=RG UNR = UNG UTR = UTG UVR = UVG RETR=RETG PR = PG LOGETD = .FALSE. GAMR = GAMG ELSEIF(FR1D.GT.VGRIL) THEN C C********** ETAT : 1-DETENTE C TAU=0.5D0*GM1G/GAMG AUX = GM1G*(UNG-VGRIL)/(CG*GP1G) + 2.D0/GP1G XX = -LOG(AUX)/TAU IF(XX .LT. 0.0D0) THEN LOGAN = .TRUE. MESERR= 'RIEMAN, subroutine riecom.eso ' GOTO 9999 ENDIF PI=PG*EXP(-XX) RR=RG*EXP(-XX/GAMG) UNR = UNG + 2.D0*CG*USGM1G * (1.D0-EXP(-TAU*XX)) UTR = UTG UVR = UVG RETR = PI*USGM1G + 0.5D0*RR*(UNR*UNR+UTR*UTR+UVR*UVR) PR = PI LOGETD = .FALSE. GAMR = GAMG C C********** ETAT : INTERMEDIAIRE 1 C RR = R1 UNR = UI UTR = UTG UVR = UVG RETR = PI*USGM1G + 0.5D0*RR*(UNR*UNR+UTR*UTR & +UVR*UVR) PR = PI LOGETD = .FALSE. GAMR = GAMG ELSEIF(FR3G.GT.VGRIL) THEN C C********** ETAT : INTERMEDIAIRE 2 C RR=R2 UNR = UI UTR = UTD UVR = UVD RETR = PI*USGM1D + 0.5D0*RR*(UNR*UNR+UTR*UTR & +UVR*UVR) PR = PI LOGETD = .TRUE. GAMR = GAMD ELSEIF(FR3D.GT.VGRIL) THEN C C********** ETAT : 3-DETENTE C TAU=0.5D0*GM1D/GAMD XX = LOG(AUX) / TAU IF(XX.LT.0) THEN LOGAN = .TRUE. MESERR = & 'RIEMAN, subroutine riecom.eso ' GOTO 9999 ENDIF PI=PI*EXP(XX) RR=R2*EXP(XX/GAMD) UTR = UTD UVR = UVD RETR = PI*USGM1D + 0.5D0*RR*(UNR*UNR+UTR*UTR & +UVR*UVR) PR = PI LOGETD = .TRUE. GAMR = GAMD ELSE C C********** ETAT : DROITE C RR = RD UNR = UND UTR = UTD UVR = UVD RETR = RETD PR = PD LOGETD = .TRUE. GAMR = GAMD ENDIF ENDIF C 9999 CONTINUE C RETURN END
© Cast3M 2003 - Tous droits réservés.
Mentions légales