kintheory_u_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_EXPLICIT_MOM_SOURCE_S                              C
!  Purpose: Determine any additional momentum source terms that are    C
!  calculated explicitly here                                          C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE calc_explicit_mom_source_s()

!-----------------------------------------------
! Modules
!-----------------------------------------------

      USE param1, only: zero

! kinetic theories
      USE run, only: kt_type_enum
      USE run, only: ia_2005

! number of solids phases
      USE physprop, only: smax

! solids source term
      USE kintheory, only: ktmom_u_s, ktmom_v_s, ktmom_w_s

      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Solids phase index
      INTEGER :: M
!-----------------------------------------------
! Additional interphase interaction terms that arise from kinetic theory
      DO m = 1, smax
         ktmom_u_s(:,m) = zero
         ktmom_v_s(:,m) = zero
         ktmom_w_s(:,m) = zero

         IF (kt_type_enum == ia_2005) THEN
            CALL calc_ia_momsource_u_s(m)
            CALL calc_ia_momsource_v_s(m)
            CALL calc_ia_momsource_w_s(m)
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_explicit_mom_source_s


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_IA_MOMSOURCE_U_S                                   C
!  Purpose: Determine source terms for U_S momentum equation arising   C
!           from IA kinetic theory constitutive relations for stress   C
!           and solids-solids drag (collisional momentum source)       C
!                                                                      C
!  Literature/Document References:                                     C
!    Idir, 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 calc_ia_momsource_u_s(M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1, only: half, zero
      USE constant, only: pi

! trace of D_s at i, j, k
      USE visc_s, only: trd_s

! number of solids phases
      USE physprop, only: mmax

! x,y,z-components of solids velocity
      USE fldvar, only: u_s, v_s, w_s
! particle diameter, bulk density, material density
      USE fldvar, only: d_p, rop_s, ro_s
! granular temperature
      USE fldvar, only: theta_m, ep_s
! dilute threshold
      USE toleranc, only: dil_ep_s

      Use kintheory

      USE geometry
      USE indices
      USE compar
      USE fun_avg
      USE functions

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! solids phase index
      INTEGER, INTENT(IN) :: M
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Temporary variable
      DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
! Indices
      INTEGER :: I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
                 IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
                 IJKB, IJKBE, IJKM, IPJKM, &
                 IJPK, IJKP, IM, IJKW
! solids phase index
      INTEGER :: L
! Viscosity values for stress calculations
      DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
                          MU_sL_pB, Mu_sL_p
! Bulk viscosity values for calculations
      DOUBLE PRECISION :: XI_sL_pE, XI_sL_pW, LAMBDA_sL_pE, LAMBDA_sL_pW
! Variables for drag calculations
      DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pE, NU_PM_pW, NU_PM_p, &
                          NU_PL_pE, NU_PL_pW, NU_PL_p, T_PM_pE, T_PM_pW, &
                          T_PL_pE, T_PL_pW, Fnu_s_p, FT_sM_p, FT_sL_p
! Average volume fraction
      DOUBLE PRECISION :: EPSA
! Source terms (Surface)
      DOUBLE PRECISION :: ssux, ssuy, ssuz, ssx, ssy, ssz, ssbv
! Source terms (Volumetric)
      DOUBLE PRECISION :: tauzz, DS1, DS2, DS3, DS4, DS1plusDS2
!-----------------------------------------------

! section largely based on tau_u_g:

      DO ijk = ijkstart3, ijkend3

! Skip walls where some values are undefined.
         IF(wall_at(ijk)) cycle

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

         i = i_of(ijk)
         ijke = east_of(ijk)
         epsa = avg_x(ep_s(ijk,m),ep_s(ijke,m),i)
         IF ( .NOT.sip_at_e(ijk) .AND. epsa>dil_ep_s) THEN

            ip = ip1(i)
            im = im1(i)
            ijkw  = west_of(ijk)
            j = j_of(ijk)
            jm = jm1(j)
            k = k_of(ijk)
            km = km1(k)
            ipjk = ip_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)
            ipjmk = jm_of(ipjk)
            ipjkm = ip_of(ijkm)

            ijkn = north_of(ijk)
            ijks = south_of(ijk)
            ijkne = east_of(ijkn)
            ijkse = east_of(ijks)

            ijkt = top_of(ijk)
            ijkb = bottom_of(ijk)
            ijkte = east_of(ijkt)
            ijkbe = east_of(ijkb)

! additional required quantities:
            ijpk = jp_of(ijk)
            ijkp = kp_of(ijk)

! initialize variable
            stress_terms = zero
            drag_terms = zero

            DO l = 1,mmax
               IF (m .ne. l) THEN

!--------------------- Sources from Stress Terms ---------------------
! Surface Forces
! standard shear stress terms (i.e. ~diffusion)
                  mu_sl_pe = mu_sl_ip(ijke,m,l)
                  mu_sl_pw = mu_sl_ip(ijk,m,l)
                  ssux = mu_sl_pe*(u_s(ipjk,l)-u_s(ijk,l))*ayz_u(ijk)*odx(ip)&
                       -mu_sl_pw*(u_s(ijk,l)-u_s(imjk,l))*ayz_u(imjk)*odx(i)

                  mu_sl_pn = avg_x_h(avg_y_h(mu_sl_ip(ijk,m,l),mu_sl_ip(ijkn,m,l), j),&
                       avg_y_h(mu_sl_ip(ijke,m,l),mu_sl_ip(ijkne,m,l), j), i)
                  mu_sl_ps = avg_x_h(avg_y_h(mu_sl_ip(ijks,m,l),mu_sl_ip(ijk,m,l), jm),&
                       avg_y_h(mu_sl_ip(ijkse,m,l),mu_sl_ip(ijke,m,l), jm), i)
                  ssuy = mu_sl_pn*(u_s(ijpk,l)-u_s(ijk,l))*axz_u(ijk)*ody_n(j)&
                       -mu_sl_ps*(u_s(ijk,l)-u_s(ijmk,l))*axz_u(ijmk)*ody_n(jm)

                  mu_sl_pt = avg_x_h(avg_z_h(mu_sl_ip(ijk,m,l),mu_sl_ip(ijkt,m,l),k),&
                       avg_z_h(mu_sl_ip(ijke,m,l),mu_sl_ip(ijkte,m,l),k),i)
                  mu_sl_pb = avg_x_h(avg_z_h(mu_sl_ip(ijkb,m,l),mu_sl_ip(ijk,m,l),km),&
                       avg_z_h(mu_sl_ip(ijkbe,m,l),mu_sl_ip(ijke,m,l),km),i)
                  ssuz = mu_sl_pt*(u_s(ijkp,l)-u_s(ijk,l))*axy_u(ijk)*odz_t(k)*ox_e(i)&
                       -mu_sl_pb*(u_s(ijk,l)-u_s(ijkm,l))*axy_u(ijkm)*odz_t(km)*ox_e(i)

! bulk viscosity term
                  xi_sl_pe = xi_sl_ip(ijke,m,l)
                  xi_sl_pw = xi_sl_ip(ijk,m,l)
                  lambda_sl_pe = -(2.d0/3.d0)*mu_sl_pe + xi_sl_pe
                  lambda_sl_pw = -(2.d0/3.d0)*mu_sl_pw + xi_sl_pw
                  ssbv = (lambda_sl_pe*trd_s(ijke,l)-lambda_sl_pw*trd_s(ijk,l))*ayz(ijk)

! off diagonal shear stress terms
                  ssx = ssux
                  ssy = mu_sl_pn*(v_s(ipjk,l)-v_s(ijk,l))*axz_u(ijk)*odx_e(i)&
                       -mu_sl_ps*(v_s(ipjmk,l)-v_s(ijmk,l))*axz_u(ijmk)*odx_e(i)
                  ssz = mu_sl_pt*(w_s(ipjk,l)-w_s(ijk,l))*axy_u(ijk)*odx_e(i)&
                       -mu_sl_pb*(w_s(ipjkm,l)-w_s(ijkm,l))*axy_u(ijkm)*odx_e(i)

! special terms for cylindrical coordinates
                  IF (cylindrical) THEN

! modify Ssz: (1/x) (d/dz) (tau_zz)
!                    integral of (1/x) (d/dz) (mu*(-w/x))
!                    (normally part of tau_u_s) - explicit
                     ssz = ssz - (mu_sl_pt*(w_s(ipjk,l)+w_s(ijk,l))*half*ox_e(i)*axy_u(ijk)&
                          -mu_sl_pb*(w_s(ipjkm,l)+w_s(ijkm,l))*half*ox_e(i)*axy_u(ijkm))

! tau_zz/x terms: (tau_zz/x)
!                    integral of -(2mu/x)*(1/x)*dw/dz
!                    (normally part of tau_u_s) - explicit
                     mu_sl_p = avg_x(mu_sl_ip(ijk,m,l),mu_sl_ip(ijke,m,l),i)
                     tauzz = -2.d0*mu_sl_p*ox_e(i)*half*(&
                          ((w_s(ipjk,l)-w_s(ipjkm,l))*ox(ip)*odz(k))+&
                          ((w_s(ijk,l)-w_s(ijkm,l))*ox(i)*odz(k)) &
                          ) * vol_u(ijk)

!                    integral of -(2mu/x)*(1/x)*u
!                    (normally part of source_u_s)
                     tauzz = tauzz + (-2.d0*mu_sl_p*ox_e(i)*ox_e(i)*&
                          u_s(ijk,l)*vol_u(ijk))
                  ELSE
                     tauzz = zero
                  ENDIF
!--------------------- End Sources from Stress Term ---------------------


!--------------------- Sources from Momentum Source Term ---------------------
! Momentum source associated with the difference in the gradients in
! number density of solids phase m and all other solids phases
                  d_pl = d_p(ijk,l)
                  m_pl = (pi/6.d0)*d_pl**3 * ro_s(ijk,l)

                  nu_pm_pe = rop_s(ijke,m)/m_pm
                  nu_pm_pw = rop_s(ijk,m)/m_pm
                  nu_pm_p = avg_x(nu_pm_pw,nu_pm_pe,i)

                  nu_pl_pe = rop_s(ijke,l)/m_pl
                  nu_pl_pw = rop_s(ijk,l)/m_pl
                  nu_pl_p = avg_x(nu_pl_pw,nu_pl_pe,i)

                  fnu_s_p = avg_x(fnu_s_ip(ijk,m,l),fnu_s_ip(ijke,m,l),i)
                  ds1 = fnu_s_p*nu_pl_p*(nu_pm_pe-nu_pm_pw)*odx_e(i)
                  ds2 = -fnu_s_p*nu_pm_p*(nu_pl_pe-nu_pl_pw)*odx_e(i)
                  ds1plusds2 = ds1 + ds2

! Momentum source associated with the gradient in granular
! temperature of species M
                  t_pm_pe = theta_m(ijke,m)
                  t_pm_pw = theta_m(ijk,m)

                  ft_sm_p = avg_x(ft_sm_ip(ijk,m,l),ft_sm_ip(ijke,m,l),i)
                  ds3 = ft_sm_p*(t_pm_pe-t_pm_pw)*odx_e(i)

! Momentum source associated with the gradient in granular
! temperature of species L
                  t_pl_pe = theta_m(ijke,l)
                  t_pl_pw = theta_m(ijk,l)

                  ft_sl_p = avg_x(ft_sl_ip(ijk,m,l),ft_sl_ip(ijke,m,l),i)
                  ds4 = ft_sl_p*(t_pl_pe-t_pl_pw)*odx_e(i)
!--------------------- End Sources from Momentum Source Term ---------------------


! Add the terms
                  stress_terms = stress_terms + ssux + ssuy + ssuz + &
                       ssbv + ssx + ssy + ssz + tauzz
                  drag_terms = drag_terms + (ds1plusds2+ds3+ds4)*vol_u(ijk)

               ELSE ! if m .ne. L
! for m = l all stress terms should already be handled in existing routines
! for m = l all drag terms should become zero
                  stress_terms = stress_terms + zero
                  drag_terms = drag_terms + zero

               ENDIF ! if m .ne. L
            ENDDO     ! over L

            ktmom_u_s(ijk,m) = stress_terms + drag_terms
         ELSE
            ktmom_u_s(ijk,m) = zero
         ENDIF     ! dilute
      ENDDO        ! over ijk

      RETURN
      END SUBROUTINE calc_ia_momsource_u_s


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: COLL_MOMENTUM_COEFF_IA                                  !
!  Purpose: Compute collisional momentum source terms betweem solids   !
!  phase M and solids phase L using Iddir Arastoopour (2005) kinetic   !
!  theory model that are not proportional to the relative velocity     !
!  between the two phases.  Specifically, terms proportional to the    !
!  gradient in number density and gradient in temperature              !
!                                                                      !
!  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                                                !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

      SUBROUTINE coll_momentum_coeff_ia(L, M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE constant
      USE drag
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE kintheory
      USE param1
      USE physprop
      USE rdf
      USE sendrecv

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Solids phase index
      INTEGER, INTENT(IN) :: L
      INTEGER, INTENT(IN) :: M
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK
! Particle diameters
      DOUBLE PRECISION :: D_PM, D_PL
! Sum of particle diameters
      DOUBLE PRECISION :: DPSUM
!
      DOUBLE PRECISION :: M_PM, M_PL, MPSUM, DPSUMo2, NU_PL, NU_PM
      DOUBLE PRECISION :: Ap_lm, Dp_lm, Bp_lm
      DOUBLE PRECISION :: R0p_lm, R3p_lm, R4p_lm, R10p_lm
      DOUBLE PRECISION :: Fnus_ip, FTsM_ip, FTsL_ip, F_common_term
!-----------------------------------------------

      DO ijk = ijkstart3, ijkend3
         IF (.NOT.wall_at(ijk)) THEN

            IF (m == l) THEN
               fnu_s_ip(ijk,m,l) = zero
               ft_sm_ip(ijk,m,l) = zero
               ft_sl_ip(ijk,m,l) = zero

            ELSE
               d_pm = d_p(ijk,m)
               d_pl = d_p(ijk,l)
               dpsum = d_pl + d_pm
               m_pm = (pi/6.d0) * d_pm**3 *ro_s(ijk,m)
               m_pl = (pi/6.d0) * d_pl**3 *ro_s(ijk,l)
               mpsum = m_pm + m_pl
               dpsumo2 = dpsum/2.d0
               nu_pm = rop_s(ijk,m)/m_pm
               nu_pl = rop_s(ijk,l)/m_pl

               IF(theta_m(ijk,m) > zero .AND. theta_m(ijk,l) > zero) THEN

                  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)

                  r0p_lm = ( 1.d0/( ap_lm**1.5 * dp_lm**2.5 ) )+ &
                            ( (15.d0*bp_lm*bp_lm)/( 2.d0* ap_lm**2.5 * dp_lm**3.5 ) )+&
                            ( (175.d0*(bp_lm**4))/( 8.d0*ap_lm**3.5 * dp_lm**4.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 ) )

                  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 ) )

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

! Momentum source associated with the difference in the gradients in
! number density of solids phase m and all other solids phases
                  fnus_ip = f_common_term*(pi*dpsumo2/12.d0)*r0p_lm*&
                         (theta_m(ijk,m)*theta_m(ijk,l))**2.5

! Momentum source associated with the gradient in granular temperature
! of solid phase M
                  ftsm_ip = f_common_term*nu_pm*nu_pl*dpsumo2*pi*&
                            (theta_m(ijk,m)**1.5 * theta_m(ijk,l)**2.5) *&
                            (  (-1.5d0/12.d0*r0p_lm)+&
                            theta_m(ijk,l)/16.d0*(  (-m_pm*r10p_lm) - &
                            ((5.d0*m_pl*m_pl*m_pm/(192.d0*mpsum*mpsum))*r3p_lm)+&
                            ((5.d0*m_pm*m_pl)/(96.d0*mpsum)*r4p_lm*bp_lm)  )  )

! Momentum source associated with the gradient in granular temperature
! of solid phase L ! no need to recompute (sof Aug 30 2006)
                  ftsl_ip = f_common_term*nu_pm*nu_pl*dpsumo2*pi*&
                           (theta_m(ijk,l)**1.5 * theta_m(ijk,m)**2.5) *&
                            (  (1.5d0/12.d0*r0p_lm)+&
                            theta_m(ijk,m)/16.d0*(  (m_pl*r10p_lm)+&
                            (5.d0*m_pm*m_pm*m_pl/(192.d0*mpsum*mpsum)*r3p_lm)+&
                            (5.d0*m_pm*m_pl/(96.d0*mpsum) *r4p_lm*bp_lm)  )  )


                  fnu_s_ip(ijk,m,l) = fnus_ip

! WARNING: the following two terms have caused some convergence problems
! earlier. Set them to ZERO for debugging in case of converegence
! issues. (sof)
                  ft_sm_ip(ijk,m,l) = ftsm_ip  ! ZERO
                  ft_sl_ip(ijk,m,l) = ftsl_ip  ! ZERO
               ELSE
                  fnu_s_ip(ijk,m,l) = zero
                  ft_sm_ip(ijk,m,l) = zero
                  ft_sl_ip(ijk,m,l) = zero
               ENDIF
            ENDIF
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE coll_momentum_coeff_ia