drag_ss.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DRAG_ss                                                 C
!  Purpose: This module computes the coefficient of drag between       C
!  two solids phases (M and L).                                        C
!                                                                      C
!  Author: M. Syamlal                                 Date: 20-MAY-96  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE drag_ss(L, M)

! Modules
!---------------------------------------------------------------------//
      use compar, only: mype
      use exit, only: mfix_exit
      use functions, only: funlm
      use run, only: granular_energy, kt_type, kt_type_enum
      use run, only: lun_1984, simonin_1996, ahmadi_1995
      use run, only: ia_2005, gd_1999, ghd_2007, gtsh_2012
      use usr_prop, only: usr_fss, calc_usr_prop
      use usr_prop, only: solidssolids_drag
      IMPLICIT NONE

! Dummy Arguments
!---------------------------------------------------------------------//
! Index of solids phases
      INTEGER, INTENT(IN) :: L, M

! Local Variables
!---------------------------------------------------------------------//
! index for storing solids-solids drag coefficients in the upper
! triangle of the matrix
      INTEGER :: LM

!---------------------------------------------------------------------//
      lm = funlm(l,m)

      IF (usr_fss(lm)) THEN
         CALL calc_usr_prop(solidssolids_drag,ll=l,lm=m)
      ELSE
         IF (granular_energy) THEN
            SELECT CASE (kt_type_enum)
               CASE(lun_1984, simonin_1996, ahmadi_1995)
                  CALL drag_ss_syam(l, m)
               CASE (ia_2005)
                  CALL drag_ss_ia(l, m)
               CASE (gd_1999, gtsh_2012)
! strictly speaking gd and gtsh are monodisperse theories and so
! do not have solids-solids drag
                  RETURN
! ghd theory is a polydisperse theory but is self contained and does
! not invoke this routine
               CASE (ghd_2007)
                  RETURN
            CASE DEFAULT
! should never hit this
               WRITE (*, '(A)') 'DRAG_SS'
               WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
            call mfix_exit(mype)
            END SELECT
         ELSE
            CALL drag_ss_syam(l, m)
         ENDIF
      ENDIF

      RETURN
      END SUBROUTINE drag_ss


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DRAG_SS_SYAM                                            C
!  Purpose: Calculate the solids-solids drag coefficient between a     C
!           continuous solids phase and discrete solids                C
!                                                                      C
!  Literature/Document References:                                     C
!     M. Syamlal. 1987. The particle-particle drag term in a           C
!        multiparticle model of fluidization. Technical Report.        C
!        DOE/MC/21353-2373. Office of Fossil Energy, Morgantown        C
!        Energy Technology Center, Morgantown, West Virginia.          C
!     Gera, D., Syamlal, M., O'Brien T.J. 2004. International Journal  C
!        of Multiphase Flow, 30, p419-428.                             C
!                                                                      C
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
      SUBROUTINE drag_ss_syam(L,M)

! Modules
!---------------------------------------------------------------------//
      use compar, only: mype, ijkstart3, ijkend3
      USE constant, only: segregation_slope_coefficient
      USE drag, only: f_ss
      USE fldvar, only: d_p, ro_s, rop_s, theta_m, p_star
      USE fldvar, only: u_s, v_s, w_s
      USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
      use functions, only: wall_at, funlm
      use functions, only: im_of, jm_of, km_of
      USE indices, only: i_of
      USE physprop, only: close_packed
      USE rdf, only: g_0
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids phase indices
      INTEGER, INTENT(IN) :: M, L

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: I, IJK, IMJK, IJMK, IJKM
! index for storing solids-solids drag coefficients in the upper
! triangle of the matrix
      INTEGER :: LM
! cell center value of U_sm, U_sl, V_sm, V_sl, W_sm, W_sl
      DOUBLE PRECISION :: USCM, USCL, VSCM, VSCL, WSCM, WSCL
! relative velocity between solids phase m and l
      DOUBLE PRECISION :: VREL
! particle diameters of phase M and phase L
      DOUBLE PRECISION :: D_pm, D_pl
! particle densities of phase M and phase L
      DOUBLE PRECISION :: RO_M, RO_L
! radial distribution function between phases M and L
      DOUBLE PRECISION :: G0_ML
! solid-solid drag coefficient
      DOUBLE PRECISION :: lDss

!---------------------------------------------------------------------//
      lm = funlm(l,m)

      DO ijk = ijkstart3, ijkend3
         IF (wall_at(ijk)) cycle
! Evaluate at all flow boundaries and fluid cellls. This is unlike
! calculation of the fluid-solid drag coefficient, which is only
! evaluated in fluid cells and pressure inflow cells

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

! calculating velocity components at i, j, k (cell center)
         uscl = avg_x_e(u_s(imjk,l),u_s(ijk,l),i)
         vscl = avg_y_n(v_s(ijmk,l),v_s(ijk,l))
         wscl = avg_z_t(w_s(ijkm,l),w_s(ijk,l))

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

! magnitude of solids-solids relative velocity
         vrel = sqrt((uscl - uscm)**2 + (vscl - vscm)**2 + &
                     (wscl - wscm)**2)

! setting aliases for easy reference
         d_pm = d_p(ijk,m)
         d_pl = d_p(ijk,l)
         ro_m = ro_s(ijk,m)
         ro_l = ro_s(ijk,l)
         g0_ml = g_0(ijk,l,m)

! evaluating the solids-solids drag coefficient
         CALL drag_ss_syam0(ldss, d_pm, d_pl, ro_m, ro_l, g0_ml, vrel)
         f_ss(ijk,lm) = ldss*rop_s(ijk,m)*rop_s(ijk,l)

! Gera: accounting for particle-particle drag due to enduring contact
! in a close-packed system.
         IF(close_packed(m) .AND. close_packed(l)) &
            f_ss(ijk,lm) = f_ss(ijk,lm) + &
               segregation_slope_coefficient*p_star(ijk)

      ENDDO    ! end do (ijk=ijkstart3,ijkend3)

      RETURN
      END SUBROUTINE drag_ss_syam


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DRAG_SS_SYAM0                                           C
!  Purpose: Created as a work-around for des/hybrid cases that need    C
!  to use this model while passing particle type information.          C
!                                                                      C
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
      SUBROUTINE drag_ss_syam0(ldss, d_pm, d_pl, ro_m, ro_l, g0_ml, vrel)

! Modules 
!---------------------------------------------------------------------//
      use param1, only: one
      use constant, only: c_e, c_f, pi
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solid-solid drag coefficient
      DOUBLE PRECISION, INTENT(OUT) :: lDss
! particle diameters of phase M and phase L
      DOUBLE PRECISION, INTENT(IN) :: D_pm, D_pl
! particle densities of phase M and phase L
      DOUBLE PRECISION, INTENT(IN) :: RO_M, RO_L
! radial distribution function between phases M and L
      DOUBLE PRECISION, INTENT(IN) :: G0_ML
! relative velocity between solids phase m and l
      DOUBLE PRECISION, INTENT(IN) :: VREL

! Local variables
!---------------------------------------------------------------------//
! Sum of particle diameters
      DOUBLE PRECISION :: DPSUM
! Intermediate calculation
      DOUBLE PRECISION :: const
!---------------------------------------------------------------------//

      dpsum = d_pl + d_pm
      const = 3.d0*(one + c_e)*(pi/2.d0 + c_f*pi*pi/8.d0)*&
              dpsum**2/(2.d0*pi*(ro_l*d_pl**3+ro_m*d_pm**3))
      ldss = const * g0_ml * vrel

      RETURN
      END SUBROUTINE drag_ss_syam0

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DRAG_SS_IA                                              !
!  Purpose: Compute collisional momentum source terms betweem solids   !
!  phase M and solids phase L using Iddir Arastoopour (2005) kinetic   !
!  theory mdoel that is proportional to the relative velocity between  !
!  between the two phases (i.e. solids-solids drag).                   !
!                                                                      !
!  Literature/Document References:                                     !
!    Iddir, Y.H., "Modeling of the multiphase mixture of particles     !
!       using the kinetic theory approach," PhD Thesis, Illinois       !
!       Institute of Technology, Chicago, Illinois, 2004               !
!    Iddir, Y.H., & H. Arastoopour, "Modeling of Multitype particle    !
!      flow using the kinetic theory approach," AIChE J., Vol 51,      !
!      no. 6, June 2005                                                !
!     Gera, D., Syamlal, M., O'Brien T.J. 2004. International Journal  !
!        of Multiphase Flow, 30, p419-428.                             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE drag_ss_ia(L, M)

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      USE constant, only: segregation_slope_coefficient
      USE constant, only: c_e, pi
      USE drag, only: f_ss
      USE fldvar, only: d_p, ro_s, rop_s, theta_m, p_star
      use functions, only: wall_at, funlm
      USE param1, only: zero
      USE physprop, only: close_packed
      USE rdf, only: g_0
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids phase indices
      INTEGER, INTENT(IN) :: M, L

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: IJK
! index for storing solids-solids drag coefficients in the upper
! triangle of the matrix
      INTEGER :: LM
! particle diameters of phase M and phase L
      DOUBLE PRECISION :: D_pm, D_pl
! particle densities of phase M and phase L
      DOUBLE PRECISION :: RO_M, RO_L
! radial distribution function between phases M and L
      DOUBLE PRECISION :: G0_ML
! granular temperature of phase M and phase L
      DOUBLE PRECISION :: thetaM, thetaL
! solid-solid drag coefficient
      DOUBLE PRECISION :: lDss
! Sum of particle diameters
      DOUBLE PRECISION :: DPSUM
! Intermediate calculation
      DOUBLE PRECISION :: M_PM, M_PL, MPSUM, DPSUMo2
      DOUBLE PRECISION :: Ap_lm, Dp_lm, Bp_lm, R2p_lm
      DOUBLE PRECISION :: F_common_term
!---------------------------------------------------------------------//
      lm = funlm(l,m)

      DO ijk = ijkstart3, ijkend3
         IF (wall_at(ijk)) cycle
! Evaluate at all flow boundaries and fluid cellls. This is unlike
! calculation of the fluid-solid drag coefficient, which is only
! evaluated in fluid cells and pressure inflow cells

! setting aliases for easy reference
         d_pm = d_p(ijk,m)
         d_pl = d_p(ijk,l)
         ro_m = ro_s(ijk,m)
         ro_l = ro_s(ijk,l)
         thetam = theta_m(ijk,m)
         thetal = theta_m(ijk,l)
         g0_ml = g_0(ijk,l,m)

         dpsum = d_pl + d_pm
         m_pm = (pi/6.d0) * d_pm**3 *ro_m
         m_pl = (pi/6.d0) * d_pl**3 *ro_l
         mpsum = m_pm + m_pl
         dpsumo2 = dpsum/2.d0
         ldss = zero

! evaluating the solids-solids drag coefficient
         IF(thetam > zero .AND. thetal > zero) THEN
            ap_lm = (m_pm*thetal+m_pl*thetam)/2.d0
            bp_lm = (m_pm*m_pl*(thetal-thetam ))/(2.d0*mpsum)
            dp_lm = (m_pl*m_pm*(m_pm*thetam+m_pl*thetal ))/&
                (2.d0*mpsum*mpsum)

            r2p_lm = ( 1.d0/( 2.d0*ap_lm**1.5 * dp_lm*dp_lm ) )+&
                     ( (3.d0*bp_lm*bp_lm)/( ap_lm**2.5 * dp_lm**3 ) )+&
                     ( (15.d0*bp_lm**4)/( 2.d0*ap_lm**3.5 * dp_lm**4 ) )

            f_common_term = (dpsumo2*dpsumo2/4.d0)*(m_pm*m_pl/mpsum)*&
                 g0_ml*(1.d0+c_e)*(m_pm*m_pl)**1.5

! Momentum source associated with relative velocity between solids
! phase m and solid phase l
! Factor of rop_sm and rop_sl included in calling routine
            ldss = f_common_term/(m_pm*m_pl)*dsqrt(pi)*r2p_lm*&
                        (thetam*thetal)**2
         ENDIF
         f_ss(ijk,lm) = ldss*rop_s(ijk,m)*rop_s(ijk,l)

! Gera: accounting for particle-particle drag due to enduring contact
! in a close-packed system.
         IF(close_packed(m) .AND. close_packed(l)) &
            f_ss(ijk,lm) = f_ss(ijk,lm) + &
               segregation_slope_coefficient*p_star(ijk)

      ENDDO    ! end do (ijk=ijkstart3,ijkend3)

      RETURN
      END SUBROUTINE drag_ss_ia