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 ! solids phase index used throughout routine... 121 ! may be inappropriate for multiple solids phases 122 M = 1 ! for solids phase 123 124 125 !!$omp parallel do private(ijk) schedule(dynamic,chunk_size) 126 DO IJK = ijkstart3, ijkend3 127 IF (FLUID_AT(IJK)) THEN 128 129 130 ! Gas viscosity (in Poise or Pa.s) 131 ! Calculating gas viscosity using Sutherland's formula with 132 ! Sutherland's constant (C) given by Vogel's equation C = 1.47*Tb. 133 ! For air C = 110 (Tb=74.82) 134 ! mu = 1.71*10-4 poise at T = 273K 135 136 IF (MU_G0 == UNDEFINED) MU_G(IJK) = to_SI*1.7D-4 * & 137 (T_G(IJK)/273.0D0)**1.5D0 * (383.D0/(T_G(IJK)+110.D0)) 138 139 MU_GT(IJK) = MU_G(IJK) 140 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 141 142 ! K_epsilon model 143 ! ---------------------------------------------------------------->>> 144 IF (K_Epsilon) THEN 145 146 C_MU = 9D-02 147 148 ! I'm not very confident about this correction in Peirano paper, but it's made 149 ! available here, uncomment to use it. sof@fluent.com --> 02/01/05 150 ! IF(SIMONIN .AND. F_GS(IJK,1) > SMALL_NUMBER) THEN 151 ! Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1) 152 ! X_21 = Ep_s(IJK,M)*RO_S(IJK,M)/(EP_g(IJK)*RO_g(IJK)) 153 ! new definition of C_mu (equation A.12, Peirano et al. (2002) Powder tech. 122,69-82) 154 ! IF( K_12(IJK)/(2.0D0*K_Turb_G(IJK)) < ONE) & 155 ! C_MU = C_MU/(ONE+ 0.314D0*X_21*Tau_12_st / Tau_1(IJK) * & 156 ! (ONE - K_12(IJK)/(2.0D0*K_Turb_G(IJK))) ) 157 ! ENDIF 158 ! Correction in Ahmadi paper (Cao and Ahmadi) 159 IF(AHMADI .AND. F_GS(IJK,1) > SMALL_NUMBER) THEN 160 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1) 161 C_MU = C_MU/(ONE+ Tau_12_st/Tau_1(IJK) * & 162 (EP_s(IJK,M)/(ONE-EP_star_array(IJK)))**3) 163 ENDIF 164 165 ! Definition of the turbulent viscosity 166 MU_GT(IJK) = MU_G(IJK) + RO_G(IJK)*C_MU*& 167 K_Turb_G(IJK)**2 / (E_Turb_G(IJK) + SMALL_NUMBER) 168 MU_GT(IJK) = MIN(MU_GMAX, MU_GT(IJK)) 169 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 170 ENDIF 171 ! ----------------------------------------------------------------<<< 172 173 ELSE 174 MU_G(IJK) = ZERO 175 MU_GT(IJK) = ZERO 176 LAMBDA_GT(IJK) = ZERO 177 ENDIF ! end if (fluid_at(ijk)) 178 179 ENDDO ! end do (ijk=ijkstart3,ijkend3) 180 181 182 ! MMS: Force constant gas viscosity at all cells including ghost cells. 183 IF (USE_MMS) THEN 184 DO IJK = ijkstart3, ijkend3 185 MU_G(IJK) = MU_G0 186 MU_GT(IJK) = MU_G(IJK) 187 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 188 ENDDO 189 END IF ! end if (USE_MMS) 190 191 192 193 ! L_scale0 model 194 ! ---------------------------------------------------------------->>> 195 !!$omp parallel do & 196 !!$omp$ schedule(dynamic,chunk_size) & 197 !!$omp$ private(IJK, I,J,K,IM,JM,KM, & 198 !!$omp& IMJK,IPJK,IJMK,IJPK,IJKM,IJKP,IMJPK,IMJMK,IMJKP, & 199 !!$omp& IMJKM,IPJKM,IPJMK,IJMKP,IJMKM,IJPKM, & 200 !!$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, & 201 !!$omp$ W_G_N,W_G_S,W_G_E,W_G_W, U_G_C,W_G_C, D_G,I2_DEVD_G ) 202 DO IJK = ijkstart3, ijkend3 203 IF ( FLUID_AT(IJK) .AND. L_SCALE(IJK)/=ZERO) THEN 204 I = I_OF(IJK) 205 J = J_OF(IJK) 206 K = K_OF(IJK) 207 IM = IM1(I) 208 JM = JM1(J) 209 KM = KM1(K) 210 211 IMJK = IM_OF(IJK) 212 IPJK = IP_OF(IJK) 213 IJMK = JM_OF(IJK) 214 IJPK = JP_OF(IJK) 215 IJKM = KM_OF(IJK) 216 IJKP = KP_OF(IJK) 217 IMJPK = IM_OF(IJPK) 218 IMJMK = IM_OF(IJMK) 219 IMJKP = IM_OF(IJKP) 220 IMJKM = IM_OF(IJKM) 221 IPJKM = IP_OF(IJKM) 222 IPJMK = IP_OF(IJMK) 223 IJMKP = JM_OF(IJKP) 224 IJMKM = JM_OF(IJKM) 225 IJPKM = JP_OF(IJKM) 226 227 ! Find fluid velocity values at faces of the cell 228 U_G_N = AVG_Y(AVG_X_E(U_G(IMJK),U_G(IJK),I),& 229 AVG_X_E(U_G(IMJPK),U_G(IJPK),I),J) !i, j+1/2, k 230 U_G_S = AVG_Y(AVG_X_E(U_G(IMJMK),U_G(IJMK),I),& 231 AVG_X_E(U_G(IMJK),U_G(IJK),I),JM) !i, j-1/2, k 232 U_G_T = AVG_Z(AVG_X_E(U_G(IMJK),U_G(IJK),I),& 233 AVG_X_E(U_G(IMJKP),U_G(IJKP),I),K) !i, j, k+1/2 234 U_G_B = AVG_Z(AVG_X_E(U_G(IMJKM),U_G(IJKM),I),& 235 AVG_X_E(U_G(IMJK),U_G(IJK),I),KM) !i, j, k-1/2 236 V_G_E = AVG_X(AVG_Y_N(V_G(IJMK),V_G(IJK)),& 237 AVG_Y_N(V_G(IPJMK),V_G(IPJK)),I) !i+1/2, j, k 238 V_G_W = AVG_X(AVG_Y_N(V_G(IMJMK),V_G(IMJK)),& 239 AVG_Y_N(V_G(IJMK),V_G(IJK)),IM) !i-1/2, j, k 240 V_G_T = AVG_Z(AVG_Y_N(V_G(IJMK),V_G(IJK)),& 241 AVG_Y_N(V_G(IJMKP),V_G(IJKP)),K) !i, j, k+1/2 242 V_G_B = AVG_Z(AVG_Y_N(V_G(IJMKM),V_G(IJKM)),& 243 AVG_Y_N(V_G(IJMK),V_G(IJK)),KM) !i, j, k-1/2 244 W_G_N = AVG_Y(AVG_Z_T(W_G(IJKM),W_G(IJK)),& 245 AVG_Z_T(W_G(IJPKM),W_G(IJPK)),J) !i, j+1/2, k 246 W_G_S = AVG_Y(AVG_Z_T(W_G(IJMKM),W_G(IJMK)),& 247 AVG_Z_T(W_G(IJKM),W_G(IJK)),JM) !i, j-1/2, k 248 W_G_E = AVG_X(AVG_Z_T(W_G(IJKM),W_G(IJK)),& 249 AVG_Z_T(W_G(IPJKM),W_G(IPJK)),I) !i+1/2, j, k 250 W_G_W = AVG_X(AVG_Z_T(W_G(IMJKM),W_G(IMJK)),& 251 AVG_Z_T(W_G(IJKM),W_G(IJK)),IM) !i-1/2, j, k 252 253 IF (CYLINDRICAL) THEN 254 U_G_C = AVG_X_E(U_G(IMJK),U_G(IJK),I) !i, j, k 255 W_G_C = AVG_Z_T(W_G(IJKM),W_G(IJK)) !i, j, k 256 ELSE 257 U_G_C = ZERO 258 W_G_C = ZERO 259 ENDIF 260 261 ! Find components of fluid phase strain rate tensor, D_g, at center of 262 ! the cell - (i,j,k) 263 D_G(1,1) = (U_G(IJK)-U_G(IMJK))*ODX(I) 264 D_G(1,2) = HALF*((U_G_N - U_G_S)*ODY(J)+(V_G_E-V_G_W)*& 265 ODX(I)) 266 D_G(1,3) = HALF*((W_G_E - W_G_W)*ODX(I)+(U_G_T-U_G_B)*& 267 (OX(I)*ODZ(K))-W_G_C*OX(I)) 268 D_G(2,1) = D_G(1,2) 269 D_G(2,2) = (V_G(IJK)-V_G(IJMK))*ODY(J) 270 D_G(2,3) = HALF*((V_G_T-V_G_B)*(OX(I)*ODZ(K))+& 271 (W_G_N-W_G_S)*ODY(J)) 272 D_G(3,1) = D_G(1,3) 273 D_G(3,2) = D_G(2,3) 274 D_G(3,3) = (W_G(IJK)-W_G(IJKM))*(OX(I)*ODZ(K)) + U_G_C*OX(I) 275 276 ! Calculate the second invariant of the deviator of D_g 277 I2_DEVD_G = ((D_G(1,1)-D_G(2,2))**2+(D_G(2,2)-D_G(3,3))**2+& 278 (D_G(3,3)-D_G(1,1))**2)/6.D0 + & 279 D_G(1,2)**2 + D_G(2,3)**2 + D_G(3,1)**2 280 281 MU_GT(IJK) = MIN(MU_GMAX, MU_G(IJK)+2.0*& 282 L_SCALE(IJK)*L_SCALE(IJK)*RO_G(IJK)*SQRT(I2_DEVD_G)) 283 LAMBDA_GT(IJK) = -F2O3*MU_GT(IJK) 284 285 ENDIF ! end if (fluid_at(ijk) and l_scale(ijk)/=0)) 286 ENDDO ! end loop (ijk=ijkstart3,ijkend3) 287 ! end calculations for L_scale0 model 288 ! ----------------------------------------------------------------<<< 289 290 RETURN 291 END SUBROUTINE CALC_MU_G 292 293