File: RELATIVE:/../../../mfix.git/model/chem/stiff_chem_rrates.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: STIFF_CHEM_RRATES(IER)                                  C
4     !                                                                      C
5     !  Purpose: Calculate reaction rates for various reactions in cell ijk C
6     !           using information from the data file                       C
7     !                                                                      C
8     !  Author: J. Musser                                  Date: 12-Feb-13  C
9     !  Reviewer:                                          Date: dd-mmm-yy  C
10     !                                                                      C
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
12           SUBROUTINE STIFF_CHEM_RRATES(lNEQ, lTime, Y, YDOT)
13     
14     ! External Module Procedures:
15     !---------------------------------------------------------------------//
16           use stiff_chem_maps, only: mapODEtoMFIX
17     
18     ! Global Parameters:
19     !---------------------------------------------------------------------//
20     ! Maximum number of solids phases
21           use param, only: DIMENSION_M
22     ! Maximum number of gas species
23           use param, only: DIMENSION_N_g
24     ! Maximum number of solid species
25           use param, only: DIMENSION_N_s
26     
27           use param1, only: ZERO
28     
29     ! Global Variables:
30     !---------------------------------------------------------------------//
31     ! Number of solids phases
32           use physprop, only: MMAX
33     ! Number of species
34           use physprop, only: NMAX
35     
36     ! Gas Phase
37     !--------------------------------------------//
38     !Temperature
39           use fldvar, only: T_g
40     ! Species Mass Fractions
41           use fldvar, only: X_g
42     ! Gas Phase apparent density
43           use fldvar, only: ROP_g
44     ! Gas phase specific heat
45           use physprop, only: C_pg
46     
47     ! Solids Phases
48     !--------------------------------------------//
49     ! Temperatures
50           use fldvar, only: T_s
51     ! Species mass fractions
52           use fldvar, only: X_s
53     ! Apparent densities
54           use fldvar, only: ROP_s
55     ! Specific heats
56           use physprop, only: C_ps
57     
58           use rxns, only: REACTION
59           use rxns, only: NO_OF_RXNS
60     
61           use run, only: UNITS
62     
63           use stiff_chem, only: ODE_DIMN_all
64           use stiff_chem, only: NEQ_DIMN
65     
66           implicit none
67     
68     ! Passed Variables: Dummy argument format required by ODEPACK.
69     !---------------------------------------------------------------------//
70     ! (1) Number of ODEs to be solve
71     ! (2) Fluid cell index
72           INTEGER, intent(in) :: lNEQ(NEQ_DIMN)
73     ! Independent variable (not used)
74           DOUBLE PRECISION, intent(in) :: lTime
75     ! Array of dependent variable initial values.
76           DOUBLE PRECISION, dimension(ODE_DIMN_all), intent(in)  :: Y
77     ! Rate of change of dependent variables.
78           DOUBLE PRECISION, dimension(ODE_DIMN_all), intent(out) :: YDOT
79     
80     ! Local Variables:
81     !---------------------------------------------------------------------//
82           INTEGER :: IJK  ! Fluid Cell Index
83           INTEGER :: H    ! Reaction loop counter
84           INTEGER :: L, M ! Global Phase index loop counters
85           INTEGER :: N    ! Global species index
86           INTEGER :: lN   ! Local reaction speices index/loop counter
87           INTEGER :: mXfr ! Global phase index for mass transfer
88           INTEGER :: Node ! Variable index
89     
90     ! User-defined reaction rates returned from USR_RATES
91           DOUBLE PRECISION :: RATES(NO_OF_RXNS)
92     ! Rate of formation (+) or consumption (-) of a species (GENERIC).
93           DOUBLE PRECISION :: lRate
94     ! Rate of formation (+) or consumption (-) of each gas phase species.
95     ! > Accumulative over all reactions.
96           DOUBLE PRECISION :: lRg(DIMENSION_N_g)
97     ! Rate of formation (+) or consumption (-) of each gas phase species.
98     ! > Single reaction.
99           DOUBLE PRECISION :: rRg(DIMENSION_N_g)
100     ! Heat of reaction assigned to gas phase for a reaction.
101           DOUBLE PRECISION :: rHORg
102     ! Cummulative heat of reaction assigned to gas phase.
103           DOUBLE PRECISION :: lHORg
104     ! Rate of formation (+) or consumption (-) of each solids phase species.
105     ! > Accumulative over all reactions.
106           DOUBLE PRECISION :: lRs(DIMENSION_M, DIMENSION_N_s)
107     ! Rate of formation (+) or consumption (-) of each solids phase species.
108     ! > Single reaction.
109           DOUBLE PRECISION :: rRs(DIMENSION_M, DIMENSION_N_s)
110     ! Heat of reaction assigned to the solids phase for a reaction.
111           DOUBLE PRECISION :: rHORs(DIMENSION_M)
112     ! Cummulative heat of reaction assigned to the solids phase.
113           DOUBLE PRECISION :: lHORs(DIMENSION_M)
114     ! Rate of interphase enthalpy transfer due to mass transfer.
115           DOUBLE PRECISION :: RxH(0:DIMENSION_M, 0:DIMENSION_M)
116     
117     ! Net rate of change of gas phase material.
118           DOUBLE PRECISION :: sumlRg
119     ! Net rate of change of solids phase material.
120           DOUBLE PRECISION :: sumlRs(DIMENSION_M)
121     
122     ! Reaction limiters.
123           DOUBLE PRECISION, parameter :: speciesLimiter = 1.0d-7
124     
125     ! External Function for comparing two numbers.
126           LOGICAL, external :: COMPARE
127           DOUBLE PRECISION, external :: CALC_H
128     
129     ! UDF for Reaction Rates:
130           external USR_RATES
131     
132     ! Initialize variables:
133           IJK    = lNEQ(2)
134           RATES  = ZERO
135           YDOT   = ZERO
136     
137           lRg = ZERO; lHORg = ZERO
138           lRs = ZERO; lHORs = ZERO
139     
140     ! Map the current ODE independent variables to MFIX variables.
141           CALL mapODEtoMFIX(NEQ_DIMN, lNEQ, ODE_DIMN_all, Y)
142     
143     ! Calculate user defined reaction rates.
144           CALL USR_RATES(IJK, RATES)
145     
146     ! Loop over reactions.
147           RXN_LP: DO H = 1, NO_OF_RXNS
148     
149     ! Skip empty reactions
150              IF(Reaction(H)%nSpecies == 0) CYCLE RXN_LP
151     
152     ! Initialize local loop arrays
153              rRg = ZERO; rHORg = ZERO
154              rRs = ZERO; rHORs = ZERO
155              RxH = ZERO
156     
157     ! Calculate the rate of formation/consumption for each species.
158     !---------------------------------------------------------------------//
159              DO lN = 1, Reaction(H)%nSpecies
160     ! Global phase index.
161                 M = Reaction(H)%Species(lN)%pMap
162     ! Global species index.
163                 N = Reaction(H)%Species(lN)%sMap
164     ! Index for interphase mass transfer. For a gas/solid reaction, the
165     ! index is stored with the gas phase. For solid/solid mass transfer
166     ! the index is stored with the source phase.
167                 mXfr = Reaction(H)%Species(lN)%mXfr
168                 lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
169     ! Gas Phase:
170                 IF(M == 0) THEN
171     ! Consumption of gas phase species.
172                    IF(lRate < ZERO) THEN
173     ! Check that there is a sufficient amount of reactant.
174                       IF(X_g(IJK,N) > speciesLimiter) THEN
175                          rRg(N) = rRg(N) + lRate
176     ! Enthalpy transfer associated with mass transfer. (gas/solid)
177                          IF(M /= mXfr) RxH(M,mXfr) =  RxH(M,mXfr) + &
178                             lRate * CALC_H(T_G(IJK),0,N)
179                       ELSE
180     ! There is an insignificant amount of reactant. Skip this reaction.
181                          CYCLE RXN_LP
182                       ENDIF
183                    ELSE
184     ! Formation of gas phase species.
185                       rRg(N) = rRg(N) + lRate
186     ! Enthalpy transfer associated with mass transfer. (gas/solid)
187                       IF(M /= mXfr) RxH(M,mXfr) =  RxH(M,mXfr) + &
188                          lRate * CALC_H(T_s(IJK,mXfr),0,N)
189                    ENDIF
190     ! Solids Phase M:
191                 ELSE
192     ! Consumption of solids phase species.
193                    IF(lRate < ZERO) THEN
194                       IF(X_s(IJK,M,N) > speciesLimiter) THEN
195                          rRs(M,N) = rRs(M,N) + lRate
196     ! Enthalpy transfer associated with mass transfer. (solid/solid) This
197     ! is only calculated from the source (reactant) material.
198                          IF(M /= mXfr) THEN
199                             IF(M < mXfr) THEN
200                                RxH(M,mXfr) =  RxH(M,mXfr) +                &
201                                   lRate * CALC_H(T_s(IJK,M),M,N) *         &
202                                   Reaction(H)%Species(lN)%xXfr
203                             ELSE
204                                RxH(mXfr,M) =  RxH(mXfr,M) -                &
205                                  lRate * CALC_H(T_s(IJK,M),M,N) *          &
206                                  Reaction(H)%Species(lN)%xXfr
207                             ENDIF
208                          ENDIF
209                       ELSE
210     ! There is an insignificant amount of reactant. Skip this reaction.
211                          CYCLE RXN_LP
212                       ENDIF
213                    ELSE
214     ! Formation of solids phase species.
215                       rRs(M,N) = rRs(M,N) + lRate
216                    ENDIF
217                 ENDIF
218              ENDDO ! Loop of species
219     
220     ! Copy the single reaction rates of formation/consumption to the total
221     ! (accumulative) rates of formation/consumption arrays. This ensures
222     ! that any reactions without sufficient reactants are not included.
223              lRg = lRg + rRg
224              lRs = lRs + rRs
225     
226     ! Calculate and store the heat of reaction.
227     !---------------------------------------------------------------------//
228     ! Automated heat of reaction calculations
229              IF(Reaction(H)%Calc_DH) THEN
230     ! Loop over reaction species.
231                 DO lN = 1, Reaction(H)%nSpecies
232     ! Global phase index.
233                    M = Reaction(H)%Species(lN)%pMap
234     ! Global species index.
235                    N = Reaction(H)%Species(lN)%sMap
236     ! Rate of formation/consumption for speices N
237                    lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
238     ! Gas phase enthalpy chnage from energy equation derivation.
239                    IF(M == 0) THEN
240                       rHORg = rHORg + lRate * CALC_H(T_g(IJK),0,N)
241     ! Solid phase enthalpy change from energy equation derivation.
242                    ELSE
243                       rHORs(M) = rHORs(M) + lRate * CALC_H(T_s(IJK,M),M,N)
244                    ENDIF
245                 ENDDO
246     
247     ! Complete the skew-symettric for enthalpy transfer with mass transfer
248                 DO M=1, MMAX
249                    DO L=0, M-1
250                       RxH(M,L) = - RxH(L,M)
251                    ENDDO
252                 ENDDO
253     ! Apply enthalpy transfer associated with mass transfer to get the
254     ! complete heat of reaction of heat phse for Reaction H.
255                 DO L=0, MMAX
256                    DO M = 0, MMAX
257                       IF(L == M) CYCLE
258                       IF(L == 0) THEN
259                          rHORg = rHORg - RxH(L,M)
260                       ELSE
261                          rHORs(L) = rHORs(L) - RxH(L,M)
262                       ENDIF
263                    ENDDO
264                 ENDDO
265     
266     ! Convert the heat of reaction to the appropriate units (if SI), and
267     ! store in the global array.
268                 IF(UNITS == 'SI') THEN
269                    lHORg = lHORg + 4.183925d3*rHORg
270                    DO M=1,MMAX
271                       lHORs(M) = lHORs(M) + 4.183925d3*rHORs(M)
272                    ENDDO
273                 ELSE
274                    lHORg = lHORg + rHORg
275                    DO M=1,MMAX
276                       lHORs(M) = lHORs(M) + rHORs(M)
277                    ENDDO
278                 ENDIF
279              ELSE
280     ! User-defined heat of reaction.
281                 lHORg = lHORg + Reaction(H)%HoR(0) * RATES(H)
282                 DO M=1, MMAX
283                    lHORs(M) = lHORs(M) + Reaction(H)%HoR(M) * RATES(H)
284                 ENDDO
285              ENDIF
286     
287           ENDDO RXN_LP ! Loop over reactions.
288     
289     
290     
291     
292     ! Calculate and store the heat of reaction.
293     !---------------------------------------------------------------------//
294     
295     ! Calculate the net change for the gas phase.
296           sumlRg = sum(lRg)
297     
298     ! Calculate the net change for solids phases.
299           sumlRs = 0.0d0
300           DO M=1, MMAX
301              IF(lNEQ(2+M) == 1) sumlRs(M) = sum(lRs(M,:))
302           ENDDO
303     
304     
305     ! Calculate and store the heat of reaction.
306     !---------------------------------------------------------------------//
307     ! Initialize counter.
308           Node = 1
309     
310     ! Density:  ROP_g
311           YDOT(Node) = sumlRg
312           Node = Node + 1
313     
314     ! Temperature: T_g
315           YDOT(Node) = -lHORg/(ROP_g(IJK)*C_pg(IJK))
316           Node = Node + 1
317     
318     ! Species mass fractions: X_g
319           DO N=1,NMAX(0)
320              YDOT(Node) = (lRg(N) - X_g(IJK,N)*sumlRg)/(ROP_g(IJK))
321              Node = Node + 1
322           ENDDO
323     
324     
325     ! Temperature: T_s
326           DO M = 1, MMAX
327              IF(ROP_s(IJK,M) > 1.0d-8) THEN
328                 YDOT(Node) = -lHORs(M)/(ROP_s(IJK,M)*C_ps(IJK,M))
329              ELSE
330                 YDOT(Node) = 0.0d0
331              ENDIF
332              Node = Node + 1
333           ENDDO
334     
335     
336           DO M = 1, MMAX
337              IF(lNEQ(2+M) == 1) THEN
338     ! Solids bulk density: ROP_s
339                 YDOT(Node) = sumlRs(M)
340                 Node = Node + 1
341     ! Species Mass Fraction: X_s
342                 DO N=1, NMAX(M)
343                    YDOT(Node) = (lRs(M,N) - X_s(IJK,M,N)*sumlRs(M)) /     &
344                       ROP_s(IJK,M)
345                    Node = Node + 1
346                 ENDDO
347              ENDIF
348           ENDDO
349     
350           RETURN
351           END SUBROUTINE STIFF_CHEM_RRATES
352