File: N:\mfix\model\rrates0.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: RRATES0(IER)                                           C
4     !  Purpose: Calculate reaction rates for various reactions in cell ijk C
5     !           using information from the data file                       C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 3-10-98    C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number: 1                                                  C
11     !  Purpose:Replaced routines with new proceedures for automated        C
12     !          reaction rate calculations.                                 C
13     !  Author: J. Musser                                  Date: 10-Oct-12  C
14     !  Reviewer:                                          Date: dd-mmm-yy  C
15     !                                                                      C
16     !  Literature/Document References:                                     C
17     !                                                                      C
18     !  Variables referenced: MMAX, IJK, T_g, T_s1, D_p, X_g, X_s, EP_g,    C
19     !            P_g, HOR_g, HOR_s                                         C
20     !                                                                      C
21     !                                                                      C
22     !  Variables modified: M, N, R_gp, R_sp, RoX_gc, RoX_sc, SUM_R_g,      C
23     !                      SUM_R_s                                         C
24     !                                                                      C
25     !  Local variables:                                                    C
26     !                                                                      C
27     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
28           SUBROUTINE RRATES0()
29     
30     ! Global domain parameters.
31     !`````````````````````````````````````````````````````````````````````//
32     ! Blanket access is provided to the following modules as they don't
33     ! contain variables which are likely to be modified 'by accident.'
34           use compar
35           use indices
36           use geometry
37     
38     ! Gas phase global variables:
39     !`````````````````````````````````````````````````````````````````````//
40     ! Species mass fractions.
41           use fldvar, only: X_g
42     ! Temperature.
43           use fldvar, only: T_g
44     ! Rate of production.
45           use rxns, only: R_gp
46     ! Rate of consumption divided by species mass fraction.
47           use rxns, only: RoX_gc
48     ! Net rate of change of gas phase mass.
49           use rxns, only: SUM_R_g
50     ! Heat of reaction source term
51           use energy, only: HOR_g
52     
53     ! Gas phase global variables:
54     !`````````````````````````````````````````````````````````````````````//
55     ! Species mass fractions.
56           use fldvar, only: X_s
57     ! Temperatures.
58           use fldvar, only: T_s
59     ! Rate of production.
60           use rxns, only: R_sp
61     ! Rate of consumption divided by species mass fraction.
62           use rxns, only: RoX_sc
63     ! Net rate of change of gas phase mass.
64           use rxns, only: SUM_R_s
65     ! Heat of reaction source term
66           use energy, only: HOR_s
67     
68     
69     ! Reaction global variables:
70     !`````````````````````````````````````````````````````````````````````//
71     ! Number of TFM reactions.
72           use rxns, only: NO_OF_RXNS
73     ! Rate of interphase mass transfer.
74           use rxns, only: R_phase
75     ! Data structure storing reaction information.
76           use rxns, only: Reaction
77     
78     
79     ! Global run-time parameters.
80     !`````````````````````````````````````````````````````````````````````//
81     ! Number of solids phases (or one if MMAX is zero)
82           use param, only: DIMENSION_M
83     ! Number of gas species.
84           use param, only: DIMENSION_N_g
85     ! Number of solids phase species.
86           use param, only: DIMENSION_N_s
87           use param1, only: dimension_lm
88           use param1, only: zero
89     ! Indicates that the energy equations are solved.
90           use run, only: ENERGY_EQ
91     ! Indicates that the species equations are solved.
92           use run, only: SPECIES_EQ
93     ! Simulation units (CGS or SI)
94           use run, only: UNITS
95     ! Number of solids phases.
96           use physprop, only: MMAX
97     ! Number of species comprising each phase.
98           use physprop, only: NMAX
99     ! Indicate that DEM is used.
100           use discretelement, only: DISCRETE_ELEMENT
101     ! Indicate that TFM/DEM hybrid is used.
102           use discretelement, only: DES_CONTINUUM_HYBRID
103     ! Small value for species mass fractions
104           use toleranc
105           use functions
106     
107           implicit none
108     
109     ! Local variables:
110     !`````````````````````````````````````````````````````````````````````//
111     ! Index/loop counters
112           INTEGER :: IJK    ! fluid cell index
113           INTEGER :: H      ! reaction loop counter
114           INTEGER :: L, M   ! global phase index loop counters
115           INTEGER :: NN      ! global species index
116           INTEGER :: lN     ! reaction speices index
117           INTEGER :: lM     ! reaction phase index
118           INTEGER :: mXfr   ! global phase index for mass transfer
119     
120     ! This is a mask for the number of solids phases. The mask is used to
121     ! avoid issues for DEM based reactions where MMAX and DES_MMAX may be
122     ! in conflict.
123           INTEGER :: lMMAX
124           LOGICAL :: noSolids
125     
126     ! User-defined reaction rates returned from USR_RATES
127           DOUBLE PRECISION RATES(NO_OF_RXNS)
128     ! Reaction and species specific reaction rate.
129           DOUBLE PRECISION lRate
130     
131     ! Single reaction formation/consumption rate of each gas species.
132           DOUBLE PRECISION :: r_Rgp(DIMENSION_N_g)
133           DOUBLE PRECISION :: r_ROXgc(DIMENSION_N_g)
134     ! Single reaction formation/consumption rate of each solids species.
135           DOUBLE PRECISION :: r_Rsp(DIMENSION_M, DIMENSION_N_s)
136           DOUBLE PRECISION :: r_ROXsc(DIMENSION_M, DIMENSION_N_s)
137     
138     ! Rate of interphase enthalpy transfer
139           DOUBLE PRECISION :: r_RxH(0:MMAX, 0:MMAX)
140     ! Local heat of reactions
141           DOUBLE PRECISION :: r_HoRg, r_HoRs(1:MMAX)
142     ! Conversion factor for HoR (CGS to SI)
143           DOUBLE PRECISION :: l_2SI
144     
145     ! Reaction limiters. If a species mass fraction is less than this
146     ! value, then the reaction is suppressed.
147           DOUBLE PRECISION :: speciesLimiter
148     
149     ! External functions:
150     !`````````````````````````````````````````````````````````````````````//
151     ! Calculates enthalpy (cal/gram)
152           DOUBLE PRECISION, EXTERNAL :: CALC_H
153     
154     ! Initialize global storage arrays to zero
155     !---------------------------------------------------------------------//
156     !     Gas                Solids
157           R_GP    = ZERO;    R_SP    = ZERO
158           ROX_GC  = ZERO;    ROX_SC  = ZERO
159           SUM_R_G = ZERO;    SUM_R_S = ZERO
160           HOR_G   = ZERO;    HOR_S   = ZERO
161     
162           R_PHASE = ZERO
163     
164     ! Set the species limiter:
165           speciesLimiter = ZERO_X_gs
166     
167     ! Set the conversion factor for heat of reaction.
168           l_2SI = merge(4.183925d3, ONE, UNITS(1:2) == 'SI')
169     
170     ! Set the number of TFM solids phases.
171     ! the following 2 lines could be replaced by assing lmmax=smax
172           noSolids = DISCRETE_ELEMENT .AND. (.NOT.DES_CONTINUUM_HYBRID)
173           lMMAX = merge(0, MMAX, noSolids)
174     
175     ! Loop over each fluid cell.
176           DO IJK = ijkstart3, ijkend3
177           IF (FLUID_AT(IJK)) THEN
178     
179           RATES(:) = ZERO
180     
181     ! Calculate user defined reaction rates.
182           CALL USR_RATES(IJK, RATES)
183     
184     ! Loop over reactions.
185           RXN_LP: DO H = 1, NO_OF_RXNS
186     
187     ! Skip empty reactions
188              IF(Reaction(H)%nSpecies == 0) CYCLE RXN_LP
189              IF(COMPARE(RATES(H),ZERO)) CYCLE RXN_LP
190     
191     ! Initialize reaction loop arrays
192              r_Rgp = ZERO; r_ROXgc = ZERO; r_HoRg = ZERO
193              r_Rsp = ZERO; r_ROXsc = ZERO; r_HoRs = ZERO
194     
195              r_RxH = ZERO
196     
197     ! Calculate the rate of formation/consumption for each species.
198     !---------------------------------------------------------------------//
199              DO lN = 1, Reaction(H)%nSpecies
200     ! Global phase index.
201                 M = Reaction(H)%Species(lN)%pMap
202     ! Global species index.
203                 NN = Reaction(H)%Species(lN)%sMap
204     ! Index for interphase mass transfer. For a gas/solid reaction, the
205     ! index is stored with the gas phase. For solid/solid mass transfer
206     ! the index is stored with the source phase.
207                 mXfr = Reaction(H)%Species(lN)%mXfr
208                 lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
209     ! Gas Phase:
210                 IF(M == 0) THEN
211     ! Consumption of gas phase species.
212                    IF(lRate < ZERO) THEN
213                       IF(X_g(IJK,NN) > speciesLimiter) THEN
214                          r_ROXgc(NN) = r_ROXgc(NN) - lRate/X_g(IJK,NN)
215     ! Enthalpy transfer associated with mass transfer. (gas/solid)
216                          IF(M /= mXfr) r_RxH(M,mXfr) = r_RxH(M,mXfr) +     &
217                             lRate * CALC_H(T_G(IJK),0,NN)
218                       ELSE
219     ! There is an insignificant amount of reactant. Skip this reaction.
220                          RATES(H) = ZERO
221                          CYCLE RXN_LP
222                       ENDIF
223                    ELSE
224     ! Formation of gas phase species.
225                       r_Rgp(NN) = r_Rgp(NN) + lRate
226     ! Enthalpy transfer associated with mass transfer. (gas/solid)
227                       IF(M /= mXfr) r_RxH(M,mXfr) = r_RxH(M,mXfr) +        &
228                          lRate * CALC_H(T_s(IJK,mXfr),0,NN)
229                    ENDIF
230     
231     
232     ! Solids Phase M:
233                 ELSE
234     ! Consumption of solids phase species.
235                    IF(lRate < ZERO) THEN
236                       IF(X_s(IJK,M,NN) > speciesLimiter) THEN
237                          r_ROXsc(M,NN) = r_ROXsc(M,NN) - lRate/X_s(IJK,M,NN)
238     ! Enthalpy transfer associated with mass transfer. (solid/solid) This
239     ! is only calculated from the source (reactant) material.
240                          IF(M /= mXfr) THEN
241                             IF(M < mXfr) THEN
242                                r_RxH(M,mXfr) =  r_RxH(M,mXfr) + lRate *    &
243                                   Reaction(H)%Species(lN)%xXfr *           &
244                                   CALC_H(T_s(IJK,M),M,NN)
245                             ELSE
246                                r_RxH(mXfr,M) =  r_RxH(mXfr,M) - lRate *    &
247                                   Reaction(H)%Species(lN)%xXfr *           &
248                                   CALC_H(T_s(IJK,M),M,NN)
249                             ENDIF
250                          ENDIF
251                       ELSE
252     ! There is an insignificant amount of reactant. Skip this reaction.
253                          RATES(H) = ZERO
254                          CYCLE RXN_LP
255                       ENDIF
256                    ELSE
257     ! Formation of solids phase species.
258                       r_Rsp(M,NN) = r_Rsp(M,NN) + lRate
259                    ENDIF
260                 ENDIF
261              ENDDO ! Loop of species
262     
263     
264     
265     ! Map the local reaction production/consumption into global variables.
266     !---------------------------------------------------------------------//
267     ! The outer reaction loop was cycled if there was insufficient amount
268     ! of species available.
269     
270     ! Gas phase production
271              R_gp(IJK,:NMAX(0)) = R_gp(IJK,:NMAX(0)) +                     &
272                 r_Rgp(:NMAX(0))
273     ! Gas phase consumption divided by species mass fraction
274              ROX_gc(IJK,:NMAX(0)) = ROX_gc(IJK,:NMAX(0)) +                 &
275                 r_ROXgc(:NMAX(0))
276     
277              DO M=1,lMMAX
278     ! Solids phase species production
279                 R_sp(IJK,M,:NMAX(M)) = R_sp(IJK,M,:NMAX(M)) +              &
280                    r_Rsp(M,:NMAX(M))
281     ! Solids phase species consumption divided by species mass fraction
282                 ROX_sc(IJK,M,:NMAX(M)) = ROX_sc(IJK,M,:NMAX(M)) +          &
283                    r_ROXsc(M,:NMAX(M))
284              ENDDO
285     
286     
287     ! Calculate and store the heat of reaction.
288     !---------------------------------------------------------------------//
289              IF(ENERGY_EQ) THEN
290     ! Automated heat of reaction calculations
291                 IF(Reaction(H)%Calc_DH) THEN
292     ! Loop over reaction species.
293                    DO lN = 1, Reaction(H)%nSpecies
294     ! Global phase index.
295                       M = Reaction(H)%Species(lN)%pMap
296     ! Global species index.
297                       NN = Reaction(H)%Species(lN)%sMap
298     ! Rate of formation/consumption for speices N
299                       lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
300     ! Gas phase enthalpy change from energy equation derivation.
301                       IF(M == 0) THEN
302                          r_HORg = r_HORg + CALC_H(T_g(IJK),0,NN) * lRate
303     ! Solid phase enthalpy change from energy equation derivation.
304                       ELSE
305                          r_HORs(M) = r_HORs(M) + CALC_H(T_s(IJK,M),M,NN)*lRate
306                       ENDIF
307                    ENDDO
308     
309     ! Complete the skew-symettric for enthalpy transfer with mass transfer
310                    DO M=1, lMMAX
311                        DO L=0, M-1
312                         r_RxH(M,L) = - r_RxH(L,M)
313                       ENDDO
314                    ENDDO
315     ! Apply enthalpy transfer associated with mass transfer to get the
316     ! complete heat of reaction of heat phse for Reaction H.
317                    DO L=0, lMMAX
318                       DO M = 0, lMMAX
319                          IF(L == M) CYCLE
320                          IF(L == 0) THEN
321                             r_HORg = r_HORg - r_RxH(L,M)
322                          ELSE
323                             r_HORs(L) = r_HORs(L) - r_RxH(L,M)
324                          ENDIF
325                       ENDDO
326                    ENDDO
327     ! Apply unit conversion and store heats of reaction in the global array.
328                    HOR_g(IJK) = HOR_g(IJK) + l_2SI*r_HORg
329                    DO M=1,lMMAX
330                       HOR_s(IJK,M) = HOR_s(IJK,M) + l_2SI*r_HORs(M)
331                    ENDDO
332                 ELSE
333     ! User-defined heat of reaction.
334                    HOR_g(IJK) = HOR_g(IJK) + Reaction(H)%HoR(0) * RATES(H)
335                    DO M=1, lMMAX
336                       HOR_s(IJK,M) = HOR_s(IJK,M) +                        &
337                          Reaction(H)%HoR(M) * RATES(H)
338                    ENDDO
339                 ENDIF
340              ENDIF  ! ENERGY_EQ
341     
342     ! Update rate of interphase mass transfer.
343     !---------------------------------------------------------------------//
344               DO LM=1, (DIMENSION_LM+DIMENSION_M-1)
345                  R_PHASE(IJK,LM) = R_PHASE(IJK,LM) +                       &
346                     RATES(H) * Reaction(H)%rPHASE(LM)
347               ENDDO
348           ENDDO RXN_LP ! Loop over reactions.
349     
350     
351     ! Calculate the toal rate of formation and consumption for each species.
352     !---------------------------------------------------------------------//
353           IF(SPECIES_EQ(0)) THEN
354              SUM_R_G(IJK) = SUM( R_gp(IJK,:NMAX(0)) -                      &
355                 ROX_gc(IJK,:NMAX(0))*X_g(IJK,:NMAX(0)))
356           ELSE
357              DO H=1, NO_OF_RXNS
358                 IF(Reaction(H)%nPhases <= 0) CYCLE
359                 DO M=1, lMMAX
360                    LM = 1 + ((M-1)*M)/2
361                    SUM_R_G(IJK) = SUM_R_G(IJK) +                           &
362                       RATES(H) * Reaction(H)%rPHASE(LM)
363                 ENDDO
364              ENDDO
365           ENDIF
366     
367           DO M=1, lMMAX
368              IF(SPECIES_EQ(M)) THEN
369                 SUM_R_S(IJK,M) = SUM( R_sp(IJK,M,:NMAX(M)) -               &
370                    RoX_sc(IJK,M,:NMAX(M))*X_s(IJK,M,:NMAX(M)))
371              ELSE
372                 DO H=1, NO_OF_RXNS
373                    IF(Reaction(H)%nPhases <= 0) CYCLE
374                    DO L=0, M-1
375                       LM = 1 + L + ((M-1)*M)/2
376                       SUM_R_S(IJK,M) = SUM_R_S(IJK,M) -                    &
377                          RATES(H) * Reaction(H)%rPHASE(LM)
378                    ENDDO
379                    DO L=M+1, lMMAX
380                       LM = 1 + M + ((L-1)*L)/2
381                       SUM_R_S(IJK,M) = SUM_R_S(IJK,M) +                    &
382                          RATES(H) * Reaction(H)%rPHASE(LM)
383                    ENDDO
384                 ENDDO
385              ENDIF
386           ENDDO
387     
388     
389           ENDIF  ! Fluid_At(IJK)
390           END DO ! IJK
391     
392           RETURN
393     
394           END SUBROUTINE RRATES0
395