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 ! Author: W. Sams/M. Syamlal Date: 18-JUL-94 C 8 ! Reviewer: Date: dd-mmm-yy C 9 ! C 10 ! Revision Number: 1 C 11 ! Purpose: MFIX 2.0 mods (previous name CALC_MU_gt) C 12 ! Author: Date: dd-mmm-yy C 13 ! Reviewer: Date: dd-mmm-yy C 14 ! C 15 ! Revision Number: 2 C 16 ! Purpose: allow SI unit C 17 ! Author: S. Dartevelle Date: 01-Jul-02 C 18 ! Reviewer: Date: dd-mmm-yy C 19 ! C 20 ! Revision Number: 3 C 21 ! Purpose: compute turbulent eddy viscosity C 22 ! Author: S. Benyahia Date: May-13-04 C 23 ! Reviewer: Date: dd-mmm-yy C 24 ! C 25 ! Literature/Document References: C 26 ! Perry, R. H., and Chilton, C. H., Chemical Engineers' Handbook, C 27 ! 5th Edition, McGraw-Hill Inc., 1973, pp. 248, eqn. 3-133. C 28 ! Arnold, J. H., Vapor viscosities and the Sutherland equation, C 29 ! Journal of Chemical Physics, 1 (2), 1933, pp. 170-176. C 30 ! Sutherland, W., The Viscosity of Gases and Molecular Force, C 31 ! Phil. Mag. 5:507-531, 1893. C 32 ! C 33 ! Cao, J. and Ahmadi, G., 1995, Gas-particle two-phase turbulent C 34 ! flow in a vertical duct. Int. J. Multiphase Flow, vol. 21, C 35 ! No. 6, pp. 1203-1228. C 36 ! C 37 ! C 38 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 39 40 SUBROUTINE CALC_MU_G() 41 42 !----------------------------------------------- 43 ! Modules 44 !----------------------------------------------- 45 USE param 46 USE param1 47 USE parallel 48 USE physprop 49 USE geometry 50 USE fldvar 51 USE visc_g 52 USE visc_s 53 USE indices 54 USE constant 55 USE toleranc 56 USE compar 57 USE drag 58 USE run !S. Dartevelle 59 USE turb 60 USE sendrecv 61 USE mms 62 USE fun_avg 63 USE functions 64 65 IMPLICIT NONE 66 !----------------------------------------------- 67 ! Local parameters 68 !----------------------------------------------- 69 DOUBLE PRECISION, PARAMETER :: F2O3 = 2.D0/3.D0 70 !----------------------------------------------- 71 ! Local variables 72 !----------------------------------------------- 73 ! Cell indices 74 INTEGER :: I, J, K, IM, JM, KM 75 INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP 76 INTEGER :: IMJPK, IMJMK, IMJKP, IMJKM, IPJKM, IPJMK 77 INTEGER :: IJMKP, IJMKM, IJPKM 78 ! Solids phase index 79 INTEGER :: M 80 ! Strain rate tensor components for mth solids phase 81 DOUBLE PRECISION :: D_g(3,3) 82 ! Second invariant of the deviator of D_g 83 DOUBLE PRECISION :: I2_devD_g 84 ! Constant in turbulent viscosity formulation 85 DOUBLE PRECISION :: C_MU 86 ! particle relaxation time 87 DOUBLE PRECISION :: Tau_12_st 88 ! U_g at the north face of the THETA cell-(i, j+1/2, k) 89 DOUBLE PRECISION :: U_g_N 90 ! U_g at the south face of the THETA cell-(i, j-1/2, k) 91 DOUBLE PRECISION :: U_g_S 92 ! U_g at the top face of the THETA cell-(i, j, k+1/2) 93 DOUBLE PRECISION :: U_g_T 94 ! U_g at the bottom face of the THETA cell-(i, j, k-1/2) 95 DOUBLE PRECISION :: U_g_B 96 ! U_g at the center of the THETA cell-(i, j, k) 97 ! Calculated for Cylindrical coordinates only. 98 DOUBLE PRECISION :: U_g_C 99 ! V_g at the east face of the THETA cell-(i+1/2, j, k) 100 DOUBLE PRECISION :: V_g_E 101 ! V_g at the west face of the THETA cell-(i-1/2, j, k) 102 DOUBLE PRECISION :: V_g_W 103 ! V_g at the top face of the THETA cell-(i, j, k+1/2) 104 DOUBLE PRECISION :: V_g_T 105 ! V_g at the bottom face of the THETA cell-(i, j, k-1/2) 106 DOUBLE PRECISION :: V_g_B 107 ! W_g at the east face of the THETA cell-(i+1/2, j, k) 108 DOUBLE PRECISION :: W_g_E 109 ! W_g at the west face of the THETA cell-(1-1/2, j, k) 110 DOUBLE PRECISION :: W_g_W 111 ! W_g at the north face of the THETA cell-(i, j+1/2, k) 112 DOUBLE PRECISION :: W_g_N 113 ! W_g at the south face of the THETA cell-(i, j-1/2, k) 114 DOUBLE PRECISION :: W_g_S 115 ! W_g at the center of the THETA cell-(i, j, k). 116 ! Calculated for Cylindrical coordinates only. 117 DOUBLE PRECISION :: W_g_C 118 !----------------------------------------------- 119 120 ! JFD: Calling SET_EP_FACTORS here to make sure EPG_IFAC is defined 121 ! when calc_mu_g is called for the firt time. 122 ! There may be a better way to do this... 123 CALL SET_EP_FACTORS 124 125 ! solids phase index used throughout routine... 126 ! may be inappropriate for multiple solids phases 127 M = 1 ! for solids phase 128 129 130 !!$omp parallel do private(ijk) schedule(dynamic,chunk_size) 131 DO IJK = ijkstart3, ijkend3 132 IF (FLUID_AT(IJK)) THEN 133 134 135 ! Gas viscosity (in Poise or Pa.s) 136 ! Calculating gas viscosity using Sutherland's formula with 137 ! Sutherland's constant (C) given by Vogel's equation C = 1.47*Tb. 138 ! For air C = 110 (Tb=74.82) 139 ! mu = 1.71*10-4 poise at T = 273K 140 141 IF (MU_G0 == UNDEFINED) MU_G(IJK) = to_SI*1.7D-4 * & 142 (T_G(IJK)/273.0D0)**1.5D0 * (383.D0/(T_G(IJK)+110.D0)) 143 144 MU_GT(IJK) = MU_G(IJK) 145 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 146 ! if ishii then multiply by void fraction otherwise multiply by 1 147 EPMU_GT(IJK)= EPG_IFAC(IJK)*MU_GT(IJK) 148 EPLAMBDA_GT(IJK) = EPG_IFAC(IJK)*LAMBDA_GT(IJK) 149 150 ! K_epsilon model 151 ! ---------------------------------------------------------------->>> 152 IF (K_Epsilon) THEN 153 154 C_MU = 9D-02 155 156 ! I'm not very confident about this correction in Peirano paper, but it's made 157 ! available here, uncomment to use it. sof@fluent.com --> 02/01/05 158 ! IF(SIMONIN .AND. F_GS(IJK,1) > SMALL_NUMBER) THEN 159 ! Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1) 160 ! X_21 = Ep_s(IJK,M)*RO_S(IJK,M)/(EP_g(IJK)*RO_g(IJK)) 161 ! new definition of C_mu (equation A.12, Peirano et al. (2002) Powder tech. 122,69-82) 162 ! IF( K_12(IJK)/(2.0D0*K_Turb_G(IJK)) < ONE) & 163 ! C_MU = C_MU/(ONE+ 0.314D0*X_21*Tau_12_st / Tau_1(IJK) * & 164 ! (ONE - K_12(IJK)/(2.0D0*K_Turb_G(IJK))) ) 165 ! ENDIF 166 ! Correction in Ahmadi paper (Cao and Ahmadi) 167 IF(AHMADI .AND. F_GS(IJK,1) > SMALL_NUMBER) THEN 168 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1) 169 C_MU = C_MU/(ONE+ Tau_12_st/Tau_1(IJK) * & 170 (EP_s(IJK,M)/(ONE-EP_star_array(IJK)))**3) 171 ENDIF 172 173 ! Definition of the turbulent viscosity 174 MU_GT(IJK) = MU_G(IJK) + RO_G(IJK)*C_MU*& 175 K_Turb_G(IJK)**2 / (E_Turb_G(IJK) + SMALL_NUMBER) 176 MU_GT(IJK) = MIN(MU_GMAX, MU_GT(IJK)) 177 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 178 ! if ishii then multiply by void fraction otherwise multiply by 1 179 EPMU_GT(IJK) = EPG_IFAC(IJK)*MU_GT(IJK) 180 EPLAMBDA_GT(IJK) = EPG_IFAC(IJK)*LAMBDA_GT(IJK) 181 ENDIF 182 ! ----------------------------------------------------------------<<< 183 184 ELSE 185 MU_G(IJK) = ZERO 186 MU_GT(IJK) = ZERO 187 LAMBDA_GT(IJK) = ZERO 188 EPMU_GT(IJK) = ZERO 189 EPLAMBDA_GT(IJK) = ZERO 190 ENDIF ! end if (fluid_at(ijk)) 191 192 ENDDO ! end do (ijk=ijkstart3,ijkend3) 193 194 195 ! MMS: Force constant gas viscosity at all cells including ghost cells. 196 IF (USE_MMS) THEN 197 DO IJK = ijkstart3, ijkend3 198 MU_G(IJK) = MU_G0 199 MU_GT(IJK) = MU_G(IJK) 200 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 201 ! if ishii then multiply by void fraction otherwise multiply by 1 202 EPMU_GT(IJK) = EPG_IFAC(IJK)*MU_GT(IJK) 203 EPLAMBDA_GT(IJK) = EPG_IFAC(IJK)*LAMBDA_GT(IJK) 204 ENDDO 205 END IF ! end if (USE_MMS) 206 207 208 209 ! L_scale0 model 210 ! ---------------------------------------------------------------->>> 211 !!$omp parallel do & 212 !!$omp$ schedule(dynamic,chunk_size) & 213 !!$omp$ private(IJK, I,J,K,IM,JM,KM, & 214 !!$omp& IMJK,IPJK,IJMK,IJPK,IJKM,IJKP,IMJPK,IMJMK,IMJKP, & 215 !!$omp& IMJKM,IPJKM,IPJMK,IJMKP,IJMKM,IJPKM, & 216 !!$omp& U_G_N,U_G_S,U_G_T,U_G_B,V_G_E,V_G_W,V_G_T,V_G_B, & 217 !!$omp$ W_G_N,W_G_S,W_G_E,W_G_W, U_G_C,W_G_C, D_G,I2_DEVD_G ) 218 DO IJK = ijkstart3, ijkend3 219 IF ( FLUID_AT(IJK) .AND. L_SCALE(IJK)/=ZERO) THEN 220 I = I_OF(IJK) 221 J = J_OF(IJK) 222 K = K_OF(IJK) 223 IM = IM1(I) 224 JM = JM1(J) 225 KM = KM1(K) 226 227 IMJK = IM_OF(IJK) 228 IPJK = IP_OF(IJK) 229 IJMK = JM_OF(IJK) 230 IJPK = JP_OF(IJK) 231 IJKM = KM_OF(IJK) 232 IJKP = KP_OF(IJK) 233 IMJPK = IM_OF(IJPK) 234 IMJMK = IM_OF(IJMK) 235 IMJKP = IM_OF(IJKP) 236 IMJKM = IM_OF(IJKM) 237 IPJKM = IP_OF(IJKM) 238 IPJMK = IP_OF(IJMK) 239 IJMKP = JM_OF(IJKP) 240 IJMKM = JM_OF(IJKM) 241 IJPKM = JP_OF(IJKM) 242 243 ! Find fluid velocity values at faces of the cell 244 U_G_N = AVG_Y(AVG_X_E(U_G(IMJK),U_G(IJK),I),& 245 AVG_X_E(U_G(IMJPK),U_G(IJPK),I),J) !i, j+1/2, k 246 U_G_S = AVG_Y(AVG_X_E(U_G(IMJMK),U_G(IJMK),I),& 247 AVG_X_E(U_G(IMJK),U_G(IJK),I),JM) !i, j-1/2, k 248 U_G_T = AVG_Z(AVG_X_E(U_G(IMJK),U_G(IJK),I),& 249 AVG_X_E(U_G(IMJKP),U_G(IJKP),I),K) !i, j, k+1/2 250 U_G_B = AVG_Z(AVG_X_E(U_G(IMJKM),U_G(IJKM),I),& 251 AVG_X_E(U_G(IMJK),U_G(IJK),I),KM) !i, j, k-1/2 252 V_G_E = AVG_X(AVG_Y_N(V_G(IJMK),V_G(IJK)),& 253 AVG_Y_N(V_G(IPJMK),V_G(IPJK)),I) !i+1/2, j, k 254 V_G_W = AVG_X(AVG_Y_N(V_G(IMJMK),V_G(IMJK)),& 255 AVG_Y_N(V_G(IJMK),V_G(IJK)),IM) !i-1/2, j, k 256 V_G_T = AVG_Z(AVG_Y_N(V_G(IJMK),V_G(IJK)),& 257 AVG_Y_N(V_G(IJMKP),V_G(IJKP)),K) !i, j, k+1/2 258 V_G_B = AVG_Z(AVG_Y_N(V_G(IJMKM),V_G(IJKM)),& 259 AVG_Y_N(V_G(IJMK),V_G(IJK)),KM) !i, j, k-1/2 260 W_G_N = AVG_Y(AVG_Z_T(W_G(IJKM),W_G(IJK)),& 261 AVG_Z_T(W_G(IJPKM),W_G(IJPK)),J) !i, j+1/2, k 262 W_G_S = AVG_Y(AVG_Z_T(W_G(IJMKM),W_G(IJMK)),& 263 AVG_Z_T(W_G(IJKM),W_G(IJK)),JM) !i, j-1/2, k 264 W_G_E = AVG_X(AVG_Z_T(W_G(IJKM),W_G(IJK)),& 265 AVG_Z_T(W_G(IPJKM),W_G(IPJK)),I) !i+1/2, j, k 266 W_G_W = AVG_X(AVG_Z_T(W_G(IMJKM),W_G(IMJK)),& 267 AVG_Z_T(W_G(IJKM),W_G(IJK)),IM) !i-1/2, j, k 268 269 IF (CYLINDRICAL) THEN 270 U_G_C = AVG_X_E(U_G(IMJK),U_G(IJK),I) !i, j, k 271 W_G_C = AVG_Z_T(W_G(IJKM),W_G(IJK)) !i, j, k 272 ELSE 273 U_G_C = ZERO 274 W_G_C = ZERO 275 ENDIF 276 277 ! Find components of fluid phase strain rate tensor, D_g, at center of 278 ! the cell - (i,j,k) 279 D_G(1,1) = (U_G(IJK)-U_G(IMJK))*ODX(I) 280 D_G(1,2) = HALF*((U_G_N - U_G_S)*ODY(J)+(V_G_E-V_G_W)*& 281 ODX(I)) 282 D_G(1,3) = HALF*((W_G_E - W_G_W)*ODX(I)+(U_G_T-U_G_B)*& 283 (OX(I)*ODZ(K))-W_G_C*OX(I)) 284 D_G(2,1) = D_G(1,2) 285 D_G(2,2) = (V_G(IJK)-V_G(IJMK))*ODY(J) 286 D_G(2,3) = HALF*((V_G_T-V_G_B)*(OX(I)*ODZ(K))+& 287 (W_G_N-W_G_S)*ODY(J)) 288 D_G(3,1) = D_G(1,3) 289 D_G(3,2) = D_G(2,3) 290 D_G(3,3) = (W_G(IJK)-W_G(IJKM))*(OX(I)*ODZ(K)) + U_G_C*OX(I) 291 292 ! Calculate the second invariant of the deviator of D_g 293 I2_DEVD_G = ((D_G(1,1)-D_G(2,2))**2+(D_G(2,2)-D_G(3,3))**2+& 294 (D_G(3,3)-D_G(1,1))**2)/6.D0 + & 295 D_G(1,2)**2 + D_G(2,3)**2 + D_G(3,1)**2 296 297 MU_GT(IJK) = MIN(MU_GMAX, MU_G(IJK)+2.0*& 298 L_SCALE(IJK)*L_SCALE(IJK)*RO_G(IJK)*SQRT(I2_DEVD_G)) 299 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 300 ! if ishii then multiply by void fraction otherwise multiply by 1 301 EPMU_GT(IJK) = EPG_IFAC(IJK)*MU_GT(IJK) 302 EPLAMBDA_GT(IJK) = EPG_IFAC(IJK)*LAMBDA_GT(IJK) 303 ENDIF ! end if (fluid_at(ijk) and l_scale(ijk)/=0)) 304 ENDDO ! end loop (ijk=ijkstart3,ijkend3) 305 ! end calculations for L_scale0 model 306 ! ----------------------------------------------------------------<<< 307 308 RETURN 309 END SUBROUTINE CALC_MU_G 310 311