adjust_theta.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: ADJUST_THETA                                            C
!  Purpose: Remove small negative values of theta caused by linear     C
!           solvers                                                    C
!                                                                      C
!  Author: M. Syamlal                                 Date: 02-APR-98  C
!                                                                      C
!  Modified: S. Benyahia                              Date: 02-AUG-06  C
!  Purpose: check for small negative numbers at walls                  C
!           (not just fluid cells)                                     C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE adjust_theta(M, IER)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1, only: zero
      USE toleranc, only: zero_ep_s
      USE constant, only: pi, to_si
! granular temperature of solids phase m
      USE fldvar, only: theta_m
! material density of solids phase m
      USE fldvar, only: ro_s
! particle diameter of solids phase m
      USE fldvar, only: d_p
! number of solids phases
      USE physprop, only: smax
! kt types
      USE run, only: kt_type
      USE run, only: kt_type_enum
      USE run, only: lun_1984
      USE run, only: simonin_1996
      USE run, only: ahmadi_1995
      USE run, only: gd_1999
      USE run, only: gtsh_2012
      USE run, only: ia_2005
      USE run, only: ghd_2007
! needed for function.inc
      USE compar
      USE exit, only: mfix_exit
      USE functions
      USE geometry
      USE indices

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! solids phase
      INTEGER, INTENT(IN) :: M
! error indicator
      INTEGER, INTENT(INOUT) :: IER
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK
! Solids phase index
      INTEGER :: L
! Particle mass and diameter for use with those kinetic theories
! that include mass of particle in definition of theta
      DOUBLE PRECISION :: M_PM, D_PM
! small value of theta_m
      DOUBLE PRECISION :: smallTheta
!-----------------------------------------------

      ier = 0
      smalltheta = (to_si)**4 * zero_ep_s

      DO ijk = ijkstart3, ijkend3
        IF ( fluid_at(ijk) ) THEN

          SELECT CASE(kt_type_enum)
            CASE (lun_1984, simonin_1996, ahmadi_1995, gd_1999, &
                  gtsh_2012)
              IF (theta_m(ijk,m) < smalltheta) &
                 theta_m(ijk,m) = smalltheta

            CASE (ia_2005)
              d_pm = d_p(ijk,m)
              m_pm = (pi/6.d0)*(d_pm**3)*ro_s(ijk,m)
              IF (theta_m(ijk,m) < smalltheta*m_pm) &
                theta_m(ijk,m) = smalltheta*m_pm

            CASE (ghd_2007)
              m_pm = zero
              DO l = 1,smax
                d_pm = d_p(ijk,l)
                m_pm = m_pm +(pi/6.d0)*(d_pm**3)*ro_s(ijk,l)
              ENDDO
              m_pm = m_pm/dble(smax)
              IF (theta_m(ijk,m) < smalltheta*m_pm) &
                theta_m(ijk,m) = smalltheta*m_pm

            CASE DEFAULT
! should never hit this
               WRITE (*, '(A)') 'ADJUST_THETA'
               WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', kt_type
               call mfix_exit(mype)
          END SELECT   ! end selection of kt_type_enum
        ENDIF   ! end if (fluid_at)
      ENDDO  ! end do ijk

      RETURN
      END SUBROUTINE adjust_theta