kintheory_v_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_IA_MOMSOURCE_V_S                                   C
!  Purpose: Determine source terms for V_S momentum equation arising   C
!           from IA kinetic theory constitutive relations for stress   C
!           and solid-solid drag                                       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_v_s(M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1, only: 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, K, IJK, IJKN, JP, IM,  KM, IJPK, IJMK,&
                 IJKE, IJKNE, IJKW, IJKNW, IMJPK, IMJK, IJKT,&
                 IJKTN, IJKB, IJKBN, IJKM, IJPKM,&
                 IPJK, IJKP
! Phase index
      INTEGER :: L
! Viscosity values
      DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
                          MU_sL_pB
! Bulk viscosity values
      DOUBLE PRECISION :: XI_sL_pN, XI_sL_pS, LAMBDA_sL_pN, LAMBDA_sL_pS
! Variables for drag calculations
      DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pN, NU_PM_pS, NU_PM_p, &
                          NU_PL_pN, NU_PL_pS, NU_PL_p, T_PM_pN, T_PM_pS, &
                          T_PL_pN, T_PL_pS, Fnu_s_p, FT_sM_p, FT_sL_p
! Average volume fraction
      DOUBLE PRECISION :: EPSA
! Source terms (Surface)
      DOUBLE PRECISION :: ssvx, ssvy, ssvz, ssx, ssy, ssz, ssbv
! Source terms (Volumetric)
      DOUBLE PRECISION :: DS1, DS2, DS3, DS4, DS1plusDS2
!-----------------------------------------------

! section largely based on tau_v_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)

          j = j_of(ijk)
          ijkn = north_of(ijk)
          epsa = avg_y(ep_s(ijk,m),ep_s(ijkn,m),j)
          IF ( .NOT.sip_at_n(ijk) .AND. epsa>dil_ep_s) THEN

               jp = jp1(j)
               i = i_of(ijk)
               im = im1(i)
               k = k_of(ijk)
               km = km1(k)
               ijpk = jp_of(ijk)
               ijmk = jm_of(ijk)
               imjk = im_of(ijk)
               ijkm = km_of(ijk)
               imjpk = im_of(ijpk)
               ijpkm = jp_of(ijkm)

               ijkw = west_of(ijk)
               ijke = east_of(ijk)
               ijkne = east_of(ijkn)
               ijknw = north_of(ijkw)

               ijkb = bottom_of(ijk)
               ijkt = top_of(ijk)
               ijktn = north_of(ijkt)
               ijkbn = north_of(ijkb)

! additional required quantities:
               ipjk = ip_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 = avg_y_h(avg_x_h(mu_sl_ip(ijk,m,l),mu_sl_ip(ijke,m,l),i),&
                              avg_x_h(mu_sl_ip(ijkn,m,l),mu_sl_ip(ijkne,m,l),i),j)
                         mu_sl_pw = avg_y_h(avg_x_h(mu_sl_ip(ijkw,m,l),mu_sl_ip(ijk,m,l),im),&
                              avg_x_h(mu_sl_ip(ijknw,m,l),mu_sl_ip(ijkn,m,l),im),j)
                         ssvx = mu_sl_pe*(v_s(ipjk,l)-v_s(ijk,l))*ayz_v(ijk)*odx_e(i)&
                              -mu_sl_pw*(v_s(ijk,l)-v_s(imjk,l))*ayz_v(imjk)*odx_e(im)

                         mu_sl_pn = mu_sl_ip(ijkn,m,l)
                         mu_sl_ps = mu_sl_ip(ijk,m,l)
                         ssvy = mu_sl_pn*(v_s(ijpk,l)-v_s(ijk,l))*ody(jp)*axz_v(ijk)&
                              -mu_sl_ps*(v_s(ijk,l)-v_s(ijmk,l))*ody(j)*axz_v(ijmk)

                         mu_sl_pt = avg_y_h(avg_z_h(mu_sl_ip(ijk,m,l),mu_sl_ip(ijkt,m,l),k),&
                              avg_z_h(mu_sl_ip(ijkn,m,l),mu_sl_ip(ijktn,m,l),k),j)
                         mu_sl_pb = avg_y_h(avg_z_h(mu_sl_ip(ijkb,m,l),mu_sl_ip(ijk,m,l),km),&
                              avg_z_h(mu_sl_ip(ijkbn,m,l),mu_sl_ip(ijkn,m,l),km),j)
                         ssvz = mu_sl_pt*(v_s(ijkp,l)-v_s(ijk,l))*axy_v(ijk)*odz_t(k)*ox(i)&
                              -mu_sl_pb*(v_s(ijk,l)-v_s(ijkm,l))*axy_v(ijkm)*odz_t(km)*ox(i)

! bulk viscosity term
                         xi_sl_pn = xi_sl_ip(ijkn,m,l)
                         xi_sl_ps = xi_sl_ip(ijk,m,l)
                         lambda_sl_pn = -(2.d0/3.d0)*mu_sl_pn + xi_sl_pn
                         lambda_sl_ps = -(2.d0/3.d0)*mu_sl_ps + xi_sl_ps
                         ssbv = (lambda_sl_pn*trd_s(ijkn,l)-lambda_sl_ps*trd_s(ijk,l))*axz(ijk)

! off diagonal shear stress terms
                         ssx = mu_sl_pe*(u_s(ijpk,l)-u_s(ijk,l))*ody_n(j)*ayz_v(ijk)&
                              -mu_sl_pw*(u_s(imjpk,l)-u_s(imjk,l))*ody_n(j)*ayz_v(imjk)
                         ssy = ssvy
                         ssz = mu_sl_pt*(w_s(ijpk,l)-w_s(ijk,l))*axy_v(ijk)*ody_n(j)&
                              -mu_sl_pb*(w_s(ijpkm,l)-w_s(ijkm,l))*axy_v(ijkm)*ody_n(j)
!--------------------- 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/6d0)* d_pl**3 *ro_s(ijk,l)

                         nu_pm_pn = rop_s(ijkn,m)/m_pm
                         nu_pm_ps = rop_s(ijk,m)/m_pm
                         nu_pm_p = avg_y(nu_pm_ps,nu_pm_pn,j)

                         nu_pl_pn = rop_s(ijkn,l)/m_pl
                         nu_pl_ps = rop_s(ijk,l)/m_pl
                         nu_pl_p = avg_y(nu_pl_ps,nu_pl_pn,j)

                         fnu_s_p = avg_y(fnu_s_ip(ijk,m,l),fnu_s_ip(ijkn,m,l),j)
                         ds1 = fnu_s_p*nu_pl_p*(nu_pm_pn-nu_pm_ps)*ody_n(j)
                         ds2 = -fnu_s_p*nu_pm_p*(nu_pl_pn-nu_pl_ps)*ody_n(j)
                         ds1plusds2 = ds1 + ds2

! Momentum source associated with the gradient in granular
! temperature of species M
                         t_pm_pn = theta_m(ijkn,m)
                         t_pm_ps = theta_m(ijk,m)

                         ft_sm_p = avg_y(ft_sm_ip(ijk,m,l),ft_sm_ip(ijkn,m,l),j)
                         ds3 = ft_sm_p*(t_pm_pn-t_pm_ps)*ody_n(j)

! Momentum source associated with the gradient in granular
! temperature of species L
                         t_pl_pn = theta_m(ijkn,l)
                         t_pl_ps = theta_m(ijk,l)

                         ft_sl_p = avg_y(ft_sl_ip(ijk,m,l),ft_sl_ip(ijkn,m,l),j)
                         ds4 = ft_sl_p*(t_pl_pn-t_pl_ps)*ody_n(j)
!--------------------- End Sources from Momentum Source Term ---------------------


! Add the terms
                         stress_terms = stress_terms + ssvx + ssvy + ssvz + &
                             ssbv + ssx + ssy + ssz
                         drag_terms = drag_terms + (ds1plusds2+ds3+ds4)*vol_v(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_v_s(ijk,m) = stress_terms + drag_terms
          ELSE
               ktmom_v_s(ijk,m) = zero
          ENDIF     ! dilute
      ENDDO        ! over ijk

      RETURN
      END SUBROUTINE calc_ia_momsource_v_s