1 : C FHLLCL SOURCE BECC 10/05/05 21:15:09 6674 2 : SUBROUTINE FHLLCL(NESP, 3 : & GAMG,ROG,PG,UNG,UTG, 4 : & GAMD,ROD,PD,UND,UTD, 5 : & YG,YD,VINF0,FLU1, 6 : & CELLT) 7 : C 8 : C************************************************************************ 9 : C 10 : C PROJET : CASTEM 2000 11 : C 12 : C NOM : FHLLCL 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 methode 19 : C de HLLC (avec gamma = const) 20 : C +Low Mach preconditioning 21 : C 22 : C Voir: * Batten et al. SIAM J. Sci. Comp. 23 : C Vol. 18 pp 1553-6, 1997 24 : C * Batten et al. JCP, 1997 25 : C Vol. 137 pp 38-- 26 : C * Luo et al., AIIA Journal, 2005 27 : C Vol. 43, 6, pp 1160-- 28 : C 29 : C N.B.: SPLIT DE FRACTIONS MASSIQUES ET UNT A LA 30 : C LARROUTUROU 31 : C 32 : C LANGUAGE : FORTRAN 77 33 : C 34 : C AUTEUR : A. BECCANTINI DRN/DMT/SEMT/LTMF 35 : C 36 : C************************************************************************ 37 : C 38 : C APPELES 39 : C 40 : C************************************************************************ 41 : C 42 : C**** Entrées: 43 : C 44 : C NESP = nombre d'especes considérées dans les Equations 45 : C d'Euler 46 : C 47 : C GAMG, GAMD = les "gamma" du gaz (gauche et droite) 48 : C 49 : C ROG, ROD = les densités 50 : C 51 : C PG, PD = les pressions 52 : C 53 : C UNG, UND = vitesses normales 54 : C 55 : C UTG, UTD = vitesses tangentielles 56 : C 57 : C YG, YD = tables des fractiones massiques 58 : C 59 : C**** Sorties: 60 : C 61 : C FLU1 = table du flux a l'interface dans le repaire 62 : C (n,t), i.e. 63 : C (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht, 64 : C rho*un*y1, ...) 65 : C 66 : C CELLT = condition de stabilité, i.e. 67 : C 68 : C dT/diamax < cellt (temps/longueur) 69 : C 70 : C************************************************************************ 71 : C 72 : C HISTORIQUE (Anomalies et modifications éventuelles) 73 : C 74 : C HISTORIQUE : Créé le 11.09.2000 75 : C 76 : C************************************************************************ 77 : C 78 : C N.B.: Toutes les variables sont DECLAREES 79 : C 80 : C It exactly captures: 81 : C 82 : C Stationary shocks (test for instance 83 : C 84 : C GAMG=1.4D0 85 : C GAMD=1.4D0 86 : C ROG=1.0D0 87 : C ROD=4.5584471D0 88 : C PG=1.0D5 89 : C PD=0.18279373D7 90 : C UNG=0.14877919D4 91 : C UND=0.32638130D3 92 : C UTG=0.0D0 93 : C UTD=0.0D0 94 : C ) and stationary contacts 95 : C 96 : IMPLICIT INTEGER(I-N) 97 : INTEGER NESP 98 : & , IESP 99 : REAL*8 GAMG,ROG,PG,UNG,UTG 100 : & ,GAMD,ROD,PD,UND,UTD 101 : & ,YG(*),YD(*),FLU1(*),CELLT,VINF0 102 : & ,HALF 103 : PARAMETER (HALF=0.5D0) 104 : REAL*8 USGGM1,USGDM1,HTG,HTD,REING,REIND,SQRG,SQRD,USOSQR 105 : & ,HTAVE,UNAVE,REIAVE,USGAVE,GAMAVE,ROAVE,PAVE,ASONAV 106 : & ,ASONG,ASOND 107 : & ,VINF,PHIG,PHID,PHIAV,XU,SG,SGS,SD,SDS 108 : & ,SM,PSTAR,RETG,RETD,OMEGA 109 : & ,ROGS,RUNGS,RETGS,RODS,RUNDS,RETDS 110 : C 111 : USGGM1=1.0D0/(GAMG-1.0D0) 112 : USGDM1=1.0D0/(GAMD-1.0D0) 113 : C 114 : C Evaluation of the Roe average state to compute the non-linear wave 115 : C speeds (According to Shuye, JCP 142, pag. 217--, 1998 116 : C 117 : REING=PG*USGGM1 118 : REIND=PD*USGDM1 119 : HTG=(HALF*UNG*UNG)+(GAMG*REING/ROG) 120 : HTD=(HALF*UND*UND)+(GAMD*REIND/ROD) 121 : SQRG=ROG**0.5D0 122 : SQRD=ROD**0.5D0 123 : USOSQR=1.0D0/(SQRG+SQRD) 124 : HTAVE=((SQRG*HTG)+(SQRD*HTD))*USOSQR 125 : UNAVE=((SQRG*UNG)+(SQRD*UND))*USOSQR 126 : USGAVE=((SQRG*USGGM1)+(SQRD*USGDM1))*USOSQR 127 : GAMAVE=(1.0D0/USGAVE)+1.0D0 128 : REIAVE=((SQRG*REING)+(SQRD*REIND))*USOSQR 129 : C 130 : ROAVE=REIAVE*GAMAVE/(HTAVE-(HALF*UNAVE*UNAVE)) 131 : PAVE=REIAVE*(GAMAVE-1.0D0) 132 : C 133 : ASONAV=(GAMAVE*PAVE/ROAVE)**0.5D0 134 : C 135 : C**** Evaluation of the speed of the Genuinely-Non-Linear Waves 136 : C for the Preconditioned System 137 : C 138 : ASONG=(GAMG*PG/ROG)**0.5D0 139 : ASOND=(GAMD*PD/ROD)**0.5D0 140 : C 141 : C**** We include in the cut-off UNG,UND,UTG,UTD in order to 142 : C avoid low diffusivity in stagnation regions 143 : C 144 : VINF=MAX(VINF0,((UNG*UNG+UTG*UTG)**0.5D0)) 145 : VINF=MAX(VINF0,((UND*UND+UTD*UTD)**0.5D0)) 146 : C 147 : C*** Calcul des nombres de Mach de référence 148 : C PHIG, PHID, PHIAV 149 : C 150 : IF(VINF.GE.ASONG)THEN 151 : PHIG = 1.D0 152 : ELSE 153 : PHIG = VINF/ASONG 154 : ENDIF 155 : IF(VINF.GE.ASOND)THEN 156 : PHID = 1.D0 157 : ELSE 158 : PHID = VINF/ASOND 159 : ENDIF 160 : IF(VINF.GE.ASONAV)THEN 161 : PHIAV = 1.D0 162 : ELSE 163 : PHIAV = VINF/ASONAV 164 : ENDIF 165 : C 166 : C*** Calcul des valeurs propres 167 : C 168 : XU = (((1.D0 - PHIG**2) * UNG)**2) + (4.D0 * ((PHIG * ASONG)**2)) 169 : SG= ((1.D0 + (PHIG**2)) * UNG) - (XU**0.5D0) 170 : SG= SG*HALF 171 : C write(*,*) SG,(UNG-ASONG) 172 : C 173 : XU = (((1.D0 - PHIAV**2) * UNAVE)**2) + 174 : & (4.D0 * ((PHIAV * ASONAV)**2)) 175 : SGS= ((1.D0 + (PHIAV**2)) * UNAVE) - (XU**0.5D0) 176 : SGS= SGS*HALF 177 : C write(*,*) SGS,(UNAVE-ASONAV) 178 : C 179 : SG=MIN(SG,SGS) 180 : C 181 : SDS= ((1.D0 + (PHIAV**2)) * UNAVE) + (XU**0.5D0) 182 : SDS= SDS*HALF 183 : C write(*,*) SDS,(UNAVE+ASONAV) 184 : C 185 : XU = (((1.D0 - PHID**2) * UND)**2) + (4.D0 * ((PHID * ASOND)**2)) 186 : SD= ((1.D0 + (PHID**2)) * UND) + (XU**0.5D0) 187 : SD= SD*HALF 188 : C write(*,*) SD,(UND+ASOND) 189 : C 190 : SD=MAX(SD,SDS) 191 : C 192 : C**** Speed of the contact discontinuity 193 : C 194 : SM=ROD*UND*(SD-UND)-ROG*UNG*(SG-UNG)+PG-PD 195 : SM=SM/((ROD*(SD-UND))-(ROG*(SG-UNG))) 196 : C 197 : C**** Interfacial pressure 198 : C 199 : PSTAR=ROG*(UNG-SG)*(UNG-SM)+PG 200 : RETG=(USGGM1*PG)+(HALF*ROG*((UNG*UNG)+(UTG*UTG))) 201 : RETD=(USGDM1*PD)+(HALF*ROD*((UND*UND)+(UTD*UTD))) 202 : C 203 : IF(SG .GT. 0)THEN 204 : FLU1(1)=ROG*UNG 205 : FLU1(2)=ROG*UNG*UNG+PG 206 : FLU1(4)=UNG*(RETG+PG) 207 : ELSEIF(SM .GT. 0)THEN 208 : OMEGA=1.0D0/(SG-SM) 209 : ROGS=ROG*(SG-UNG) 210 : ROGS=ROGS*OMEGA 211 : RUNGS=((SG-UNG)*(ROG*UNG))+(PSTAR - PG) 212 : RUNGS=RUNGS*OMEGA 213 : RETGS=((SG-UNG)*RETG)+(PSTAR*SM)-(PG*UNG) 214 : RETGS=RETGS*OMEGA 215 : C 216 : FLU1(1)=ROGS*SM 217 : FLU1(2)=RUNGS*SM+PSTAR 218 : FLU1(4)=(RETGS+PSTAR)*SM 219 : C 220 : ELSEIF(SD .GT. 0)THEN 221 : OMEGA=1.0D0/(SD-SM) 222 : RODS=ROD*(SD-UND) 223 : RODS=RODS*OMEGA 224 : RUNDS=((SD-UND)*(ROD*UND))+(PSTAR - PD) 225 : RUNDS=RUNDS*OMEGA 226 : RETDS=((SD-UND)*RETD)+(PSTAR*SM)-(PD*UND) 227 : RETDS=RETDS*OMEGA 228 : C 229 : FLU1(1)=RODS*SM 230 : FLU1(2)=RUNDS*SM+PSTAR 231 : FLU1(4)=(RETDS+PSTAR)*SM 232 : ELSE 233 : FLU1(1)=ROD*UND 234 : FLU1(2)=ROD*UND*UND+PD 235 : FLU1(4)=UND*(RETD+PD) 236 : ENDIF 237 : C 238 : FLU1(3)=(FLU1(1)*(UTG+UTD))+(ABS(FLU1(1))*(UTG-UTD)) 239 : FLU1(3)=HALF*FLU1(3) 240 : C 241 : DO IESP=1,NESP,1 242 : FLU1(4+IESP)=(FLU1(1)*(YG(IESP)+YD(IESP)))+ 243 : & (ABS(FLU1(1))*(YG(IESP)-YD(IESP))) 244 : FLU1(4+IESP)=HALF*FLU1(4+IESP) 245 : ENDDO 246 : C 247 : CELLT=1.0D0/(MAX(SG,SD)) 248 : C 249 : RETURN 250 : END 251 : 252 : 253 : 254 : 255 : 256 : 257 :
© Cast3M 2003 - All rights reserved.
Disclaimer