calc_ps_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_PIC                                             !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Calculate the particle stress.                             !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_ps_pic


! Global Variables:
!---------------------------------------------------------------------//
! Flag to use Snider's particle stress model
      use mfix_pic, only: mppic_solid_stress_snider
! Particle stress
      use mfix_pic, only: pic_p_s

! Module procedures:
!---------------------------------------------------------------------//
      use sendrecv, only: send_recv

      IMPLICIT NONE

!......................................................................!


      IF(mppic_solid_stress_snider) THEN
         CALL calc_ps_pic_snider
      ELSE
         CALL calc_ps_pic_garg
      ENDIF

      CALL send_recv(pic_p_s,1)

      RETURN
      END SUBROUTINE calc_ps_pic



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_PIC                                             !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Evaluate the particle stress model of Snider.              !
!                                                                      !
!  REF: D.M. Snider, "Three-Dimensional Multiphase Particle-in-Cell    !
!     Model for Dense Particle Flows," Journal of Computational        !
!     Physics, Vol. 170, No. 2, pp. 523-549, 2001.                     !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_ps_pic_snider

! Global Variables:
!---------------------------------------------------------------------//
! Model parameters for Snider particle stress model
      use mfix_pic, only: fric_exp_pic
      use mfix_pic, only: psfac_fric_pic
      use mfix_pic, only: fric_non_sing_fac
! Calculated particle stress
      use mfix_pic, only: pic_p_s
! Fluid phase volume fraction
      use fldvar, only: ep_g
! Fluid volume fraction at close-pack
      use constant, only: ep_star
! Domain bounds
      use compar, only: ijkstart3, ijkend3
! Double precision parameters
      use param1, only: one

! Module procedures:
!---------------------------------------------------------------------//
      use functions, only: fluid_at
      use functions, only: west_of, east_of
      use functions, only: south_of, north_of
      use functions, only: bottom_of, top_of

      IMPLICIT NONE

! Local Variables:
!---------------------------------------------------------------------//
! Loop counter
      INTEGER :: IJK
! Volume fraction of cell, modified for wall cells.
      DOUBLE PRECISION :: lEPg
!......................................................................!

      DO ijk = ijkstart3, ijkend3

         IF(fluid_at(ijk)) THEN
            lepg = ep_g(ijk)
         ELSE

! Set the volume fraction in the wall to close pack.
            lepg = ep_star
 
! Use the lowest value across all adjacent fluid cells. This is to keep
! cells below close pack from pushing parcels through the walls.
!            lIJK = EAST_OF(IJK)
!            IF(FLUID_AT(lIJK))  lEPg = min(lEPg, EP_G(lIJK))
!!            lIJK = WEST_OF(IJK)
!            IF(FLUID_AT(lIJK))  lEPg = min(lEPg, EP_G(lIJK))
!            lIJK = NORTH_OF(IJK)
!            IF(FLUID_AT(lIJK))  lEPg = min(lEPg, EP_G(lIJK))
!!            lIJK = SOUTH_OF(IJK)
!            IF(FLUID_AT(lIJK))  lEPg = min(lEPg, EP_G(lIJK))
!            IF(DO_K) THEN
!               lIJK = TOP_OF(IJK)
!               IF(FLUID_AT(lIJK))  lEPg = min(lEPg, EP_G(lIJK))
!               lIJK = BOTTOM_OF(IJK)
!               IF(FLUID_AT(lIJK))  lEPg = min(lEPg, EP_G(lIJK))
!            ENDIF
         ENDIF

! Particle stress :: Snider (Eq 33)
         pic_p_s(ijk,1) = psfac_fric_pic *((one - lepg)**fric_exp_pic)/&
            max(lepg - ep_star, fric_non_sing_fac*lepg)
      ENDDO



      RETURN
      END SUBROUTINE calc_ps_pic_snider


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_PIC                                             !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Evaluate the particle stress as a coloring function:       !
!     X=0.0 :: cells below packing limit                               !
!     X=EPg :: cells above packing limit                               !
!     X=1.0 :: wall cells (far avobe max packing)                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_ps_pic_garg


! Global Variables:
!---------------------------------------------------------------------//
! Calculated particle stress
      use mfix_pic, only: pic_p_s
! Resulting particle stress force
      use mfix_pic, only: ps_force_pic
! Fluid phase volume fraction
      use fldvar, only: ep_g
! Fluid volume fraction at close-pack
      use constant, only: ep_star
! Domain bounds
      use compar, only: ijkstart3, ijkend3
! Double precision parameters
      use param1, only: zero, one

! Module procedures:
!---------------------------------------------------------------------//
      use functions, only: fluid_at


      IMPLICIT NONE


! Local Variables:
!---------------------------------------------------------------------//
! Loop counter
      INTEGER :: IJK
!......................................................................!


! The Garg model uses a coloring function approach. 
      DO ijk = ijkstart3, ijkend3
         ps_force_pic(:,ijk) = zero
         IF(fluid_at(ijk)) THEN
            IF(ep_g(ijk) < ep_star) THEN
               pic_p_s(ijk,1) = (one - ep_g(ijk))
            ELSE
               pic_p_s(ijk,1) = zero
            ENDIF
         ELSE
            pic_p_s(ijk,1) = one
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE calc_ps_pic_garg