calc_mu_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_MU_g                                               C
!  Purpose: Calculate the effective viscosity for a turbulent flow,    C
!           which is the sum of molecular and eddy viscosities         C
!                                                                      C
!                                                                      C
!  Comments: This routine is called even if mu_g0 is defined           C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_mu_g()

! Modules
!---------------------------------------------------------------------//
      use constant, only: l_scale0
      use param1, only: undefined, zero
      use physprop, only: mu_g0
      use run, only: k_epsilon
! invoke user defined quantity
      USE usr_prop, only: usr_mug, calc_usr_prop
      USE usr_prop, only: gas_viscosity
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! Cell indices
      INTEGER :: IJK
!---------------------------------------------------------------------//

      IF (usr_mug) THEN
         CALL calc_usr_prop(gas_viscosity,lm=0)
      ELSEIF (mu_g0 == undefined) THEN
! this is a necessary check as calc_mu_g is called at least once for
! initialization and for other possible reasons (turbulence, ishii, etc)
         CALL calc_default_mug
      ENDIF

! adjust viscosity for tubulence
      IF (k_epsilon) THEN
         CALL calc_k_epsilon_mu
      ELSEIF (l_scale0 /= zero) THEN
         CALL calc_lscale_mu
      ENDIF

      CALL set_epmug_values

      RETURN
      END SUBROUTINE calc_mu_g


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_EPMUG_VALUES                                        C
!  Purpose: This routine sets the internal variables epmu_g and        C
!  eplambda_g that are used in the stress calculations. If the         C
!  keyword Ishii is invoked then these quantities represent the        C
!  viscosity and second viscosity multiplied by the volume fraction    C
!  otherwise they are simply viscosity/second viscosity (i.e. are      C
!  multipled by one).                                                  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_epmug_values

! Modules
!---------------------------------------------------------------------//
      USE compar, only: ijkstart3, ijkend3
      USE fldvar, only: epg_ifac
      USE functions, only: fluid_at
      use mms, only: use_mms
      USE param1, only: zero
      USE visc_g, only: mu_gt, epmu_gt, lambda_gt, eplambda_gt
! Local variables
!---------------------------------------------------------------------//
! cell index
      INTEGER :: IJK
!---------------------------------------------------------------------//

!      EPMU_GT(:) = EPG_IFAC(:)*MU_gt(:)
!      EPLAMBDA_GT(:) = EPG_IFAC(:)*LAMBDA_gt(:)

! Assign and update ep_g*mu_gt and ep_g*lambda_gt. This is needed even
! for a constant viscosity case
      DO ijk = ijkstart3, ijkend3
! MMS: Force constant gas viscosity at all cells including ghost cells.
         IF (fluid_at(ijk) .OR. use_mms) THEN

! if ishii then multiply by void fraction otherwise multiply by 1
            epmu_gt(ijk)= epg_ifac(ijk)*mu_gt(ijk)
            eplambda_gt(ijk) = epg_ifac(ijk)*lambda_gt(ijk)
         ELSE
            epmu_gt(ijk) = zero
            eplambda_gt(ijk) = zero
         ENDIF   ! end if (fluid_at(ijk) .or. use_mms)
      ENDDO   ! end do (ijk=ijkstart3,ijkend3)

      RETURN
      END SUBROUTINE set_epmug_values


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Compute a default value of gas viscosity where gas is      C
!  assumed to be air                                                   C
!  Author: W. Sams/M. Syamlal                         Date: 18-JUL-94  C
!                                                                      C
!  Literature/Document References:                                     C
!     Perry, R. H., and Chilton, C. H., Chemical Engineers' Handbook,  C
!        5th Edition, McGraw-Hill Inc., 1973, pp. 248, eqn. 3-133.     C
!     Arnold, J. H., Vapor viscosities and the Sutherland equation,    C
!        Journal of Chemical Physics, 1 (2), 1933, pp. 170-176.        C
!     Sutherland, W., The Viscosity of Gases and Molecular Force,      C
!        Phil. Mag. 5:507-531, 1893.                                   C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_default_mug

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use constant, only: to_si
      use fldvar, only: t_g
      use functions, only: fluid_at
      use param1, only: zero
      use physprop, only: mu_g
      use visc_g, only: mu_gt
      use visc_g, only: lambda_gt
      IMPLICIT NONE

! Local parameters
!---------------------------------------------------------------------//
      DOUBLE PRECISION, PARAMETER :: F2O3 = 2.d0/3.d0

! Local variables
!---------------------------------------------------------------------//
! Cell indices
      INTEGER :: IJK
!---------------------------------------------------------------------//

!!$omp parallel do private(ijk) schedule(dynamic,chunk_size)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN

! Gas viscosity   (in Poise or Pa.s)
! Calculating gas viscosity using Sutherland's formula with
! Sutherland's constant (C) given by Vogel's equation C = 1.47*Tb.
! For air  C = 110 (Tb=74.82)
!         mu = 1.71*10-4 poise at T = 273K
            mu_g(ijk) = to_si*1.7d-4 * &
               (t_g(ijk)/273.0d0)**1.5d0 * (383.d0/(t_g(ijk)+110.d0))

! assign values to quantities used throughout code
            mu_gt(ijk) = mu_g(ijk)
            lambda_gt(ijk) = -f2o3*mu_gt(ijk)
         ELSE
            mu_g(ijk) = zero
            mu_gt(ijk) = zero
            lambda_gt(ijk) = zero
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_default_mug



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: compute turbulent eddy viscosity                           C
!  Author: S. Benyahia                                Date: May-13-04  C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!   Cao, J. and Ahmadi, G., 1995, Gas-particle two-phase turbulent     C
!      flow in a vertical duct. Int. J. Multiphase Flow, vol. 21,      C
!      No. 6, pp. 1203-1228.                                           C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_k_epsilon_mu

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use constant, only: mu_gmax
      use drag, only: f_gs
      use fldvar, only: k_turb_g, e_turb_g, ro_g
      use fldvar, only: ep_s, ro_s
      use functions, only: fluid_at
      use param1, only: zero, one, small_number
      use physprop, only: mu_g
      use run, only: kt_type_enum, ahmadi_1995
      use turb, only: tau_1
      use visc_g, only: mu_gt, lambda_gt
      use visc_s, only: ep_star_array
      IMPLICIT NONE

! Local parameters
!---------------------------------------------------------------------//
      DOUBLE PRECISION, PARAMETER :: F2O3 = 2.d0/3.d0

! Local variables
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER :: M
! Constant in turbulent viscosity formulation
      DOUBLE PRECISION :: C_MU
! particle relaxation time
      DOUBLE PRECISION :: Tau_12_st
! cell index
      INTEGER :: IJK
!---------------------------------------------------------------------//

      c_mu = 9d-02

      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN

! Correction in Ahmadi paper (Cao and Ahmadi)
            IF(kt_type_enum == ahmadi_1995 .AND.&
               f_gs(ijk,1) > small_number) THEN
! solids phase index used throughout routine...
               m = 1 ! for solids phase
               tau_12_st = ep_s(ijk,m)*ro_s(ijk,m)/f_gs(ijk,1)
               c_mu = c_mu/(one+ tau_12_st/tau_1(ijk) * &
                  (ep_s(ijk,m)/(one-ep_star_array(ijk)))**3)
            ENDIF

! I'm not very confident about this correction in Peirano paper,
! but it's made available here, uncomment to use it.
! sof@fluent.com --> 02/01/05
!            IF(KT_TYPE_ENUM==SIMONIN_1996 .AND.&
!               F_GS(IJK,1) > SMALL_NUMBER) THEN
! solids phase index used throughout routine...
!               M = 1 ! for solids phase
!               Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1)
!               X_21 = Ep_s(IJK,M)*RO_S(IJK,M)/(EP_g(IJK)*RO_g(IJK))
! new definition of C_mu (equation A.12, Peirano et al. (2002),
! Powder tech. 122,69-82)
!               IF( K_12(IJK)/(2.0D0*K_Turb_G(IJK)) < ONE) &
!                  C_MU = C_MU/(ONE+ 0.314D0*X_21*Tau_12_st / Tau_1(IJK) * &
!                         (ONE - K_12(IJK)/(2.0D0*K_Turb_G(IJK))) )
!            ENDIF

! Definition of the turbulent viscosity
            mu_gt(ijk) = mu_g(ijk) + ro_g(ijk)*c_mu*&
               k_turb_g(ijk)**2 / (e_turb_g(ijk) + small_number)

            mu_gt(ijk) = min(mu_gmax, mu_gt(ijk))
            lambda_gt(ijk) = -f2o3*mu_gt(ijk)
         ELSE
            mu_gt(ijk) = zero
            lambda_gt(ijk) = zero
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_k_epsilon_mu


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: compute l_scale0 model                                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_lscale_mu

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use constant, only: mu_gmax
      use fldvar, only: ro_g
      use functions, only: fluid_at
      use param1, only: zero
      use physprop, only: mu_g
      use visc_g, only: mu_gt, lambda_gt
      use visc_g, only: l_scale
      IMPLICIT NONE

! Local parameters
!---------------------------------------------------------------------//
      DOUBLE PRECISION, PARAMETER :: F2O3 = 2.d0/3.d0

! Local variables
!---------------------------------------------------------------------//
      DOUBLE PRECISION :: D_g(3,3)
! Gas velocity gradient
      DOUBLE PRECISION :: DelV_g(3,3)
! Second invariant of the deviator of D_g
      DOUBLE PRECISION :: I2_devD_g
! cell index
      INTEGER :: IJK
!---------------------------------------------------------------------//

      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
! Calculate the rate of strain tensor D_g
            CALL calc_deriv_vel_gas(ijk, delv_g, d_g)

! Calculate the second invariant of the deviator of D_g
            i2_devd_g = ((d_g(1,1)-d_g(2,2))**2+(d_g(2,2)-d_g(3,3))**2+&
                         (d_g(3,3)-d_g(1,1))**2)/6.d0 + &
                        d_g(1,2)**2 + d_g(2,3)**2 + d_g(3,1)**2

            mu_gt(ijk) =  mu_g(ijk)+2.0*l_scale(ijk)*l_scale(ijk)*&
                                     ro_g(ijk)*sqrt(i2_devd_g)

            mu_gt(ijk) = min(mu_gmax, mu_gt(ijk))
            lambda_gt(ijk) = -f2o3*mu_gt(ijk)
         ELSE
            mu_gt(ijk) = zero
            lambda_gt(ijk) = zero
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_lscale_mu