File: RELATIVE:/../../../mfix.git/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 :: N      ! 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           noSolids = DISCRETE_ELEMENT .AND. (.NOT.DES_CONTINUUM_HYBRID)
172           lMMAX = merge(0, MMAX, noSolids)
173     
174     ! Loop over each fluid cell.
175           DO IJK = ijkstart3, ijkend3
176           IF (FLUID_AT(IJK)) THEN
177     
178           RATES(:) = ZERO
179     
180     ! Calculate user defined reaction rates.
181           CALL USR_RATES(IJK, RATES)
182     
183     ! Loop over reactions.
184           RXN_LP: DO H = 1, NO_OF_RXNS
185     
186     ! Skip empty reactions
187              IF(Reaction(H)%nSpecies == 0) CYCLE RXN_LP
188              IF(COMPARE(RATES(H),ZERO)) CYCLE RXN_LP
189     
190     ! Initialize reaction loop arrays
191              r_Rgp = ZERO; r_ROXgc = ZERO; r_HoRg = ZERO
192              r_Rsp = ZERO; r_ROXsc = ZERO; r_HoRs = ZERO
193     
194              r_RxH = ZERO
195     
196     ! Calculate the rate of formation/consumption for each species.
197     !---------------------------------------------------------------------//
198              DO lN = 1, Reaction(H)%nSpecies
199     ! Global phase index.
200                 M = Reaction(H)%Species(lN)%pMap
201     ! Global species index.
202                 N = Reaction(H)%Species(lN)%sMap
203     ! Index for interphase mass transfer. For a gas/solid reaction, the
204     ! index is stored with the gas phase. For solid/solid mass transfer
205     ! the index is stored with the source phase.
206                 mXfr = Reaction(H)%Species(lN)%mXfr
207                 lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
208     ! Gas Phase:
209                 IF(M == 0) THEN
210     ! Consumption of gas phase species.
211                    IF(lRate < ZERO) THEN
212                       IF(X_g(IJK,N) > speciesLimiter) THEN
213                          r_ROXgc(N) = r_ROXgc(N) - lRate/X_g(IJK,N)
214     ! Enthalpy transfer associated with mass transfer. (gas/solid)
215                          IF(M /= mXfr) r_RxH(M,mXfr) = r_RxH(M,mXfr) +     &
216                             lRate * CALC_H(T_G(IJK),0,N)
217                       ELSE
218     ! There is an insignificant amount of reactant. Skip this reaction.
219                          RATES(H) = ZERO
220                          CYCLE RXN_LP
221                       ENDIF
222                    ELSE
223     ! Formation of gas phase species.
224                       r_Rgp(N) = r_Rgp(N) + lRate
225     ! Enthalpy transfer associated with mass transfer. (gas/solid)
226                       IF(M /= mXfr) r_RxH(M,mXfr) = r_RxH(M,mXfr) +        &
227                          lRate * CALC_H(T_s(IJK,mXfr),0,N)
228                    ENDIF
229     
230     
231     ! Solids Phase M:
232                 ELSE
233     ! Consumption of solids phase species.
234                    IF(lRate < ZERO) THEN
235                       IF(X_s(IJK,M,N) > speciesLimiter) THEN
236                          r_ROXsc(M,N) = r_ROXsc(M,N) - lRate/X_s(IJK,M,N)
237     ! Enthalpy transfer associated with mass transfer. (solid/solid) This
238     ! is only calculated from the source (reactant) material.
239                          IF(M /= mXfr) THEN
240                             IF(M < mXfr) THEN
241                                r_RxH(M,mXfr) =  r_RxH(M,mXfr) + lRate *    &
242                                   Reaction(H)%Species(lN)%xXfr *           &
243                                   CALC_H(T_s(IJK,M),M,N)
244                             ELSE
245                                r_RxH(mXfr,M) =  r_RxH(mXfr,M) - lRate *    &
246                                   Reaction(H)%Species(lN)%xXfr *           &
247                                   CALC_H(T_s(IJK,M),M,N)
248                             ENDIF
249                          ENDIF
250                       ELSE
251     ! There is an insignificant amount of reactant. Skip this reaction.
252                          RATES(H) = ZERO
253                          CYCLE RXN_LP
254                       ENDIF
255                    ELSE
256     ! Formation of solids phase species.
257                       r_Rsp(M,N) = r_Rsp(M,N) + lRate
258                    ENDIF
259                 ENDIF
260              ENDDO ! Loop of species
261     
262     
263     
264     ! Map the local reaction production/consumption into global variables.
265     !---------------------------------------------------------------------//
266     ! The outer reaction loop was cycled if there was insufficient amount
267     ! of species available.
268     
269     ! Gas phase production
270              R_gp(IJK,:NMAX(0)) = R_gp(IJK,:NMAX(0)) +                     &
271                 r_Rgp(:NMAX(0))
272     ! Gas phase consumption divided by species mass fraction
273              ROX_gc(IJK,:NMAX(0)) = ROX_gc(IJK,:NMAX(0)) +                 &
274                 r_ROXgc(:NMAX(0))
275     
276              DO M=1,lMMAX
277     ! Solids phase species production
278                 R_sp(IJK,M,:NMAX(M)) = R_sp(IJK,M,:NMAX(M)) +              &
279                    r_Rsp(M,:NMAX(M))
280     ! Solids phase species consumption divided by species mass fraction
281                 ROX_sc(IJK,M,:NMAX(M)) = ROX_sc(IJK,M,:NMAX(M)) +          &
282                    r_ROXsc(M,:NMAX(M))
283              ENDDO
284     
285     
286     ! Calculate and store the heat of reaction.
287     !---------------------------------------------------------------------//
288              IF(ENERGY_EQ) THEN
289     ! Automated heat of reaction calculations
290                 IF(Reaction(H)%Calc_DH) THEN
291     ! Loop over reaction species.
292                    DO lN = 1, Reaction(H)%nSpecies
293     ! Global phase index.
294                       M = Reaction(H)%Species(lN)%pMap
295     ! Global species index.
296                       N = Reaction(H)%Species(lN)%sMap
297     ! Rate of formation/consumption for speices N
298                       lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
299     ! Gas phase enthalpy change from energy equation derivation.
300                       IF(M == 0) THEN
301                          r_HORg = r_HORg + CALC_H(T_g(IJK),0,N) * lRate
302     ! Solid phase enthalpy change from energy equation derivation.
303                       ELSE
304                          r_HORs(M) = r_HORs(M) + CALC_H(T_s(IJK,M),M,N)*lRate
305                       ENDIF
306                    ENDDO
307     
308     ! Complete the skew-symettric for enthalpy transfer with mass transfer
309                    DO M=1, lMMAX
310                        DO L=0, M-1
311                         r_RxH(M,L) = - r_RxH(L,M)
312                       ENDDO
313                    ENDDO
314     ! Apply enthalpy transfer associated with mass transfer to get the
315     ! complete heat of reaction of heat phse for Reaction H.
316                    DO L=0, lMMAX
317                       DO M = 0, lMMAX
318                          IF(L == M) CYCLE
319                          IF(L == 0) THEN
320                             r_HORg = r_HORg - r_RxH(L,M)
321                          ELSE
322                             r_HORs(L) = r_HORs(L) - r_RxH(L,M)
323                          ENDIF
324                       ENDDO
325                    ENDDO
326     ! Apply unit conversion and store heats of reaction in the global array.
327                    HOR_g(IJK) = HOR_g(IJK) + l_2SI*r_HORg
328                    DO M=1,lMMAX
329                       HOR_s(IJK,M) = HOR_s(IJK,M) + l_2SI*r_HORs(M)
330                    ENDDO
331                 ELSE
332     ! User-defined heat of reaction.
333                    HOR_g(IJK) = HOR_g(IJK) + Reaction(H)%HoR(0) * RATES(H)
334                    DO M=1, lMMAX
335                       HOR_s(IJK,M) = HOR_s(IJK,M) +                        &
336                          Reaction(H)%HoR(M) * RATES(H)
337                    ENDDO
338                 ENDIF
339              ENDIF  ! ENERGY_EQ
340     
341     ! Update rate of interphase mass transfer.
342     !---------------------------------------------------------------------//
343               DO LM=1, (DIMENSION_LM+DIMENSION_M-1)
344                  R_PHASE(IJK,LM) = R_PHASE(IJK,LM) +                       &
345                     RATES(H) * Reaction(H)%rPHASE(LM)
346               ENDDO
347           ENDDO RXN_LP ! Loop over reactions.
348     
349     
350     ! Calculate the toal rate of formation and consumption for each species.
351     !---------------------------------------------------------------------//
352           IF(SPECIES_EQ(0)) THEN
353              SUM_R_G(IJK) = SUM( R_gp(IJK,:NMAX(0)) -                      &
354                 ROX_gc(IJK,:NMAX(0))*X_g(IJK,:NMAX(0)))
355           ELSE
356              DO H=1, NO_OF_RXNS
357                 IF(Reaction(H)%nPhases <= 0) CYCLE
358                 DO M=1, lMMAX
359                    LM = 1 + ((M-1)*M)/2
360                    SUM_R_G(IJK) = SUM_R_G(IJK) +                           &
361                       RATES(H) * Reaction(H)%rPHASE(LM)
362                 ENDDO
363              ENDDO
364           ENDIF
365     
366           DO M=1, lMMAX
367              IF(SPECIES_EQ(M)) THEN
368                 SUM_R_S(IJK,M) = SUM( R_sp(IJK,M,:NMAX(M)) -               &
369                    RoX_sc(IJK,M,:NMAX(M))*X_s(IJK,M,:NMAX(M)))
370              ELSE
371                 DO H=1, NO_OF_RXNS
372                    IF(Reaction(H)%nPhases <= 0) CYCLE
373                    DO L=0, M-1
374                       LM = 1 + L + ((M-1)*M)/2
375                       SUM_R_S(IJK,M) = SUM_R_S(IJK,M) -                    &
376                          RATES(H) * Reaction(H)%rPHASE(LM)
377                    ENDDO
378                    DO L=M+1, lMMAX
379                       LM = 1 + M + ((L-1)*L)/2
380                       SUM_R_S(IJK,M) = SUM_R_S(IJK,M) +                    &
381                          RATES(H) * Reaction(H)%rPHASE(LM)
382                    ENDDO
383                 ENDDO
384              ENDIF
385           ENDDO
386     
387     
388           ENDIF  ! Fluid_At(IJK)
389           END DO ! IJK
390     
391           RETURN
392     
393           END SUBROUTINE RRATES0
394