1 : C FROEBM SOURCE KLOCZKO 05/06/14 21:15:08 5111 2 : C FROEBM SOURCE BECC 02/06/18 21:15:55 4365 3 : SUBROUTINE FROEBM(NESP, 4 : & GAMG,ROG,PG,UNG,UTG, 5 : & GAMD,ROD,PD,UND,UTD, 6 : & YG,YD,VINF,FLU, 7 : & CELLT) 8 : C 9 : C 10 : C PROJET : CASTEM 2000 11 : C 12 : C NOM : FROEBM 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 VINF = vitesse de cut-off 48 : C 49 : C**** Sorties: 50 : C 51 : C FLU = table du flux a l'interface dans le repaire 52 : C (n,t), i.e. 53 : C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht, 54 : C rho*un*y1, ...) 55 : C 56 : C CELLT = condition de stabilité, i.e. 57 : C 58 : C dT/diamax < cellt 59 : C 60 : C************************************************************************ 61 : C 62 : C HISTORIQUE (Anomalies et modifications éventuelles) 63 : C 64 : C HISTORIQUE : 29.04.05 65 : C 66 : C 67 : C 68 : C************************************************************************ 69 : C 70 : C N.B.: Toutes les variables sont DECLAREES 71 : C 72 : C 73 : C IMPLICIT NONE 74 : C 75 : C*** Déclaration des indices de boucles et du nombre d'espèces 76 : INTEGER NESP,I1,I2,I3 77 : C 78 : C*** Déclaration des variables primitives 79 : REAL*8 ROG,PG,UNG,UTG,RETG,GAMG, 80 : & ROD,PD,UND,UTD,RETD,GAMD 81 : C 82 : C*** Déclaration des variables de cacul des incréments conservatifs 83 : REAL*8 XP(4) 84 : C 85 : C*** Déclaration des variables de l'état moyen de Roe 86 : REAL*8 SQROG,SQROD,ROM,ROS, 87 : & UNM,UTM, 88 : & HG,HD,HM, 89 : & CM,GAMM,GAM1 90 : C 91 : C*** Déclaration de la condition de stabilité 92 : REAL*8 CELLT 93 : C 94 : C*** Délaration des variables de calcul du Mach de référence 95 : REAL*8 CG,CD,PHIM,VINF 96 : C 97 : C*** Déclaration des variables de calcul des valeurs propres 98 : REAL*8 VP(4),XU 99 : C 100 : C*** Déclaration des variables de calcul des matrices de passage 101 : REAL*8 R,S,TU,Q2, 102 : & MPI(4,4),MP(4,4) 103 : C 104 : C*** Déclaration des variables de calcul de la dissipation 105 : REAL*8 VDIS(4),DS(4) 106 : C 107 : C*** Déclaration des variables de calcul du flux numérique 108 : REAL*8 FC(4),FLU(*) 109 : C 110 : C*** Déclaration de la partie multiespèce 111 : REAL*8 YG(*),YD(*),FLUM 112 : C 113 : C 114 : C 115 : C*** DEBUT DES CALCULS *** 116 : C 117 : C 118 : C**** We include in the cut-off UNG,UND,UTG,UTD in order to 119 : C avoid low diffusivity in stagnation regions 120 : C 121 : VINF=MAX(VINF,((UNG*UNG+UTG*UTG)**0.5D0)) 122 : VINF=MAX(VINF,((UND*UND+UTD*UTD)**0.5D0)) 123 : C 124 : C**** Calcul des Deltas XP 125 : C 126 : RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG) 127 : & + (0.5D0*ROG*((UNG*UNG)+(UTG*UTG))) 128 : RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD) 129 : & + (0.5D0*ROD*((UND*UND)+(UTD*UTD))) 130 : C 131 : XP(1) = ROG - ROD 132 : XP(2) = (ROG * UNG) - (ROD * UND) 133 : XP(3) = (ROG * UTG) - (ROD * UTD) 134 : XP(4) = RETG - RETD 135 : C 136 : C*** Calcul de l'état moyen de Roe 137 : C 138 : GAMM = SQRT(GAMG) * SQRT(GAMD) 139 : GAM1 = GAMM - 1.D0 140 : C 141 : SQROG = SQRT(ROG) 142 : SQROD = SQRT(ROD) 143 : ROM = SQROG * SQROD 144 : ROS = SQROG + SQROD 145 : C 146 : UNM = ((SQROG * UNG) + (SQROD * UND)) / ROS 147 : UTM = ((SQROG * UTG) + (SQROD * UTD)) / ROS 148 : C 149 : HG = (1.D0 / ROG) * (RETG + PG) 150 : HD = (1.D0 / ROD) * (RETD + PD) 151 : HM = ((SQROG * HG) + (SQROD * HD)) / ROS 152 : C 153 : CM = SQRT(GAM1 * (HM - 0.5D0 * (UNM**2 + UTM**2))) 154 : C 155 : C*** Calcul de la condition de stabilité 156 : C 157 : CELLT = 1.D0 / (CM + ABS(UNM)) 158 : C 159 : C*** Calcul du nombre de Mach de référence PHIM 160 : C 161 : CG = SQRT((GAMG * PG) / ROG) 162 : CD = SQRT((GAMD * PD) / ROD) 163 : IF(VINF.GE.MAX(CD,CG))THEN 164 : PHIM = 1.D0 165 : ELSE 166 : PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD)) 167 : ENDIF 168 : C 169 : C*** Calcul des valeurs propres 170 : C 171 : XU = ((1.D0 - PHIM**2) * UNM)**2 + 4.D0 * (PHIM * CM)**2 172 : C 173 : VP(1) = UNM 174 : VP(2) = UNM 175 : VP(3) = 0.5D0 * ((1.D0 + PHIM**2) * UNM + SQRT(XU)) 176 : VP(4) = 0.5D0 * ((1.D0 + PHIM**2) * UNM - SQRT(XU)) 177 : C 178 : C*** Calcul des matrices de passages 179 : C 180 : R = VP(3) - UNM * PHIM**2 181 : S = VP(4) - UNM * PHIM**2 182 : TU = 0.5D0 * (VP(4) - VP(3)) 183 : Q2 = 0.5D0 * (UNM**2 + UTM**2) 184 : C 185 : MPI(1,1) = 1.D0 - GAM1 * Q2 / (CM**2) 186 : MPI(1,2) = GAM1 * UNM / (CM**2) 187 : MPI(1,3) = GAM1 * UTM / (CM**2) 188 : MPI(1,4) = -GAM1 / (CM**2) 189 : MPI(2,1) = UTM 190 : MPI(2,2) = 0.D0 191 : MPI(2,3) = -1.D0 192 : MPI(2,4) = 0.D0 193 : MPI(3,1) = (S * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU 194 : MPI(3,2) = -(S * UNM * GAM1 + (PHIM*CM)**2) / TU 195 : MPI(3,3) = -(S * UTM * GAM1) / TU 196 : MPI(3,4) = S * GAM1 /TU 197 : MPI(4,1) = -(R * Q2 * GAM1 + UNM * (PHIM*CM)**2) / TU 198 : MPI(4,2) = (R * UNM * GAM1 + (PHIM*CM)**2) / TU 199 : MPI(4,3) = (R * UTM * GAM1) / TU 200 : MPI(4,4) = -R * GAM1 / TU 201 : C 202 : MP(1,1) = 1.D0 203 : MP(1,2) = 0.D0 204 : MP(1,3) = 0.5D0 / (PHIM*CM)**2 205 : MP(1,4) = 0.5D0 / (PHIM*CM)**2 206 : MP(2,1) = UNM 207 : MP(2,2) = 0.D0 208 : MP(2,3) = 0.5D0 * (UNM + R) / (PHIM*CM)**2 209 : MP(2,4) = 0.5D0 * (UNM + S) / (PHIM*CM)**2 210 : MP(3,1) = UTM 211 : MP(3,2) = -1.D0 212 : MP(3,3) = 0.5D0 * UTM / (PHIM*CM)**2 213 : MP(3,4) = 0.5D0 * UTM / (PHIM*CM)**2 214 : MP(4,1) = Q2 215 : MP(4,2) = -UTM 216 : MP(4,3) = 0.5D0 * (HM + R * UNM) / (PHIM*CM)**2 217 : MP(4,4) = 0.5D0 * (HM + S * UNM) / (PHIM*CM)**2 218 : C 219 : C*** Calcul de la dissipation 220 : C 221 : DO I1 = 1,4 222 : VDIS(I1) = 0.5D0 * ABS(VP(I1)) 223 : ENDDO 224 : C 225 : DO I1 = 1,4 226 : DS(I1) = 0.D0 227 : DO I2 = 1,4 228 : DO I3 = 1,4 229 : DS(I1) = DS(I1) + MP(I1,I2) * VDIS(I2) * MPI(I2,I3) * XP(I3) 230 : ENDDO 231 : ENDDO 232 : ENDDO 233 : C 234 : C*** Calcul du flux convectif 235 : C 236 : FC(1) = 0.5D0 * ((ROG * UNG) + (ROD * UND)) 237 : FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD)) 238 : FC(3) = 0.5D0 * ((ROG * UNG * UTG) + (ROD * UND * UTD)) 239 : FC(4) = 0.5D0 * ((UNG * (RETG + PG)) + (UND * (RETD + PD))) 240 : C 241 : C*** Calcul du flux numérique 242 : C 243 : DO I1 = 1,4 244 : FLU(I1) = FC(I1) + DS(I1) 245 : ENDDO 246 : C 247 : C*** Partie multiespèce 248 : C 249 : FLUM = FLU(1) 250 : IF(FLUM .GT. 0)THEN 251 : DO I1 = 1, NESP, 1 252 : FLU(4+I1)=FLUM * YG(I1) 253 : ENDDO 254 : ELSE 255 : DO I1 = 1, NESP, 1 256 : FLU(4+I1)=FLUM * YD(I1) 257 : ENDDO 258 : ENDIF 259 : C 260 : RETURN 261 : END 262 : 263 :
© Cast3M 2003 - All rights reserved.
Disclaimer