1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CALC_MU_g C 4 ! Purpose: Calculate the effective viscosity for a turbulent flow, C 5 ! which is the sum of molecular and eddy viscosities C 6 ! C 7 ! C 8 ! Comments: This routine is called even if mu_g0 is defined C 9 ! C 10 ! C 11 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 12 SUBROUTINE CALC_MU_G() 13 14 ! Modules 15 !---------------------------------------------------------------------// 16 use constant, only: l_scale0 17 use param1, only: undefined, zero 18 use physprop, only: mu_g0 19 use run, only: k_epsilon 20 ! invoke user defined quantity 21 USE usr_prop, only: usr_mug, calc_usr_prop 22 USE usr_prop, only: gas_viscosity 23 IMPLICIT NONE 24 25 ! Local variables 26 !---------------------------------------------------------------------// 27 ! Cell indices 28 INTEGER :: IJK 29 !---------------------------------------------------------------------// 30 31 IF (USR_MUg) THEN 32 CALL CALC_USR_PROP(Gas_Viscosity,lm=0) 33 ELSEIF (MU_G0 == UNDEFINED) THEN 34 ! this is a necessary check as calc_mu_g is called at least once for 35 ! initialization and for other possible reasons (turbulence, ishii, etc) 36 CALL CALC_DEFAULT_MUg 37 ENDIF 38 39 ! adjust viscosity for tubulence 40 IF (K_Epsilon) THEN 41 CALL CALC_K_EPSILON_MU 42 ELSEIF (L_SCALE0 /= ZERO) THEN 43 CALL CALC_LSCALE_MU 44 ENDIF 45 46 CALL SET_EPMUG_VALUES 47 48 RETURN 49 END SUBROUTINE CALC_MU_G 50 51 52 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 53 ! C 54 ! Subroutine: SET_EPMUG_VALUES C 55 ! Purpose: This routine sets the internal variables epmu_g and C 56 ! eplambda_g that are used in the stress calculations. If the C 57 ! keyword Ishii is invoked then these quantities represent the C 58 ! viscosity and second viscosity multiplied by the volume fraction C 59 ! otherwise they are simply viscosity/second viscosity (i.e. are C 60 ! multipled by one). C 61 ! C 62 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 63 SUBROUTINE SET_EPMUG_VALUES 64 65 ! Modules 66 !---------------------------------------------------------------------// 67 USE compar, only: ijkstart3, ijkend3 68 USE fldvar, only: epg_ifac 69 USE functions, only: fluid_at 70 use mms, only: use_mms 71 USE param1, only: zero 72 USE visc_g, only: mu_gt, epmu_gt, lambda_gt, eplambda_gt 73 ! Local variables 74 !---------------------------------------------------------------------// 75 ! cell index 76 INTEGER :: IJK 77 !---------------------------------------------------------------------// 78 79 ! EPMU_GT(:) = EPG_IFAC(:)*MU_gt(:) 80 ! EPLAMBDA_GT(:) = EPG_IFAC(:)*LAMBDA_gt(:) 81 82 ! Assign and update ep_g*mu_gt and ep_g*lambda_gt. This is needed even 83 ! for a constant viscosity case 84 DO IJK = ijkstart3, ijkend3 85 ! MMS: Force constant gas viscosity at all cells including ghost cells. 86 IF (FLUID_AT(IJK) .OR. USE_MMS) THEN 87 88 ! if ishii then multiply by void fraction otherwise multiply by 1 89 EPMU_GT(IJK)= EPG_IFAC(IJK)*MU_GT(IJK) 90 EPLAMBDA_GT(IJK) = EPG_IFAC(IJK)*LAMBDA_GT(IJK) 91 ELSE 92 EPMU_GT(IJK) = ZERO 93 EPLAMBDA_GT(IJK) = ZERO 94 ENDIF ! end if (fluid_at(ijk) .or. use_mms) 95 ENDDO ! end do (ijk=ijkstart3,ijkend3) 96 97 RETURN 98 END SUBROUTINE SET_EPMUG_VALUES 99 100 101 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 102 ! C 103 ! Purpose: Compute a default value of gas viscosity where gas is C 104 ! assumed to be air C 105 ! Author: W. Sams/M. Syamlal Date: 18-JUL-94 C 106 ! C 107 ! Literature/Document References: C 108 ! Perry, R. H., and Chilton, C. H., Chemical Engineers' Handbook, C 109 ! 5th Edition, McGraw-Hill Inc., 1973, pp. 248, eqn. 3-133. C 110 ! Arnold, J. H., Vapor viscosities and the Sutherland equation, C 111 ! Journal of Chemical Physics, 1 (2), 1933, pp. 170-176. C 112 ! Sutherland, W., The Viscosity of Gases and Molecular Force, C 113 ! Phil. Mag. 5:507-531, 1893. C 114 ! C 115 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 116 SUBROUTINE CALC_DEFAULT_MUG 117 118 ! Modules 119 !---------------------------------------------------------------------// 120 use compar, only: ijkstart3, ijkend3 121 use constant, only: to_si 122 use fldvar, only: T_g 123 use functions, only: fluid_at 124 use param1, only: zero 125 use physprop, only: mu_g 126 use visc_g, only: mu_gt 127 use visc_g, only: lambda_gt 128 IMPLICIT NONE 129 130 ! Local parameters 131 !---------------------------------------------------------------------// 132 DOUBLE PRECISION, PARAMETER :: F2O3 = 2.D0/3.D0 133 134 ! Local variables 135 !---------------------------------------------------------------------// 136 ! Cell indices 137 INTEGER :: IJK 138 !---------------------------------------------------------------------// 139 140 !!$omp parallel do private(ijk) schedule(dynamic,chunk_size) 141 DO IJK = ijkstart3, ijkend3 142 IF (FLUID_AT(IJK)) THEN 143 144 ! Gas viscosity (in Poise or Pa.s) 145 ! Calculating gas viscosity using Sutherland's formula with 146 ! Sutherland's constant (C) given by Vogel's equation C = 1.47*Tb. 147 ! For air C = 110 (Tb=74.82) 148 ! mu = 1.71*10-4 poise at T = 273K 149 MU_G(IJK) = to_SI*1.7D-4 * & 150 (T_G(IJK)/273.0D0)**1.5D0 * (383.D0/(T_G(IJK)+110.D0)) 151 152 ! assign values to quantities used throughout code 153 MU_GT(IJK) = MU_G(IJK) 154 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 155 ELSE 156 MU_G(IJK) = ZERO 157 MU_GT(IJK) = ZERO 158 LAMBDA_GT(IJK) = ZERO 159 ENDIF 160 ENDDO 161 162 RETURN 163 END SUBROUTINE CALC_DEFAULT_MUG 164 165 166 167 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 168 ! C 169 ! Purpose: compute turbulent eddy viscosity C 170 ! Author: S. Benyahia Date: May-13-04 C 171 ! C 172 ! C 173 ! Literature/Document References: C 174 ! Cao, J. and Ahmadi, G., 1995, Gas-particle two-phase turbulent C 175 ! flow in a vertical duct. Int. J. Multiphase Flow, vol. 21, C 176 ! No. 6, pp. 1203-1228. C 177 ! C 178 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 179 SUBROUTINE CALC_K_EPSILON_MU 180 181 ! Modules 182 !---------------------------------------------------------------------// 183 use compar, only: ijkstart3, ijkend3 184 use constant, only: mu_gmax 185 use drag, only: f_gs 186 use fldvar, only: k_turb_g, e_turb_g, ro_g 187 use fldvar, only: ep_s, ro_s 188 use functions, only: fluid_at 189 use param1, only: zero, one, small_number 190 use physprop, only: mu_g 191 use run, only: kt_type_enum, ahmadi_1995 192 use turb, only: tau_1 193 use visc_g, only: mu_gt, lambda_gt 194 use visc_s, only: ep_star_array 195 IMPLICIT NONE 196 197 ! Local parameters 198 !---------------------------------------------------------------------// 199 DOUBLE PRECISION, PARAMETER :: F2O3 = 2.D0/3.D0 200 201 ! Local variables 202 !---------------------------------------------------------------------// 203 ! Solids phase index 204 INTEGER :: M 205 ! Constant in turbulent viscosity formulation 206 DOUBLE PRECISION :: C_MU 207 ! particle relaxation time 208 DOUBLE PRECISION :: Tau_12_st 209 ! cell index 210 INTEGER :: IJK 211 !---------------------------------------------------------------------// 212 213 C_MU = 9D-02 214 215 DO IJK = ijkstart3, ijkend3 216 IF (FLUID_AT(IJK)) THEN 217 218 ! Correction in Ahmadi paper (Cao and Ahmadi) 219 IF(KT_TYPE_ENUM == AHMADI_1995 .AND.& 220 F_GS(IJK,1) > SMALL_NUMBER) THEN 221 ! solids phase index used throughout routine... 222 M = 1 ! for solids phase 223 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1) 224 C_MU = C_MU/(ONE+ Tau_12_st/Tau_1(IJK) * & 225 (EP_s(IJK,M)/(ONE-EP_star_array(IJK)))**3) 226 ENDIF 227 228 ! I'm not very confident about this correction in Peirano paper, 229 ! but it's made available here, uncomment to use it. 230 ! sof@fluent.com --> 02/01/05 231 ! IF(KT_TYPE_ENUM==SIMONIN_1996 .AND.& 232 ! F_GS(IJK,1) > SMALL_NUMBER) THEN 233 ! solids phase index used throughout routine... 234 ! M = 1 ! for solids phase 235 ! Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1) 236 ! X_21 = Ep_s(IJK,M)*RO_S(IJK,M)/(EP_g(IJK)*RO_g(IJK)) 237 ! new definition of C_mu (equation A.12, Peirano et al. (2002), 238 ! Powder tech. 122,69-82) 239 ! IF( K_12(IJK)/(2.0D0*K_Turb_G(IJK)) < ONE) & 240 ! C_MU = C_MU/(ONE+ 0.314D0*X_21*Tau_12_st / Tau_1(IJK) * & 241 ! (ONE - K_12(IJK)/(2.0D0*K_Turb_G(IJK))) ) 242 ! ENDIF 243 244 ! Definition of the turbulent viscosity 245 MU_GT(IJK) = MU_G(IJK) + RO_G(IJK)*C_MU*& 246 K_Turb_G(IJK)**2 / (E_Turb_G(IJK) + SMALL_NUMBER) 247 248 MU_GT(IJK) = MIN(MU_GMAX, MU_GT(IJK)) 249 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 250 ELSE 251 MU_GT(IJK) = ZERO 252 LAMBDA_GT(IJK) = ZERO 253 ENDIF 254 ENDDO 255 256 RETURN 257 END SUBROUTINE CALC_K_EPSILON_MU 258 259 260 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 261 ! C 262 ! Purpose: compute l_scale0 model C 263 ! C 264 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 265 SUBROUTINE CALC_LSCALE_MU 266 267 ! Modules 268 !---------------------------------------------------------------------// 269 use compar, only: ijkstart3, ijkend3 270 use constant, only: mu_gmax 271 use fldvar, only: ro_g 272 use functions, only: fluid_at 273 use param1, only: zero 274 use physprop, only: mu_g 275 use visc_g, only: mu_gt, lambda_gt 276 use visc_g, only: l_scale 277 IMPLICIT NONE 278 279 ! Local parameters 280 !---------------------------------------------------------------------// 281 DOUBLE PRECISION, PARAMETER :: F2O3 = 2.D0/3.D0 282 283 ! Local variables 284 !---------------------------------------------------------------------// 285 DOUBLE PRECISION :: D_g(3,3) 286 ! Gas velocity gradient 287 DOUBLE PRECISION :: DelV_g(3,3) 288 ! Second invariant of the deviator of D_g 289 DOUBLE PRECISION :: I2_devD_g 290 ! cell index 291 INTEGER :: IJK 292 !---------------------------------------------------------------------// 293 294 DO IJK = ijkstart3, ijkend3 295 IF (FLUID_AT(IJK)) THEN 296 ! Calculate the rate of strain tensor D_g 297 CALL CALC_DERIV_VEL_GAS(ijk, DelV_G, D_G) 298 299 ! Calculate the second invariant of the deviator of D_g 300 I2_DEVD_G = ((D_G(1,1)-D_G(2,2))**2+(D_G(2,2)-D_G(3,3))**2+& 301 (D_G(3,3)-D_G(1,1))**2)/6.D0 + & 302 D_G(1,2)**2 + D_G(2,3)**2 + D_G(3,1)**2 303 304 MU_GT(IJK) = MU_G(IJK)+2.0*L_SCALE(IJK)*L_SCALE(IJK)*& 305 RO_G(IJK)*SQRT(I2_DEVD_G) 306 307 MU_GT(IJK) = MIN(MU_GMAX, MU_GT(IJK)) 308 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 309 ELSE 310 MU_GT(IJK) = ZERO 311 LAMBDA_GT(IJK) = ZERO 312 ENDIF 313 ENDDO 314 315 RETURN 316 END SUBROUTINE CALC_LSCALE_MU 317