calc_u_friction.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_Gw_Hw_Cw                                           C
!                                                                      C
!  Purpose: Calculate Gw, Hw, and Cw                                   C
!                                                                      C
!  Author: A. Srivastava & K. Agrawal, Princeton Univ. Date: 10-APR-98 C
!  Reviewer:                                           Date:           C
!                                                                      C
!  Modified: Sofiane Benyahia, Fluent Inc.             Date: 03-FEB-05 C
!  Purpose: Include conductivity defined by Simonin and Ahmadi         C
!           Also included Jenkins small frictional limit               C
!                                                                      C
!  Literature/Document References:                                     C
!     Johnson, P. C., and Jackson, R., Frictional-collisional          C
!        constitutive relations for granular materials, with           C
!        application to plane shearing, Journal of Fluid Mechanics,    C
!        1987, 176, pp. 67-93                                          C
!     Jenkins, J. T., and Louge, M. Y., On the flux of fluctuating     C
!        energy in a collisional grain flow at a flat frictional wall, C
!        Physics of Fluids, 1997, 9(10), pp. 2835-2840                 C
!                                                                      C
!     See calc_mu_s.f for references on frictional theory models       C
!     See calc_mu_s.f for references on Ahmadi and Simonin models      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE calc_gw_hw_cw(g0, EPS, EPG, ep_star_avg, &
            g0ep_avg, th, mu_g_avg, ro_g_avg, ros_avg, &
            dp_avg, k_12_avg, tau_12_avg, tau_1_avg, vrel, vslip,&
            del_u, s_s, s_dd, vel, w_vel, m, gw, hw, cw)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE constant
      USE physprop
      USE fldvar
      USE bc
      USE run
      USE mpi_utility
      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Radial distribution function of solids phase M with each
! other solids phase
      DOUBLE PRECISION, INTENT(IN) :: g0
! Average solids volume fraction of each solids phase
      DOUBLE PRECISION, INTENT(IN) :: EPS
! Average solids and gas volume fraction
      DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
! Sum of eps*G_0
      DOUBLE PRECISION, INTENT(IN) :: g0EP_avg
! Average theta_m
      DOUBLE PRECISION, INTENT(INOUT) :: TH
! Average gas viscosity
      DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
! Average gas density
      DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
! Average solids density
      DOUBLE PRECISION, INTENT(IN) :: ROS_avg
! Average particle diameter of each solids phase
      DOUBLE PRECISION, INTENT(IN) :: DP_avg
! Average cross-correlation K_12 and lagrangian integral time-scale
      DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
! Magnitude of slip velocity between two phases
      DOUBLE PRECISION, INTENT(IN) :: VREL
! Slip velocity between wall and particles
      DOUBLE PRECISION, INTENT(IN) :: VSLIP
! Relevant solids velocity at wall
      DOUBLE PRECISION, INTENT(IN) :: VEL
! Relevant wall velocity
      DOUBLE PRECISION, INTENT(IN) :: W_VEL
! del.u
      DOUBLE PRECISION, INTENT(IN) :: DEL_U
! S:S
      DOUBLE PRECISION, INTENT(IN) :: S_S
! S_dd (d can be x,y or z)
      DOUBLE PRECISION, INTENT(IN) :: S_dd
! Solids phase index
      INTEGER, INTENT(IN) :: M
! Wall momentum coefficients:
! 1st term on LHS
      DOUBLE PRECISION, INTENT(INOUT) :: Gw
! 2nd term on LHS
      DOUBLE PRECISION, INTENT(INOUT) :: Hw
! all terms appearing on RHS
      DOUBLE PRECISION, INTENT(INOUT) :: Cw
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! Part of wall momentum coefficient in 1st term of LHS
      DOUBLE PRECISION :: Mu_s
! Term appearing in wall momenutm coefficient of 2nd term LHS
      DOUBLE PRECISION :: F_2
! Frictional Pressure Pf
      DOUBLE PRECISION :: Pf
! Critical pressure Pc
      DOUBLE PRECISION :: Pc
! Chi (appears in frictional boundary condition)
      DOUBLE PRECISION :: Chi, N_Pff
! Viscosity
      DOUBLE PRECISION :: Mu
! Bulk viscosity
      DOUBLE PRECISION :: Mu_b
! Viscosity corrected for interstitial fluid effects
      DOUBLE PRECISION :: Mu_star
! Reynolds number based on slip velocity
      DOUBLE PRECISION :: Re_g
! Friction Factor in drag coefficient
      DOUBLE PRECISION :: C_d
! Drag Coefficient
      DOUBLE PRECISION :: Beta, DgA
! Square root of S:S or the form suggested by Savage
      DOUBLE PRECISION :: ZETA
! Constants in Simonin model
      DOUBLE PRECISION :: Sigma_c, Tau_2_c, Tau_12_st, Nu_t
      DOUBLE PRECISION :: Tau_2, zeta_c_2, MU_2_T_Kin, Mu_2_Col
      DOUBLE PRECISION :: Tmp_Ahmadi_Const
! Other local terms
      DOUBLE PRECISION :: dpc_dphi
! Local variable for rdf
      DOUBLE PRECISION :: g0_loc
!-----------------------------------------------

! This is done here similar to bc_theta to avoid small negative values of
! Theta coming most probably from linear solver
      IF(th .LE. zero)THEN
        th = 1d-8
        IF (mype.eq.pe_io) THEN
           WRITE(*,*) &
              'Warning: Negative granular temp at wall set to 1e-8'
!          CALL WRITE_ERROR('THETA_HW_CW', LINE, 1)
        ENDIF
      ENDIF

      g0_loc = g0

! modify F_2 if Jenkins BC is used (sof)
      IF(jenkins) THEN

        IF (vslip == zero) THEN
! if solids velocity field is initialized to zero, use free slip bc
           f_2 = zero

        ELSE
           IF(ahmadi) THEN
! Ahmadi model uses different solids pressure model
! the coefficient mu in Jenkins paper is defined as tan_Phi_w, that's how
! I understand it from soil mechanic papers, i.e., G.I. Tardos, powder
! Tech. 92 (1997), 61-74. See his equation (1). Define Phi_w in mfix.dat!
! here F_2 divided by VSLIP to use the same bc as Johnson&Jackson

              f_2 = tan_phi_w*ros_avg*eps* &
                 ((one + 4.0d0*g0ep_avg) + half*(one -c_e*c_e))*th/vslip

           ELSE
! Simonin or granular models use same solids pressure
              f_2 = tan_phi_w*ros_avg*eps*(1d0+ 4.d0 * eta *g0ep_avg)*th/vslip
           ENDIF   ! end if(Ahmadi)

        ENDIF ! endif(vslip==0)

      ELSE    ! if(.not.jenkins)

         f_2 = (phip*dsqrt(3d0*th)*pi*ros_avg*eps*g0_loc)&
              /(6d0*(one-ep_star_avg))

      ENDIF   ! end if(Jenkins)/else


      mu = (5d0*dsqrt(pi*th)*dp_avg*ros_avg)/96d0
      mu_b = (256d0*mu*eps*g0ep_avg)/(5d0*pi)

! This is from Wen-Yu correlation, you can put here your own single particle drag
      re_g = epg*ro_g_avg*dp_avg*vrel/mu_g_avg
      IF (re_g.lt.1000d0) THEN
         c_d = (24.d0/(re_g+small_number))*(one + 0.15d0 * re_g**0.687d0)
      ELSE
         c_d = 0.44d0
      ENDIF
      dga = 0.75d0*c_d*ro_g_avg*epg*vrel/(dp_avg*epg**(2.65d0))
      IF(vrel == zero) dga = large_number
      beta = switch*eps*dga

! SWITCH enables us to turn on/off the modification to the
! particulate phase viscosity. If we want to simulate gas-particle
! flow then SWITCH=1 to incorporate the effect of drag on the
! particle viscosity. If we want to simulate granular flow
! without the effects of an interstitial gas, SWITCH=0.
      IF(switch == zero .OR. ro_g_avg == zero)THEN
         mu_star = mu
      ELSEIF(th .LT. small_number)THEN
         mu_star = zero
      ELSE
         mu_star = ros_avg*eps* g0_loc*th* mu/ &
            (ros_avg*g0ep_avg*th + 2.0d0*switch*dga/ros_avg* mu)
      ENDIF

      mu_s = ((2d0+alpha)/3d0)*((mu_star/(eta*(2d0-eta)*&
                   g0_loc))*(one+1.6d0*eta*g0ep_avg&
                   )*(one+1.6d0*eta*(3d0*eta-2d0)*&
                   g0ep_avg)+(0.6d0*mu_b*eta))

! particle relaxation time
      tau_12_st = ros_avg/(dga+small_number)

      IF(simonin) THEN !see calc_mu_s for explanation of these definitions

         sigma_c = (one+ c_e)*(3.d0-c_e)/5.d0
         tau_2_c = dp_avg/(6.d0*eps*g0_loc*dsqrt(16.d0*(th+small_number)/pi))
         zeta_c_2= 2.d0/5.d0*(one+ c_e)*(3.d0*c_e-one)
         nu_t =  tau_12_avg/tau_12_st
         tau_2 = one/(2.d0/tau_12_st+sigma_c/tau_2_c)
         mu_2_t_kin = (2.0d0/3.0d0*k_12_avg*nu_t + th * &
                     (one+ zeta_c_2*eps*g0_loc))*tau_2
         mu_2_col = 8.d0/5.d0*eps*g0_loc*eta* (mu_2_t_kin+ &
                   dp_avg*dsqrt(th/pi))
         mu_s = eps*ros_avg*(mu_2_t_kin + mu_2_col)
      ELSEIF(ahmadi) THEN
         IF(eps < (one-ep_star_avg)) THEN
            tmp_ahmadi_const = &
            one/(one+ tau_1_avg/tau_12_st * (one-eps/(one-ep_star_avg))**3)
         ELSE
            tmp_ahmadi_const = one
         ENDIF
         mu_s = tmp_ahmadi_const &
             *0.1045d0*(one/g0_loc+3.2d0*eps+12.1824d0*g0_loc*eps*eps)  &
             *dp_avg*ros_avg* dsqrt(th)
      ENDIF

! Calculating frictional terms
      IF ((one-epg)<= eps_f_min) THEN
         pf = zero
         chi = zero
         zeta = 1d0
      ELSE
         IF (savage.EQ.1) THEN    !form of Savage
            zeta = ((48d0*eta*(1d0-eta)*ros_avg*eps*eps*g0_loc*&
                    (th**1.5d0))/&
                    (sqrt_pi*dp_avg*2d0*mu_s))**0.5d0
         ELSEIF (savage.EQ.0)  THEN !S:S form
            zeta = dsqrt(s_s)
         ELSE
            zeta = dsqrt(s_s + (th/(dp_avg*dp_avg)))
         ENDIF

         IF (epg < ep_star_avg) THEN
            dpc_dphi = (to_si*fr)*((delta**5)*(2d0*(one-ep_star_avg-delta) - &
               2d0*eps_f_min)+((one-ep_star_avg-delta)-eps_f_min)&
               *(5*delta**4))/(delta**10)

            pc = (to_si*fr)*(((one-ep_star_avg-delta) - eps_f_min)**n_pc)/(delta**d_pc)
!            Pc=  1d25*(((ONE-EPG)-(ONE-ep_star_avg))**10d0)  ! this is old Pc
         ELSE
            pc = fr*(((one-epg) - eps_f_min)**n_pc)/ &
               (((one-ep_star_avg) - (one-epg) + small_number)**d_pc)
         ENDIF

         IF (del_u .GE. zero) THEN
            n_pff = dsqrt(3d0)/(2d0*sin_phi) !dilatation
         ELSE
            n_pff = n_pf !compaction
         ENDIF

         IF ((del_u/(zeta*n_pff*dsqrt(2d0)*sin_phi)) .GT. 1d0) THEN
            pf = zero
         ELSEIF( del_u == zero ) THEN
            pf = pc
         ELSE
            pf = pc*(1d0 - (del_u/(zeta*n_pff*dsqrt(2d0)*sin_phi)))**&
               (n_pff-1d0)
         ENDIF

         chi =  dsqrt(2d0)*pf*sin_phi*(n_pff - (n_pff-1d0)*&
            ((pf/pc)**(1d0/(n_pff-1d0))))

         IF (chi< zero) THEN
            pf = pc*((n_pff/(n_pff-1d0))**(n_pff-1d0))
            chi = zero
         ENDIF

! by writing Pf & Chi in the following form, we ensure distribution
! of stresses amoung all solids phases (sof, Oct 24 2005)

         pf = pf * eps/(one-epg)
         chi = chi * eps/(one-epg)

      ENDIF

! Calculating gw, hw, cw

      gw = one   ! we write this in the exact same form as non-frictional JJ BC

      IF(vslip == zero .OR. zeta == zero) THEN
        hw = zero
      ELSE
        hw = f_2 + pf*tan_phi_w - chi*s_dd*tan_phi_w/zeta
        hw = hw / (mu_s + chi/(2d0*zeta))  ! this is because Gw is set to one.
      ENDIF
      cw = hw * w_vel

      RETURN
      END SUBROUTINE calc_gw_hw_cw