stiff_chem_rrates.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: STIFF_CHEM_RRATES(IER)                                  C
!                                                                      C
!  Purpose: Calculate reaction rates for various reactions in cell ijk C
!           using information from the data file                       C
!                                                                      C
!  Author: J. Musser                                  Date: 12-Feb-13  C
!  Reviewer:                                          Date: dd-mmm-yy  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE stiff_chem_rrates(lNEQ, lTime, Y, YDOT)

! External Module Procedures:
!---------------------------------------------------------------------//
      use stiff_chem_maps, only: mapodetomfix

! Global Parameters:
!---------------------------------------------------------------------//
! Maximum number of solids phases
      use param, only: dimension_m
! Maximum number of gas species
      use param, only: dimension_n_g
! Maximum number of solid species
      use param, only: dimension_n_s

      use param1, only: zero

! Global Variables:
!---------------------------------------------------------------------//
! Number of solids phases
      use physprop, only: mmax
! Number of species
      use physprop, only: nmax

! Gas Phase
!--------------------------------------------//
!Temperature
      use fldvar, only: t_g
! Species Mass Fractions
      use fldvar, only: x_g
! Gas Phase apparent density
      use fldvar, only: rop_g
! Gas phase specific heat
      use physprop, only: c_pg

! Solids Phases
!--------------------------------------------//
! Temperatures
      use fldvar, only: t_s
! Species mass fractions
      use fldvar, only: x_s
! Apparent densities
      use fldvar, only: rop_s
! Specific heats
      use physprop, only: c_ps

      use rxns, only: reaction
      use rxns, only: no_of_rxns

      use run, only: units

      use stiff_chem, only: ode_dimn_all
      use stiff_chem, only: neq_dimn

      implicit none

! Passed Variables: Dummy argument format required by ODEPACK.
!---------------------------------------------------------------------//
! (1) Number of ODEs to be solve
! (2) Fluid cell index
      INTEGER, intent(in) :: lNEQ(neq_dimn)
! Independent variable (not used)
      DOUBLE PRECISION, intent(in) :: lTime
! Array of dependent variable initial values.
      DOUBLE PRECISION, dimension(ODE_DIMN_all), intent(in)  :: Y
! Rate of change of dependent variables.
      DOUBLE PRECISION, dimension(ODE_DIMN_all), intent(out) :: YDOT

! Local Variables:
!---------------------------------------------------------------------//
      INTEGER :: IJK  ! Fluid Cell Index
      INTEGER :: H    ! Reaction loop counter
      INTEGER :: L, M ! Global Phase index loop counters
      INTEGER :: N    ! Global species index
      INTEGER :: lN   ! Local reaction speices index/loop counter
      INTEGER :: mXfr ! Global phase index for mass transfer
      INTEGER :: Node ! Variable index

! User-defined reaction rates returned from USR_RATES
      DOUBLE PRECISION :: RATES(no_of_rxns)
! Rate of formation (+) or consumption (-) of a species (GENERIC).
      DOUBLE PRECISION :: lRate
! Rate of formation (+) or consumption (-) of each gas phase species.
! > Accumulative over all reactions.
      DOUBLE PRECISION :: lRg(dimension_n_g)
! Rate of formation (+) or consumption (-) of each gas phase species.
! > Single reaction.
      DOUBLE PRECISION :: rRg(dimension_n_g)
! Heat of reaction assigned to gas phase for a reaction.
      DOUBLE PRECISION :: rHORg
! Cummulative heat of reaction assigned to gas phase.
      DOUBLE PRECISION :: lHORg
! Rate of formation (+) or consumption (-) of each solids phase species.
! > Accumulative over all reactions.
      DOUBLE PRECISION :: lRs(dimension_m, dimension_n_s)
! Rate of formation (+) or consumption (-) of each solids phase species.
! > Single reaction.
      DOUBLE PRECISION :: rRs(dimension_m, dimension_n_s)
! Heat of reaction assigned to the solids phase for a reaction.
      DOUBLE PRECISION :: rHORs(dimension_m)
! Cummulative heat of reaction assigned to the solids phase.
      DOUBLE PRECISION :: lHORs(dimension_m)
! Rate of interphase enthalpy transfer due to mass transfer.
      DOUBLE PRECISION :: RxH(0:dimension_m, 0:dimension_m)

! Net rate of change of gas phase material.
      DOUBLE PRECISION :: sumlRg
! Net rate of change of solids phase material.
      DOUBLE PRECISION :: sumlRs(dimension_m)

! Reaction limiters.
      DOUBLE PRECISION, parameter :: speciesLimiter = 1.0d-7

! External Function for comparing two numbers.
      LOGICAL, external :: COMPARE
      DOUBLE PRECISION, external :: CALC_H

! UDF for Reaction Rates:
      external usr_rates

! Initialize variables:
      ijk    = lneq(2)
      rates  = zero
      ydot   = zero

      lrg = zero; lhorg = zero
      lrs = zero; lhors = zero

! Map the current ODE independent variables to MFIX variables.
      CALL mapodetomfix(neq_dimn, lneq, ode_dimn_all, y)

! 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

! Initialize local loop arrays
         rrg = zero; rhorg = zero
         rrs = zero; rhors = zero
         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.
            n = 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
! Check that there is a sufficient amount of reactant.
                  IF(x_g(ijk,n) > specieslimiter) THEN
                     rrg(n) = rrg(n) + lrate
! Enthalpy transfer associated with mass transfer. (gas/solid)
                     IF(m /= mxfr) rxh(m,mxfr) =  rxh(m,mxfr) + &
                        lrate * calc_h(t_g(ijk),0,n)
                  ELSE
! There is an insignificant amount of reactant. Skip this reaction.
                     cycle rxn_lp
                  ENDIF
               ELSE
! Formation of gas phase species.
                  rrg(n) = rrg(n) + lrate
! Enthalpy transfer associated with mass transfer. (gas/solid)
                  IF(m /= mxfr) rxh(m,mxfr) =  rxh(m,mxfr) + &
                     lrate * calc_h(t_s(ijk,mxfr),0,n)
               ENDIF
! Solids Phase M:
            ELSE
! Consumption of solids phase species.
               IF(lrate < zero) THEN
                  IF(x_s(ijk,m,n) > specieslimiter) THEN
                     rrs(m,n) = rrs(m,n) + lrate
! 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
                           rxh(m,mxfr) =  rxh(m,mxfr) +                &
                              lrate * calc_h(t_s(ijk,m),m,n) *         &
                              reaction(h)%Species(ln)%xXfr
                        ELSE
                           rxh(mxfr,m) =  rxh(mxfr,m) -                &
                             lrate * calc_h(t_s(ijk,m),m,n) *          &
                             reaction(h)%Species(ln)%xXfr
                        ENDIF
                     ENDIF
                  ELSE
! There is an insignificant amount of reactant. Skip this reaction.
                     cycle rxn_lp
                  ENDIF
               ELSE
! Formation of solids phase species.
                  rrs(m,n) = rrs(m,n) + lrate
               ENDIF
            ENDIF
         ENDDO ! Loop of species

! Copy the single reaction rates of formation/consumption to the total
! (accumulative) rates of formation/consumption arrays. This ensures
! that any reactions without sufficient reactants are not included.
         lrg = lrg + rrg
         lrs = lrs + rrs

! Calculate and store the heat of reaction.
!---------------------------------------------------------------------//
! 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.
               n = reaction(h)%Species(ln)%sMap
! Rate of formation/consumption for speices N
               lrate = rates(h) * reaction(h)%Species(ln)%MWxStoich
! Gas phase enthalpy chnage from energy equation derivation.
               IF(m == 0) THEN
                  rhorg = rhorg + lrate * calc_h(t_g(ijk),0,n)
! Solid phase enthalpy change from energy equation derivation.
               ELSE
                  rhors(m) = rhors(m) + lrate * calc_h(t_s(ijk,m),m,n)
               ENDIF
            ENDDO

! Complete the skew-symettric for enthalpy transfer with mass transfer
            DO m=1, mmax
               DO l=0, m-1
                  rxh(m,l) = - 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, mmax
               DO m = 0, mmax
                  IF(l == m) cycle
                  IF(l == 0) THEN
                     rhorg = rhorg - rxh(l,m)
                  ELSE
                     rhors(l) = rhors(l) - rxh(l,m)
                  ENDIF
               ENDDO
            ENDDO

! Convert the heat of reaction to the appropriate units (if SI), and
! store in the global array.
            IF(units == 'SI') THEN
               lhorg = lhorg + 4.183925d3*rhorg
               DO m=1,mmax
                  lhors(m) = lhors(m) + 4.183925d3*rhors(m)
               ENDDO
            ELSE
               lhorg = lhorg + rhorg
               DO m=1,mmax
                  lhors(m) = lhors(m) + rhors(m)
               ENDDO
            ENDIF
         ELSE
! User-defined heat of reaction.
            lhorg = lhorg + reaction(h)%HoR(0) * rates(h)
            DO m=1, mmax
               lhors(m) = lhors(m) + reaction(h)%HoR(m) * rates(h)
            ENDDO
         ENDIF

      ENDDO rxn_lp ! Loop over reactions.




! Calculate and store the heat of reaction.
!---------------------------------------------------------------------//

! Calculate the net change for the gas phase.
      sumlrg = sum(lrg)

! Calculate the net change for solids phases.
      sumlrs = 0.0d0
      DO m=1, mmax
         IF(lneq(2+m) == 1) sumlrs(m) = sum(lrs(m,:))
      ENDDO


! Calculate and store the heat of reaction.
!---------------------------------------------------------------------//
! Initialize counter.
      node = 1

! Density:  ROP_g
      ydot(node) = sumlrg
      node = node + 1

! Temperature: T_g
      ydot(node) = -lhorg/(rop_g(ijk)*c_pg(ijk))
      node = node + 1

! Species mass fractions: X_g
      DO n=1,nmax(0)
         ydot(node) = (lrg(n) - x_g(ijk,n)*sumlrg)/(rop_g(ijk))
         node = node + 1
      ENDDO


! Temperature: T_s
      DO m = 1, mmax
         IF(rop_s(ijk,m) > 1.0d-8) THEN
            ydot(node) = -lhors(m)/(rop_s(ijk,m)*c_ps(ijk,m))
         ELSE
            ydot(node) = 0.0d0
         ENDIF
         node = node + 1
      ENDDO


      DO m = 1, mmax
         IF(lneq(2+m) == 1) THEN
! Solids bulk density: ROP_s
            ydot(node) = sumlrs(m)
            node = node + 1
! Species Mass Fraction: X_s
            DO n=1, nmax(m)
               ydot(node) = (lrs(m,n) - x_s(ijk,m,n)*sumlrs(m)) /     &
                  rop_s(ijk,m)
               node = node + 1
            ENDDO
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE stiff_chem_rrates