1 : C FROEB3 SOURCE KLOCZKO 05/06/14 21:15:06 5111 2 : C FROEB3 SOURCE BECC 02/06/18 21:15:55 4365 3 : SUBROUTINE FROEB3(NESP, 4 : & GAMG,ROG,PG,UNG,UTG,UVG, 5 : & GAMD,ROD,PD,UND,UTD,UVD, 6 : & YG,YD,VINF,FLU, 7 : & CELLT) 8 : C 9 : C 10 : C PROJET : CASTEM 2000 11 : C 12 : C NOM : FROEB3 13 : C 14 : C DESCRIPTION : Formulation Volumes Finis pour les Equations 15 : C d'Euler Multi-Espèces relatives à un mélange 16 : C de gaz parfaits. 17 : C 18 : C Calcul du flux aux interfaces avec la méthode 19 : C de Roe-Turkel (cf. Thèse Cécile Viozat) adaptée 20 : C aux écoulements à faible nombre de Mach. 21 : C 22 : C LANGUAGE : FORTRAN 77 23 : C 24 : C AUTEUR : T. KLOCZKO DM2S/SFME/LTMF 25 : C 26 : C************************************************************************ 27 : C 28 : C APPELES : 29 : C 30 : C************************************************************************ 31 : C 32 : C**** Entrées: 33 : C 34 : C NESP = nombre d'especes considérées dans les Equations 35 : C d'Euler 36 : C 37 : C GAMG, GAMD = les "gamma" du gaz (gauche et droite) 38 : C 39 : C ROG, ROD = les densités 40 : C 41 : C PG, PD = les pressions 42 : C 43 : C UNG, UND = vitesses normales 44 : C 45 : C UTG, UTD = vitesses tangentielles 46 : C 47 : C UVG, UVD = vitesses tangentielles 48 : C 49 : C VINF = vitesse de cut-off 50 : C 51 : C**** Sorties: 52 : C 53 : C FLU = table du flux a l'interface dans le repère 54 : C (n,t), i.e. 55 : C (rho*un, rho*un*un + p, rho*un*ut, rho*un*uv, 56 : C rho*un*ht, rho*un*y1, ...) 57 : C 58 : C CELLT = condition de stabilité, i.e. 59 : C 60 : C dT/diamax < cellt 61 : C 62 : C************************************************************************ 63 : C 64 : C HISTORIQUE (Anomalies et modifications éventuelles) 65 : C 66 : C HISTORIQUE : Créé le 02.05.05 67 : C 68 : C 69 : C 70 : C************************************************************************ 71 : C 72 : C N.B.: Toutes les variables sont DECLAREES 73 : C 74 : C 75 : C IMPLICIT NONE 76 : C 77 : C*** Déclaration des indices de boucles et du nombre d'espèces 78 : INTEGER NESP,I1,I2,I3 79 : C 80 : C*** Déclaration des variables primitives 81 : REAL*8 ROG,PG,UNG,UTG,UVG,RETG,GAMG, 82 : & ROD,PD,UND,UTD,UVD,RETD,GAMD 83 : C 84 : C*** Déclaration des variables de cacul des incréments conservatifs 85 : REAL*8 XP(5) 86 : C 87 : C*** Déclaration des variables de l'état moyen de Roe 88 : REAL*8 SQROG,SQROD,ROM,ROS, 89 : & UNM,UTM,UVM, 90 : & HG,HD,HM, 91 : & CM,GAMM,GAM1 92 : C 93 : C*** Déclaration de la condition de stabilité 94 : REAL*8 CELLT 95 : C 96 : C*** Délaration des variables de calcul du Mach de référence 97 : REAL*8 CG,CD,PHIM,VINF 98 : C 99 : C*** Déclaration des variables de calcul des valeurs propres 100 : REAL*8 VP(5),XU 101 : C 102 : C*** Déclaration des variables de calcul des matrices de passage 103 : REAL*8 R,S,TU,Q2, 104 : & MPI(5,5),MP(5,5) 105 : C 106 : C*** Déclaration des variables de calcul de la dissipation 107 : REAL*8 VDIS(5),DS(5) 108 : C 109 : C*** Déclaration des variables de calcul du flux numérique 110 : REAL*8 FC(5),FLU(*) 111 : C 112 : C*** Déclaration de la partie multiespèce 113 : REAL*8 YG(*),YD(*),FLUM 114 : C 115 : C 116 : C 117 : C*** DEBUT DES CALCULS *** 118 : C 119 : C 120 : C**** We include in the cut-off UNG,UND,UTG,UTD in order to 121 : C avoid low diffusivity in stagnation regions 122 : C 123 : VINF=MAX(VINF,((UNG**2+UTG**2+UVG**2)**0.5D0)) 124 : VINF=MAX(VINF,((UND**2+UTD**2+UVD**2)**0.5D0)) 125 : C 126 : C**** Calcul des Deltas XP 127 : C 128 : RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG) 129 : & + (0.5D0*ROG*(UNG**2+UTG**2+UVG**2)) 130 : RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD) 131 : & + (0.5D0*ROD*(UND**2+UTD**2+UVD**2)) 132 : C 133 : XP(1) = ROG - ROD 134 : XP(2) = (ROG * UNG) - (ROD * UND) 135 : XP(3) = (ROG * UTG) - (ROD * UTD) 136 : XP(4) = (ROG * UVG) - (ROD * UVD) 137 : XP(5) = RETG - RETD 138 : C 139 : C*** Calcul de l'état moyen de Roe 140 : C 141 : GAMM = SQRT(GAMG) * SQRT(GAMD) 142 : GAM1 = GAMM - 1.D0 143 : C 144 : SQROG = SQRT(ROG) 145 : SQROD = SQRT(ROD) 146 : ROM = SQROG * SQROD 147 : ROS = SQROG + SQROD 148 : C 149 : UNM = ((SQROG * UNG) + (SQROD * UND)) / ROS 150 : UTM = ((SQROG * UTG) + (SQROD * UTD)) / ROS 151 : UVM = ((SQROG * UVG) + (SQROD * UVD)) / ROS 152 : C 153 : HG = (1.D0 / ROG) * (RETG + PG) 154 : HD = (1.D0 / ROD) * (RETD + PD) 155 : HM = ((SQROG * HG) + (SQROD * HD)) / ROS 156 : C 157 : CM = SQRT(GAM1 * (HM - 0.5D0 * (UNM**2 + UTM**2+UVM**2))) 158 : C 159 : C*** Calcul de la condition de stabilité 160 : C 161 : CELLT = 1.D0 / (CM + ABS(UNM)) 162 : C 163 : C*** Calcul du nombre de Mach de référence PHIM 164 : C 165 : CG = SQRT((GAMG * PG) / ROG) 166 : CD = SQRT((GAMD * PD) / ROD) 167 : IF(VINF.GE.MAX(CD,CG))THEN 168 : PHIM = 1.D0 169 : ELSE 170 : PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD)) 171 : ENDIF 172 : C 173 : C*** Calcul des valeurs propres 174 : C 175 : XU = ((1.D0 - PHIM**2) * UNM)**2 + 4.D0 * (PHIM * CM)**2 176 : C 177 : VP(1) = UNM 178 : VP(2) = UNM 179 : VP(3) = UNM 180 : VP(4) = 0.5D0 * ((1.D0 + PHIM**2) * UNM + SQRT(XU)) 181 : VP(5) = 0.5D0 * ((1.D0 + PHIM**2) * UNM - SQRT(XU)) 182 : C 183 : C*** Calcul des matrices de passages 184 : C 185 : R = VP(4) - UNM * PHIM**2 186 : S = VP(5) - UNM * PHIM**2 187 : TU = 0.5D0 * (VP(5) - VP(4)) 188 : Q2 = 0.5D0 * (UNM**2 + UTM**2+UVM**2) 189 : C 190 : MPI(1,1) = 1.D0 - GAM1 * Q2 / (CM**2) 191 : MPI(1,2) = GAM1 * UNM / (CM**2) 192 : MPI(1,3) = GAM1 * UTM / (CM**2) 193 : MPI(1,4) = GAM1 * UVM / (CM**2) 194 : MPI(1,5) = -GAM1 / (CM**2) 195 : MPI(2,1) = UVM 196 : MPI(2,2) = 0.D0 197 : MPI(2,3) = 0.D0 198 : MPI(2,4) = -1.D0 199 : MPI(2,5) = 0.D0 200 : MPI(3,1) = UTM 201 : MPI(3,2) = 0.D0 202 : MPI(3,3) = -1.D0 203 : MPI(3,4) = 0.D0 204 : MPI(3,5) = 0.D0 205 : MPI(4,1) = (S * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU 206 : MPI(4,2) = -(S * UNM * GAM1 + (PHIM*CM)**2) / TU 207 : MPI(4,3) = -(S * UTM * GAM1) / TU 208 : MPI(4,4) = -(S * UVM * GAM1) / TU 209 : MPI(4,5) = S * GAM1 /TU 210 : MPI(5,1) = -(R * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU 211 : MPI(5,2) = (R * UNM * GAM1 + (PHIM*CM)**2) / TU 212 : MPI(5,3) = (R * UTM * GAM1) / TU 213 : MPI(5,4) = (R * UVM * GAM1) / TU 214 : MPI(5,5) = -R * GAM1 / TU 215 : C 216 : MP(1,1) = 1.D0 217 : MP(1,2) = 0.D0 218 : MP(1,3) = 0.D0 219 : MP(1,4) = 0.5D0 / (PHIM*CM)**2 220 : MP(1,5) = 0.5D0 / (PHIM*CM)**2 221 : MP(2,1) = UNM 222 : MP(2,2) = 0.D0 223 : MP(2,3) = 0.D0 224 : MP(2,4) = 0.5D0 * (UNM + R) / (PHIM*CM)**2 225 : MP(2,5) = 0.5D0 * (UNM + S) / (PHIM*CM)**2 226 : MP(3,1) = UTM 227 : MP(3,2) = 0.D0 228 : MP(3,3) = -1.D0 229 : MP(3,4) = 0.5D0 * UTM / (PHIM*CM)**2 230 : MP(3,5) = 0.5D0 * UTM / (PHIM*CM)**2 231 : MP(4,1) = UVM 232 : MP(4,2) = -1.D0 233 : MP(4,3) = 0.D0 234 : MP(4,4) = 0.5D0 * UVM / (PHIM*CM)**2 235 : MP(4,5) = 0.5D0 * UVM / (PHIM*CM)**2 236 : MP(5,1) = Q2 237 : MP(5,2) = -UVM 238 : MP(5,3) = -UTM 239 : MP(5,4) = 0.5D0 * (HM + R * UNM) / (PHIM*CM)**2 240 : MP(5,5) = 0.5D0 * (HM + S * UNM) / (PHIM*CM)**2 241 : C 242 : C*** Calcul de la dissipation 243 : C 244 : DO I1 = 1,5 245 : VDIS(I1) = 0.5D0 * ABS(VP(I1)) 246 : ENDDO 247 : C 248 : DO I1 = 1,5 249 : DS(I1) = 0.D0 250 : DO I2 = 1,5 251 : DO I3 = 1,5 252 : DS(I1) = DS(I1) + MP(I1,I2) * VDIS(I2) * MPI(I2,I3) * XP(I3) 253 : ENDDO 254 : ENDDO 255 : ENDDO 256 : C 257 : C*** Calcul du flux convectif 258 : C 259 : FC(1) = 0.5D0 * ((ROG * UNG) + (ROD * UND)) 260 : FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD)) 261 : FC(3) = 0.5D0 * ((ROG * UNG * UTG) + (ROD * UND * UTD)) 262 : FC(4) = 0.5D0 * ((ROG * UNG * UVG) + (ROD * UND * UVD)) 263 : FC(5) = 0.5D0 * ((UNG * (RETG + PG)) + (UND * (RETD + PD))) 264 : C 265 : C*** Calcul du flux numérique 266 : C 267 : DO I1 = 1,5 268 : FLU(I1) = FC(I1) + DS(I1) 269 : ENDDO 270 : C 271 : C*** Partie multiespèce 272 : C 273 : FLUM = FLU(1) 274 : IF(FLUM .GT. 0)THEN 275 : DO I1 = 1, NESP, 1 276 : FLU(5+I1)=FLUM * YG(I1) 277 : ENDDO 278 : ELSE 279 : DO I1 = 1, NESP, 1 280 : FLU(5+I1)=FLUM * YD(I1) 281 : ENDDO 282 : ENDIF 283 : C 284 : RETURN 285 : END 286 : 287 :
© Cast3M 2003 - All rights reserved.
Disclaimer