1 : C FLHU21 SOURCE CHAT 05/01/13 00:03:38 5004 2 : SUBROUTINE FLHU21(NESP, 3 : & GAMG,ROG,PG,UNG,UTG,UVG, 4 : & GAMD,ROD,PD,UND,UTD,UVD, 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,UVG,MG,HTG 101 : & ,GM1D,AD2,AD,UND,UTD,UVD,MD,HTD 102 : & ,UTG2,UTD2,UVG2,UVD2 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 : UVG2= UVG*UVG 142 : MG = UNG / AG 143 : HTG = AG2 /GM1G + 0.5D0 * (UNG*UNG + UTG2 + UVG2) 144 : AMG = ABS(MG) 145 : CELLT = (2.0D0*GAMG + (AMG*(3.0D0-GAMG)))/(GAMG+3.0D0) 146 : & /((AMG+1.0D0)*AG) 147 : C 147 : CALL FLUPV2(NESP,GAMG,ROG,MG,AG,UTG,UVG,HTG,YG,FLU1) 148 : C 149 : C 150 : C**** Onde de "droite" a "gauche". 151 : C 152 : GM1D = GAMD - 1.0D0 153 : AD2 = GAMD * PD/ ROD 154 : AD = SQRT(AD2) 155 : UTD2 = UTD*UTD 156 : UVD2 = UVD*UVD 157 : MD = UND / AD 158 : HTD = AD2 / GM1D +0.5D0 * (UND*UND + UTD2 + UVD2) 159 : AMD = ABS(MD) 160 : CELLT0 = (2.0D0*GAMD + (AMD*(3.0D0-GAMD)))/(GAMD+3.0D0) 161 : & /((AMD+1.0D0)*AD) 162 : CELLT = MIN(CELLT,CELLT0) 163 : C 163 : CALL FLUMV2(NESP,GAMD,ROD,MD,AD,UTD,UVD,HTD,YD,FLU2) 164 : C 165 : DO I1 = 1, NESP+5 166 : FLU1(I1) = FLU1(I1) + FLU2(I1) 167 : ENDDO 168 : C 169 : C**** CALCUL des etats gauche et droite sur le champ linearment 170 : C dégénéré en imposant la 'reverse order decomposition' 171 : C d'Osher, i.e. la decomposition d'Osher 'physique' 172 : C 173 : C 174 : C**** Osher reversal composition: 175 : C 176 : C Vide -> Les invariants de Riemann ne se coupent pas 177 : C -> RETURN 178 : C 179 : COEFFG = 2.0D0 / GM1G 180 : COEFFD = 2.0D0 / GM1D 181 : UTILG = UNG + COEFFG * AG 182 : UTILD = UND - COEFFD * AD 183 : DU = UTILG - UTILD 184 : IF(DU .GT. ZERO)THEN 185 : C 186 : C 187 : C******* Calcul des intersections des Invariants de Riemann 188 : C 189 : C******* Protection des PARAMETER NMAX, ZERO 190 : C N.B. Les subroutines pouvent les changer 191 : C 192 : NMAX0 = NMAX 193 : ZERO0 = ZERO 194 : C 194 : CALL RACHUS(NMAX0,ZERO0, 195 : & GAMG,COEFFG,AG,UTILG,PG, 196 : & GAMD,COEFFD,AD,UTILD,PD, 197 : & USTAR,PSTAR,LOGNC,MESERR,LOGAN) 198 : C 199 : C******* Pas de convergence de la mèthode de Newton-Rapson en RACHUS 200 : C 201 : IF(LOGNC) GOTO 9999 202 : C 203 : C******* Calcul des etats a Gauche et a Droite du champ linearment dégénéré 204 : C 205 : CELLG = (PSTAR/PG)**(1.0D0/GAMG) 206 : CELLD = (PSTAR/PD)**(1.0D0/GAMD) 207 : ASTARG=AG*CELLG**(1.0D0/COEFFG) 208 : ASTARD=AD*CELLD**(1.0D0/COEFFD) 209 : RSTARG=ROG*CELLG 210 : RSTARD=ROD*CELLD 211 : MSTARG = USTAR / ASTARG 212 : MSTARD = USTAR / ASTARD 213 : USTAR2 = USTAR*USTAR 214 : HTSG = (ASTARG * ASTARG) / GM1G + 0.5D0*(USTAR2+UTG2+UVG2) 215 : HTSD = (ASTARD * ASTARD) / GM1D + 0.5D0*(USTAR2+UTD2+UVD2) 216 : C 217 : C******* CFL 218 : C 219 : AMG = ABS(MSTARG) 220 : CELLT0 = (2.0D0*GAMG + (AMG*(3.0D0-GAMG)))/(GAMG+3.0D0) 221 : & /((AMG+1.0D0)*ASTARG) 222 : CELLT = MIN(CELLT,CELLT0) 223 : AMD = ABS(MSTARD) 224 : CELLT0 = (2.0D0*GAMD + (AMD*(3.0D0-GAMD)))/(GAMD+3.0D0) 225 : & /((AMD+1.0D0)*ASTARD) 226 : CELLT = MIN(CELLT,CELLT0) 227 : C 228 : C******* Partie antidiffusive du flux 229 : C 230 : IF(USTAR .GT. 0.0D0)THEN 230 : CALL FLUMV2(NESP, 231 : & GAMD,RSTARD,MSTARD,ASTARD,UTD,UVD,HTSD,YD,FLU2) 232 : C 233 : DO I1 = 1, NESP+5 234 : FLU1(I1) = FLU1(I1) - FLU2(I1) 235 : ENDDO 236 : C 236 : CALL FLUMV2(NESP, 237 : & GAMG,RSTARG,MSTARG,ASTARG,UTG,UVG,HTSG,YG,FLU2) 238 : C 239 : DO I1 = 1, NESP+5 240 : FLU1(I1) = FLU1(I1) + FLU2(I1) 241 : ENDDO 242 : ELSE 242 : CALL FLUPV2(NESP, 243 : & GAMD,RSTARD,MSTARD,ASTARD,UTD,UVD,HTSD,YD,FLU2) 244 : C 245 : DO I1 = 1, NESP+5 246 : FLU1(I1) = FLU1(I1) + FLU2(I1) 247 : ENDDO 248 : C 248 : CALL FLUPV2(NESP, 249 : & GAMG,RSTARG,MSTARG,ASTARG,UTG,UVG,HTSG,YG,FLU2) 250 : C 251 : DO I1 = 1, NESP+5 252 : FLU1(I1) = FLU1(I1) - FLU2(I1) 253 : ENDDO 254 : ENDIF 255 : C 256 : C******* Correction de Larrouturou pour la TVD des fractions massiques. 257 : C pour la vitesse tangentielle 258 : C 259 : C N.B. Pour le solver de Riemann exact cette propriete est 260 : C valide 261 : C Dans le cas de vide ce n'est pas necessaire: 262 : C VanLeer a la TVD des fractions massiques. 263 : C 264 : LOGI = LARCO 265 : IF(LOGI)THEN 266 : CELL = FLU1(1) 267 : IF(CELL .GT. 0.0D0)THEN 268 : FLU1(3) = CELL * UTG 269 : FLU1(4) = CELL * UVG 270 : DO I1 = 1, NESP 271 : FLU1(I1+5) = CELL * YG(I1) 272 : ENDDO 273 : ELSE 274 : FLU1(3) = CELL * UTD 275 : FLU1(4) = CELL * UVD 276 : DO I1 = 1, NESP 277 : FLU1(I1+5) = CELL * YD(I1) 278 : ENDDO 279 : ENDIF 280 : ENDIF 281 : ENDIF 282 : C 283 : 9999 CONTINUE 284 : RETURN 285 : END 286 : 287 : 288 : 289 : 290 : 291 : 292 : 293 : 294 : 295 :
© Cast3M 2003 - All rights reserved.
Disclaimer