kintheory_w_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_IA_MOMSOURCE_W_S                                   C
!  Purpose: Determine source terms for W_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_w_s(M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1, only: half, zero, undefined
      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 :: IJK, J, I, IM, IJKP, IMJK, IJKN, IJKNT, IJKS,&
                 IJKST, IJMKP, IJMK, IJKE, IJKTE, IJKW, IJKTW,&
                 IMJKP, K, IJKT, JM, KP, IJKM, IPJK,&
                 IJPK
! 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, MU_sL_p
! Bulk viscosity values
      DOUBLE PRECISION :: XI_sL_pT, XI_sL_pB, LAMBDA_sL_pT, LAMBDA_sL_pB
! Variables for drag calculations
      DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pT, NU_PM_pB, NU_PM_p, &
                          NU_PL_pT, NU_PL_pB, NU_PL_p, T_PM_pT, T_PM_pB,&
                          T_PL_pT, T_PL_pB, Fnu_s_p, FT_sM_p, FT_sL_p
! Average velocity gradients
      DOUBLE PRECISION :: duodz
! Average volume fraction
      DOUBLE PRECISION :: EPSA
! Source terms (Surface)
      DOUBLE PRECISION :: sswx, sswy, sswz, ssx, ssy, ssz, ssbv, tauxz_x
! Source terms (Volumetric)
      DOUBLE PRECISION :: tauxz_ox, DS1, DS2, DS3, DS4, DS1plusDS2
!-----------------------------------------------

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

          k = k_of(ijk)
          ijkt = top_of(ijk)
          epsa = avg_z(ep_s(ijk,m),ep_s(ijkt,m),k)
          IF ( .NOT.sip_at_t(ijk) .AND. epsa>dil_ep_s) THEN

               j = j_of(ijk)
               i = i_of(ijk)
               kp = kp1(k)
               im = im1(i)
               jm = jm1(j)
               imjk = im_of(ijk)
               ijmk = jm_of(ijk)
               ijkm = km_of(ijk)
               ijkp = kp_of(ijk)
               ijmkp = jm_of(ijkp)
               imjkp = kp_of(imjk)

               ijkn = north_of(ijk)
               ijks = south_of(ijk)
               ijknt = top_of(ijkn)
               ijkst = top_of(ijks)

               ijke = east_of(ijk)
               ijkw = west_of(ijk)
               ijkte = east_of(ijkt)
               ijktw = west_of(ijkt)

! additional required quantities:
               ipjk = ip_of(ijk)
               ijpk = jp_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_z_h(avg_x_h(mu_sl_ip(ijk,m,l),mu_sl_ip(ijke,m,l),i),&
                              avg_x_h(mu_sl_ip(ijkt,m,l),mu_sl_ip(ijkte,m,l),i),k)
                         mu_sl_pw = avg_z_h(avg_x_h(mu_sl_ip(ijkw,m,l),mu_sl_ip(ijk,m,l),im),&
                              avg_x_h(mu_sl_ip(ijktw,m,l),mu_sl_ip(ijkt,m,l),im),k)
                         sswx = mu_sl_pe*(w_s(ipjk,l)-w_s(ijk,l))*ayz_w(ijk)*odx_e(i)&
                              -mu_sl_pw*(w_s(ijk,l)-w_s(imjk,l))*ayz_w(imjk)*odx_e(im)

                         mu_sl_pn = avg_z_h(avg_y_h(mu_sl_ip(ijk,m,l),mu_sl_ip(ijkn,m,l), j),&
                              avg_y_h(mu_sl_ip(ijkt,m,l),mu_sl_ip(ijknt,m,l), j), k)
                         mu_sl_ps = avg_z_h(avg_y_h(mu_sl_ip(ijks,m,l),mu_sl_ip(ijk,m,l), jm),&
                              avg_y_h(mu_sl_ip(ijkst,m,l),mu_sl_ip(ijkt,m,l), jm), k)
                         sswy = mu_sl_pn*(w_s(ijpk,l)-w_s(ijk,l))*axz_w(ijk)*ody_n(j)&
                              -mu_sl_ps*(w_s(ijk,l)-w_s(ijmk,l))*axz_w(ijkm)*ody_n(jm)

                         mu_sl_pt = mu_sl_ip(ijkt,m,l)
                         mu_sl_pb = mu_sl_ip(ijk,m,l)
                         sswz = mu_sl_pt*(w_s(ijkp,l)-w_s(ijk,l))*axy_w(ijk)*odz(kp)*ox(i)&
                              -mu_sl_pb*(w_s(ijk,l)-w_s(ijkm,l))*axy_w(ijkm)*odz(k)*ox(i)

! bulk viscosity term
                         xi_sl_pt = xi_sl_ip(ijkt,m,l)
                         xi_sl_pb = xi_sl_ip(ijk,m,l)
                         lambda_sl_pt = -(2d0/3d0)*mu_sl_pt + xi_sl_pt
                         lambda_sl_pb = -(2d0/3d0)*mu_sl_pb + xi_sl_pb
                         ssbv = (lambda_sl_pt*trd_s(ijkt,l)-lambda_sl_pb*trd_s(ijk,l))*axy(ijk)

! off diagonal shear stress terms
                         ssx = mu_sl_pe*(u_s(ijkp,l)-u_s(ijk,l))*ox_e(i)*ayz_w(ijk)*odz_t(k)&
                              -mu_sl_pw*(u_s(imjkp,l)-u_s(imjk,l))*(dy(j)*half*(dz(k)+&
                              dz(kp)))*odz_t(k)
                         !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity

                         ssy = mu_sl_pn*(v_s(ijkp,l)-v_s(ijk,l))*axz_w(ijk)*odz_t(k)*ox(i)&
                              -mu_sl_ps*(v_s(ijmkp,l)-v_s(ijmk,l))*axz_w(ijmk)*odz_t(k)*ox(i)

                         ssz = sswz

! special terms for cylindrical coordinates
                         IF (cylindrical) THEN
! tau_zz term: modify Ssz: (1/x) (d/dz) (tau_zz)
!                             integral of (1/x) (d/dz) (2mu*(u/x))
!                             (normally part of tau_w_s) - explicit
                              ssz = ssz +(&
                                   (mu_sl_pt*(u_s(ijkp,l)+u_s(imjkp,l))*ox(i)*axy_w(ijk))-&
                                   (mu_sl_pb*(u_s(ijk,l)+u_s(imjk,l))*ox(i)*axy_w(ijkm)))

                              mu_sl_p = avg_z(mu_sl_ip(ijk,m,l),mu_sl_ip(ijkt,m,l),k)
! tau_xz/x terms: (tau_xz/x)
!                             integral of (1/x)*(mu/x)*du/dz
!                             (normally part of tau_w_s) - explicit
                              IF (ox_e(im) /= undefined) THEN
                                   duodz = (u_s(imjkp,l)-u_s(imjk,l))*ox_e(im)*odz_t(k)
                              ELSE
                                   duodz = zero
                              ENDIF

                              tauxz_ox = mu_sl_p*ox(i)*half*(&
                                   (ox_e(i)*(u_s(ijkp,l)-u_s(ijk,l))*odz_t(k))+&
                                   duodz)

!                             integral of (1/x)*mu*dw/dx
!                             (normally part of source_w_s)
                              tauxz_ox = tauxz_ox + (mu_sl_p*ox(i)*half*&
                                   ( (w_s(ipjk,l)-w_s(ijk,l))*odx_e(i) + &
                                   (w_s(ijk,l)-w_s(imjk,l))*odx_e(im) ) )

!                             integral of (1/x)*(-mu/x)*w
!                             (normally part of source_w_s)
                              tauxz_ox = tauxz_ox - (ox(i)*ox(i)*mu_sl_p*w_s(ijk,l))

!                             multiply all tau_xz/x terms by volume:
                              tauxz_ox = tauxz_ox*vol_w(ijk)

! x*tau_xz term: (1/x) (d/dz) (x*tau_xz)
!                             integral of (1/x)*d( (-x)*(mu/x)*w )/dx
!                             (normally part of source_w_s)
                              tauxz_x = -( ( mu_sl_pe*ox_e(i)*half*(w_s(ipjk,l)+&
                                   w_s(ijk,l))*ayz_w(ijk) )-( mu_sl_pw*half*&
                                   (w_s(ijk,l)+w_s(imjk,l))*&
                                   (dy(j)*half*(dz(k)+dz(kp)) ) ) )
                                   !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
                         ELSE
                              tauxz_ox = zero
                              tauxz_x = 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/6d0)*(d_pl**3.)*ro_s(ijk,l)

                         nu_pm_pt = rop_s(ijkt,m)/m_pm
                         nu_pm_pb = rop_s(ijk,m)/m_pm
                         nu_pm_p = avg_z(nu_pm_pb,nu_pm_pt,k)

                         nu_pl_pt = rop_s(ijkt,l)/m_pl
                         nu_pl_pb = rop_s(ijk,l)/m_pl
                         nu_pl_p = avg_z(nu_pl_pb,nu_pl_pt,k)

                         fnu_s_p = avg_z(fnu_s_ip(ijk,m,l),fnu_s_ip(ijkt,m,l),k)
                         ds1 = fnu_s_p*nu_pl_p*(nu_pm_pt-nu_pm_pb)*ox(i)*odz_t(k)
                         ds2 = -fnu_s_p*nu_pm_p*(nu_pl_pt-nu_pl_pb)*ox(i)*odz_t(k)
                         ds1plusds2 = ds1 + ds2

! Momentum source associated with the gradient in granular
! temperature of species M
                         t_pm_pt = theta_m(ijkt,m)
                         t_pm_pb = theta_m(ijk,m)

                         ft_sm_p = avg_z(ft_sm_ip(ijk,m,l),ft_sm_ip(ijkt,m,l),k)
                         ds3 = ft_sm_p*(t_pm_pt-t_pm_pb)*ox(i)*odz_t(k)

! Momentum source associated with the gradient in granular
! temperature of species L
                         t_pl_pt = theta_m(ijkt,l)
                         t_pl_pb = theta_m(ijk,l)

                         ft_sl_p = avg_z(ft_sl_ip(ijk,m,l),ft_sl_ip(ijkt,m,l),k)
                         ds4 = ft_sl_p*(t_pl_pt-t_pl_pb)*ox(i)*odz_t(k)
!--------------------- End Sources from Momentum Source Term ---------------------


! Add the terms
                    stress_terms = stress_terms + sswx + sswy + sswz + &
                        ssbv + ssx + ssy + ssz + tauxz_ox + tauxz_x
                    drag_terms = drag_terms + (ds1plusds2+ds3+ds4)*vol_w(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_w_s(ijk,m) = stress_terms + drag_terms
          ELSE
               ktmom_w_s(ijk,m) = zero
          ENDIF     ! dilute
      ENDDO     ! over ijk

      RETURN
      END SUBROUTINE calc_ia_momsource_w_s