1 : C FLHUS1 SOURCE CHAT 05/01/13 00:03:45 5004 2 : SUBROUTINE FLHUS1(NESP, 3 : & GAMG,ROG,PG,UNG,UTG, 4 : & GAMD,ROD,PD,UND,UTD, 5 : & YG,YD,FLU1,FLU2, 6 : & CELLT, 7 : & LOGNC,MESERR,LOGAN) 8 : C************************************************************************ 9 : C 10 : C PROJET : CASTEM 2000 11 : C 12 : C NOM : FLHUS1 13 : C 14 : C DESCRIPTION : Formulation Volumes Finis pour les Equations 15 : C d'Euler Multi-Especes relatives à un melange 16 : C de gaz ideals. 17 : C 18 : C Calcul du flux aux interfaces avec la mèthode 19 : C HUS de Coquel - Liou. 20 : C 21 : C Hybridation de van Leer "FVS" 22 : C avec Osher "FDS" in "reverse order decomposition." 23 : C 24 : C (voir: 25 : C 1) Beccantini, Paillere, 26 : C "Upwind Flux Splitting Schemes..." 27 : C RAPPORT DMT 97//268 28 : C 2) Coquel, Liou 29 : C "Hybrid Upwind Splitting (HUS) by a Field-by- 30 : C Field Decomposition " 31 : C 1993. NASA TM 106843) 32 : C 33 : C LANGUAGE : 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 FLHUS1 ------ FLUPVL 42 : C | 43 : C --------- FLUMVL 44 : C | 45 : C --------- RACHUS 46 : C 47 : C 48 : C************************************************************************ 49 : C 50 : C**** Entrées: 51 : C 52 : C NESP = nombre d'especes considérées dans les Equations 53 : C d'Euler 54 : C 55 : C GAMG, GAMD = les "gamma" du gaz (gauche et droite) 56 : C 57 : C ROG, ROD = les densités 58 : C 59 : C PG, PD = les pressions 60 : C 61 : C UNG, UND = vitesses normales 62 : C 63 : C UTG, UTD = vitesses tangentielles 64 : C 65 : C YG, YD = tables des fractiones massiques 66 : C 67 : C**** Sorties: 68 : C 69 : C FLU1 = table du flux a l'interface, i.e. 70 : C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht, 71 : C rho*un*y1, ...) 72 : C 73 : C IFLU2 = table de travaille, utilizé ici mais definie 74 : C avant (pour gagner du temps CPU). 75 : C 76 : C CELLT = condition de stabilité, i.e. 77 : C 78 : C dT/diamax < cellt 79 : C 80 : C LOGNC, MESSER = paramètres qui segnalent des eventuels problèmes 81 : C de convergence 82 : C 83 : C LOGAN = .TRUE. -> anomalie detectée 84 : C 85 : C************************************************************************ 86 : C 87 : C HISTORIQUE (Anomalies et modifications éventuelles) 88 : C 89 : C HISTORIQUE : Créé le 6.1.98 90 : C 91 : C************************************************************************ 92 : C 93 : C N.B.: Toutes les variables sont DECLAREES 94 : C 95 : 96 : IMPLICIT INTEGER(I-N) 97 : INTEGER NESP, I1, NMAX, NMAX0 98 : REAL*8 GAMG,ROG,PG 99 : & ,GAMD,ROD,PD 100 : & ,GM1G,AG2,AG,UNG,UTG,MG,HTG 101 : & ,GM1D,AD2,AD,UND,UTD,MD,HTD 102 : & ,UTG2,UTD2 103 : & ,ZERO, ZERO0 104 : & ,DU 105 : & ,COEFFG, COEFFD, UTILG, UTILD 106 : & ,USTAR,RSTARG,RSTARD,PSTAR,CELLG,CELLD 107 : & ,ASTARG,ASTARD,MSTARG,MSTARD,USTAR2,HTSG,HTSD 108 : & ,CELL 109 : & ,AMG,AMD,CELLT0,CELLT 110 : & ,YG(*),YD(*),FLU1(*),FLU2(*) 111 : CHARACTER*(40) MESERR 112 : LOGICAL LOGNC, LARCO, LOGI, LOGAN 113 : C 114 : C**** ZERO = tolerance d'egalite pour REAL*8 115 : C NMAX = numero max d'iterations in RACHUS 116 : C pour calculer l'etats d'intersection 117 : C des invariants de Riemann 118 : C 119 : C LARCO = Logical; si .TRUE. on calcule le flux 120 : C relative aux fractiones massiques avec 121 : C la correction de Larrouturou 122 : C (voir J. Comp. Phys., 95, 1991) 123 : C 124 : PARAMETER(NMAX=1000,ZERO = 1.0D-12,LARCO = .TRUE.) 125 : C 126 : LOGNC = .FALSE. 127 : MESERR = ' ' 128 : C 129 : C**** YG, YD, FLU1, FLU2 déjà defini 130 : C 131 : C 132 : C**** Calcul du flux avec le FVS de van Leer 133 : C 134 : C 135 : C**** Onde de "gauche" a "droite" 136 : C 137 : GM1G = GAMG - 1.0D0 138 : AG2 = GAMG * PG/ ROG 139 : AG = SQRT(AG2) 140 : UTG2 = UTG*UTG 141 : MG = UNG / AG 142 : HTG = AG2 /GM1G + 0.5D0 * (UNG*UNG + UTG2) 143 : AMG = ABS(MG) 144 : CELLT = (2.0D0*GAMG + (AMG*(3.0D0-GAMG)))/(GAMG+3.0D0) 145 : & /((AMG+1.0D0)*AG) 146 : C 146 : CALL FLUPVL(NESP,GAMG,ROG,MG,AG,UTG,HTG,YG,FLU1) 147 : C 148 : C 149 : C**** Onde de "droite" a "gauche". 150 : C 151 : GM1D = GAMD - 1.0D0 152 : AD2 = GAMD * PD/ ROD 153 : AD = SQRT(AD2) 154 : UTD2 = UTD*UTD 155 : MD = UND / AD 156 : HTD = AD2 / GM1D +0.5D0 * (UND*UND + UTD2) 157 : AMD = ABS(MD) 158 : CELLT0 = (2.0D0*GAMD + (AMD*(3.0D0-GAMD)))/(GAMD+3.0D0) 159 : & /((AMD+1.0D0)*AD) 160 : CELLT = MIN(CELLT,CELLT0) 161 : C 161 : CALL FLUMVL(NESP,GAMD,ROD,MD,AD,UTD,HTD,YD,FLU2) 162 : C 163 : DO I1 = 1, NESP+4 164 : FLU1(I1) = FLU1(I1) + FLU2(I1) 165 : ENDDO 166 : C 167 : C**** CALCUL des etats gauche et droite sur le champ linearment 168 : C dégénéré en imposant la 'reverse order decomposition' 169 : C d'Osher, i.e. la decomposition d'Osher 'physique' 170 : C 171 : C 172 : C**** Osher reversal composition: 173 : C 174 : C Vide -> Les invariants de Riemann ne se coupent pas 175 : C -> RETURN 176 : C 177 : COEFFG = 2.0D0 / GM1G 178 : COEFFD = 2.0D0 / GM1D 179 : UTILG = UNG + COEFFG * AG 180 : UTILD = UND - COEFFD * AD 181 : DU = UTILG - UTILD 182 : IF(DU .GT. ZERO)THEN 183 : C 184 : C 185 : C******* Calcul des intersections des Invariants de Riemann 186 : C 187 : C******* Protection des PARAMETER NMAX, ZERO 188 : C N.B. Les subroutines pouvent les changer 189 : C 190 : NMAX0 = NMAX 191 : ZERO0 = ZERO 192 : C 192 : CALL RACHUS(NMAX0,ZERO0, 193 : & GAMG,COEFFG,AG,UTILG,PG, 194 : & GAMD,COEFFD,AD,UTILD,PD, 195 : & USTAR,PSTAR,LOGNC,MESERR,LOGAN) 196 : C 197 : C******* Pas de convergence de la mèthode de Newton-Rapson en RACHUS 198 : C 199 : IF(LOGNC) GOTO 9999 200 : C 201 : C******* Calcul des etats a Gauche et a Droite du champ linearment dégénéré 202 : C 203 : CELLG = (PSTAR/PG)**(1.0D0/GAMG) 204 : CELLD = (PSTAR/PD)**(1.0D0/GAMD) 205 : ASTARG=AG*CELLG**(1.0D0/COEFFG) 206 : ASTARD=AD*CELLD**(1.0D0/COEFFD) 207 : RSTARG=ROG*CELLG 208 : RSTARD=ROD*CELLD 209 : MSTARG = USTAR / ASTARG 210 : MSTARD = USTAR / ASTARD 211 : USTAR2 = USTAR*USTAR 212 : HTSG = (ASTARG * ASTARG) / GM1G + 0.5D0 * (USTAR2 + UTG2) 213 : HTSD = (ASTARD * ASTARD) / GM1D + 0.5D0 * (USTAR2 + UTD2) 214 : C 215 : C******* CFL 216 : C 217 : AMG = ABS(MSTARG) 218 : CELLT0 = (2.0D0*GAMG + (AMG*(3.0D0-GAMG)))/(GAMG+3.0D0) 219 : & /((AMG+1.0D0)*ASTARG) 220 : CELLT = MIN(CELLT,CELLT0) 221 : AMD = ABS(MSTARD) 222 : CELLT0 = (2.0D0*GAMD + (AMD*(3.0D0-GAMD)))/(GAMD+3.0D0) 223 : & /((AMD+1.0D0)*ASTARD) 224 : CELLT = MIN(CELLT,CELLT0) 225 : C 226 : C******* Partie antidiffusive du flux 227 : C 228 : IF(USTAR .GT. 0.0D0)THEN 228 : CALL FLUMVL(NESP, 229 : & GAMD,RSTARD,MSTARD,ASTARD,UTD,HTSD,YD,FLU2) 230 : C 231 : DO I1 = 1, NESP+4 232 : FLU1(I1) = FLU1(I1) - FLU2(I1) 233 : ENDDO 234 : C 234 : CALL FLUMVL(NESP, 235 : & GAMG,RSTARG,MSTARG,ASTARG,UTG,HTSG,YG,FLU2) 236 : C 237 : DO I1 = 1, NESP+4 238 : FLU1(I1) = FLU1(I1) + FLU2(I1) 239 : ENDDO 240 : ELSE 240 : CALL FLUPVL(NESP, 241 : & GAMD,RSTARD,MSTARD,ASTARD,UTD,HTSD,YD,FLU2) 242 : C 243 : DO I1 = 1, NESP+4 244 : FLU1(I1) = FLU1(I1) + FLU2(I1) 245 : ENDDO 246 : C 246 : CALL FLUPVL(NESP, 247 : & GAMG,RSTARG,MSTARG,ASTARG,UTG,HTSG,YG,FLU2) 248 : C 249 : DO I1 = 1, NESP+4 250 : FLU1(I1) = FLU1(I1) - FLU2(I1) 251 : ENDDO 252 : ENDIF 253 : C 254 : C******* Correction de Larrouturou pour la TVD des fractions massiques. 255 : C pour la vitesse tangentielle 256 : C 257 : C N.B. Pour le solver de Riemann exact cette propriete est 258 : C valide 259 : C Dans le cas de vide ce n'est pas necessaire: 260 : C VanLeer a la TVD des fractions massiques. 261 : C 262 : LOGI = LARCO 263 : IF(LOGI)THEN 264 : CELL = FLU1(1) 265 : IF(CELL .GT. 0.0D0)THEN 266 : FLU1(3) = CELL * UTG 267 : DO I1 = 1, NESP 268 : FLU1(I1+4) = CELL * YG(I1) 269 : ENDDO 270 : ELSE 271 : FLU1(3) = CELL * UTD 272 : DO I1 = 1, NESP 273 : FLU1(I1+4) = CELL * YD(I1) 274 : ENDDO 275 : ENDIF 276 : ENDIF 277 : ENDIF 278 : C 279 : 9999 CONTINUE 280 : RETURN 281 : END 282 : 283 : 284 : 285 : 286 : 287 : 288 : 289 : 290 : 291 :
© Cast3M 2003 - All rights reserved.
Disclaimer