source_granular_energy.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_GRANULAR_ENERGY                                  C
!  Purpose: Calculate the source terms in the granular energy equation C
!  Note: The center coefficient and source vector are negative.        C
!  The off-diagonal coefficients are positive.                         C
!                                                                      C
!  Author: Kapil Agrawal, Princeton University        Date: 04-FEB-98  C
!  Reviewer: M. Syamlal                               Date:            C
!                                                                      C
!  Revision Number:1                                                   C
!  Purpose: Add Simonin and Ahmadi models                              C
!  Author: Sofiane Benyahia, Fluent Inc.               Date: 02-01-05  C
!                                                                      C
!  Literature/Document References:                                     C
!     Lun, C.K.K., S.B. Savage, D.J. Jeffrey, and N. Chepurniy,        C
!        Kinetic theories for granular flow - inelastic particles in   C
!        Couette-flow and slightly inelastic particles in a general    C
!        flow field. Journal of Fluid Mechanics, 1984. 140(MAR):       C
!        p. 223-256                                                    C
!     Gidaspow, D., Multiphase flow and fluidziation, 1994, Academic   C
!        Press Inc., California, Chapter 9.                            C
!     Koch, D. L., and Sangani, A. S., Particle pressure and marginal  C
!        stability limits for a homogeneous monodisperse gas-fluidized C
!        bed: kinetic theory and numerical simulations, Journal of     C
!        Fluid Mechanics, 1999, 400, 229-263.                          C
!                                                                      C
!     Simonin, O., 1996. Combustion and turbulence in two-phase flows, C
!        Von Karman institute for fluid dynamics, lecture series,      C
!        1996-02                                                       C
!     Balzer, G., Simonin, O., Boelle, A., and Lavieville, J., 1996,   C
!        A unifying modelling approach for the numerical prediction    C
!        of dilute and dense gas-solid two phase flow. CFB5, 5th int.  C
!        conf. on circulating fluidized beds, Beijing, China.          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
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE source_granular_energy(SOURCELHS, &
                    sourcerhs, ijk, m)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE constant
      USE drag
      USE fldvar
      USE fun_avg
      USE functions
      USE geometry
      USE indices
      USE kintheory
      USE mms
      USE parallel
      USE param
      USE param1
      USE physprop
      USE rdf
      USE run
      USE solids_pressure
      USE toleranc
      USE trace
      USE turb
      USE visc_g
      USE visc_s

      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Source terms to be kept on rhs (source vector)
      DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
! Source terms to be kept on lhs (center coefficient)
      DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
! Solid phase index
      INTEGER, INTENT(IN) :: M
! Indices
      INTEGER, INTENT(IN) :: IJK
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM, &
                 IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
! Solid phase index
      INTEGER :: MM
! Particle relaxation time
      DOUBLE PRECISION :: Tau_12_st
! Sum of eps*G_0
      DOUBLE PRECISION :: SUM_EpsGo
! Slip velocity
      DOUBLE PRECISION :: VSLIP
! coefficients in heat flux term
      DOUBLE PRECISION :: Knu_e, Knu_w, Knu_n, Knu_s, &
                          Knu_t, Knu_b
! Source terms to be kept on rhs
      DOUBLE PRECISION :: S1_rhs, S2_rhs, S3_rhs, &
                          S5a_rhs, S5b_rhs, S5c_rhs, S5_rhs, &
                          S7_rhs
! Source terms to be kept on lhs
      DOUBLE PRECISION :: S1_lhs, S3_lhs, S4_lhs, S6_lhs
!-----------------------------------------------

      i = i_of(ijk)
      j = j_of(ijk)
      k = k_of(ijk)
      im = im1(i)
      jm = jm1(j)
      km = km1(k)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijke = east_of(ijk)
      ijkw = west_of(ijk)
      ijkn = north_of(ijk)
      ijks = south_of(ijk)
      ijkt = top_of(ijk)
      ijkb = bottom_of(ijk)

! initialize summation variables
      s1_rhs = zero
      s1_lhs = zero
      s2_rhs = zero
      s3_rhs = zero
      s3_lhs = zero
      s4_lhs = zero
      s6_lhs = zero
      s7_rhs = zero

! Changes needed for multitype particles, sof June 16 2005
! Sum of eps*G_0 is used instead of Eps*G_0
      sum_epsgo = zero
      DO mm = 1, mmax
         sum_epsgo =  sum_epsgo+ep_s(ijk,mm)*g_0(ijk,mm,mm)
      ENDDO

! Production by shear: (S:grad(vi))
! Pi_s*tr(Di)     (Lun et al. 1984)
      s1_rhs = p_s_c(ijk,m)*zmax(( -trd_s_c(ijk,m) ))
      s1_lhs = p_s_c(ijk,m)*zmax(( trd_s_c(ijk,m) ))

! Production by shear: (S:grad(vi))
! Mu_s*tr(Di^2)     (Lun et al. 1984)
      s2_rhs = 2.d0*mu_s_c(ijk,m)*trd_s2(ijk,m)

! Production by shear: (S:grad(vi))
! Lambda_s*tr(Di)^2     (Lun et al. 1984)
      s3_rhs = (trd_s_c(ijk,m)**2)*zmax( lambda_s_c(ijk,m) )
      s3_lhs = (trd_s_c(ijk,m)**2)*zmax( -lambda_s_c(ijk,m) )

! Energy dissipation by collisions
      s4_lhs = (48.d0/dsqrt(pi))*eta*(one-eta)*rop_s(ijk,m)*&
         sum_epsgo*dsqrt(theta_m(ijk,m))/d_p(ijk,m)

! Need to revisit this part about granular energy MMS source terms.
      IF (use_mms) s4_lhs = zero

! Energy dissipation by viscous dampening
! Gidaspow (1994)  : addition due to role of interstitial fluid
      IF(kt_type_enum == simonin_1996 .OR. &
         kt_type_enum == ahmadi_1995) THEN
         s6_lhs = 3.d0 * f_gs(ijk,m)
      ELSE IF(switch > zero .AND. ro_g0 /= zero) THEN
         s6_lhs = switch * 3.d0 * f_gs(ijk,m)
      ENDIF

! Energy production due to gas-particle slip
! Koch & Sangani (1999) : addition due to role of interstitial fluid
      IF(kt_type_enum == simonin_1996) THEN
         s7_rhs = f_gs(ijk,m)*k_12(ijk)
      ELSEIF(kt_type_enum == ahmadi_1995) THEN
! note specific reference to F_GS of solids phase 1!
         IF(ep_s(ijk,m) > dil_ep_s .AND. f_gs(ijk,1) > small_number) THEN
            tau_12_st = ep_s(ijk,m)*ro_s(ijk,m)/f_gs(ijk,1)
            s7_rhs = 2.d0*f_gs(ijk,m) * (one/(one+tau_12_st/  &
               (tau_1(ijk)+small_number)))*k_turb_g(ijk)
         ELSE
            s7_rhs = zero
         ENDIF
      ELSEIF(switch > zero .AND. ro_g0 /= zero) THEN
         vslip = dsqrt( (u_s(ijk,m)-u_g(ijk))**2 + &
            (v_s(ijk,m)-v_g(ijk))**2 + (w_s(ijk,m)-w_g(ijk))**2 )
         s7_rhs = switch*81.d0*ep_s(ijk,m)*(mu_g(ijk)*vslip)**2 /&
            (g_0(ijk,m,m)*d_p(ijk,m)**3 * ro_s(ijk,m)*&
            dsqrt(pi)*dsqrt( theta_m(ijk,m)+small_number ) )
! Need to revisit this part about granular energy MMS source terms.
         IF (use_mms) s7_rhs = zero

      ENDIF


! The following lines are essentially commented out since Kphi_s has
! been set to zero in subroutine calc_mu_s.  To activate the feature
! activate the following lines and the lines in calc_mu_s.
! Part of Heat Flux: div (q)
! Kphi_s*grad(eps)     (Lun et al. 1984)
      IF (.false.) THEN
         knu_e = avg_x_h(kphi_s(ijk,m), kphi_s(ijke,m),i)
         knu_w = avg_x_h(kphi_s(ijkw,m),kphi_s(ijk,m),im)
         knu_n = avg_y_h(kphi_s(ijk,m), kphi_s(ijkn,m),j)
         knu_s = avg_y_h(kphi_s(ijks,m),kphi_s(ijk,m),jm)
         knu_t = avg_z_h(kphi_s(ijk,m), kphi_s(ijkt,m),k)
         knu_b = avg_z_h(kphi_s(ijkb,m),kphi_s(ijk,m),km)

         s5a_rhs = knu_e*(ep_s(ijke,m)-ep_s(ijk,m))*&
            odx_e(i)*ayz(ijk) - &
            knu_w*(ep_s(ijk,m)-ep_s(ijkw,m))*&
            odx_e(im)*ayz(imjk)
         s5b_rhs = knu_n*(ep_s(ijkn,m)-ep_s(ijk,m))*&
            ody_n(j)*axz(ijk) - &
            knu_s*(ep_s(ijk,m)-ep_s(ijks,m))*&
            ody_n(jm)*axz(ijmk)
         s5c_rhs = knu_t*(ep_s(ijkt,m)-ep_s(ijk,m))*&
            ox(i)*odz_t(k)*axy(ijk) - &
            knu_b*(ep_s(ijk,m)-ep_s(ijkb,m))*&
            ox(i)*odz_t(km)*axy(ijkm)
         s5_rhs = s5a_rhs + s5b_rhs + s5c_rhs
      ENDIF
      s5_rhs = zero


      sourcerhs = (s1_rhs + s2_rhs + s3_rhs + s7_rhs)*vol(ijk) + &
           s5_rhs

      sourcelhs = ( ((s1_lhs + s3_lhs)/(theta_m(ijk,m)+small_number)) + &
          s4_lhs + s6_lhs) * vol(ijk)

      RETURN
      END SUBROUTINE source_granular_energy




!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_IA_GRANULAR_ENERGY                               C
!  Purpose: Calculate the source terms of the Mth solids phase         C
!           granular energy equation                                   C
!                                                                      C
!  Author: Janine E. Galvin, Univeristy of Colorado                    C
!                                                                      C
!  Literature/Document References:                                     C
!    Iddir, Y.H., Modeling of the multiphase mixture of particles      C
!       using the kinetic theory approach, PhD Thesis, Illinois        C
!       Institute of Technology, Chicago, Illinois, 2004               C
!    Iddir, Y.H., & H. Arastoopour, Modeling of multitype particle     C
!       flow using the kinetic theory approach, AIChE J., Vol 51,      C
!       No 6, June 2005                                                C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE source_granular_energy_ia(SOURCELHS, &
                    sourcerhs, ijk, m)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE physprop
      USE run
      USE drag
      USE geometry
      USE fldvar
      USE visc_g
      USE visc_s
      USE trace
      USE turb
      USE indices
      USE constant
      USE toleranc
      USE residual
      use kintheory
      USE compar
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Source terms to be kept on rhs (source vector)
      DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
! Source terms to be kept on lhs (center coefficient)
      DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
! Solid phase index
      INTEGER, INTENT(IN) :: M
! Indices
      INTEGER, INTENT(IN) :: IJK
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
                 IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
! phase index
      INTEGER :: L, LM
! velocities
      DOUBLE PRECISION :: UsM_e, UsM_w, VsM_n, VsM_s, WsM_t, WsM_b,&
                          UsL_e, UsL_w, VsL_n, VsL_s, WsL_t, WsL_b,&
                          UsM_p, VsM_p, WsM_P, UsL_p, VsL_p, WsL_p
! number densities
      DOUBLE PRECISION :: NU_PM_E, NU_PM_W, NU_PM_N, NU_PM_S, NU_PM_T,&
                          NU_PM_B, NU_PM_p,&
                          NU_PL_E, NU_PL_W, NU_PL_N, NU_PL_S, NU_PL_T,&
                          NU_PL_B, NU_PL_p
! temperature of species L
      DOUBLE PRECISION :: T_PL_E, T_PL_W, T_PL_N, T_PL_S, T_PL_T,&
                          T_PL_B, T_PL_p
! particle characteristics
      DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL
! coefficients in heat flux term
      DOUBLE PRECISION :: Knu_sL_e, Knu_sL_w, Knu_sL_n, Knu_sL_s, Knu_sL_t,&
                          Knu_sL_b, Knu_sM_e, Knu_sM_w, Knu_sM_n, Knu_sM_s,&
                          Knu_sM_t, Knu_sM_b,&
                          Kvel_s_e, Kvel_s_w, Kvel_s_n, Kvel_s_s, Kvel_s_t,&
                          Kvel_s_b,&
                          Kth_sL_e, Kth_sL_w, Kth_sL_n, Kth_sL_s, Kth_sL_t,&
                          Kth_sL_b
! Source terms to be kept on rhs
      DOUBLE PRECISION :: S10_rhs, S15_rhs, S16_rhs,&
                          S11_sum_rhs, S12_sum_rhs,&
                          S13_sum_rhs, S17_sum_rhs, S18_sum_rhs
! Source terms
      DOUBLE PRECISION :: S14a_sum, S14b_sum, S14c_sum, &
                          S9_sum, s21a_sum, s21b_sum, s21c_sum
! Source terms to be kept on lhs
      DOUBLE PRECISION :: S10_lhs, S16_lhs,&
                          S11_sum_lhs, S12_sum_lhs, S13_sum_lhs,&
                          S17_sum_lhs, S18_sum_lhs, S20_sum_lhs
!-----------------------------------------------

      i = i_of(ijk)
      j = j_of(ijk)
      k = k_of(ijk)
      im = im1(i)
      jm = jm1(j)
      km = km1(k)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijke = east_of(ijk)
      ijkw = west_of(ijk)
      ijkn = north_of(ijk)
      ijks = south_of(ijk)
      ijkt = top_of(ijk)
      ijkb = bottom_of(ijk)

! initialize summation variables
      s9_sum = zero
      s14a_sum = zero
      s14b_sum = zero
      s14c_sum = zero
      s21a_sum = zero
      s21b_sum = zero
      s21c_sum = zero

      s11_sum_rhs = zero
      s12_sum_rhs = zero
      s13_sum_rhs = zero
      s17_sum_rhs = zero
      s18_sum_rhs = zero

      s11_sum_lhs = zero
      s12_sum_lhs = zero
      s13_sum_lhs = zero
      s17_sum_lhs = zero
      s18_sum_lhs = zero
      s20_sum_lhs = zero

      usm_e = u_s(ijk,m)
      usm_w = u_s(imjk,m)
      vsm_n = v_s(ijk,m)
      vsm_s = v_s(ijmk,m)
      wsm_t = w_s(ijk,m)
      wsm_b = w_s(ijkm,m)
      usm_p = avg_x_e(u_s(imjk,m),u_s(ijk,m),i)
      vsm_p = avg_y_n(v_s(ijmk,m),v_s(ijk,m) )
      wsm_p = avg_z_t(w_s(ijkm,m),w_s(ijk,m) )

      d_pm = d_p(ijk,m)
      m_pm = (pi/6.d0)*d_pm**3 * ro_s(ijk,m)
      nu_pm_p = rop_s(ijk,m)/m_pm
      nu_pm_e = rop_s(ijke,m)/m_pm
      nu_pm_w = rop_s(ijkw,m)/m_pm
      nu_pm_n = rop_s(ijkn,m)/m_pm
      nu_pm_s = rop_s(ijks,m)/m_pm
      nu_pm_t = rop_s(ijkt,m)/m_pm
      nu_pm_b = rop_s(ijkb,m)/m_pm

! Production by shear: (S:grad(vi))
! Pi_s*tr(Di)
      s10_lhs = p_s_c(ijk,m) * zmax(trd_s_c(ijk,m))
      s10_rhs = p_s_c(ijk,m) * zmax(-trd_s_c(ijk,m))


! Production by shear: (S:grad(vi))
! Mu_s*tr(Di^2)
      s15_rhs = 2.d0*mu_s_c(ijk,m)*trd_s2(ijk,m)

! Production by shear: (S:grad(vi))
! Lambda_s*tr(Di)^2
      s16_lhs = (trd_s_c(ijk,m)**2)*zmax( -lambda_s_c(ijk,m) )
      s16_rhs = (trd_s_c(ijk,m)**2)*zmax(  lambda_s_c(ijk,m) )

      DO l = 1, mmax
          d_pl = d_p(ijk,l)
          m_pl = (pi/6.d0)*d_pl**3 * ro_s(ijk,l)
          nu_pl_p = rop_s(ijk,l)/m_pl
          nu_pl_e = rop_s(ijke,l)/m_pl
          nu_pl_w = rop_s(ijkw,l)/m_pl
          nu_pl_n = rop_s(ijkn,l)/m_pl
          nu_pl_s = rop_s(ijks,l)/m_pl
          nu_pl_t = rop_s(ijkt,l)/m_pl
          nu_pl_b = rop_s(ijkb,l)/m_pl

          t_pl_p = theta_m(ijk,l)
          t_pl_e = theta_m(ijke,l)
          t_pl_w = theta_m(ijkw,l)
          t_pl_n = theta_m(ijkn,l)
          t_pl_s = theta_m(ijks,l)
          t_pl_t = theta_m(ijkt,l)
          t_pl_b = theta_m(ijkb,l)

          usl_e = u_s(ijk,l)
          usl_w = u_s(imjk,l)
          vsl_n = v_s(ijk,l)
          vsl_s = v_s(ijmk,l)
          wsl_t = w_s(ijk,l)
          wsl_b = w_s(ijkm,l)
          usl_p = avg_x_e(u_s(imjk,l),u_s(ijk,l),i)
          vsl_p = avg_y_n(v_s(ijmk,l),v_s(ijk,l) )
          wsl_p = avg_z_t(w_s(ijkm,l),w_s(ijk,l) )

! Energy dissipation by collisions: Sum(Nip)
! SUM( EDT_s_ip )
          s20_sum_lhs = s20_sum_lhs + edt_s_ip(ijk,m,l)

! Energy dissipation by collisions: SUM(Nip)
! SUM( EDvel_sL_ip* div(vp) ) !Modified by sof to include trace of V_s_L
          s11_sum_lhs = s11_sum_lhs + zmax(-edvel_sl_ip(ijk,m,l)* &
               trd_s_c(ijk,l) ) * vol(ijk)
          s11_sum_rhs = s11_sum_rhs + zmax( edvel_sl_ip(ijk,m,l)* &
               trd_s_c(ijk,l) ) * vol(ijk)

! Energy dissipation by collisions: Sum(Nip)
! SUM( EDvel_sM_ip* div(vi) ) !Modified by sof to include trace of V_s_M
          s12_sum_lhs = s12_sum_lhs + zmax(-edvel_sm_ip(ijk,m,l)*&
               trd_s_c(ijk,m) ) * vol(ijk)
          s12_sum_rhs = s12_sum_rhs + zmax( edvel_sm_ip(ijk,m,l)*&
               trd_s_c(ijk,m) ) * vol(ijk)

          IF (m .NE. l) THEN
               lm = funlm(l,m)

! Production by shear: (S:grad(vi))
! SUM(2*Mu_sL_ip*tr(Dk*Di) )
               s17_sum_lhs = s17_sum_lhs + 2.d0*mu_sl_ip(ijk,m,l)*&
                    zmax( - trd_s2_ip(ijk,m,l) )
               s17_sum_rhs = s17_sum_rhs + 2.d0*mu_sl_ip(ijk,m,l)*&
                    zmax( trd_s2_ip(ijk,m,l) )

! These two terms can be treated explicitly here by uncommenting the following
! two lines. They are currently treated by PEA algorithm in solve_granular_energy
!
!               S13_sum_lhs = S13_sum_lhs + ED_ss_ip(IJK,LM)*Theta_m(IJK,M)
!               S13_sum_rhs = S13_sum_rhs + ED_ss_ip(IJK,LM)*Theta_m(IJK,L)
!
!
! Production by shear: (S:grad(vi))
! SUM( (Xi_sL_ip-(2/3)*Mu_sL_ip)*tr(Dk)tr(Di) )
               s18_sum_lhs = s18_sum_lhs + zmax(-1.d0* &
                    (xi_sl_ip(ijk,m,l)-(2.d0/3.d0)*mu_sl_ip(ijk,m,l))*&
                    trd_s_c(ijk,m)*trd_s_c(ijk,l) )
               s18_sum_rhs = s18_sum_rhs + zmax(&
                    (xi_sl_ip(ijk,m,l)-(2.d0/3.d0)*mu_sl_ip(ijk,m,l))*&
                    trd_s_c(ijk,m)*trd_s_c(ijk,l) )

! Part of Heat Flux: div (q)
! Kth_sL_ip*[grad(Tp)]
! Note for L=M S21 terms cancel with similar term arising from grad(Ti)
               kth_sl_e = avg_x_s(kth_sl_ip(ijk,m,l), kth_sl_ip(ijke,m,l),i)
               kth_sl_w = avg_x_s(kth_sl_ip(ijkw,m,l),kth_sl_ip(ijk,m,l), im)
               kth_sl_n = avg_y_s(kth_sl_ip(ijk,m,l), kth_sl_ip(ijkn,m,l),j)
               kth_sl_s = avg_y_s(kth_sl_ip(ijks,m,l),kth_sl_ip(ijk,m,l), jm)
               kth_sl_t = avg_z_s(kth_sl_ip(ijk,m,l), kth_sl_ip(ijkt,m,l),k)
               kth_sl_b = avg_z_s(kth_sl_ip(ijkb,m,l),kth_sl_ip(ijk,m,l), km)

               s21a_sum = s21a_sum + ( (kth_sl_e*(t_pl_e-t_pl_p) )*&
                    odx_e(i)*ayz(ijk) - (kth_sl_w*(t_pl_p-t_pl_w) )*odx_e(im)*&
                    ayz(imjk) )
               s21b_sum = s21b_sum + ( (kth_sl_n*(t_pl_n-t_pl_p) )*&
                    ody_n(j)*axz(ijk) - (kth_sl_s*(t_pl_p-t_pl_s) )*ody_n(jm)*&
                    axz(ijmk) )
               s21c_sum = s21c_sum + ( (kth_sl_t*(t_pl_t-t_pl_p) )*&
                    odz_t(k)*ox(i)*axy(ijk) - (kth_sl_b*(t_pl_p-t_pl_b) )*&
                    odz_t(km)*ox(i)*axy(ijkm) )

! Part of Heat Flux: div (q)
! Knu_s_ip*[ni*grad(np)-np*grad(ni)]
! Note S14 terms should evaluate to zero for particles from the same phase
               knu_sl_e = avg_x_s(knu_sl_ip(ijk,m,l), knu_sl_ip(ijke,m,l),i)
               knu_sm_e = avg_x_s(knu_sm_ip(ijk,m,l), knu_sm_ip(ijke,m,l),i)
               knu_sl_w = avg_x_s(knu_sl_ip(ijkw,m,l),knu_sl_ip(ijk,m,l), im)
               knu_sm_w = avg_x_s(knu_sm_ip(ijkw,m,l),knu_sm_ip(ijk,m,l), im)
               knu_sl_n = avg_y_s(knu_sl_ip(ijk,m,l), knu_sl_ip(ijkn,m,l),j)
               knu_sm_n = avg_y_s(knu_sm_ip(ijk,m,l), knu_sm_ip(ijkn,m,l),j)
               knu_sl_s = avg_y_s(knu_sl_ip(ijks,m,l),knu_sl_ip(ijk,m,l), jm)
               knu_sm_s = avg_y_s(knu_sm_ip(ijks,m,l),knu_sm_ip(ijk,m,l), jm)
               knu_sl_t = avg_z_s(knu_sl_ip(ijk,m,l), knu_sl_ip(ijkt,m,l),k)
               knu_sm_t = avg_z_s(knu_sm_ip(ijk,m,l), knu_sm_ip(ijkt,m,l),k)
               knu_sl_b = avg_z_s(knu_sl_ip(ijkb,m,l),knu_sl_ip(ijk,m,l), km)
               knu_sm_b = avg_z_s(knu_sm_ip(ijkb,m,l),knu_sm_ip(ijk,m,l), km)

               s14a_sum = s14a_sum + ( (knu_sl_e*(nu_pl_e-nu_pl_p) - &
                    knu_sm_e*(nu_pm_e-nu_pm_p) )*odx_e(i)*ayz(ijk) - (knu_sl_w*&
                    (nu_pl_p-nu_pl_w) - knu_sm_w*(nu_pm_p-nu_pm_w) )*odx_e(im)*&
                    ayz(imjk) )
               s14b_sum = s14b_sum + ( (knu_sl_n*(nu_pl_n-nu_pl_p) - &
                    knu_sm_n*(nu_pm_n-nu_pm_p) )*ody_n(j)*axz(ijk) - (knu_sl_s*&
                    (nu_pl_p-nu_pl_s) - knu_sm_s*(nu_pm_p-nu_pm_s) )*ody_n(jm)*&
                    axz(ijmk) )
               s14c_sum = s14c_sum + ( (knu_sl_t*(nu_pl_t-nu_pl_p) - &
                    knu_sm_t*(nu_pm_t-nu_pm_p) )*odz_t(k)*ox(i)*axy(ijk) - &
                    (knu_sl_b*(nu_pl_p-nu_pl_b) - knu_sm_b*(nu_pm_p-nu_pm_b) )*&
                    odz_t(km)*ox(i)*axy(ijkm) )

! Part of Heat Flux: div (q)
! Kvel_s_ip*[vi-vp]
! Note S9 terms should evaluate to zero for particles from the same phase
               kvel_s_e = avg_x_h(kvel_s_ip(ijk,m,l), kvel_s_ip(ijke,m,l),i)
               kvel_s_w = avg_x_h(kvel_s_ip(ijkw,m,l),kvel_s_ip(ijk,m,l), im)
               kvel_s_n = avg_y_h(kvel_s_ip(ijk,m,l), kvel_s_ip(ijkn,m,l),j)
               kvel_s_s = avg_y_h(kvel_s_ip(ijks,m,l),kvel_s_ip(ijk,m,l), jm)
               kvel_s_t = avg_z_h(kvel_s_ip(ijk,m,l), kvel_s_ip(ijkt,m,l),k)
               kvel_s_b = avg_z_h(kvel_s_ip(ijkb,m,l),kvel_s_ip(ijk,m,l), km)

               s9_sum = s9_sum + ( kvel_s_e*(usm_e-usl_e)*ayz(ijk) - &
                    kvel_s_w*(usm_w-usl_w)*ayz(imjk) + kvel_s_n*(vsm_n-vsl_n)*axz(ijk)-&
                    kvel_s_s*(vsm_s-vsl_s)*axz(ijmk) + kvel_s_t*(wsm_t-wsl_t)*axy(ijk)-&
                    kvel_s_b*(wsm_b-wsl_b)*axy(ijkm) )

          ENDIF    ! (IF M.NE.L)

      ENDDO

!  WARNING: The terms due to granular temperature gradients S21 (a,b,c) have caused
!           some converegence issues, remove them from LHS and RHS for debugging (sof).

      sourcelhs = ( (s11_sum_lhs+s12_sum_lhs)+&
          (s10_lhs+s16_lhs+s17_sum_lhs+&
          s18_sum_lhs-s20_sum_lhs+s13_sum_lhs)*vol(ijk) + &
          zmax(s21a_sum+s21b_sum+s21c_sum)+ &
          zmax(s14a_sum+s14b_sum+s14c_sum)+ zmax(s9_sum) ) / &
          theta_m(ijk,m)

      sourcerhs = ( s10_rhs+s15_rhs+s16_rhs+s17_sum_rhs+s18_sum_rhs+&
          s13_sum_rhs) * vol(ijk) + s11_sum_rhs+s12_sum_rhs+ &
          zmax(- (s14a_sum+s14b_sum+s14c_sum) ) + zmax(-s9_sum) + &
          zmax(- (s21a_sum+s21b_sum+s21c_sum) )


      RETURN
      END SUBROUTINE source_granular_energy_ia


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_GRANULAR_ENERGY_GD                               C
!  Purpose: Calculate the source terms of the Mth pahse granular       C
!     energy equation                                                  C
!                                                                      C
!  Author: Janine E. Galvin                                            C
!                                                                      C
!  Literature/Document References:                                     C
!    Garzo, V., and Dufty, J., Homogeneous cooling state for a         C
!      granular mixture, Physical Review E, 1999, Vol 60 (5), 5706-    C
!      5713                                                            C
!    Garzo, Tenneti, Subramaniam, Hrenya, J. Fluid Mech., 2012, 712,   C
!      pp 129-404                                                      C
!                                                                      C
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC

      SUBROUTINE source_granular_energy_gd(SOURCELHS, &
                    sourcerhs, ijk, m)

!-----------------------------------------------
!  Modules
!-----------------------------------------------
      USE compar
      USE constant
      USE drag
      USE fldvar
      USE fun_avg
      USE functions
      USE geometry
      USE indices
      USE parallel
      USE param
      USE param1
      USE physprop
      USE rdf
      USE residual
      USE run
      USE toleranc
      USE trace
      USE turb
      USE visc_g
      USE visc_s
      use kintheory
      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Source terms to be kept on rhs (source vector)
      DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
! Source terms to be kept on lhs (center coefficient)
      DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
! Solid phase index
      INTEGER, INTENT(IN) :: M
! Indices
      INTEGER, INTENT(IN) :: IJK
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
                 IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
! number densities
      DOUBLE PRECISION :: NU_PM_E, NU_PM_W, NU_PM_N, NU_PM_S, NU_PM_T,&
                          NU_PM_B, NU_PM_p
! particle characteristics
      DOUBLE PRECISION :: M_PM, D_PM
! coefficients in heat flux term
      DOUBLE PRECISION :: Knu_sM_e, Knu_sM_w, Knu_sM_n, Knu_sM_s,&
                          Knu_sM_t, Knu_sM_b
! Source terms to be kept on rhs
      DOUBLE PRECISION :: S8_rhs, S9_rhs, S10_rhs, S11_rhs, S12_rhs,&
                          S13a_rhs, S13b_rhs, S14a_rhs, S14b_rhs, S15a_rhs,&
                          S15b_rhs
! Source terms to be kept on lhs
      DOUBLE PRECISION :: S8_lhs, S10_lhs, S11_lhs, S12_lhs, S13a_lhs,&
                          S13b_lhs, S14a_lhs, S14b_lhs, S15a_lhs, S15b_lhs
! Source terms from interstitial effects
      DOUBLE PRECISION :: VSLIP, Kslip, Tslip_rhs, Tslip_lhs, Tvis_lhs
      DOUBLE PRECISION :: ep_sM
! local var. for GTSH theory
      DOUBLE PRECISION :: chi, Sgama_lhs, Spsi_rhs
!-----------------------------------------------

      i = i_of(ijk)
      j = j_of(ijk)
      k = k_of(ijk)
      im = im1(i)
      jm = jm1(j)
      km = km1(k)
      imjk = im_of(ijk)
      ijmk = jm_of(ijk)
      ijkm = km_of(ijk)
      ijke = east_of(ijk)
      ijkw = west_of(ijk)
      ijkn = north_of(ijk)
      ijks = south_of(ijk)
      ijkt = top_of(ijk)
      ijkb = bottom_of(ijk)

      s8_lhs = zero
      s10_lhs = zero
      s11_lhs = zero
      s12_lhs = zero
      s13a_lhs = zero
      s13b_lhs = zero
      s14a_lhs = zero
      s14b_lhs = zero
      s15a_lhs = zero
      s15b_lhs = zero
      s8_rhs = zero
      s9_rhs = zero
      s10_rhs = zero
      s11_rhs = zero
      s12_rhs = zero
      s13a_rhs = zero
      s13b_rhs = zero
      s14a_rhs = zero
      s14b_rhs = zero
      s15a_rhs = zero
      s15b_rhs = zero
      sgama_lhs = zero  ! for GTSH theory
      spsi_rhs = zero   ! for GTSH theory

      d_pm = d_p(ijk,m)
      m_pm = (pi/6.d0)*d_pm**3 * ro_s(ijk,m)
      ep_sm = ep_s(ijk,m)
      chi = g_0(ijk,m,m)

      nu_pm_p = rop_s(ijk,m)/m_pm
      nu_pm_e = rop_s(ijke,m)/m_pm
      nu_pm_w = rop_s(ijkw,m)/m_pm
      nu_pm_n = rop_s(ijkn,m)/m_pm
      nu_pm_s = rop_s(ijks,m)/m_pm
      nu_pm_t = rop_s(ijkt,m)/m_pm
      nu_pm_b = rop_s(ijkb,m)/m_pm

! Production by shear: (S:grad(v))
! P_s*tr(D)
      s8_lhs = p_s_c(ijk,m) * zmax(trd_s_c(ijk,m))
      s8_rhs = p_s_c(ijk,m) * zmax(-trd_s_c(ijk,m))

! Production by shear: (S:grad(v))
! Mu_s*tr(D^2)
      s9_rhs = 2.d0*mu_s_c(ijk,m)*trd_s2(ijk,m)

! Production by shear: (S:grad(v))
! Lambda_s*tr(D)^2
      s10_lhs = (trd_s_c(ijk,m)**2)*zmax( -lambda_s_c(ijk,m) )
      s10_rhs = (trd_s_c(ijk,m)**2)*zmax(  lambda_s_c(ijk,m) )

! Energy dissipation by collisions: (3/2)*n*kboltz*T*zeta0
! linearized (3/2)*rop_s*T*zeta0
      s11_lhs = (3.d0/2.d0)*edt_s_ip(ijk,m,m)
      s11_rhs = (1.d0/2.d0)*edt_s_ip(ijk,m,m)*theta_m(ijk,m)

! Energy dissipation by collisions: (3/2)*n*kboltz*T*zeta1
! (3/2)*rop_s*T*zeta1
      s12_lhs = zmax( edvel_sm_ip(ijk,m,m) * trd_s_c(ijk,m) )
      s12_rhs = zmax( -edvel_sm_ip(ijk,m,m) * trd_s_c(ijk,m) )*theta_m(ijk,m)

! for GTSH theory, the dissipation terms above need to be multiplied
! by 3/2 rop_s(ijk,m)
! GTSH theory has two additional terms as source (Psi) and sink (gama)
! in eq (4.9) of GTSH JFM paper
      IF(kt_type_enum == gtsh_2012) THEN
        s11_lhs = 1.5d0*rop_s(ijk,m) * s11_lhs
        s11_rhs = 1.5d0*rop_s(ijk,m) * s11_rhs
        s12_lhs = 1.5d0*rop_s(ijk,m) * s12_lhs
        s12_rhs = 1.5d0*rop_s(ijk,m) * s12_rhs
        sgama_lhs = 3d0*nu_pm_p * g_gtsh(ep_sm, chi, ijk, m)
        spsi_rhs = 1.5d0*rop_s(ijk,m) * xsi_gtsh(ijk)
      ENDIF

! Part of Heat Flux: div (q)
! Knu_s_ip*grad(nu)
      knu_sm_e = avg_x_s(kphi_s(ijk,m), kphi_s(ijke,m),i)
      knu_sm_w = avg_x_s(kphi_s(ijkw,m),kphi_s(ijk,m), im)
      knu_sm_n = avg_y_s(kphi_s(ijk,m), kphi_s(ijkn,m),j)
      knu_sm_s = avg_y_s(kphi_s(ijks,m),kphi_s(ijk,m), jm)
      knu_sm_t = avg_z_s(kphi_s(ijk,m), kphi_s(ijkt,m),k)
      knu_sm_b = avg_z_s(kphi_s(ijkb,m),kphi_s(ijk,m), km)

      s13a_rhs = knu_sm_e*zmax( (nu_pm_e-nu_pm_p) )*odx_e(i)*ayz(ijk)
      s13a_lhs = knu_sm_e*zmax( -(nu_pm_e-nu_pm_p) )*odx_e(i)*ayz(ijk)

      s13b_rhs = knu_sm_w*zmax( -(nu_pm_p-nu_pm_w) )*odx_e(im)*ayz(imjk)
      s13b_lhs = knu_sm_w*zmax( (nu_pm_p-nu_pm_w) )*odx_e(im)*ayz(imjk)

      s14a_rhs = knu_sm_n*zmax( (nu_pm_n-nu_pm_p) )*ody_n(j)*axz(ijk)
      s14a_lhs = knu_sm_n*zmax( -(nu_pm_n-nu_pm_p) )*ody_n(j)*axz(ijk)

      s14b_rhs = knu_sm_s*zmax( -(nu_pm_p-nu_pm_s) )*ody_n(jm)*axz(ijmk)
      s14b_lhs = knu_sm_s*zmax( (nu_pm_p-nu_pm_s) )*ody_n(jm)*axz(ijmk)

      s15a_rhs = knu_sm_t*zmax( (nu_pm_t-nu_pm_p) )*odz_t(k)*ox(i)*axy(ijk)
      s15a_lhs = knu_sm_t*zmax( -(nu_pm_t-nu_pm_p) )*odz_t(k)*ox(i)*axy(ijk)

      s15b_rhs = knu_sm_b*zmax( -(nu_pm_p-nu_pm_b) )*odz_t(km)*ox(i)*axy(ijkm)
      s15b_lhs = knu_sm_b*zmax( (nu_pm_p-nu_pm_b) )*odz_t(km)*ox(i)*axy(ijkm)

      sourcelhs = ( (s8_lhs+s10_lhs)*vol(ijk) + &
          s13a_lhs+s13b_lhs+s14a_lhs+s14b_lhs+s15a_lhs+s15b_lhs)/theta_m(ijk,m)&
          + (s11_lhs + s12_lhs + sgama_lhs)*vol(ijk)


      sourcerhs = ( s8_rhs+s9_rhs+s10_rhs+s11_rhs+s12_rhs+spsi_rhs) * vol(ijk) + &
          s13a_rhs+s13b_rhs+s14a_rhs+s14b_rhs+s15a_rhs+s15b_rhs


! this is only done for GD_99, do not add these terms to GTSH
      IF(switch > zero .AND. ro_g0 /= zero .AND. &
         (kt_type_enum == gd_1999)) THEN

         vslip = (u_s(ijk,m)-u_g(ijk))**2 + (v_s(ijk,m)-v_g(ijk))**2 +&
            (w_s(ijk,m)-w_g(ijk))**2
         vslip = dsqrt(vslip)

! production by gas-particle slip: Koch & Sangani (1999)
         kslip = switch*81.d0*ep_sm*(mu_g(ijk)*vslip)**2.d0 / &
            (chi*d_p(ijk,m)**3.d0*ro_s(ijk,m)*dsqrt(pi))

         tslip_rhs = 1.5d0*kslip/( (theta_m(ijk,m)+small_number)**0.50)*vol(ijk)
         tslip_lhs = 0.5d0*kslip/( (theta_m(ijk,m)+small_number)**1.50)*vol(ijk)

! dissipation by viscous damping: Gidaspow (1994)
         tvis_lhs = switch*3d0*f_gs(ijk,m)*vol(ijk)

         sourcelhs = sourcelhs + tslip_lhs + tvis_lhs
         sourcerhs = sourcerhs + tslip_rhs
      ENDIF

      RETURN
      END SUBROUTINE source_granular_energy_gd