File: /nfs/home/0/users/jenkins/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(IER)
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     ! Indicates that the energy equations are solved.
88           use run, only: ENERGY_EQ
89     ! Indicates that the species equations are solved.
90           use run, only: SPECIES_EQ
91     ! Simulation units (CGS or SI)
92           use run, only: UNITS
93     ! Number of solids phases.
94           use physprop, only: MMAX
95     ! Number of species comprising each phase.
96           use physprop, only: NMAX
97     ! Indicate that DEM is used.
98           use discretelement, only: DISCRETE_ELEMENT
99     ! Indicate that TFM/DEM hybrid is used.
100           use discretelement, only: DES_CONTINUUM_HYBRID
101     ! Small value for species mass fractions
102           use toleranc
103           use functions
104     
105           implicit none
106     
107     ! Dummy Arguments:
108     !`````````````````````````````````````````````````````````````````````//
109     ! Error index
110           INTEGER :: IER
111     
112     
113     ! Local variables:
114     !`````````````````````````````````````````````````````````````````````//
115     ! Index/loop counters
116           INTEGER :: IJK    ! fluid cell index
117           INTEGER :: H      ! reaction loop counter
118           INTEGER :: L, M   ! global phase index loop counters
119           INTEGER :: N      ! global species index
120           INTEGER :: lN     ! reaction speices index
121           INTEGER :: lM     ! reaction phase index
122           INTEGER :: mXfr   ! global phase index for mass transfer
123     
124     ! This is a mask for the number of solids phases. The mask is used to
125     ! avoid issues for DEM based reactions where MMAX and DES_MMAX may be
126     ! in conflict.
127           INTEGER :: lMMAX
128           LOGICAL :: noSolids
129     
130     ! User-defined reaction rates returned from USR_RATES
131           DOUBLE PRECISION RATES(NO_OF_RXNS)
132     ! Reaction and species specific reaction rate.
133           DOUBLE PRECISION lRate
134     
135     ! Single reaction formation/consumption rate of each gas species.
136           DOUBLE PRECISION :: r_Rgp(DIMENSION_N_g)
137           DOUBLE PRECISION :: r_ROXgc(DIMENSION_N_g)
138     ! Single reaction formation/consumption rate of each solids species.
139           DOUBLE PRECISION :: r_Rsp(DIMENSION_M, DIMENSION_N_s)
140           DOUBLE PRECISION :: r_ROXsc(DIMENSION_M, DIMENSION_N_s)
141     
142     ! Rate of interphase enthalpy transfer
143           DOUBLE PRECISION :: r_RxH(0:MMAX, 0:MMAX)
144     ! Local heat of reactions
145           DOUBLE PRECISION :: r_HoRg, r_HoRs(1:MMAX)
146     ! Conversion factor for HoR (CGS to SI)
147           DOUBLE PRECISION :: l_2SI
148     
149     ! Reaction limiters. If a species mass fraction is less than this
150     ! value, then the reaction is suppressed.
151           DOUBLE PRECISION :: speciesLimiter
152     
153     ! External functions:
154     !`````````````````````````````````````````````````````````````````````//
155     ! Calculates enthalpy (cal/gram)
156           DOUBLE PRECISION, EXTERNAL :: CALC_H
157     
158     ! Initialize global storage arrays to zero
159     !---------------------------------------------------------------------//
160     !     Gas                Solids
161           R_GP    = ZERO;    R_SP    = ZERO
162           ROX_GC  = ZERO;    ROX_SC  = ZERO
163           SUM_R_G = ZERO;    SUM_R_S = ZERO
164           HOR_G   = ZERO;    HOR_S   = ZERO
165     
166           R_PHASE = ZERO
167     
168     ! Set the species limiter:
169           speciesLimiter = ZERO_X_gs
170     
171     ! Set the conversion factor for heat of reaction.
172           l_2SI = merge(4.183925d3, ONE, UNITS(1:2) == 'SI')
173     
174     ! Set the number of TFM solids phases.
175           noSolids = DISCRETE_ELEMENT .AND. (.NOT.DES_CONTINUUM_HYBRID)
176           lMMAX = merge(0, MMAX, noSolids)
177     
178     ! Loop over each fluid cell.
179           DO IJK = ijkstart3, ijkend3
180           IF (FLUID_AT(IJK)) THEN
181     
182           RATES(:) = ZERO
183     
184     ! Calculate user defined reaction rates.
185           CALL USR_RATES(IJK, RATES)
186     
187     ! Loop over reactions.
188           RXN_LP: DO H = 1, NO_OF_RXNS
189     
190     ! Skip empty reactions
191              IF(Reaction(H)%nSpecies == 0) CYCLE RXN_LP
192              IF(COMPARE(RATES(H),ZERO)) CYCLE RXN_LP
193     
194     ! Initialize reaction loop arrays
195              r_Rgp = ZERO; r_ROXgc = ZERO; r_HoRg = ZERO
196              r_Rsp = ZERO; r_ROXsc = ZERO; r_HoRs = ZERO
197     
198              r_RxH = ZERO
199     
200     ! Calculate the rate of formation/consumption for each species.
201     !---------------------------------------------------------------------//
202              DO lN = 1, Reaction(H)%nSpecies
203     ! Global phase index.
204                 M = Reaction(H)%Species(lN)%pMap
205     ! Global species index.
206                 N = Reaction(H)%Species(lN)%sMap
207     ! Index for interphase mass transfer. For a gas/solid reaction, the
208     ! index is stored with the gas phase. For solid/solid mass transfer
209     ! the index is stored with the source phase.
210                 mXfr = Reaction(H)%Species(lN)%mXfr
211                 lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
212     ! Gas Phase:
213                 IF(M == 0) THEN
214     ! Consumption of gas phase species.
215                    IF(lRate < ZERO) THEN
216                       IF(X_g(IJK,N) > speciesLimiter) THEN
217                          r_ROXgc(N) = r_ROXgc(N) - lRate/X_g(IJK,N)
218     ! Enthalpy transfer associated with mass transfer. (gas/solid)
219                          IF(M /= mXfr) r_RxH(M,mXfr) = r_RxH(M,mXfr) +     &
220                             lRate * CALC_H(T_G(IJK),0,N)
221                       ELSE
222     ! There is an insignificant amount of reactant. Skip this reaction.
223                          RATES(H) = ZERO
224                          CYCLE RXN_LP
225                       ENDIF
226                    ELSE
227     ! Formation of gas phase species.
228                       r_Rgp(N) = r_Rgp(N) + lRate
229     ! Enthalpy transfer associated with mass transfer. (gas/solid)
230                       IF(M /= mXfr) r_RxH(M,mXfr) = r_RxH(M,mXfr) +        &
231                          lRate * CALC_H(T_s(IJK,mXfr),0,N)
232                    ENDIF
233     
234     
235     ! Solids Phase M:
236                 ELSE
237     ! Consumption of solids phase species.
238                    IF(lRate < ZERO) THEN
239                       IF(X_s(IJK,M,N) > speciesLimiter) THEN
240                          r_ROXsc(M,N) = r_ROXsc(M,N) - lRate/X_s(IJK,M,N)
241     ! Enthalpy transfer associated with mass transfer. (solid/solid) This
242     ! is only calculated from the source (reactant) material.
243                          IF(M /= mXfr) THEN
244                             IF(M < mXfr) THEN
245                                r_RxH(M,mXfr) =  r_RxH(M,mXfr) + lRate *    &
246                                   Reaction(H)%Species(lN)%xXfr *           &
247                                   CALC_H(T_s(IJK,M),M,N)
248                             ELSE
249                                r_RxH(mXfr,M) =  r_RxH(mXfr,M) - lRate *    &
250                                   Reaction(H)%Species(lN)%xXfr *           &
251                                   CALC_H(T_s(IJK,M),M,N)
252                             ENDIF
253                          ENDIF
254                       ELSE
255     ! There is an insignificant amount of reactant. Skip this reaction.
256                          RATES(H) = ZERO
257                          CYCLE RXN_LP
258                       ENDIF
259                    ELSE
260     ! Formation of solids phase species.
261                       r_Rsp(M,N) = r_Rsp(M,N) + lRate
262                    ENDIF
263                 ENDIF
264              ENDDO ! Loop of species
265     
266     
267     
268     ! Map the local reaction production/consumption into global variables.
269     !---------------------------------------------------------------------//
270     ! The outer reaction loop was cycled if there was insufficient amount
271     ! of species available.
272     
273     ! Gas phase production
274              R_gp(IJK,:NMAX(0)) = R_gp(IJK,:NMAX(0)) +                     &
275                 r_Rgp(:NMAX(0))
276     ! Gas phase consumption divided by species mass fraction
277              ROX_gc(IJK,:NMAX(0)) = ROX_gc(IJK,:NMAX(0)) +                 &
278                 r_ROXgc(:NMAX(0))
279     
280              DO M=1,lMMAX
281     ! Solids phase species production
282                 R_sp(IJK,M,:NMAX(M)) = R_sp(IJK,M,:NMAX(M)) +              &
283                    r_Rsp(M,:NMAX(M))
284     ! Solids phase species consumption divided by species mass fraction
285                 ROX_sc(IJK,M,:NMAX(M)) = ROX_sc(IJK,M,:NMAX(M)) +          &
286                    r_ROXsc(M,:NMAX(M))
287              ENDDO
288     
289     
290     ! Calculate and store the heat of reaction.
291     !---------------------------------------------------------------------//
292              IF(ENERGY_EQ) THEN
293     ! Automated heat of reaction calculations
294                 IF(Reaction(H)%Calc_DH) THEN
295     ! Loop over reaction species.
296                    DO lN = 1, Reaction(H)%nSpecies
297     ! Global phase index.
298                       M = Reaction(H)%Species(lN)%pMap
299     ! Global species index.
300                       N = Reaction(H)%Species(lN)%sMap
301     ! Rate of formation/consumption for speices N
302                       lRate = RATES(H) * Reaction(H)%Species(lN)%MWxStoich
303     ! Gas phase enthalpy change from energy equation derivation.
304                       IF(M == 0) THEN
305                          r_HORg = r_HORg + CALC_H(T_g(IJK),0,N) * lRate
306     ! Solid phase enthalpy change from energy equation derivation.
307                       ELSE
308                          r_HORs(M) = r_HORs(M) + CALC_H(T_s(IJK,M),M,N)*lRate
309                       ENDIF
310                    ENDDO
311     
312     ! Complete the skew-symettric for enthalpy transfer with mass transfer
313                    DO M=1, lMMAX
314                        DO L=0, M-1
315                         r_RxH(M,L) = - r_RxH(L,M)
316                       ENDDO
317                    ENDDO
318     ! Apply enthalpy transfer associated with mass transfer to get the
319     ! complete heat of reaction of heat phse for Reaction H.
320                    DO L=0, lMMAX
321                       DO M = 0, lMMAX
322                          IF(L == M) CYCLE
323                          IF(L == 0) THEN
324                             r_HORg = r_HORg - r_RxH(L,M)
325                          ELSE
326                             r_HORs(L) = r_HORs(L) - r_RxH(L,M)
327                          ENDIF
328                       ENDDO
329                    ENDDO
330     ! Apply unit conversion and store heats of reaction in the global array.
331                    HOR_g(IJK) = HOR_g(IJK) + l_2SI*r_HORg
332                    DO M=1,lMMAX
333                       HOR_s(IJK,M) = HOR_s(IJK,M) + l_2SI*r_HORs(M)
334                    ENDDO
335                 ELSE
336     ! User-defined heat of reaction.
337                    HOR_g(IJK) = HOR_g(IJK) + Reaction(H)%HoR(0) * RATES(H)
338                    DO M=1, lMMAX
339                       HOR_s(IJK,M) = HOR_s(IJK,M) +                        &
340                          Reaction(H)%HoR(M) * RATES(H)
341                    ENDDO
342                 ENDIF
343              ENDIF  ! ENERGY_EQ
344     
345     ! Update rate of interphase mass transfer.
346     !---------------------------------------------------------------------//
347               DO LM=1, (DIMENSION_LM+DIMENSION_M-1)
348                  R_PHASE(IJK,LM) = R_PHASE(IJK,LM) +                       &
349                     RATES(H) * Reaction(H)%rPHASE(LM)
350               ENDDO
351           ENDDO RXN_LP ! Loop over reactions.
352     
353     
354     ! Calculate the toal rate of formation and consumption for each species.
355     !---------------------------------------------------------------------//
356           IF(SPECIES_EQ(0)) THEN
357              SUM_R_G(IJK) = SUM( R_gp(IJK,:NMAX(0)) -                      &
358                 ROX_gc(IJK,:NMAX(0))*X_g(IJK,:NMAX(0)))
359           ELSE
360              DO H=1, NO_OF_RXNS
361                 IF(Reaction(H)%nPhases <= 0) CYCLE
362                 DO M=1, lMMAX
363                    LM = 1 + ((M-1)*M)/2
364                    SUM_R_G(IJK) = SUM_R_G(IJK) +                           &
365                       RATES(H) * Reaction(H)%rPHASE(LM)
366                 ENDDO
367              ENDDO
368           ENDIF
369     
370           DO M=1, lMMAX
371              IF(SPECIES_EQ(M)) THEN
372                 SUM_R_S(IJK,M) = SUM( R_sp(IJK,M,:NMAX(M)) -               &
373                    RoX_sc(IJK,M,:NMAX(M))*X_s(IJK,M,:NMAX(M)))
374              ELSE
375                 DO H=1, NO_OF_RXNS
376                    IF(Reaction(H)%nPhases <= 0) CYCLE
377                    DO L=0, M-1
378                       LM = 1 + L + ((M-1)*M)/2
379                       SUM_R_S(IJK,M) = SUM_R_S(IJK,M) -                    &
380                          RATES(H) * Reaction(H)%rPHASE(LM)
381                    ENDDO
382                    DO L=M+1, lMMAX
383                       LM = 1 + M + ((L-1)*L)/2
384                       SUM_R_S(IJK,M) = SUM_R_S(IJK,M) +                    &
385                          RATES(H) * Reaction(H)%rPHASE(LM)
386                    ENDDO
387                 ENDDO
388              ENDIF
389           ENDDO
390     
391     
392           ENDIF  ! Fluid_At(IJK)
393           END DO ! IJK
394     
395           RETURN
396     
397           END SUBROUTINE RRATES0
398