rrates0.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: RRATES0(IER)                                           C
!  Purpose: Calculate reaction rates for various reactions in cell ijk C
!           using information from the data file                       C
!                                                                      C
!  Author: M. Syamlal                                 Date: 3-10-98    C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose:Replaced routines with new proceedures for automated        C
!          reaction rate calculations.                                 C
!  Author: J. Musser                                  Date: 10-Oct-12  C
!  Reviewer:                                          Date: dd-mmm-yy  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced: MMAX, IJK, T_g, T_s1, D_p, X_g, X_s, EP_g,    C
!            P_g, HOR_g, HOR_s                                         C
!                                                                      C
!                                                                      C
!  Variables modified: M, N, R_gp, R_sp, RoX_gc, RoX_sc, SUM_R_g,      C
!                      SUM_R_s                                         C
!                                                                      C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE rrates0()

! Global domain parameters.
!`````````````````````````````````````````````````````````````````````//
! Blanket access is provided to the following modules as they don't
! contain variables which are likely to be modified 'by accident.'
      use compar
      use indices
      use geometry

! Gas phase global variables:
!`````````````````````````````````````````````````````````````````````//
! Species mass fractions.
      use fldvar, only: x_g
! Temperature.
      use fldvar, only: t_g
! Rate of production.
      use rxns, only: r_gp
! Rate of consumption divided by species mass fraction.
      use rxns, only: rox_gc
! Net rate of change of gas phase mass.
      use rxns, only: sum_r_g
! Heat of reaction source term
      use energy, only: hor_g

! Gas phase global variables:
!`````````````````````````````````````````````````````````````````````//
! Species mass fractions.
      use fldvar, only: x_s
! Temperatures.
      use fldvar, only: t_s
! Rate of production.
      use rxns, only: r_sp
! Rate of consumption divided by species mass fraction.
      use rxns, only: rox_sc
! Net rate of change of gas phase mass.
      use rxns, only: sum_r_s
! Heat of reaction source term
      use energy, only: hor_s


! Reaction global variables:
!`````````````````````````````````````````````````````````````````````//
! Number of TFM reactions.
      use rxns, only: no_of_rxns
! Rate of interphase mass transfer.
      use rxns, only: r_phase
! Data structure storing reaction information.
      use rxns, only: reaction


! Global run-time parameters.
!`````````````````````````````````````````````````````````````````````//
! Number of solids phases (or one if MMAX is zero)
      use param, only: dimension_m
! Number of gas species.
      use param, only: dimension_n_g
! Number of solids phase species.
      use param, only: dimension_n_s
      use param1, only: dimension_lm
      use param1, only: zero
! Indicates that the energy equations are solved.
      use run, only: energy_eq
! Indicates that the species equations are solved.
      use run, only: species_eq
! Simulation units (CGS or SI)
      use run, only: units
! Number of solids phases.
      use physprop, only: mmax
! Number of species comprising each phase.
      use physprop, only: nmax
! Indicate that DEM is used.
      use discretelement, only: discrete_element
! Indicate that TFM/DEM hybrid is used.
      use discretelement, only: des_continuum_hybrid
! Small value for species mass fractions
      use toleranc
      use functions

      implicit none

! Local variables:
!`````````````````````````````````````````````````````````````````````//
! Index/loop counters
      INTEGER :: IJK    ! fluid cell index
      INTEGER :: H      ! reaction loop counter
      INTEGER :: L, M   ! global phase index loop counters
      INTEGER :: NN      ! global species index
      INTEGER :: lN     ! reaction speices index
      INTEGER :: lM     ! reaction phase index
      INTEGER :: mXfr   ! global phase index for mass transfer

! This is a mask for the number of solids phases. The mask is used to
! avoid issues for DEM based reactions where MMAX and DES_MMAX may be
! in conflict.
      INTEGER :: lMMAX
      LOGICAL :: noSolids

! User-defined reaction rates returned from USR_RATES
      DOUBLE PRECISION RATES(no_of_rxns)
! Reaction and species specific reaction rate.
      DOUBLE PRECISION lRate

! Single reaction formation/consumption rate of each gas species.
      DOUBLE PRECISION :: r_Rgp(dimension_n_g)
      DOUBLE PRECISION :: r_ROXgc(dimension_n_g)
! Single reaction formation/consumption rate of each solids species.
      DOUBLE PRECISION :: r_Rsp(dimension_m, dimension_n_s)
      DOUBLE PRECISION :: r_ROXsc(dimension_m, dimension_n_s)

! Rate of interphase enthalpy transfer
      DOUBLE PRECISION :: r_RxH(0:mmax, 0:mmax)
! Local heat of reactions
      DOUBLE PRECISION :: r_HoRg, r_HoRs(1:mmax)
! Conversion factor for HoR (CGS to SI)
      DOUBLE PRECISION :: l_2SI

! Reaction limiters. If a species mass fraction is less than this
! value, then the reaction is suppressed.
      DOUBLE PRECISION :: speciesLimiter

! External functions:
!`````````````````````````````````````````````````````````````````````//
! Calculates enthalpy (cal/gram)
      DOUBLE PRECISION, EXTERNAL :: CALC_H

! Initialize global storage arrays to zero
!---------------------------------------------------------------------//
!     Gas                Solids
      r_gp    = zero;    r_sp    = zero
      rox_gc  = zero;    rox_sc  = zero
      sum_r_g = zero;    sum_r_s = zero
      hor_g   = zero;    hor_s   = zero

      r_phase = zero

! Set the species limiter:
      specieslimiter = zero_x_gs

! Set the conversion factor for heat of reaction.
      l_2si = merge(4.183925d3, one, units(1:2) == 'SI')

! Set the number of TFM solids phases.
! the following 2 lines could be replaced by assing lmmax=smax
      nosolids = discrete_element .AND. (.NOT.des_continuum_hybrid)
      lmmax = merge(0, mmax, nosolids)

! Loop over each fluid cell.
      DO ijk = ijkstart3, ijkend3
      IF (fluid_at(ijk)) THEN

      rates(:) = zero

! Calculate user defined reaction rates.
      CALL usr_rates(ijk, rates)

! Loop over reactions.
      rxn_lp: DO h = 1, no_of_rxns

! Skip empty reactions
         IF(reaction(h)%nSpecies == 0) cycle rxn_lp
         IF(compare(rates(h),zero)) cycle rxn_lp

! Initialize reaction loop arrays
         r_rgp = zero; r_roxgc = zero; r_horg = zero
         r_rsp = zero; r_roxsc = zero; r_hors = zero

         r_rxh = zero

! Calculate the rate of formation/consumption for each species.
!---------------------------------------------------------------------//
         DO ln = 1, reaction(h)%nSpecies
! Global phase index.
            m = reaction(h)%Species(ln)%pMap
! Global species index.
            nn = reaction(h)%Species(ln)%sMap
! Index for interphase mass transfer. For a gas/solid reaction, the
! index is stored with the gas phase. For solid/solid mass transfer
! the index is stored with the source phase.
            mxfr = reaction(h)%Species(ln)%mXfr
            lrate = rates(h) * reaction(h)%Species(ln)%MWxStoich
! Gas Phase:
            IF(m == 0) THEN
! Consumption of gas phase species.
               IF(lrate < zero) THEN
                  IF(x_g(ijk,nn) > specieslimiter) THEN
                     r_roxgc(nn) = r_roxgc(nn) - lrate/x_g(ijk,nn)
! Enthalpy transfer associated with mass transfer. (gas/solid)
                     IF(m /= mxfr) r_rxh(m,mxfr) = r_rxh(m,mxfr) +     &
                        lrate * calc_h(t_g(ijk),0,nn)
                  ELSE
! There is an insignificant amount of reactant. Skip this reaction.
                     rates(h) = zero
                     cycle rxn_lp
                  ENDIF
               ELSE
! Formation of gas phase species.
                  r_rgp(nn) = r_rgp(nn) + lrate
! Enthalpy transfer associated with mass transfer. (gas/solid)
                  IF(m /= mxfr) r_rxh(m,mxfr) = r_rxh(m,mxfr) +        &
                     lrate * calc_h(t_s(ijk,mxfr),0,nn)
               ENDIF


! Solids Phase M:
            ELSE
! Consumption of solids phase species.
               IF(lrate < zero) THEN
                  IF(x_s(ijk,m,nn) > specieslimiter) THEN
                     r_roxsc(m,nn) = r_roxsc(m,nn) - lrate/x_s(ijk,m,nn)
! Enthalpy transfer associated with mass transfer. (solid/solid) This
! is only calculated from the source (reactant) material.
                     IF(m /= mxfr) THEN
                        IF(m < mxfr) THEN
                           r_rxh(m,mxfr) =  r_rxh(m,mxfr) + lrate *    &
                              reaction(h)%Species(ln)%xXfr *           &
                              calc_h(t_s(ijk,m),m,nn)
                        ELSE
                           r_rxh(mxfr,m) =  r_rxh(mxfr,m) - lrate *    &
                              reaction(h)%Species(ln)%xXfr *           &
                              calc_h(t_s(ijk,m),m,nn)
                        ENDIF
                     ENDIF
                  ELSE
! There is an insignificant amount of reactant. Skip this reaction.
                     rates(h) = zero
                     cycle rxn_lp
                  ENDIF
               ELSE
! Formation of solids phase species.
                  r_rsp(m,nn) = r_rsp(m,nn) + lrate
               ENDIF
            ENDIF
         ENDDO ! Loop of species



! Map the local reaction production/consumption into global variables.
!---------------------------------------------------------------------//
! The outer reaction loop was cycled if there was insufficient amount
! of species available.

! Gas phase production
         r_gp(ijk,:nmax(0)) = r_gp(ijk,:nmax(0)) +                     &
            r_rgp(:nmax(0))
! Gas phase consumption divided by species mass fraction
         rox_gc(ijk,:nmax(0)) = rox_gc(ijk,:nmax(0)) +                 &
            r_roxgc(:nmax(0))

         DO m=1,lmmax
! Solids phase species production
            r_sp(ijk,m,:nmax(m)) = r_sp(ijk,m,:nmax(m)) +              &
               r_rsp(m,:nmax(m))
! Solids phase species consumption divided by species mass fraction
            rox_sc(ijk,m,:nmax(m)) = rox_sc(ijk,m,:nmax(m)) +          &
               r_roxsc(m,:nmax(m))
         ENDDO


! Calculate and store the heat of reaction.
!---------------------------------------------------------------------//
         IF(energy_eq) THEN
! Automated heat of reaction calculations
            IF(reaction(h)%Calc_DH) THEN
! Loop over reaction species.
               DO ln = 1, reaction(h)%nSpecies
! Global phase index.
                  m = reaction(h)%Species(ln)%pMap
! Global species index.
                  nn = reaction(h)%Species(ln)%sMap
! Rate of formation/consumption for speices N
                  lrate = rates(h) * reaction(h)%Species(ln)%MWxStoich
! Gas phase enthalpy change from energy equation derivation.
                  IF(m == 0) THEN
                     r_horg = r_horg + calc_h(t_g(ijk),0,nn) * lrate
! Solid phase enthalpy change from energy equation derivation.
                  ELSE
                     r_hors(m) = r_hors(m) + calc_h(t_s(ijk,m),m,nn)*lrate
                  ENDIF
               ENDDO

! Complete the skew-symettric for enthalpy transfer with mass transfer
               DO m=1, lmmax
                   DO l=0, m-1
                    r_rxh(m,l) = - r_rxh(l,m)
                  ENDDO
               ENDDO
! Apply enthalpy transfer associated with mass transfer to get the
! complete heat of reaction of heat phse for Reaction H.
               DO l=0, lmmax
                  DO m = 0, lmmax
                     IF(l == m) cycle
                     IF(l == 0) THEN
                        r_horg = r_horg - r_rxh(l,m)
                     ELSE
                        r_hors(l) = r_hors(l) - r_rxh(l,m)
                     ENDIF
                  ENDDO
               ENDDO
! Apply unit conversion and store heats of reaction in the global array.
               hor_g(ijk) = hor_g(ijk) + l_2si*r_horg
               DO m=1,lmmax
                  hor_s(ijk,m) = hor_s(ijk,m) + l_2si*r_hors(m)
               ENDDO
            ELSE
! User-defined heat of reaction.
               hor_g(ijk) = hor_g(ijk) + reaction(h)%HoR(0) * rates(h)
               DO m=1, lmmax
                  hor_s(ijk,m) = hor_s(ijk,m) +                        &
                     reaction(h)%HoR(m) * rates(h)
               ENDDO
            ENDIF
         ENDIF  ! ENERGY_EQ

! Update rate of interphase mass transfer.
!---------------------------------------------------------------------//
          DO lm=1, (dimension_lm+dimension_m-1)
             r_phase(ijk,lm) = r_phase(ijk,lm) +                       &
                rates(h) * reaction(h)%rPHASE(lm)
          ENDDO
      ENDDO rxn_lp ! Loop over reactions.


! Calculate the toal rate of formation and consumption for each species.
!---------------------------------------------------------------------//
      IF(species_eq(0)) THEN
         sum_r_g(ijk) = sum( r_gp(ijk,:nmax(0)) -                      &
            rox_gc(ijk,:nmax(0))*x_g(ijk,:nmax(0)))
      ELSE
         DO h=1, no_of_rxns
            IF(reaction(h)%nPhases <= 0) cycle
            DO m=1, lmmax
               lm = 1 + ((m-1)*m)/2
               sum_r_g(ijk) = sum_r_g(ijk) +                           &
                  rates(h) * reaction(h)%rPHASE(lm)
            ENDDO
         ENDDO
      ENDIF

      DO m=1, lmmax
         IF(species_eq(m)) THEN
            sum_r_s(ijk,m) = sum( r_sp(ijk,m,:nmax(m)) -               &
               rox_sc(ijk,m,:nmax(m))*x_s(ijk,m,:nmax(m)))
         ELSE
            DO h=1, no_of_rxns
               IF(reaction(h)%nPhases <= 0) cycle
               DO l=0, m-1
                  lm = 1 + l + ((m-1)*m)/2
                  sum_r_s(ijk,m) = sum_r_s(ijk,m) -                    &
                     rates(h) * reaction(h)%rPHASE(lm)
               ENDDO
               DO l=m+1, lmmax
                  lm = 1 + m + ((l-1)*l)/2
                  sum_r_s(ijk,m) = sum_r_s(ijk,m) +                    &
                     rates(h) * reaction(h)%rPHASE(lm)
               ENDDO
            ENDDO
         ENDIF
      ENDDO


      ENDIF  ! Fluid_At(IJK)
      END DO ! IJK

      RETURN

      END SUBROUTINE rrates0