rxns_gs_des1.f

Go to the documentation of this file.

#include "version.inc"
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: CONV_GS_DES1                                           !
!  Author: J.Musser: 16-Jun-10                                         !
!                                                                      !
!  Purpose: This routine is called from the DISCRETE side to calculate !
!  the gas-particle convective heat transfer.                          !
!                                                                      !
!  Comments: Explicitly coupled simulations use a stored convective    !
!  heat transfer coefficient. Otherwise, the convective heat transfer  !
!  coeff is calculated every time step and the total interphase energy !
!  transfered is 'stored' and used explictly in the gas phase. The     !
!  latter conserves all energy
!                                                                      !
!  REF: Zhou, Yu, and Zulli, "Particle scale study of heat transfer in !
!       packed and bubbling fluidized beds," AIChE Journal, Vol. 55,   !
!       no 4, pp 868-884, 2009.                                        !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE rxns_gs_des1

      use constant, only: pi
! Flag: The fluid and discrete solids are explicitly coupled.
      use discretelement, only: des_explicitly_coupled
      use discretelement, only: dtsolid
      use discretelement, only: pijk
      use discretelement, only: max_pip

      use particle_filter, only: des_interp_on

      use particle_filter, only: filter_cell
      use particle_filter, only: filter_weight
      use particle_filter, only: filter_size

      use des_rxns, only: des_r_gp, des_r_gc
      use des_rxns, only: des_r_phase, des_sum_r_g
      use des_rxns, only: des_hor_g

      use physprop, only: nmax

      use functions, only: fluid_at
      use functions, only: is_normal

      use param1, only: zero
      use des_rxns, only: des_r_s
      use des_thermo, only: rxns_qs

      use param1, only: dimension_lm
      use param, only: dimension_m

      IMPLICIT NONE

      INTEGER :: IJK, LC, NP


! Local gas phase values.
      DOUBLE PRECISION :: lRgp(nmax(0)) ! Rate of species production
      DOUBLE PRECISION :: lRgc(nmax(0)) ! Rate of species consumption
      DOUBLE PRECISION :: lHoRg         ! Heat of reaction
      DOUBLE PRECISION :: lSUMRg        ! lSUMRg

! Interphase mass transfer
      DOUBLE PRECISION :: lRPhase(dimension_lm+dimension_m-1)

      DOUBLE PRECISION :: WEIGHT

      DO np=1,max_pip
         IF(.NOT.is_normal(np)) cycle

! Avoid convection calculations in cells without fluid (cut-cell)
         IF(.NOT.fluid_at(pijk(np,4))) THEN
            des_r_s(np,:) = zero
            rxns_qs(np) = zero

! No additional calculations are needed for explictly coupled
         ELSEIF(.NOT.des_explicitly_coupled) THEN

! Calculate the heat transfer coefficient.
            CALL calc_rrates_des(np, lrgp, lrgc, lrphase, lhorg, lsumrg)

! Integrate over solids time step and store in global array. This
! needs updated when interpolation is reintroduced into thermo code.
!---------------------------------------------------------------------//
! Store the gas phase source terms.
            IF(des_interp_on) THEN
               DO lc=1,filter_size
                  ijk = filter_cell(lc,np)
                  weight = dtsolid*filter_weight(lc,np)

                  des_r_gp(ijk,:) = des_r_gp(ijk,:)+lrgp*weight
                  des_r_gc(ijk,:) = des_r_gc(ijk,:)+lrgc*weight
                  des_r_phase(ijk,:) = des_r_phase(ijk,:)+lrphase*weight
                  des_hor_g(ijk) = des_hor_g(ijk) + lhorg*weight
                  des_sum_r_g(ijk) = des_sum_r_g(ijk) + lsumrg*weight
               ENDDO
            ELSE
               ijk = pijk(np,4)
               weight = dtsolid

               des_r_gp(ijk,:) = des_r_gp(ijk,:) + lrgp*weight
               des_r_gc(ijk,:) = des_r_gc(ijk,:) + lrgc*weight
               des_r_phase(ijk,:) = des_r_phase(ijk,:) + lrphase*weight
               des_hor_g(ijk) = des_hor_g(ijk) + lhorg*weight
               des_sum_r_g(ijk) = des_sum_r_g(ijk) + lsumrg*weight
            ENDIF
         ENDIF
      ENDDO

! Note that MPI sync is managed at the end of des_time_march for
! non-explicitly coupled cases that use interpolation.

      RETURN
      END SUBROUTINE rxns_gs_des1


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: RXNS_GS_GAS1                                            !
!  Author: J.Musser                                   Date: 21-NOV-14  !
!                                                                      !
!  Purpose: This routine is called from the CONTINUUM. It calculates   !
!  the scalar cell center drag force acting on the fluid using         !
!  interpolated values for the gas velocity and volume fraction. The   !
!  The resulting sources are interpolated back to the fluid grid.      !
!                                                                      !
!  NOTE: The loop over particles includes ghost particles so that MPI  !
!  communications are needed to distribute overlapping force between   !
!  neighboring grid cells. This is possible because only cells "owned" !
!  by the current process will have non-zero weights.                  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE rxns_gs_gas1

! Size of particle array on this process.
      use discretelement, only: max_pip
! Flag to use interpolation
      use particle_filter, only: des_interp_on
! Interpolation cells and weights
      use particle_filter, only: filter_cell,filter_weight,filter_size
! IJK of fluid cell containing particles center
      use discretelement, only: pijk
! Gas phase mass, species, and energy equation sources
      use des_rxns, only: des_r_gp, des_r_gc, des_sum_r_g
      use des_rxns, only: des_r_phase, des_hor_g
! Number of species for each phase
      use physprop, only: nmax
! Funtion for identifying fluid cells and normal particles.
      use functions, only: fluid_at
      use functions, only: is_normal
! MPI function for collecting interpolated data from ghost cells.
      use sendrecvnode, only: des_collect_gdata
! MPI wrapper for halo exchange.
      use sendrecv, only: send_recv
! Fluid time step size.
      use run, only: dt
! Flag to use stiff chemistry solver
      use stiff_chem, only: stiff_chemistry
! Routine to sync stiff solver results
      use stiff_chem, only: finalize_stiff_solver

! Global Parameters:
!---------------------------------------------------------------------//
! Double precision values.
      use param1, only: zero, one
      use param1, only: dimension_lm
      use param, only: dimension_m

      IMPLICIT NONE

! Loop counters: Particle, fluid cell, neighbor cells
      INTEGER :: NP, IJK, LC
! Loop bound for filter
      DOUBLE PRECISION :: GAMMAxTp
! Local gas phase values.
      DOUBLE PRECISION :: lRgp(nmax(0)) ! Rate of species production
      DOUBLE PRECISION :: lRgc(nmax(0)) ! Rate of species consumption
      DOUBLE PRECISION :: lHoRg         ! Heat of reaction
      DOUBLE PRECISION :: lSUMRg        ! lSUMRg

! Interphase mass transfer
      DOUBLE PRECISION :: lRPhase(dimension_lm+dimension_m-1)

      DOUBLE PRECISION :: WEIGHT

! Initialize fluid cell values.
      des_r_gp = zero
      des_r_gc = zero
      des_r_phase = zero
      des_hor_g = zero
      des_sum_r_g = zero

! Calculate the gas phase forces acting on each particle.

      DO np=1,max_pip

! Only calculate chemical reactions for normal particles that are inside
! fluid cells. Ex: Skip ghost particles or particles that are in a cut-
! cell dead space.
         IF(.NOT.is_normal(np)) cycle
         IF(.NOT.fluid_at(pijk(np,4))) cycle

         ijk = pijk(np,4)

! Calculate the rates of species formation/consumption.
         CALL calc_rrates_des(np, lrgp, lrgc, lrphase, lhorg, lsumrg)

! Directly integrate fluid cell reaction sources
         IF(stiff_chemistry) THEN
            CALL des_stiff_chem(ijk, dt, lrgp, lrgc, lhorg, lsumrg)

! Store the gas phase source terms.
         ELSEIF(des_interp_on) THEN
            DO lc=1,filter_size
               ijk = filter_cell(lc,np)
               weight = filter_weight(lc,np)

               des_r_gp(ijk,:) = des_r_gp(ijk,:)+lrgp*weight
               des_r_gc(ijk,:) = des_r_gc(ijk,:)+lrgc*weight
               des_r_phase(ijk,:) = des_r_phase(ijk,:)+lrphase*weight
               des_hor_g(ijk) = des_hor_g(ijk) + lhorg*weight
               des_sum_r_g(ijk) = des_sum_r_g(ijk) + lsumrg*weight
            ENDDO
         ELSE
            des_r_gp(ijk,:) = des_r_gp(ijk,:) + lrgp
            des_r_gc(ijk,:) = des_r_gc(ijk,:) + lrgc
            des_r_phase(ijk,:) = des_r_phase(ijk,:) + lrphase
            des_hor_g(ijk) = des_hor_g(ijk) + lhorg
            des_sum_r_g(ijk) = des_sum_r_g(ijk) + lsumrg
         ENDIF
      ENDDO

! Make the necessary send/recv calls for field vars
      IF(stiff_chemistry) THEN
         CALL finalize_stiff_solver
      ELSE
! Add in data stored in ghost cells from interpolation. This call must
! preceed the SEND_RECV to avoid overwriting ghost cell data.
         IF(des_interp_on) THEN
            CALL des_collect_gdata(des_r_gp)
            CALL des_collect_gdata(des_r_gc)
            CALL des_collect_gdata(des_r_phase)
            CALL des_collect_gdata(des_hor_g)
            CALL des_collect_gdata(des_sum_r_g)
         ENDIF

! Update the species mass sources in ghost layers.
         CALL send_recv(des_r_gp, 2)
         CALL send_recv(des_r_gc, 2)
         CALL send_recv(des_r_phase, 2)
         CALL send_recv(des_hor_g, 2)
         CALL send_recv(des_sum_r_g, 2)
      ENDIF

      RETURN
      END SUBROUTINE rxns_gs_gas1


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DES_STIFF_CHEM                                          !
!  Author: J.Musser                                   Date:  7-FEB-16  !
!                                                                      !
!  Purpose: This routine updates a fluid cell by directly integrating  !
!  the source terms due to reaction. Note that this routine is called  !
!  for each particle to prevent over consumption of reactants.         !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE des_stiff_chem(pIJK, pDT, pRgp, pRgc, pHoRg, pSUMRg)

      USE constant, only: gas_const
      USE fldvar, only: ep_g, ro_g, rop_g, t_g, x_g, p_g
      USE geometry, only: vol
      USE param1, only: zero, small_number, one
      USE physprop, only: c_pg, mw_g, mw_mix_g
      use toleranc, only: zero_x_gs
      use physprop, only: nmax

      use compar, only: mype
      IMPLICIT NONE

! Passed variables
!---------------------------------------------------------------------//
      INTEGER, INTENT(IN) :: pIJK   ! fluid cell index
      DOUBLE PRECISION, INTENT(IN) :: pDT           ! fluid time step
      DOUBLE PRECISION, INTENT(IN) :: pRgp(nmax(0)) ! Rate of production
      DOUBLE PRECISION, INTENT(IN) :: pRgc(nmax(0)) ! Rate of consumption
      DOUBLE PRECISION, INTENT(IN) :: pHoRg         ! Heat of reaction
      DOUBLE PRECISION, INTENT(IN) :: pSUMRg        ! Net mass change

! Local variables
!---------------------------------------------------------------------//
      DOUBLE PRECISION :: lRg(nmax(0)), sumXg
      DOUBLE PRECISION :: lDToVOL
      INTEGER :: N
!......................................................................!

! Time step divided by cell volume.
      ldtovol = pdt/vol(pijk)
! Net change in species mass
      lrg = prgp - prgc

! Gas phase density: ROP_g
      rop_g(pijk) = rop_g(pijk) + ldtovol*psumrg

! Temperature: T_g
      t_g(pijk) = t_g(pijk) - ldtovol*(phorg/(rop_g(pijk)*c_pg(pijk)))

! Species mass fractions: X_g
      DO n=1,nmax(0)
         x_g(pijk,n) = x_g(pijk,n) + ldtovol * &
            (lrg(n) - x_g(pijk,n)*psumrg)/(rop_g(pijk))
      ENDDO

! Cleanup over/under shoots.
      x_g(pijk,:) = min(max(x_g(pijk,:), zero), one)
      sumxg = sum(x_g(pijk,:))
      x_g(pijk,:) = x_g(pijk,:)/sumxg

! Gas phase bulk density is updated within the stiff solver (lVar(1)).
! Now that the gas phase volume fraction is updated, the gas phase
! density can be backed out. RO_g * EP_g = ROP_g
      IF(ep_g(pijk) > small_number) THEN
         ro_g(pijk) = rop_g(pijk) / ep_g(pijk)
      ELSE
! This case shouldn't happen, however HUGE is used to aid in tracking
! errors should this somehow become and issue.
         ro_g(pijk) = huge(0.0)
      ENDIF

! Calculate the mixture molecular weight.
      mw_mix_g(pijk) = sum(x_g(pijk,1:nmax(0))/mw_g(1:nmax(0)))
      mw_mix_g(pijk) = one/mw_mix_g(pijk)

! Calculate the gas phase pressure.
      p_g(pijk) = (ro_g(pijk)*gas_const*t_g(pijk))/mw_mix_g(pijk)

      RETURN
      END SUBROUTINE des_stiff_chem