!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv! ! ! ! Module name: USR_RATES ! ! ! ! Purpose: ! ! ! ! Author: M.Denison Date: 26-Oct-18 ! ! ! ! ! !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^! SUBROUTINE USR_RATES(IJK, RATES) USE compar USE constant USE energy USE fldvar USE fun_avg USE functions USE funits USE geometry USE indices USE parallel USE param USE param1 USE physprop USE run USE rxns USE sendrecv USE toleranc ! USE usr IMPLICIT NONE INTEGER, INTENT(IN) :: IJK INTEGER ir DOUBLE PRECISION, DIMENSION(NO_OF_RXNS), INTENT(OUT) :: RATES ! Reaction specific variables: !`````````````````````````````````````````````````````````````````````// ! Bounded phase temperatures (K) DOUBLE PRECISION xTg ! gas phase DOUBLE PRECISION xTs ! solids phase ! Gas phase pressure in atm ! DOUBLE PRECISION Pg_KPa DOUBLE PRECISION Pg_atm ! Gas phase concentrations (mol/m^3) DOUBLE PRECISION c_CO, c_O2, c_CO2, c_H2, c_H2O, c_CH4 ! Partial pressures (atm) DOUBLE PRECISION p_O2, p_H2O, p_H2, p_CO2, p_CO, p_CH4 ! Total solids surface area per volume (1/m) DOUBLE PRECISION Sa ! Diffusion coefficient for SiH2 (m^2/sec) ! DOUBLE PRECISION DIFF_SiH2 ! SiH2 gas constant (KPa.m^3)/(kg-SiH2.K) ! DOUBLE PRECISION R_SiH2 ! Calculate film diffusion resistance. (kg-SiH2/(KPa.sec.m^2) ! DOUBLE PRECISION K_f ! Apparent first order heterogeneous reation rate constant (m/sec) ! DOUBLE PRECISION k_so ! Constant of hydroggen inhibition of heterogeneous reaction (1/KPa) ! DOUBLE PRECISION K_H ! Rate constants DOUBLE PRECISION k13, k14, k15, k16, kox, kdry, kdevol, k1 DOUBLE PRECISION k6, k7, k8, k9, k10, k11 ! fraction of Carbon converted to CO2 DOUBLE PRECISION psi ! Temporary values for reversible reactions DOUBLE PRECISION FWD, RVS, NET ! INTEGER M ! Prevent reactions when mass fraction falls below c_Limiter. DOUBLE PRECISION, parameter :: c_Limiter = 1.0d-7 !`````````````````````````````````````````````````````````````````````// INCLUDE 'species.inc' !`````````````````````````````````````````````````````````````````````// ! Calculate the bounded temperatures. xTg = max(min(T_g(IJK), TMAX), 300.0d0) xTs = max(min(T_s(IJK,1), TMAX), 300.0d0) ! Convert gas pressure to atm ! Pg_KPa = 1.0d-3 * P_g(IJK) ! KPa Pg_atm = P_g(IJK) / 101325.0 ! atm ! Initialize partial pressures and concentrations. c_CO = ZERO; c_O2 = ZERO; c_CO2 = ZERO; c_H2 = ZERO; c_H2O = ZERO; c_CH4=ZERO; ! Calculate the partial pressure and molar concentration of CO IF(X_g(IJK,CO) > c_Limiter) then p_CO = Pg_atm * X_g(IJK,CO) * (MW_MIX_g(IJK) / MW_g(CO)) c_CO = RO_g(IJK) * X_g(IJK,CO) / Mw_g(CO) ENDIF ! Calculate the partial pressure and molar concentration of O2 IF(X_g(IJK,O2) > c_Limiter) then p_O2 = Pg_atm * X_g(IJK,O2) * (MW_MIX_g(IJK) / MW_g(O2)) c_O2 = RO_g(IJK) * X_g(IJK,O2) / Mw_g(O2) ENDIF ! Calculate the partial pressure and molar concentration of CO2 IF(X_g(IJK,CO2) > c_Limiter) then p_CO2 = Pg_atm * X_g(IJK,CO2) * (MW_MIX_g(IJK) / MW_g(CO2)) c_CO2 = RO_g(IJK) * X_g(IJK,CO2) / Mw_g(CO2) ENDIF ! Calculate the partial pressure and molar concentration of H2 IF(X_g(IJK,H2_REF_ELEMENT) > c_Limiter) then p_H2 = Pg_atm * X_g(IJK,H2_REF_ELEMENT) * (MW_MIX_g(IJK) / MW_g(H2_REF_ELEMENT)) c_H2 = RO_g(IJK)*X_g(IJK,H2_REF_ELEMENT)/Mw_g(H2_REF_ELEMENT) ENDIF ! Calculate the partial pressure and molar concentration of H2O IF(X_g(IJK,H2O) > c_Limiter) then p_H2O = Pg_atm * X_g(IJK,H2O) * (MW_MIX_g(IJK) / MW_g(H2O)) c_H2O = RO_g(IJK) * X_g(IJK,H2O) / Mw_g(H2O) ENDIF ! Calculate the partial pressure and molar concentration of CH4 IF(X_g(IJK,CH4) > c_Limiter) then p_CH4 = Pg_atm * X_g(IJK,CH4) * (MW_MIX_g(IJK) / MW_g(CH4)) c_CH4 = RO_g(IJK) * X_g(IJK,CH4) / Mw_g(CH4) ENDIF !**********************************************************************! ! Homogeneous Reactions ! !**********************************************************************! ! CO + 1/2 O2 --> CO2 (kg-mole/m^3.s) !---------------------------------------------------------------------// IF(xTg < 1150.0) then FWD = 2.61d12 * exp(-2.2929d4/xTg) * EP_g(IJK) * c_CO * c_H2O**0.5 * c_O2**0.25 ELSE FWD = 6.52d6 * exp(-8.0352d3/xTg) * EP_g(IJK) * c_CO * c_H2O**0.5 * c_O2**0.25 ENDIF RATES(CO_oxidation) = FWD ! CO + H2O <--> CO2 + H2 (kg-mole/m^3.s) !---------------------------------------------------------------------// FWD = 7.68d10 * exp(-3.664d4/xTg) * EP_g(IJK) * sqrt(c_CO)*c_H2O RVS = 6.4d9 * exp(-3.926d4/xTg) * EP_g(IJK) * sqrt(c_H2) * c_CO2 NET = FWD - RVS IF(NET >= ZERO) THEN RATES(WGS_forward) = NET RATES(WGS_reverse) = ZERO ELSE RATES(WGS_forward) = ZERO RATES(WGS_reverse) = -NET ENDIF ! H2 + 1/2 O2 --> H2O (kg-mole/m^3.s) !---------------------------------------------------------------------// FWD = 1.8d13 * exp(-1.7614d4/xTg) * EP_g(IJK) * c_H2 * sqrt(c_O2) RATES(H2_oxidation) = FWD ! CH4 + 2 O2 --> CO2 + 2 H2O (kg-mole/m^3.s) !---------------------------------------------------------------------// ! FWD = 1.3d8 * exp(-2.4367d4/xTg) * EP_g(IJK) * c_O2**1.3 * c_CH4**(-0.3) FWD = 1.3d8 * exp(-2.4367d4/xTg) * EP_g(IJK) * c_O2**1.3 * c_CH4**(0.3) ! print*,EP_g(IJK),c_O2,c_CH4 RATES(CH4_oxidation) = FWD ! C_g + 1/2 O2 --> CO2 (kg-mole/m^3.s) !---------------------------------------------------------------------// FWD = 1.0d10 * exp(-1.7614d4/xTg) * EP_g(IJK) * c_O2 RATES(C_oxidation) = FWD !**********************************************************************! ! Heterogeneous Reactions ! !**********************************************************************! IF(EP_s(IJK,1) > ZERO_EP_s) THEN ! Cu2O + 1/2 O2 --> 2 CuO (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(Cu2Ooxidation)%Species(CU2O)%pmap k13 = 6.3145d13 * exp(-4.1497d4/xTg) * EP_g(IJK) * ROP_s(IJK,1) * X_s(IJK,1,CU2O) k14 = 5.6d3 * exp(-8179.0/xTg) *EP_g(IJK) *ROP_s(IJK,1) * X_s(IJK,1,CU2O) ! Reaction rate: RATES(Cu2Ooxidation) = (k13 - k14 * p_O2**1.3)/MW_s(1,CU2O) ! CuO --> 1/2 Cu2O+ 1/4 O2 (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(CuOdecomp)%Species(CUO)%pmap k15 = 3.0 * exp(-2350.0/xTg) * EP_g(IJK) * ROP_s(IJK,1) * X_s(IJK,1,CUO) k16 = 1.62d8 * exp(-2.795d4/xTg) *EP_g(IJK) *ROP_s(IJK,1) * X_s(IJK,1,CUO) ! Reaction rate: RATES(CuOdecomp) = (k15 * p_O2 - k16)/MW_s(1,CUO) ENDIF IF(EP_s(IJK,2) > ZERO_EP_s) THEN ! Drying H2O_L_1 --> H2O (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(Drying)%Species(H2O_L_1)%pmap kdry = 1.3d7 * exp(-20133.0/xTg) * ROP_s(IJK,2) * X_s(IJK,2,H2O_L_1) ! Reaction rate: RATES(Drying) = kdry/MW_s(2,H2O_L_1) ! Coal --> Volatiles (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(VolatileRelease)%Species(Coal)%pmap kdevol = 4.28d14 * exp(-27477.3/xTg) * ROP_s(IJK,2) * X_s(IJK,2,Coal) ! Reaction rate: RATES(VolatileRelease) = 0.536068 * kdevol / MW_s(2,Coal) ! Coal --> Char (kg-mole/m^3.s) !---------------------------------------------------------------------// ! Reaction rate: RATES(char_production) = 0.463932 * kdevol / MW_s(2,Coal) ! CGR_REF_ELEMENT + O2 --> CO2 (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(Char_C_oxidation)%Species(CGR_REF_ELEMENT)%pmap ! Total surface area of solids in computational cell (m^2/m^3) Sa = 6.0d0 * EP_s(IJK,2) / D_p0(2) ! (1/m) ! Sa = 6.0d0 * EP_s(IJK,M) / 0.0001 ! (1/m) kox = 1.45d6 * exp(-10049.3/xTg) * EP_g(IJK) * Sa / MW_s(2,Coal) ! fraction of carbon converted to CO2 k1 = 200.0 * exp(-4529.4/xTg) psi = 1.0 / (1.0 + k1) ! Reaction rate: RATES(Char_C_oxidation) = 0.51926 * psi * kox * p_O2 ! print*,RATES(Char_C_oxidation),psi,kox,p_O2,Sa,EP_s(IJK,2),MW_s(2,Coal) ! CGR_REF_ELEMENT + 1/2 O2 --> CO (kg-mole/m^3.s) !---------------------------------------------------------------------// ! Reaction rate: RATES(Char_to_CO_oxidation) = 0.51926 * (1.0 - psi) * kox * p_O2 ! Other char oxidation RATES(Char_other_oxidation) = 0.48074 * kox * p_O2 ! CGR_REF_ELEMENT + H2O --> CO + H2_REF_ELEMENT (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(steam_gasification)%Species(CGR_REF_ELEMENT)%pmap k6 = 7.0d1 * exp(-15096.6/xTg) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT) ! Reaction rate: RATES(steam_gasification) = k6 * p_H2O /MW_s(2,CGR_REF_ELEMENT) ! CO + H2_REF_ELEMENT --> CGR_REF_ELEMENT + H2O (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(steam_gasification_reverse)%Species(CGR_REF_ELEMENT)%pmap k7 = 7.0d1 * exp(1233.4/xTg - 17.29) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT) ! Reaction rate: RATES(steam_gasification_reverse) = k7 * p_H2 * p_H2O /MW_s(2,CGR_REF_ELEMENT) ! Reaction rate: RATES(Char_to_CO_oxidation) = 0.51926 * (1.0 - psi) * kox * p_O2/MW_s(2,CGR_REF_ELEMENT) ! CGR_REF_ELEMENT + CO2 --> 2 CO (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(CO2_gasification)%Species(CGR_REF_ELEMENT)%pmap k8 = 7.0d1 * exp(-15096.6/xTg) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT) ! Reaction rate: RATES(CO2_gasification) = k8 * p_CO2 /MW_s(2,CGR_REF_ELEMENT) ! 2 CO --> CGR_REF_ELEMENT + CO2 (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(CO2_gasification_reverse)%Species(CGR_REF_ELEMENT)%pmap k9 = 7.0d1 * exp(5185.4/xTg - 20.92) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT) ! Reaction rate: RATES(CO2_gasification_reverse) = k9 * p_CO * p_CO /MW_s(2,CGR_REF_ELEMENT) ! CGR_REF_ELEMENT + 2 H2 --> CH4 (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(Char_methanation)%Species(CGR_REF_ELEMENT)%pmap k10 = exp(-8078.0/xTg-7.087) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT) ! Reaction rate: RATES(Char_methanation) = k10 * p_CO2 /MW_s(2,CGR_REF_ELEMENT) ! CH4 --> CGR_REF_ELEMENT + 2 H2 (kg-mole/m^3.s) !---------------------------------------------------------------------// ! M = Reaction(Char_methanation_reverse)%Species(CGR_REF_ELEMENT)%pmap k11 = exp(-13578.0/xTg-0.372) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT) ! Reaction rate: RATES(Char_methanation_reverse) = k11 * sqrt(p_CH4) / MW_s(2,CGR_REF_ELEMENT) END IF ! do ir=1,20 ! print*,'RATES ',ir,RATES(ir) ! enddo ! Save reaction rates for visualization ! IF(nrr>=RX1F) ReactionRates(IJK,RX1F) = RATES(RX1F) ! IF(nrr>=RX1R) ReactionRates(IJK,RX1R) = RATES(RX1R) ! IF(nrr>=RX2F) ReactionRates(IJK,RX2F) = RATES(RX2F) ! IF(nrr>=RX2R) ReactionRates(IJK,RX2R) = RATES(RX2R) ! IF(nrr>=RX3) ReactionRates(IJK,RX3) = RATES(RX3) ! IF(nrr>=RX4) ReactionRates(IJK,RX4) = RATES(RX4) RETURN END SUBROUTINE USR_RATES