drag_gp_des.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DES_DRAG_GP                                             C
!  Purpose: Calculate the gas-particle drag coefficient using          C
!           the gas velocity interpolated to the particle position     C
!           and the particle velocity.                                 C
!           Invoked from des_drag_gs and calc_des_drag_gs              C
!                                                                      C
!  Comments: The BVK drag model and all drag models with the           C
!            polydisperse correction factor (i.e., suffix _PCF)        C
!            require an average particle diameter. This has been       C
!            loosely defined for discrete particles based on their     C
!            solids phase                                              C
!                                                                      C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE des_drag_gp(NP, PARTICLE_VEL, FLUID_VEL, EPg)

! Modules
!---------------------------------------------------------------------//
      USE compar
      USE constant
      USE discretelement
      USE drag
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE machine, only: start_log, end_log
      USE param
      USE param1
      USE physprop
      USE run
      USE sendrecv
      USE ur_facs
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! particle number id.
      INTEGER , INTENT(IN) :: NP
! particle velocity
      DOUBLE PRECISION, INTENT(IN) :: PARTICLE_VEL(3)
! fluid velocity interpolated to particle position
      DOUBLE PRECISION, INTENT(IN) :: FLUID_VEL(3)
! Gas phase volume fraction.
      DOUBLE PREcISION, INTENT(IN) :: EPg

! Local variables
!---------------------------------------------------------------------//
! indices, associated with current particle
      INTEGER :: IJK
! solids phase index, associated with current particle
      INTEGER :: M
! Slip velocity and its magnitude
      DOUBLE PRECISION :: VSLP(3), VREL
! gas laminar viscosity redefined here to set viscosity at pressure
! boundaries
      DOUBLE PRECISION :: Mu
! drag coefficient
      DOUBLE PRECISION :: DgA
! current value of F_gs (i.e., without underrelaxation)
      DOUBLE PRECISION F_gstmp
! indices of solids phases (continuous, discrete)
      INTEGER :: lM
! correction factors for implementing polydisperse drag model
! proposed by van der Hoef et al. (2005)
      DOUBLE PRECISION :: F_cor, tSUM, tfac
! average particle diameter in polydisperse systems
      DOUBLE PRECISION :: DPA
! diameter ratio in polydisperse systems
      DOUBLE PRECISION :: Y_i
! total solids volume fraction and phase volume fraction
      DOUBLE PRECISION :: PHIS, phism
! aliases for gas density, gas bulk density,
      DOUBLE PRECISION :: ROg, ROPg
! particle diameter, particle density
      DOUBLE PRECISION :: DPM, ROd
!......................................................................!

! values based on current particle
      ijk = pijk(np,4)
! solids phase index of current particle
      m = pijk(np,5)
! Gas material and bulk densities
      rog = ro_g(ijk)
      ropg = ro_g(ijk) * epg
! Laminar viscosity.
      mu = mu_g(ijk)
! Slip velocity and its magnitude
      vslp = fluid_vel - particle_vel
      vrel = sqrt(dot_product(vslp, vslp))
! assign variables for short dummy arguments
      dpm = 2.0d0*des_radius(np)
      rod = ro_sol(np)
! Total solids volume fraction.
      phis = one - epg

! determine the drag coefficient
      SELECT CASE(drag_type_enum)

      CASE (syam_obrien)
         CALL drag_syam_obrien(dga,epg,mu,rog,vrel,dpm)

      CASE (gidaspow)
         CALL drag_gidaspow(dga,epg,mu,rog,ropg,vrel,dpm)

      CASE (gidaspow_blend)
         CALL drag_gidaspow_blend(dga,epg,mu,rog,ropg,vrel,dpm)

      CASE (wen_yu)
         CALL drag_wen_yu(dga,epg,mu,ropg,vrel,dpm)

      CASE (koch_hill)
         CALL drag_koch_hill(dga,epg,mu,ropg,vrel,dpm,dpm,phis)

      CASE (user_drag)
         CALL drag_usr(ijk,np,dga,epg,mu,rog,vrel,dpm,rod, &
            fluid_vel(1), fluid_vel(2), fluid_vel(3))

      CASE DEFAULT

! calculate the average particle diameter and particle ratio:
! the loop over discrete solids could be replaced using a loop
! over all particles in the given ijk cell but this should be
! done outside the current particle loop to minimize repetitive
! computational expense
         tsum = zero
         DO lm = mmax+1,des_mmax+mmax
            phism = ep_s(ijk,lm)     ! PVOL(NP)/VOL(IJK)
            IF (phis .GT. zero) THEN
               tfac = phism/phis
               tsum = tsum + tfac/dpm
            ELSE
               tsum = tsum + one/dpm
            ENDIF
         ENDDO
         IF(des_continuum_hybrid) THEN
            DO lm = 1,mmax
               phism = ep_s(ijk,lm)
               IF(phis > zero) THEN
                  tfac = phism/phis
                  tsum = tsum + tfac/d_p(ijk,lm)
                ELSE
                  tsum = tsum + one/d_p(ijk,lm)
                ENDIF
            ENDDO
         ENDIF
         dpa = one / tsum
         y_i = dpm / dpa

         SELECT CASE(drag_type_enum)
         CASE (gidaspow_pcf)
            CALL drag_gidaspow(dga,epg,mu,rog,ropg,vrel,dpa)
         CASE (gidaspow_blend_pcf)
            CALL drag_gidaspow_blend(dga,epg,mu,rog,ropg,vrel,dpa)
         CASE (wen_yu_pcf)
            CALL drag_wen_yu(dga,epg,mu,ropg,vrel,dpa)
         CASE (koch_hill_pcf)
            CALL drag_koch_hill(dga,epg,mu,ropg,vrel,dpm,dpa,phis)
         CASE (bvk)
            CALL drag_bvk(dga,epg,mu,ropg,vrel,dpm,dpa,phis)

         CASE DEFAULT
            CALL start_log
            IF(dmp_log) WRITE (*, '(A,A)') &
               'Unknown DRAG_TYPE: ', drag_type
            WRITE (unit_log, '(A,A)') 'Unknown DRAG_TYPE: ', drag_type
            CALL end_log
            CALL mfix_exit(mype)
         END SELECT   ! end selection of drag_type

! Modify drag coefficient to account for possible corrections and for
! differences between Model B and Model A; see erratum Beetstra (2007)
         IF(model_b) THEN
            IF (m == 1) THEN
               f_cor = (epg*y_i + phis*y_i**2)
            ELSE
               f_cor = (epg*y_i + phis*y_i**2 + &
                  0.064d0*epg*y_i**3)
            ENDIF
         ELSE
            f_cor = y_i
         ENDIF
         dga = one/(y_i*y_i) * dga * f_cor

      END SELECT   ! end selection of drag_type



! Calculate the drag coefficient (Model B coeff = Model A coeff/EP_g)
      IF(model_b) THEN
         f_gstmp = dga * pvol(np)/ep_g(ijk)
      ELSE
         f_gstmp = dga * pvol(np)
      ENDIF

! Determine drag force coefficient accounting for any under relaxation
! f_gp() =  single particle drag excluding vector(v_g - v_p)
      f_gp(np) = (one - ur_f_gs) * f_gp(np) + ur_f_gs * f_gstmp

      RETURN
      END SUBROUTINE des_drag_gp