calc_vol_fr.f

Go to the documentation of this file.

! TO DO:
! 1. Check the formulation based on MCP.
! 2. The pressure correction should be based on sum of close-packed
!    solids?

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_VOL_FR                                             C
!  Purpose: Calculate volume fractions of phases used in pressure      C
!           corrections.                                               C
!                                                                      C
!  Notes: see mark_phase_4_cor for more details                        C
!                                                                      C
!  Author: M. Syamlal                                 Date: 5-JUL-96   C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_vol_fr(P_STAR, RO_G, ROP_G, EP_G, ROP_S, IER)


! Modules
!---------------------------------------------------------------------//
      USE compar
      USE constant
      USE discretelement
      USE functions
      USE geometry
      USE indices
      USE mpi_utility
      USE parallel
      USE param
      USE param1
      USE pgcor
      USE physprop
      USE pscor
      USE run, only: ghd_2007, kt_type_enum
      USE sendrecv
      USE solids_pressure
      USE visc_s
      use fldvar, only: ro_s, ep_s

      IMPLICIT NONE

! Dummy Arguments
!---------------------------------------------------------------------//
! Solids pressure
      DOUBLE PRECISION, INTENT(IN) :: P_star(dimension_3)
! Gas density
      DOUBLE PRECISION, INTENT(INOUT) :: RO_g(dimension_3)
! Gas bulk density
      DOUBLE PRECISION, INTENT(INOUT) :: ROP_g(dimension_3)
! Gas volume fraction
      DOUBLE PRECISION, INTENT(INOUT) :: EP_g(dimension_3)
! solids bulk densities
      DOUBLE PRECISION, INTENT(INOUT) :: ROP_s(dimension_3, dimension_m)
! error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! volume of particle type M for GHD theory
      DOUBLE PRECISION :: VOL_M
! volume fraction of close-packed region
      DOUBLE PRECISION :: EPcp
! volume fraction of solids phase
      DOUBLE PRECISION :: EPS
! sum of volume fractions
      DOUBLE PRECISION :: SUMVF
! Whichever phase is given by MF becomes the phase whose volume fraction
! is corrected based on all other phases present.  Generally MF gets
! defined as 0 so that the gas phase void fraction becomes corrected
! based on the summation of the volume fractions of all other phases
      INTEGER :: MF
! Whichever phase is given by MCPl becomes the phase whose volume
! fraction is corrected based on value of maximum close packing and all
! other solids phase that can close pack.  This is basically a local
! variable for MCP but only in cells where close packed conditions
! exist.
      INTEGER :: MCPl
! Index of solids phase
      INTEGER :: M
! Indices
      INTEGER :: IJK

! Arrays for storing errors:
! 110 - Negative volume fraciton
! 11x - Unclassified
      INTEGER :: Err_l(0:numpes-1)  ! local
      INTEGER :: Err_g(0:numpes-1)  ! global

      LOGICAL, PARAMETER :: REPORT_NEG_VOLFRAC = .true.

! Flag to write log header
      LOGICAL :: wHeader
!---------------------------------------------------------------------//


! Initialize:
      wheader = .true.
! Initialize error flag.
      err_l = 0

!!$omp  parallel do private(MCPl, EPCP, SUMVF, MF, M) &
!!$omp&  schedule(static)

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

! calculate the volume fraction of the solids phase based on the volume
! fractions of all other solids phases and on the value of maximum
! packing when that solids phase continuity was not solved for the
! indicated cell.
!----------------------------------------------------------------->>>
            IF (phase_4_p_s(ijk) /= undefined_i) THEN
! for the current cell check if a solids phase continuity was skipped.
! this value will either be undefined (no cell was skipped in any
! solids continuity equations) or be a solids phase index (indicated
! solids continuity was skipped) for a solids phase that does close pack
! (i.e., close_packed=T) and the cell exhibits close packing conditions.
! note: this branch is never entered given the existing version of
! mark_phase_4_cor.
               mcpl   = phase_4_p_s(ijk)

! finding the value of maximum close packing based on the expression for
! plastic pressure
               epcp = 1. - inv_h(p_star(ijk),ep_g_blend_end(ijk))
! summing the solids volume fraction of all continuum solids phases
! that are marked as close_packed except for the solids phase which
! is also marked by phase_4_p_s (skipping its own continuity)
               sumvf = zero
               DO m = 1, mmax+des_mmax
                  IF (close_packed(m) .AND. m/=mcpl) sumvf = sumvf + ep_s(ijk,m)
               ENDDO
               rop_s(ijk,mcpl) = (epcp - sumvf)*ro_s(ijk,mcpl)
            ENDIF
!-----------------------------------------------------------------<<<


! calculate the volume fraction of the 'solids' phase based on the
! volume fractions of all other phases (including gas) if the gas
! continuity was solved rather than that phases own continuity.
! if the gas continuity was not solved then calculate the void
! fraction of the gas phase based on all other solids phases.
!----------------------------------------------------------------->>>
! for the current cell check if the gas phase continuity was solved for
! the gas phase while the solids phase continuity for the indicated
! solids phase was skipped. this value will either be 0 (gas phase
! continuity was not solved) or be a solids phase index (gas continuity
! was solved while the indicated solids phase continuity was skipped)
! for a solids phase that does not close pack (i.e., close_packed=F) and
! is in greater concentration than the gas phase.
! Note: MF will always be set to 0 here given the existing version of
! mark_phase_4_cor
            mf = phase_4_p_g(ijk)

            sumvf = zero
            IF (0 /= mf) THEN
! if gas continuity was solved rather than the solids phase, then
! include the gas phase void fraction in the summation here.
               ep_g(ijk) = rop_g(ijk)/ro_g(ijk)
               sumvf = sumvf + ep_g(ijk)
            ENDIF

! evaluate mixture bulk density rop_s(mmax) for GHD theory
            IF(kt_type_enum == ghd_2007) THEN
              rop_s(ijk,mmax) = zero  ! mixture density
              DO m = 1, smax
                 IF (m /= mf) THEN
                   rop_s(ijk,mmax) = rop_s(ijk,mmax) + ro_s(ijk,m)*ep_s(ijk,m)
                 ENDIF
              ENDDO
            ENDIF

! summing the solids volume fraction of all continuum solids phases
! except for the continuum solids phase which was marked by phase_4_p_g
! (skipping its continuity while solving gas continuity)
            DO m = 1, des_mmax+mmax
               IF(kt_type_enum == ghd_2007 .AND. m == mmax) cycle
               IF (m /= mf) sumvf = sumvf + ep_s(ijk,m)
            ENDDO

            IF (0 == mf) THEN
! if no gas phase continuity was solved in the current cell then correct
! the void fraction of the gas phase based on the total solids volume
! fraction of all solids phases
               ep_g(ijk) = one - sumvf

! Set error flag for negative volume fraction.
               IF (ep_g(ijk) < zero) THEN
                  err_l(mype) = 110
                  IF(report_neg_volfrac) CALL epgerr_log(ijk, wheader)
               ENDIF

               rop_g(ijk) = ep_g(ijk)*ro_g(ijk)
            ELSE
! else correct the volume fraction of the solids phase that was marked
               rop_s(ijk,mf) = (one - sumvf)*ro_s(ijk,mf)
            ENDIF
!-----------------------------------------------------------------<<<

         ENDIF    ! end if (fluid_at(ijk))
      ENDDO   ! end do (ijk=ijkstart3,ijkend3)

      CALL send_recv(ep_g, 2)
      CALL send_recv(rop_g, 2)
      CALL send_recv(rop_s, 2)

      CALL global_all_sum(err_l, err_g)
      ier = maxval(err_g)


      RETURN
      END SUBROUTINE calc_vol_fr

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_EP_FACTORS                                          C
!  Purpose: Set multiplication factors to ep_g to control the form     C
!  of the governing equations                                          C
!                                                                      C
!  References:                                                         C
!  Anderson and Jackson, A fluid mechanical description of fluidized   C
!     beds, Ind. Eng. Chem. Fundam., 6(4) 1967, 527-539.               C
!  Ishii, Thermo-fluid dynamic theory of two-phase flow, des Etudes    C
!     et Recherches D'Electricite de France, Eyrolles, Paris, France   C
!     1975.                                                            C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_ep_factors

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use fldvar, only: ep_g, ep_s
      use fldvar, only: epg_ifac, eps_ifac, epg_jfac
      use functions, only: wall_at
      use param1, only: one, undefined
      use physprop, only: mmax, mu_s0
      use run, only: jackson, ishii
      IMPLICIT NONE
! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: ijk, m
!---------------------------------------------------------------------//
      epg_jfac(:) = one
      epg_ifac(:) = one
      eps_ifac(:,:) = one

      if (jackson) then
         do ijk = ijkstart3, ijkend3
            if (wall_at(ijk)) cycle
            epg_jfac(ijk) = ep_g(ijk)
         enddo
      endif
      if (ishii) then
         do ijk = ijkstart3, ijkend3
            if (wall_at(ijk)) cycle
            epg_ifac(ijk) = ep_g(ijk)
            do m = 1, mmax
! This seems more appropriate for non-granular systems
! So currently only incorporate this for 'constant' viscosity cases
               IF (mu_s0(m) /= undefined) THEN
                  eps_ifac(ijk,m) = ep_s(ijk,m)
               ENDIF
            enddo
         enddo
      endif

      RETURN
      END SUBROUTINE set_ep_factors

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: EPgErr_LOG                                              !
!  Purpose: Record information about the location and conditions that  !
!           resulted in a negative gas phase volume fraction.          !
!                                                                      !
!  Author: J. Musser                                  Date: 16-APR-15  !
!  Reviewer:                                          Date:            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE epgerr_log(IJK, tHeader)

! Simulation time
      use run, only: time
! Gas phase pressure.
      use fldvar, only: ep_g, ep_s, rop_s, ro_s, u_s, v_s, w_s
      use cutcell

      use physprop, only: smax
      use compar
      use indices
      use functions
      use error_manager


      IMPLICIT NONE

      INTEGER, intent(in) :: IJK
      LOGICAL, intent(inout) :: tHeader

      LOGICAL :: lExists
      CHARACTER(LEN=255) :: lFile
      INTEGER, parameter :: lUnit = 4868
      LOGICAL, save :: fHeader = .true.

      INTEGER :: M

      lfile = '';
      if(numpes > 1) then
         write(lfile,"('EPgErr_',I4.4,'.log')") mype
      else
         write(lfile,"('EPgErr.log')")
      endif
      inquire(file=trim(lfile),exist=lexists)
      if(lexists) then
         open(lunit,file=trim(adjustl(lfile)),                         &
            status='old', position='append')
      else
         open(lunit,file=trim(adjustl(lfile)), status='new')
      endif

      if(fheader) then
         write(lunit,1000)
         fheader = .false.
      endif

      if(theader) then
         write(lunit,"(/2x,'Simulation time: ',g12.5)") time
         theader = .false.
      endif

      write(lunit,1001) ijk, i_of(ijk), j_of(ijk), k_of(ijk)
      write(lunit,"(6x,A,1X,g12.5)",advance='NO') 'EP_g:', ep_g(ijk)
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('EP_s',m)), ep_s(ijk,m)
      ENDDO
      write(lunit,"(' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('RO_s',m)), ro_s(ijk,m)
      ENDDO
      write(lunit,"(' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('ROPs',m)), rop_s(ijk,m)
      ENDDO
      write(lunit,"(24x,' ')")


      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('U_se',m)), u_s(ijk,m)
      ENDDO
      write(lunit,"(24x,' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('V_sn',m)), v_s(ijk,m)
      ENDDO
      write(lunit,"(24x,' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('W_st',m)), w_s(ijk,m)
      ENDDO
      write(lunit,"(24x,' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('U_sw',m)), u_s(west_of(ijk),m)
      ENDDO
      write(lunit,"(24x,' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('V_ss',m)), v_s(south_of(ijk),m)
      ENDDO
      write(lunit,"(24x,' ')")

      write(lunit,"(24x)", advance='NO')
      DO m=1,smax
         write(lunit,"(2x,A,1X,g12.5)", advance='NO') &
            trim(ivar('W_sb',m)), w_s(bottom_of(ijk),m)
      ENDDO
      write(lunit,"(24x,' ')")


      if(cartesian_grid) then
         write(lunit,"(6x,A,1X,L1)",advance='NO') 'Cut Cell:', cut_cell_at(ijk)
         write(lunit,"(6x,A,1X,L1)") 'Small Cell:', small_cell_at(ijk)
         write(lunit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
            xg_e(i_of(ijk)), yg_n(j_of(ijk)), zg_t(k_of(ijk))
      endif

      close(lunit)

      RETURN

 1000 FORMAT('One or more cells have reported a negative gas volume ', &
         'fraction (EP_g).',/)

 1001 FORMAT(/4x,'IJK: ',i8,7x,'I: ',i4,'  J: ',i4,'  K: ',i4)

      END SUBROUTINE epgerr_log