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