1 : C FDV32F SOURCE BECC 05/12/02 21:15:16 5258 2 : SUBROUTINE fDV32F(INDMET, 3 : & alphaL, uvnL, uvtL, uvvL, ulnL, ultL, ulvL, 4 : & pL, TvL, TlL, rvL, rlL, 5 : & alphaR, uvnR, uvtR, uvvR, ulnR, ultR, ulvR, 6 : & pR, TvR, TlR, rvR, rlR, 7 : & F) 8 : C*********************************************************************** 9 : C NOM : fDV32F 10 : C DESCRIPTION : Calculation of the Two phase AUSMDV flux 11 : C (3D two fluid model) 12 : C 13 : C LANGAGE : ESOPE 14 : C AUTEUR : José R. Garcia Cascales (CEA/DEN/DM2S/SFME/LTMF) 15 : C mél : fd1@semt2.smts.cea.fr 16 : C*********************************************************************** 17 : C APPELES : 18 : C APPELES (E/S) : 19 : C inputs : INDMET, variable to say if we want 20 : C preconditioned version or not 21 : C primitive variable at left side 22 : C AL, UVNL, UVTL, UVVL, ULNL, ULTL, ULVL, 23 : C PL, TVL, TLL, RVG, RLG, 24 : C primitive variables at right side 25 : C AR, UVNR, UVTR, UVVR, ULNR, ULTR, ULVR, 26 : C PR, TVR, TLR, RVR, RLR 27 : C 28 : C output : F, 3D AUSMDV or "preconditioned" AUSMDV flux 29 : C 30 : C APPELES (BLAS) : 31 : C APPELES (CALCUL) : 32 : C APPELE PAR : 33 : C*********************************************************************** 34 : C SYNTAXE GIBIANE : 35 : C ENTREES : 36 : C ENTREES/SORTIES : 37 : C SORTIES : 38 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 39 : C*********************************************************************** 40 : C VERSION : v1, 13/03/2002, version initiale 41 : C HISTORIQUE : 42 : C HISTORIQUE : 43 : C*********************************************************************** 44 : * 45 : real*8 F(10) 46 : real*8 alphaL, alphaR, uvnL, uvnR, uvtL, uvtR,uvvL, uvvR, 47 : & ulnL, ulnR, ultL, ultR, ulvL,ulvR, 48 : & pL, pR, evL, evR, elL, 49 : & elR, hvL, hvR, hlL, hlR, TvL, TvR, TlL, TlR 50 : integer INDMET 51 : C variables for AUSMDV scheme 52 : real*8 rvL, rvR, rlL, rlR, avL, avR, alL, alR 53 : real*8 alm, avm, MlL, MlR, MvL, MvR, avpm, alpm, 54 : & dotmv, dotml, Mvplus, Mvmin, Mlplus, Mlmin, 55 : & Pvplus, Pvmin, Plplus, Plmin, 56 : & fvL, fvR, flL, flR, wvplus, wvmin, wlplus, wlmin 57 : 58 : real*8 gamv, cpv, gaml, cpl, pil 59 : PARAMETER(gamv = 1.4D0, cpv = 1008.D0, 60 : & gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8) 61 : 62 : evL = cpv*TvL/gamv 63 : evR = cpv*TvR/gamv 64 : elL = cpl*TlL/gaml + pil/rlL 65 : elR = cpl*TlR/gaml + pil/rlR 66 : 67 : hvL = evL + pL/rvL 68 : hvR = evR + pR/rvR 69 : hlL = elL + pL/rlL 70 : hlR = elR + pR/rlR 71 : 72 : C Speed of sound at left and rigth sides 73 : avL = dsqrt((gamv - 1.D0)*Cpv*TvL) 74 : avR = dsqrt((gamv - 1.D0)*Cpv*TvR) 75 : avm = dsqrt(avL*avR) 76 : alL = dsqrt((gaml - 1.D0)*Cpl*TlL) 77 : alR = dsqrt((gaml - 1.D0)*Cpl*TlR) 78 : alm = dsqrt(alL*alR) 79 : 80 : C Match number at both sides 81 : MlL = ulnL/alm 82 : MlR = ulnR/alm 83 : MvL = uvnL/avm 84 : MvR = uvnR/avm 85 : 86 : C AUSMDV scheme 87 : C fvL = pL/(alphaL*rvL) 88 : C fvR = pR/(alphaR*rvR) 89 : C flL = pL/((1.D0 - alphaL)*rlL) 90 : C flR = pR/((1.D0 - alphaR)*rlR) 91 : 92 : fvL = 1.D0/(alphaL*rvL) 93 : fvR = 1.D0/(alphaR*rvR) 94 : flL = 1.D0/((1.D0 - alphaL)*rlL) 95 : flR = 1.D0/((1.D0 - alphaR)*rlR) 96 : 97 : C fvL = pL/(rvL) 98 : C fvR = pR/(rvR) 99 : C flL = pL/(rlL) 100 : C flR = pR/(rlR) 101 : 102 : C fvL = 1.D0/(rvL) 103 : C fvR = 1.D0/(rvR) 104 : C flL = 1.D0/(rlL) 105 : C flR = 1.D0/(rlR) 106 : 107 : wvplus = 2.D0*fvL/(fvL + fvR) 108 : wvmin = 2.D0*fvR/(fvL + fvR) 109 : wlplus = 2.D0*flL/(flL + flR) 110 : wlmin = 2.D0*flR/(flL + flR) 111 : 111 : call mpADV(wvplus, wvmin, MvL, MvR, Mvplus, Mvmin, Pvplus, Pvmin) 111 : call mpADV(wlplus, wlmin, MlL, MlR, Mlplus, Mlmin, Plplus, Plmin) 112 : 113 : C AUSM schemes mass fluxes 114 : dotmv = avm*(alphaL*rvL*Mvplus + alphaR*rvR*Mvmin) 115 : dotml = alm*((1.D0 - alphaL)*rlL*Mlplus + (1.D0 - alphaR)*rlR 116 : $ *Mlmin) 117 : if(INDMET .EQ. 4)then 118 : dotml = dotml + 0.01d0*(pL - pR) 119 : end if 120 : C AUSM schemes pressure fluxes 121 : avpm = Pvplus*alphaL*pL + Pvmin*alphaR*pR 122 : alpm = Plplus*(1.D0 - alphaL)*pL + Plmin*(1.D0 - alphaR) 123 : $ *pR 124 : 125 : C Definition of convective part of the flux 126 : 127 : if (dotmv .GT. 0.D0) then 128 : 129 : F(1) = dotmv*1.D0 130 : F(3) = dotmv*uvnL 131 : F(5) = dotmv*uvtL 132 : F(7) = dotmv*uvvL 133 : F(9) = dotmv*(hvL + (uvnL**2 + uvtL**2 + uvvL**2)/2.D0) 134 : else 135 : F(1) = dotmv*1.D0 136 : F(3) = dotmv*uvnR 137 : F(5) = dotmv*uvtR 138 : F(7) = dotmv*uvvR 139 : F(9) = dotmv*(hvR + (uvnR**2 + uvtR**2 + uvvR**2)/2.D0) 140 : end if 141 : 142 : if (dotml .GT. 0.D0) then 143 : F(2) = dotml*1.D0 144 : F(4) = dotml*ulnL 145 : F(6) = dotml*ultL 146 : F(8) = dotml*ulvL 147 : F(10) = dotml*(hlL + (ulnL**2 + ultL**2 + ulvL**2)/2.D0) 148 : else 149 : F(2) = dotml*1.D0 150 : F(4) = dotml*ulnR 151 : F(6) = dotml*ultR 152 : F(8) = dotml*ulvR 153 : F(10) = dotml*(hlR + (ulnR**2 + ultR**2 + ulvR**2)/2.D0) 154 : end if 155 : 156 : C Addition of the pressure term 157 : 158 : F(3) = F(3) + avpm 159 : F(4) = F(4) + alpm 160 : 161 : END 162 : 163 :
© Cast3M 2003 - All rights reserved.
Disclaimer