kintheory_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module: kintheory                                                   C
!  Purpose: Common block containing constants, variables, functions    C
!  used by various kinetic theory models                               C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Garzo, V., Tenneti, S., Subramaniam, S., and Hrenya, C. M.,         C
!      "Enskog kinetic theory for monodisperse gas-solid flows", JFM,  C
!      Vol. 712, 2012, pp. 129-168                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      MODULE kintheory


! coefficient terms needed for stress (in addition to
! mu_s and lambda_s which are defined in visc_s_mod, allocated
! in allocate_arrays and initialized in set_constprop; and
! P_s which is defined in fldvar_mod, allocated in
! allocate_arrays, and initialized in init_fvars)
!     stress term with gradient in particle M velocity
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: mu_sm_ip
!     stress term with gradient in particle L velocity
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: mu_sl_ip
!     stress term with trace in particle M velocity
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: xi_sm_ip
!     stress term with trace in particle L velocity
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: xi_sl_ip


! coefficient terms needed for momentum source (in addition to
! F_SS which is defined in drag_mod, allocated in allocate_arrays
! and initialized in set_constprop)
!     momentum source term with gradient in number density
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: fnu_s_ip
!     momentum source term with gradient in mixture temperature or
!     with the gradient in temperature of species M
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: ft_sm_ip
!     momentum source term with gradient in temperature of species L
      DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: ft_sl_ip


! coefficient terms needed for heat flux (in addition to
! kth_s which is defined in set_constprop, allocated
! allocate_arrays and initialized in init_fvars)
!     heat flux term with gradient in granular temperature of species L
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: kth_sl_ip
!     heat flux term with gradient in number density of species M
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: knu_sm_ip
!     heat flux term with gradient in number density of species L
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: knu_sl_ip
!     heat flux term with velocity difference
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: kvel_s_ip


! coefficient terms needed for energy dissipation
!     energy dissipation with difference in species granular
!     temperature: transfer between solid solid phases
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: ed_ss_ip
!     energy dissipation term
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: edt_s_ip
!     energy dissipation with divergence of velocity of species M
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: edvel_sm_ip
!     energy dissipation with divergence of velocity of species L
      DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: edvel_sl_ip
!     coefficient A2, xsi used in multiple places in GTSH theory
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: a2_gtsh
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xsi_gtsh

! Solids source terms needed for Iddir & Arastoopour (2005)
! kinetic theory model
      DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE ::  ktmom_u_s
      DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE ::  ktmom_v_s
      DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE ::  ktmom_w_s

! parameter in the theory of GTSH that is related to length scale
! of lubrication effects. For details see GTSH, 2012.
      DOUBLE PRECISION, PARAMETER :: epm = 0.01d0



      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Calculate the magnitude of the gas-solids relative         C
!  velocity at i, j, k                                                 C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION kt_rvel(IJK, m)

! Modules
!---------------------------------------------------------------------//
      use fldvar, only: u_s, v_s, w_s
      use fldvar, only: u_g, v_g, w_g
      use run, only: shear
      use vshear, only: vsh
      use indices, only: i_of
      use functions, only: im_of, jm_of, km_of
      use functions, only: fluid_at
      use fun_avg, only: avg_x_e, avg_y_n, avg_z_t
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! ijk index
      INTEGER, INTENT(IN) :: ijk
! solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! cell indices
      INTEGER :: I, IMJK, IJMK, IJKM
! Cell center value of solids and gas velocities
      DOUBLE PRECISION :: USCM, VSCM, WSCM, &
                          UGC, VGC, WGC
! y-component of velocity with 'shear' applied
      DOUBLE PRECISION :: vs_j, vs_jm
!---------------------------------------------------------------------//
      i = i_of(ijk)
      imjk  = im_of(ijk)
      ijmk  = jm_of(ijk)
      ijkm  = km_of(ijk)

! Awkward here but captures what was done in early version of
! calc_mu_s. Otherwise any call to this routine must be done with
! 'shear' already applied to v_s.
      vs_j = v_s(ijk,m)
      vs_jm = v_s(ijmk,m)
      IF (shear) THEN
         vs_j = v_s(ijk,m)+vsh(ijk)
         IF(fluid_at(ijmk)) THEN
            vs_jm = v_s(ijmk,m)+vsh(ijmk)
         ELSE
            vs_jm = v_s(ijmk,m)
         ENDIF
      ENDIF

! Start calculation for relative velocity
! Calculate velocity components at i, j, k
      ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
      vgc = avg_y_n(v_g(ijmk),v_g(ijk))
      wgc = avg_z_t(w_g(ijkm),w_g(ijk))

      uscm = avg_x_e(u_s(imjk,m),u_s(ijk,m),i)
      vscm = avg_y_n(vs_jm, vs_j)
      wscm = avg_z_t(w_s(ijkm,m),w_s(ijk,m))

! Magnitude of gas-solids relative velocity
      kt_rvel = sqrt((ugc - uscm)**2 + (vgc - vscm)**2 + &
                     (wgc - wscm)**2)

      RETURN
      END FUNCTION kt_rvel

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION kt_cos_theta(ijk, M)

! Modules
!---------------------------------------------------------------------//
      USE param1, only: zero, small_number
      USE fldvar, only: u_g, v_g, w_g
      USE fldvar, only: u_s, v_s, w_s
      USE fldvar, only: ep_s
      USE run, only: shear
      USE vshear, only: vsh
      USE toleranc, only: zero_ep_s
      USE indices, only: i_of
      USE functions, only: im_of, jm_of, km_of
      USE functions, only: fluid_at
      USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! ijk index
      INTEGER, INTENT(IN) :: ijk
! solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! cell indices
      INTEGER :: I, IMJK, IJMK, IJKM
! Cell center value of solids and gas velocities
      DOUBLE PRECISION :: USCM, VSCM, WSCM, &
                          UGC, VGC, WGC
      DOUBLE PRECISION :: vs_j, vs_jm, speed, rvel
!---------------------------------------------------------------------//

      i = i_of(ijk)
      imjk  = im_of(ijk)
      ijmk  = jm_of(ijk)
      ijkm  = km_of(ijk)

! Awkward here but captures what was done in early version of
! calc_mu_s. Otherwise any call to this routine must be done with
! 'shear' already applied to v_s.
      vs_j = v_s(ijk,m)
      vs_jm = v_s(ijmk,m)
      IF (shear) THEN
         vs_j = v_s(ijk,m)+vsh(ijk)
         IF(fluid_at(ijmk)) THEN
            vs_jm = v_s(ijmk,m)+vsh(ijmk)
         ELSE
            vs_jm = v_s(ijmk,m)
         ENDIF
      ENDIF

      ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
      vgc = avg_y_n(v_g(ijmk),v_g(ijk))
      wgc = avg_z_t(w_g(ijkm),w_g(ijk))

      uscm = avg_x_e(u_s(imjk,m),u_s(ijk,m),i)
      vscm = avg_y_n(vs_jm, vs_j)
      wscm = avg_z_t(w_s(ijkm,m),w_s(ijk,m))

      rvel =  sqrt((ugc - uscm)**2 + (vgc - vscm)**2 + &
                   (wgc - wscm)**2)

      speed = sqrt(uscm**2+vscm**2+wscm**2)

! motion viewed by the particles (crossing trajectory effect)
      IF(speed > small_number .AND. rvel > small_number .AND. &
         ep_s(ijk,m) > zero_ep_s) THEN
         kt_cos_theta = ( (ugc-uscm)*uscm + (vgc-vscm)*vscm + &
                       (wgc-wscm)*wscm )/(rvel*speed)
      ELSE
         kt_cos_theta = zero
      ENDIF

      RETURN
      END FUNCTION  kt_cos_theta


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Calculate a single particle drag coefficient/ep_s, which   C
!  is used in evaluting limiting values of certain granular kinetic    C
!  theory terms                                                        C
!                                                                      C
!  Comments: This is currently based on wen-yu but in the the future   C
!  we may want to make this consistent with drag_gs (that is, replace  C
!  this function call and assign appropriate variable within drag_gs   C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION kt_dga(IJK, m)

! Modules
!---------------------------------------------------------------------//
      use param1, only: zero, one, small_number, large_number
      use fldvar, only: rop_g
      use fldvar, only: d_p
      use physprop, only: mu_g
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! ijk index
      INTEGER, INTENT(IN) :: IJK
! solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! single particle drag coefficient, reynolds number
      DOUBLE PRECISION :: C_d, Re
! local value for relative velocity
      DOUBLE PRECISION :: rvel
!---------------------------------------------------------------------//

! initialization
      kt_dga = zero
      rvel = zero
      rvel = kt_rvel(ijk, m)

! Defining single particle drag coefficient dgA based on Wen-Yu
! correlation
      re = d_p(ijk,m)*rvel*rop_g(ijk)/(mu_g(ijk) + small_number)
      IF(re .LE. 1000.d0)THEN
         c_d = (24.d0/(re+small_number)) * &
            (one + 0.15d0 * re**0.687d0)
      ELSE
         c_d = 0.44d0
      ENDIF

! dga_s is local to this routine as is lrvel
      kt_dga = 0.75d0 * c_d * rvel * rop_g(ijk) / d_p(ijk,m)

! set value for 1st iteration and 1st time step
      IF(rvel == zero) kt_dga = large_number

      RETURN
      END FUNCTION kt_dga

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_IA_ENERGY_DISSIPATION_SS                           C
!                                                                      C
!  Purpose: Implement kinetic theory of Iddir & Arastoopour (2005)     C
!     for calculation of source terms in granular energy equation      C
!                                                                      C
!  Author: Janine E. Galvin, Univeristy of Colorado                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_ia_energy_dissipation_ss(M)

! Modules
!---------------------------------------------------------------------//
      USE constant, only: pi
      USE constant, only: c_e
      USE fldvar, only: ro_s, rop_s, d_p
      USE fldvar, only: theta_m
      USE physprop, only: mmax
      USE rdf, only: g_0
      USE param1, only: zero

      USE functions, only: fluid_at, funlm
      USE compar, only: ijkstart3, ijkend3
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! Index
      INTEGER :: IJK
! Solids phase index
      INTEGER :: L
! Index for storing solids-solids drag coefficients
! in the upper triangle of the matrix
      INTEGER :: LM
! variables for IA theory
      DOUBLE PRECISION :: ED_common_term
      DOUBLE PRECISION :: EDvel_sL, EDvel_sM
      DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, &
                          D_PL, DPSUMo2
      DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, R10p_lm, R3p_lm, &
                          R4p_lm, R5p_lm, Bp_lm
!---------------------------------------------------------------------//

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

             d_pm = d_p(ijk,m)
             m_pm = (pi/6.d0)*(d_pm**3)*ro_s(ijk,m)
             nu_pm = rop_s(ijk,m)/m_pm

             DO l = 1, mmax
                lm = funlm(l,m)
                d_pl = d_p(ijk,l)
                m_pl = (pi/6.d0)*(d_pl**3)*ro_s(ijk,l)

                mpsum = m_pm + m_pl
                dpsumo2 = (d_pm+d_pl)/2.d0
                nu_pl = rop_s(ijk,l)/m_pl

                ed_common_term = (3.d0/4.d0)*(dpsumo2*dpsumo2)*&
                   (1.d0+c_e)*g_0(ijk,m,l)*nu_pm*nu_pl*(m_pm*&
                   m_pl/mpsum)*((m_pm*m_pl)**1.5)

                IF (m .eq. l) THEN
                   ap_lm = mpsum/(2.d0)
                   dp_lm = m_pl*m_pm/(2.d0*mpsum)
                   r1p_lm = 1.d0/( (ap_lm**1.5)*(dp_lm**3) )
                   r3p_lm = 1.d0/( (ap_lm**1.5)*(dp_lm**3.5) )

! Dissipation associated with the difference in temperature
! (e.g. interphase transfer term).  For the case case (L=M)
! the term ED_s * (TL-TM) cancels. Therefore, explicity set
! the term to zero for this case.
                   ed_ss_ip(ijk,lm) = zero

! Dissipation associated with temperature
                   edt_s_ip(ijk,m,l) = -ed_common_term* (1.d0-c_e)*&
                      (m_pl/mpsum)*(dsqrt(pi)/6.d0)*r1p_lm*&
                      (theta_m(ijk,m)**1.5)

! Dissipation associated with divergence of
! velocity of solid phase L: do not explicity include
! terms that cancel when EDvel_sL is summed with EDvel_sM
                   edvel_sl = ed_common_term*((1.d0-c_e)*(m_pl/mpsum)*&
                      (dpsumo2*pi/48.d0)*(m_pm*m_pl/mpsum)*r3p_lm)
                      edvel_sl_ip(ijk,m,l) = edvel_sl*(theta_m(ijk,l))

! Dissipation associated with divergence of
! velocity of solid phase M: do not explicity include
! terms that cancel when EDvel_sL is summed with EDvel_sM
! commented by sof, no need to re-do computation
                   !EDvel_sM = ED_common_term*((1.d0-C_E)*(M_PL/MPSUM)*&
                   !   (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm)
                   !EDvel_sM_ip(IJK,M,L) = EDvel_sM*(Theta_m(IJK,M)) ! M=L

                   edvel_sm_ip(ijk,m,l) =  edvel_sl_ip(ijk,m,l)

                ELSE

                   ap_lm = (m_pm*theta_m(ijk,l)+m_pl*theta_m(ijk,m))/&
                      2.d0
                   bp_lm = (m_pm*m_pl*(theta_m(ijk,l)-&
                      theta_m(ijk,m) ))/(2.d0*mpsum)
                   dp_lm = (m_pl*m_pm*(m_pm*theta_m(ijk,m)+m_pl*&
                      theta_m(ijk,l) ))/(2.d0*mpsum*mpsum)

                   r1p_lm = (1.d0/((ap_lm**1.5)*(dp_lm**3)))+ &
                      ((9.d0*bp_lm*bp_lm)/(ap_lm**2.5 * dp_lm**4))+&
                      ((30.d0*bp_lm**4)/(2.d0*ap_lm**3.5 * dp_lm**5))

                   r3p_lm = (1.d0/((ap_lm**1.5)*(dp_lm**3.5)))+&
                      ((21.d0*bp_lm*bp_lm)/(2.d0 * ap_lm**2.5 * dp_lm**4.5))+&
                      ((315.d0*bp_lm**4)/(8.d0 * ap_lm**3.5 *dp_lm**5.5))

                   r4p_lm = (3.d0/( ap_lm**2.5 * dp_lm**3.5))+&
                      ((35.d0*bp_lm*bp_lm)/(2.d0 * ap_lm**3.5 * dp_lm**4.5))+&
                      ((441.d0*bp_lm**4)/(8.d0 * ap_lm**4.5 * dp_lm**5.5))

                   r5p_lm = (1.d0/(ap_lm**2.5 * dp_lm**3))+ &
                      ((5.d0*bp_lm*bp_lm)/(ap_lm**3.5 * dp_lm**4))+&
                      ((14.d0*bp_lm**4)/( ap_lm**4.5 * dp_lm**5))

                   r10p_lm = (1.d0/(ap_lm**2.5 * dp_lm**2.5))+&
                      ((25.d0*bp_lm*bp_lm)/(2.d0* ap_lm**3.5 * dp_lm**3.5))+&
                      ((1225.d0*bp_lm**4)/(24.d0* ap_lm**4.5 * dp_lm**4.5))

! Dissipation associated with the difference in temperature
! (e.g. interphase transfer term). to solved using PEA
                   ed_ss_ip(ijk,lm) = ed_common_term*dsqrt(pi)*(m_pm*m_pl/&
                      (2.d0*mpsum))*r5p_lm*( (theta_m(ijk,m)*theta_m(ijk,l))**3 )

! Dissipation associated with temperature
                   edt_s_ip(ijk,m,l) = -ed_common_term* (1.d0-c_e)*(m_pl/mpsum)*&
                      (dsqrt(pi)/6.d0)*r1p_lm*( (theta_m(ijk,m)*theta_m(ijk,l))**3 )

! Dissipation associated with divergence of
! velocity of solid phase L
                   edvel_sl = ed_common_term*( ((3.d0*dpsumo2*pi/40.d0)*m_pl*&
                      r10p_lm)+( (dpsumo2*pi/4.d0)*(m_pm*m_pl/mpsum)*bp_lm*&
                      r4p_lm)-( (1.d0-c_e)*(m_pl/mpsum)*(dpsumo2*pi/16.d0)*&
                      m_pl*bp_lm*r4p_lm)+( (1.d0-c_e)*(m_pl/mpsum)*&
                      (dpsumo2*pi/48.d0)*(m_pm*m_pl/mpsum)*r3p_lm))

                   edvel_sl_ip(ijk,m,l) = edvel_sl*( theta_m(ijk,m)**3.5 *&
                      theta_m(ijk,l)**2.5 )

! Dissipation associated with divergence of
! velocity of solid phase M
                   edvel_sm = ed_common_term*( (-(3.d0*dpsumo2*pi/40.d0)*m_pm*&
                      r10p_lm)+( (dpsumo2*pi/4.d0)*(m_pm*m_pl/mpsum)*bp_lm*&
                      r4p_lm)+( (1.d0-c_e)*(m_pl/mpsum)*(dpsumo2*pi/16.d0)*&
                       m_pm*bp_lm*r4p_lm)+( (1.d0-c_e)*(m_pl/mpsum)*&
                      (dpsumo2*pi/48.d0)*(m_pm*m_pl/mpsum)*r3p_lm))

                   edvel_sm_ip(ijk,m,l) = edvel_sm*( theta_m(ijk,m)**2.5 *&
                      theta_m(ijk,l)**3.5 )

                ENDIF   ! end if/else (m.eq.l)

             ENDDO   ! end do (l=1,mmax)

         ENDIF   ! end if(fluid_at)
      ENDDO   ! end do ijk

      RETURN
      END SUBROUTINE calc_ia_energy_dissipation_ss



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_GD_99_ENERGY_DISSIPATION_SS                        C
!                                                                      C
!  Purpose: Implement kinetic theory of Garzo & Dufty (1999)           C
!  for calculation of source terms in granular energy equation         C
!                                                                      C
!  Author: Janine E. Galvin                                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_gd_99_energy_dissipation_ss(M)

! Modules
!---------------------------------------------------------------------//
      USE constant, only: pi
      USE constant, only: c_e
      USE fldvar, only: rop_s, d_p
      USE fldvar, only: theta_m, ep_s
      USE rdf, only: g_0

      USE functions, only: fluid_at
      USE compar, only: ijkstart3, ijkend3
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK
! variables for GD model
      DOUBLE PRECISION :: press_star, c_star, zeta0_star, &
                          nu_gamma_star, &
                          lambda_num, cd_num, zeta1
      DOUBLE PRECISION :: D_PM, EP_SM
      DOUBLE PRECISION :: nu0, Chi
!---------------------------------------------------------------------//

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

! Note: k_boltz = M_PM

! local aliases
             chi = g_0(ijk,m,m)
             ep_sm = ep_s(ijk,m)
             d_pm = d_p(ijk,m)

!            nu0=p_k/eta0
             nu0 = (96.d0/5.d0)*(ep_sm/d_pm)*dsqrt(theta_m(ijk,m)/pi)

             press_star = 1.d0 + 2.d0*(1.d0+c_e)*ep_sm*chi

             c_star = 32.0d0*(1.0d0 - c_e)*(1.d0 - 2.0d0*c_e*c_e) &
                / (81.d0 - 17.d0*c_e + 30.d0*c_e*c_e*(1.0d0-c_e))

             zeta0_star = (5.d0/12.d0)*chi*(1.d0 - c_e*c_e) &
                * (1.d0 + (3.d0/32.d0)*c_star)

             nu_gamma_star = ((1.d0+c_e)/48.d0)*chi*(128.d0-96.d0*c_e + &
                15.d0*c_e*c_e-15.d0*c_e*c_e*c_e+ (c_star/64.d0) * (15.d0* &
                c_e*c_e*c_e-15.d0*c_e*c_e+498.d0*c_e-434.d0))

             lambda_num = (3.d0/8.d0)*( (1.d0-c_e)*(5.d0*c_e*c_e+4.d0*c_e-1.d0)+ &
                (c_star/12.d0)*(159.d0*c_e+3.d0*c_e*c_e-19.d0*c_e- &
                15.d0*c_e*c_e*c_e) ) * (1.d0+c_e)

! does not include factor of 1/nu0.  the factor 1/nu0 in cD will cancel
! with the factor of nu0 used to obtain zeta1 from zeta1_star
             cd_num = ( (4.d0/15.d0)*lambda_num*ep_sm*chi + &
                (press_star-1.d0)*(2.d0/3.d0-c_e)*c_star ) / &
                ( 0.5d0*zeta0_star+nu_gamma_star + (5.d0*c_star/64.d0) * &
                (1.d0+(3.d0*c_star/64.d0))*chi*(1.d0-c_e*c_e))

! does not include factor of 1/nu0 in the first term.  the factor 1/nu0
! will cancel with the factor of nu0 used to obtain zeta1 from zeta1_star
! zeta1 = nu0*zeta1_star
             zeta1 = -(1.d0-c_e)*(press_star-1.d0) + (5.d0/32.d0) * &
                (1.d0-c_e*c_e)*(1.d0+(3.d0*c_star/64.d0))*chi*cd_num

! in the energy equation this term is multiplied by 3/2*n*kboltz*T
! leave multiplication of theta for source routine
             edvel_sm_ip(ijk,m,m) = (3.d0/2.d0)*rop_s(ijk,m)*zeta1

! in the energy equation this term is multiplied by 3/2*n*kboltz*T
! leave multiplication of theta for source routine
             edt_s_ip(ijk,m,m) = (3.d0/2.d0)*rop_s(ijk,m)*nu0*zeta0_star

         ENDIF   ! end if (fluid_at)
      ENDDO   ! end do (ijk=ijkstart3,ijkend3)

      RETURN
      END SUBROUTINE calc_gd_99_energy_dissipation_ss


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_GTSH_ENERGY_DISSIPATION_SS                         C
!                                                                      C
!  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
!  Hrenya (2012) for calculation of source terms in granular           C
!  energy equation                                                     C
!                                                                      C
!  Author: Sofiane Benyahia                                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_gtsh_energy_dissipation_ss(M)

! Modules
!---------------------------------------------------------------------//
      USE constant, only: pi
      USE constant, only: c_e
      USE fldvar, only: ep_s
      USE fldvar, only: ro_g, rop_g
      USE fldvar, only: ro_s, d_p
      USE fldvar, only: theta_m
      USE physprop, only: mu_g
      USE rdf, only: g_0
      USE param1, only: zero, one, small_number

      USE functions, only: fluid_at
      USE compar, only: ijkstart3, ijkend3
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK
! variables
      DOUBLE PRECISION :: D_PM, EP_SM, V_p, N_p, M_p
      DOUBLE PRECISION :: nu0, Chi
      DOUBLE PRECISION :: VREL
      DOUBLE PRECISION :: Re_m, Re_T
      DOUBLE PRECISION :: zeta_star, mu2_0, mu4_0, mu4_1
      DOUBLE PRECISION :: omega, nu_j, rho_10, rho_11

!---------------------------------------------------------------------//

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

! Local aliases
            chi = g_0(ijk,m,m)
            ep_sm = ep_s(ijk,m)
            d_pm = d_p(ijk,m)
            v_p = pi*d_pm**3/6d0
            n_p = ep_sm/v_p
            m_p = v_p * ro_s(ijk,m)

            nu0 = (96.d0/5.d0)*(ep_sm/d_pm)*dsqrt(theta_m(ijk,m)/pi)

! First calculate the Re number (Re_m, Re_T) in eq. 3.1 in GTSH theory
            vrel = kt_rvel(ijk, m)
! Note: rop_g = ro_g * (1-phi)
            re_m = d_pm*vrel*rop_g(ijk)/mu_g(ijk)
            re_t = ro_g(ijk)*d_pm*dsqrt(theta_m(ijk,m)) / mu_g(ijk)

! Now calculate and store xsi (used in many spots) in eq. 8.2 in
! GTSH theory.
            xsi_gtsh(ijk) = one/6d0*d_pm*vrel**2*(3d0*pi*mu_g(ijk)*&
               d_pm/m_p)**2 / dsqrt(pi*theta_m(ijk,m)) * &
               s_star(ep_sm, chi)

! eq. (6.22) GTSH theory
            mu2_0 = dsqrt(2d0*pi) * chi * (one-c_e**2)

! eq. (6.23) GTSH theory
            mu4_0 = (4.5d0+c_e**2) * mu2_0

! eq. (6.24) GTSH theory
!            mu4_1 = (6.46875d0+0.3125d0*C_E**2 + 2d0/(one-C_E)) * mu2_0
! this is done to avoid /0 in case c_e = 1.0
            mu4_1 = (6.46875d0+0.9375d0*c_e**2)*mu2_0 + &
               2d0*dsqrt(2d0*pi)* chi*(one+c_e)

            a2_gtsh(ijk) = zero ! for EP_SM = zero
            if(ep_sm> small_number) then ! avoid singularity
! Now calculate zeta_star in eq. 8.10 in GTSH theory
               zeta_star = 4.5d0*dsqrt(2d0*pi)*&
                  (ro_g(ijk)/ro_s(ijk,m))**2*re_m**2 * &
                  s_star(ep_sm,chi) / (ep_sm*(one-ep_sm)**2 * re_t**4)

! Now calculate important parameter A2. This is used in many transport
! coefficients so that it's best to store it in an array instead of
! having many function calls.
               a2_gtsh(ijk) = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
                  (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
            endif

! Now calculate first term rho_0 in cooling rate rho, eq (6.26) in GTSH
! theory.
! note that theta_(ijk,m) does't contain mass m, so that T/m = theta_m.
! note that edt_s_ip will need to be multiplied by 3/2 rop_s(ijk,m) in
! source_granular_energy to have the same meaning as in GD_99 theory.
            edt_s_ip(ijk,m,m) = 4d0/3d0*dsqrt(pi)*(one-c_e**2)*chi* &
               (one+0.1875d0*a2_gtsh(ijk))*n_p*d_pm**2*dsqrt(theta_m(ijk,m))

! Calculate the second term in cooling rate, eq (7.25) in GTSH theory.
! First calculate eq (7.28, 7.29) of GTSH theory.
            omega = (one+c_e)*nu0*((one-c_e**2)*(5d0*c_e-one) - &
               a2_gtsh(ijk)/ 6d0 * (15d0*c_e**3-3d0*c_e**2+81d0*c_e-61d0))
            nu_j = (one+c_e)/192d0*chi*nu0* &
               (241d0-177d0*c_e+30d0*c_e**2-30d0*c_e**3)

! Now calculate eq (7.26, 7.27) of GTSH theory. ! corrected by W. Fullmer
            rho_10 = 2d0*chi*ep_sm*(c_e**2-one)
            rho_11 = 25d0/1024d0*ep_sm*chi**2*(one-c_e**2)* &
               (one+3d0/128d0*a2_gtsh(ijk)) * (omega/10d0 - &
               (one+c_e)*nu0*(one/3d0-c_e)*a2_gtsh(ijk)/2d0) / &
               (nu_j + g_gtsh(ep_sm, chi, ijk, m)/m_p + 1.5d0* &
               xsi_gtsh(ijk)/theta_m(ijk,m) -1.5d0*edt_s_ip(ijk,m,m))

! note that EDvel_sM_ip will later need to be multiplied by 3/2 rop_s(ijk,m).
            edvel_sm_ip(ijk,m,m) = rho_10 + rho_11

         ENDIF   ! end if (fluid_at)
      ENDDO   ! end do (ijk=ijkstart3,ijkend3)

      RETURN
      END SUBROUTINE calc_gtsh_energy_dissipation_ss


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Functions: G_tsh, K_phi, R_d, S_star                                C
!                                                                      C
!  Purpose: Implement kinetic theory of Garzo, Tenneti, Subramaniam    C
!  Hrenya (2012) for calculation of source terms in granular           C
!  energy equation                                                     C
!                                                                      C
!  Author: Sofiane Benyahia                                            C
!                                                                      C
!  Comments:                                                           C
!     these functions are needed only for gtsh theory                  C
!     Function gamma, eq. (8.1) in GTSH theory                         C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION g_gtsh (EP_SM, Chi, IJK, M)

! Modules
!---------------------------------------------------------------------//
      USE constant, only: pi
      USE param1, only: one
      USE physprop, only: mu_g
      USE fldvar, only: ro_g
      USE fldvar, only: d_p, theta_m
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids volume fraction of phase m at ijk
      DOUBLE PRECISION, INTENT(IN) :: EP_SM
! radial distribution function of phase M at ijk
      DOUBLE PRECISION, INTENT(IN) :: Chi
! Index
      INTEGER, INTENT(IN) :: IJK
! Solids phase index. Note M should be equal to 1 since theory
! valid for only mmax = 1.
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
      DOUBLE PRECISION :: Re_T
      DOUBLE PRECISION :: Rdiss, RdissP

!---------------------------------------------------------------------//

      if(ep_sm <= 0.1d0) then
         rdissp = one+3d0*dsqrt(ep_sm/2d0)
      else
         rdissp = &
         one + 3d0*dsqrt(ep_sm/2d0) + 135d0/64d0*ep_sm*dlog(ep_sm) + &
         11.26d0*ep_sm*(one-5.1d0*ep_sm+16.57d0*ep_sm**2-21.77d0*    &
         ep_sm**3) - ep_sm*chi*dlog(epm)
      endif

      re_t = ro_g(ijk)*d_p(ijk,m)*dsqrt(theta_m(ijk,m)) / mu_g(ijk)

      rdiss = rdissp + re_t * k_phi(ep_sm)

      g_gtsh = 3d0*pi*mu_g(ijk)*d_p(ijk,m)*rdiss  ! eq. (8.1) in GTSH theory

      RETURN
      END FUNCTION g_gtsh

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Comments:                                                           C
!     Function K(phi), eq. (2.7) in Yi/Zenk/Mitrano/Hrenya JFM (2013)  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION k_phi (phi)
      IMPLICIT NONE

! Dummy Arguments
!---------------------------------------------------------------------//
! solids volume fraction
      DOUBLE PRECISION, INTENT(IN) :: phi
!---------------------------------------------------------------------//

      k_phi = (0.096d0 + 0.142d0*phi**0.212d0) / (1d0-phi)**4.454d0
      k_phi = 0.0d0 ! set to zero for compatibility with GTSH JFM (2012)
      RETURN
      END FUNCTION k_phi

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Comments:                                                           C
!     Function R_d, eq. (8.6) in GTSH theory                           C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION r_d (phi)
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids volume fraction
      DOUBLE PRECISION, INTENT(IN) :: phi
!---------------------------------------------------------------------//

      r_d = 1.0d0  ! this avoids singularity at phi = 0.0
      if((phi > 1d-15) .and. (phi <= 0.4d0)) then
        r_d = (1d0+3d0*dsqrt(phi/2d0)+135d0/64d0*phi*dlog(phi)+17.14d0*phi) / &
              (1d0+0.681d0*phi-8.48*phi**2+8.16d0*phi**3)
      elseif(phi > 0.4d0) then
        r_d = 10d0*phi/(1d0-phi)**3 + 0.7d0
      endif
      RETURN
      END FUNCTION r_d

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Comments:                                                           C
!     Function S_Star(phi), eq. (8.3, 8.5) in GTSH theory              C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      DOUBLE PRECISION FUNCTION s_star (phi, Chi)
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids volume fraction
      DOUBLE PRECISION, INTENT(IN) :: phi
! radial distribution function
      DOUBLE PRECISION, INTENT(IN) :: Chi

!---------------------------------------------------------------------//

      s_star = 1.0d0
      if(phi >= 0.1d0) &
        s_star = r_d(phi)**2/(chi*(1d0+3.5d0*dsqrt(phi)+5.9*phi))
      RETURN
      END FUNCTION s_star

      END MODULE kintheory