1 : C FRUSB3 SOURCE KLOCZKO 05/06/14 21:15:14 5111 2 : C FRUSB3 SOURCE CHAT 05/01/13 00:09:18 5004 3 : SUBROUTINE FRUSB3(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 PROJET : CASTEM 2000 10 : C 11 : C NOM : FRUSB3 12 : C 13 : C DESCRIPTION : Formulation Volumes Finis pour les Equations 14 : C d'Euler Multi-Especes relatives à un melange 15 : C de gaz parfaits. 16 : C 17 : C Calcul du flux aux interfaces avec la methode 18 : C de Rusanov preconditionée 19 : C 20 : C LANGUAGE : FORTRAN 77 21 : C 22 : C AUTEUR : T. KLOCZKO DM2S/SFME/LTMF 23 : C 24 : C************************************************************************ 25 : C 26 : C APPELES : 27 : C 28 : C************************************************************************ 29 : C 30 : C**** Entrées: 31 : C 32 : C NESP = nombre d'especes considérées dans les Equations 33 : C d'Euler 34 : C 35 : C GAMG, GAMD = les "gamma" du gaz (gauche et droite) 36 : C 37 : C ROG, ROD = les densités 38 : C 39 : C PG, PD = les pressions 40 : C 41 : C UNG, UND = vitesses normales 42 : C 43 : C UTG, UTD = vitesses tangentielles 44 : C 45 : C UVG, UVD = 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*ut, 54 : C rho*un*ht, 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 : Créé le 02.05.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 : IMPLICIT INTEGER(I-N) 74 : INTEGER NESP,I1,I2 75 : REAL*8 GAMG,ROG,PG,UNG,UTG,UVG,HG,RETG,CG 76 : & ,GAMD,ROD,PD,UND,UTD,UVD,HD,RETD,CD 77 : & ,YG(*),YD(*),VINF 78 : & ,FLU(*),CELLT, FLUM 79 : & ,UNF,UTF,UVF,ROF,PF,GAMF,HF,VPMAX,PHIM 80 : & ,AMAT(5,5),XP(5),FC(5),DS(5) 81 : & ,GAM1,XU,CF,Q2,COEF 82 : C 83 : C**** Etat moyen 84 : C 85 : UNF = 0.5D0 * (UNG + UND) 86 : UTF = 0.5D0 * (UTG + UTD) 87 : UVF = 0.5D0 * (UVG + UVD) 88 : ROF = 0.5D0 * (ROG + ROD) 89 : PF = 0.5D0 * (PG + PD) 90 : GAMF = 0.5D0 * (GAMG + GAMD) 91 : GAM1 = GAMF - 1.D0 92 : CF = (GAMF * PF / ROF)**0.5D0 93 : C 94 : C**** We include in the cut-off UNG,UND,UTG,UTD in order to 95 : C avoid low diffusivity in stagnation regions 96 : C 97 : VINF = MAX(VINF,((UNG**2 + UTG**2 + UVG**2)**0.5D0)) 98 : VINF = MAX(VINF,((UND**2 + UTD**2 + UVD**2)**0.5D0)) 99 : C 100 : C*** Calcul de la condition de stabilité 101 : C 102 : CELLT = 1.0D0 / (CF + ABS(UNF)) 103 : C 104 : C**** DELTA 105 : C 106 : RETG = ((1.0D0 / (GAMG - 1.0D0)) * PG) 107 : & + (0.5D0 * ROG * (UNG**2 + UTG**2 + UVG**2)) 108 : RETD = ((1.0D0 / (GAMD - 1.0D0)) * PD) 109 : & + (0.5D0 * ROD * (UND**2 + UTD**2 + UVD**2)) 110 : C 111 : HG = (1.D0 / ROG) * (RETG + PG) 112 : HD = (1.D0 / ROD) * (RETD + PD) 113 : HF = 0.5D0 * (HG + HD) 114 : C 115 : XP(1) = ROG - ROD 116 : XP(2) = (ROG * UNG) - (ROD * UND) 117 : XP(3) = (ROG * UTG) - (ROD * UTD) 118 : XP(4) = (ROG * UVG) - (ROD * UVD) 119 : XP(5) = RETG - RETD 120 : C 121 : C*** Calcul du nombre de Mach de référence PHIM 122 : C 123 : CG = SQRT((GAMG * PG) / ROG) 124 : CD = SQRT((GAMD * PD) / ROD) 125 : IF(VINF.GE.MAX(CD,CG))THEN 126 : PHIM = 1.D0 127 : ELSE 128 : PHIM = 0.5D0 * VINF * ((1.D0 / CG) + (1.D0 / CD)) 129 : ENDIF 130 : C 131 : C*** Calcul des valeurs propres 132 : C 133 : XU = ((1.D0 - PHIM**2) * UNF)**2 + 4.D0 * (PHIM * CF)**2 134 : C 135 : VPMAX = 0.5D0 * ((1.D0 + PHIM**2) * ABS(UNF) + SQRT(XU)) 136 : C 137 : C 138 : C*** Calcul de la matrice de préconditionnement 139 : C 140 : Q2 = 0.5D0 * (UNF**2 + UTF**2 + UVF**2) 141 : COEF = (1.0D0 / PHIM**2 - 1.0D0) * GAM1 / CF**2 142 : C 143 : AMAT(1,1) = COEF * 1.D0 * Q2 + 1.D0 144 : AMAT(1,2) = COEF * 1.D0 * (-UNF) 145 : AMAT(1,3) = COEF * 1.D0 * (-UTF) 146 : AMAT(1,4) = COEF * 1.D0 * (-UVF) 147 : AMAT(1,5) = COEF * 1.D0 * 1.D0 148 : C 149 : AMAT(2,1) = COEF * UNF * Q2 150 : AMAT(2,2) = COEF * UNF * (-UNF) + 1.D0 151 : AMAT(2,3) = COEF * UNF * (-UTF) 152 : AMAT(2,4) = COEF * UNF * (-UVF) 153 : AMAT(2,5) = COEF * UNF * 1.D0 154 : C 155 : AMAT(3,1) = COEF * UTF * Q2 156 : AMAT(3,2) = COEF * UTF * (-UNF) 157 : AMAT(3,3) = COEF * UTF * (-UTF) + 1.D0 158 : AMAT(3,4) = COEF * UTF * (-UVF) 159 : AMAT(3,5) = COEF * UTF * 1.D0 160 : C 161 : AMAT(4,1) = COEF * UVF * Q2 162 : AMAT(4,2) = COEF * UVF * (-UNF) 163 : AMAT(4,3) = COEF * UVF * (-UTF) 164 : AMAT(4,4) = COEF * UVF * (-UVF) + 1.D0 165 : AMAT(4,5) = COEF * UVF * 1.D0 166 : C 167 : AMAT(5,1) = COEF * HF * Q2 168 : AMAT(5,2) = COEF * HF * (-UNF) 169 : AMAT(5,3) = COEF * HF * (-UTF) 170 : AMAT(5,4) = COEF * HF * (-UVF) 171 : AMAT(5,5) = COEF * HF * 1.D0 + 1.D0 172 : C 173 : C*** Calcul de la dissipation 174 : C 175 : DO I1 = 1,5 176 : DS(I1) = 0.D0 177 : DO I2 = 1,5 178 : DS(I1) = DS(I1) + 0.5D0 * VPMAX * AMAT(I1,I2) * XP(I2) 179 : ENDDO 180 : ENDDO 181 : C 182 : C*** Calcul du flux convectif 183 : C 184 : FC(1) = 0.5D0 * ((ROG * UNG) + (ROD * UND)) 185 : FC(2) = 0.5D0 * ((ROG * UNG * UNG + PG) + (ROD * UND * UND + PD)) 186 : FC(3) = 0.5D0 * ((ROG * UNG * UTG) + (ROD * UND * UTD)) 187 : FC(4) = 0.5D0 * ((ROG * UNG * UVG) + (ROD * UND * UVD)) 188 : FC(5) = 0.5D0 * ((UNG * (RETG + PG)) + (UND * (RETD + PD))) 189 : C 190 : C*** Calcul du flux numérique 191 : C 192 : DO I1 = 1,5 193 : FLU(I1) = FC(I1) + DS(I1) 194 : ENDDO 195 : C 196 : C*** Partie multi-espèces 197 : C 198 : FLUM = FLU(1) 199 : IF(FLUM .GT. 0)THEN 200 : DO I1 = 1, NESP, 1 201 : FLU(5+I1)=FLUM * YG(I1) 202 : ENDDO 203 : ELSE 204 : DO I1 = 1, NESP, 1 205 : FLU(5+I1)=FLUM * YD(I1) 206 : ENDDO 207 : ENDIF 208 : C 209 : RETURN 210 : END 211 : 212 : 213 :
© Cast3M 2003 - All rights reserved.
Disclaimer