1 : C FLURIE SOURCE CHAT 05/01/13 00:04:53 5004 2 : SUBROUTINE FLURIE(NESP, 3 : & GAMG,ROG,PG,UNG,UTG, 4 : & GAMD,ROD,PD,UND,UTD, 5 : & YG,YD,FLUX, 6 : & CELLT, 7 : & LOGNC,LOGAN,MESERR) 8 : C 9 : C************************************************************************ 10 : C 11 : C PROJET : CASTEM 2000 12 : C 13 : C NOM : FLURIE 14 : C 15 : C DESCRIPTION : Formulation Volumes Finis pour les Equations 16 : C d'Euler Multi-Especes relatives à un melange 17 : C de gaz ideals. 18 : C 19 : C Calcul du flux aux interfaces avec la mèthode 20 : C de Godunov (Solveur de Riemann exact). 21 : C 22 : C Parametrisation de Smoller 23 : C 24 : C (voir: 25 : C 1) J. SMOLLER, "Shock Waves and Reaction 26 : C Diffusion Equations", Springer Verlag, 1983; 27 : C 28 : C 2) Beccantini, Paillere, 29 : C "Upwind Flux Splitting Schemes..." 30 : C RAPPORT DMT 97//268 31 : C ) 32 : C 33 : C LANGAGE : FORTRAN 77 34 : C 35 : C AUTEUR : A. BECCANTINI DRN/DMT/SEMT/TTMF 36 : C 37 : C************************************************************************ 38 : C 39 : C APPELES 40 : C 41 : C FLURIE 42 : C | 43 : C | 44 : C ---- RIECOM ---- RACC ---- WNVXC ---- VLH1 45 : C | 46 : C | 47 : C -------- VLH1 48 : C | 49 : C | 50 : C -------- VLF1 51 : C | 52 : C | 53 : C -------- VLF3 54 : C 55 : C************************************************************************ 56 : C 57 : C**** Entrées: 58 : C 59 : C NESP = nombre d'especes considérées dans les Equations 60 : C d'Euler 61 : C 62 : C GAMG, GAMD = les "gamma" du gaz (gauche et droite) 63 : C 64 : C ROG, ROD = les densités 65 : C 66 : C PG, PD = les pressions 67 : C 68 : C UNG, UND = vitesses normales 69 : C 70 : C UTG, UTD = vitesses tangentielles 71 : C 72 : C YG, YD = tables des fractiones massiques 73 : C 74 : C**** Sorties: 75 : C 76 : C FLUX = table du flux a l'interface, i.e. 77 : C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht, 78 : C rho*un*y1, ...) 79 : C 80 : C CELLT = condition de stabilité, i.e. 81 : C 82 : C dT/diamax < cellt 83 : C 84 : C LOGNC = si .TRUE. -> no convergence 85 : C 86 : C LOGAN = si .TRUE. -> une anomalie a été detectée 87 : C 88 : C MESERR = message d'erreur 89 : C 90 : C 91 : C************************************************************************ 92 : C 93 : C HISTORIQUE (Anomalies et modifications éventuelles) 94 : C 95 : C HISTORIQUE : Créé le 7.1.98 96 : C 97 : C************************************************************************ 98 : C 99 : C N.B.: Toutes les variables sont DECLAREES 100 : C 101 : IMPLICIT INTEGER(I-N) 102 : INTEGER NESP, NMAX, NMAX0, I1 103 : REAL*8 ZERO, ZERO0, EPSI, EPSI0 104 : & ,GAMG,ROG,PG,UNG,UTG 105 : & ,GAMD,ROD,PD,UND,UTD 106 : & ,VGRIL 107 : & ,ROR,UNR,UTR,ROETR,PR 108 : & ,VREL, FLUM, RHTR 109 : & ,CELLT, AR, GAMR 110 : & ,YG(*),YD(*),FLUX(*) 111 : C 112 : CHARACTER*(40) MESERR 113 : LOGICAL LOGNC, LOGAN, LOGETD, LOGLAR 114 : PARAMETER(NMAX=1000,ZERO = 1.0D-14,EPSI=1D-12) 115 : C 116 : C**** ZERO = tolérance d'egalite pour REAL*8 117 : C 118 : C NMAX = nombre max d'itérations in racc.eso 119 : C pour calculer l'etats d'intersection 120 : C invariants de Riemann - conditions de 121 : C Rankine-Hugoniot. 122 : C 123 : C EPSI = erreur tolérée in racc.eso 124 : C 125 : C 126 : C**** Initialisation de LOGNC, LOGAN,MESERR ne doit pas etre faite ici, 127 : C mais avant, i.e. 128 : C 129 : C LOGNC = .FALSE. 130 : C LOGAN = .FALSE. 131 : C MESERR(1:40) = ' ' 132 : C 133 : C 134 : C**** Les fractions massiques 135 : C 136 : C YG, YD déjà definis 137 : C 138 : C**** Le flux à l'interfaces sont dans le repaire (n,t) 139 : C 140 : C 141 : C 142 : C**** VGRIL = pour ALE 143 : C 144 : VGRIL = 0.0D0 145 : C 146 : C**** On protege le "PARAMETER" 147 : C 148 : NMAX0 = NMAX 149 : ZERO0 = ZERO 150 : EPSI0 = EPSI 151 : C 152 : C**** Calcul de l'etat a l'interface 153 : C 153 : CALL RIECOM(NMAX0,EPSI0,ZERO0, 154 : & VGRIL, 155 : & GAMG,ROG,PG,UNG,UTG, 156 : & GAMD,ROD,PD,UND,UTD, 157 : & ROR,UNR,UTR,ROETR,PR,GAMR,LOGETD, 158 : & LOGAN,LOGNC,MESERR) 159 : C 160 : C**** Valeurs à l' interface (x/t=VGRIL) 161 : C 162 : C ROR = rho 163 : C UNR = = un 164 : C UTR = = ut 165 : C ROETR = rho*et 166 : C PR = p 167 : C LOGETD = .FALSE. -> on est à gauche de la discontinuité de contact. 168 : C .TRUE. -> on est à droite de la discontinuité de contact. 169 : C 170 : IF(LOGNC .OR. LOGAN) GOTO 9999 171 : C 172 : C**** Les flux 173 : C 174 : VREL = UNR - VGRIL 175 : FLUM = ROR * VREL 176 : RHTR = ROETR + PR 177 : FLUX(1) = FLUM 178 : FLUX(2) = FLUM * UNR + PR 179 : FLUX(3) = FLUM * UTR 180 : FLUX(4) = VREL * RHTR 181 : C 182 : C**** TEST 183 : C 184 : C On peut demontrer que le sign de la vitesse 185 : C sur la discontinuité de contact est le meme 186 : C que le sign de ROR*(UNR-VGRIL) 187 : C (voir Larrouturou, 188 : C J. Comp. Phys., 95, 1991) 189 : C 190 : C Par contre, si ROR = 0, ceci n'est plus vrai, 191 : C mais dans ce cas RHO*U*Y = 0 192 : C 193 : C 194 : C FLUM < 0 -> On est à droite 195 : C 196 : LOGLAR = FLUM .LT. 0.0D0 197 : IF((ABS(FLUM) .GT. 1.0D-8) .AND. (LOGLAR .NEQV. LOGETD)) THEN 198 : MESERR = 'RIEMANN, flurie.eso (Larrouturou) ' 199 : LOGAN = .TRUE. 200 : GOTO 9999 201 : ENDIF 202 : C 203 : C**** Fin test 204 : C 205 : C 206 : IF(NESP .GT. 0)THEN 207 : IF(LOGETD)THEN 208 : DO I1 = 1, NESP 209 : FLUX(I1+4) = FLUM * YD(I1) 210 : ENDDO 211 : ELSE 212 : DO I1 = 1, NESP 213 : FLUX(I1+4) = FLUM * YG(I1) 214 : ENDDO 215 : ENDIF 216 : ENDIF 217 : C 218 : C**** CFL 219 : C 220 : C PPw AR = SQRT (GAMR * PR / ROR) 221 : C PPw CELLT = 1.0D0 / (ABS(VREL)+AR) 222 : C+PPw Dans l'espace, on n'entend pas crier... 223 : IF(PR.EQ.0.D0.AND.ROR.EQ.0.D0)THEN 224 : AR=0.D0 225 : ELSE 226 : AR = SQRT (GAMR * PR / ROR) 227 : ENDIF 228 : IF(AR.EQ.0.D0.AND.VREL.EQ.0.D0)THEN 229 : CELLT=1.D0/ZERO 230 : ELSE 231 : CELLT = 1.0D0 / (ABS(VREL)+AR) 232 : ENDIF 233 : C+PPw 234 : C 235 : 9999 CONTINUE 236 : C 237 : RETURN 238 : END 239 : 240 : 241 : 242 : 243 :
© Cast3M 2003 - All rights reserved.
Disclaimer