File: N:\mfix\model\des\rxns_gs_des1.f

1     #include "version.inc"
2     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
3     !                                                                      !
4     !  Module name: CONV_GS_DES1                                           !
5     !  Author: J.Musser: 16-Jun-10                                         !
6     !                                                                      !
7     !  Purpose: This routine is called from the DISCRETE side to calculate !
8     !  the gas-particle convective heat transfer.                          !
9     !                                                                      !
10     !  Comments: Explicitly coupled simulations use a stored convective    !
11     !  heat transfer coefficient. Otherwise, the convective heat transfer  !
12     !  coeff is calculated every time step and the total interphase energy !
13     !  transfered is 'stored' and used explictly in the gas phase. The     !
14     !  latter conserves all energy
15     !                                                                      !
16     !  REF: Zhou, Yu, and Zulli, "Particle scale study of heat transfer in !
17     !       packed and bubbling fluidized beds," AIChE Journal, Vol. 55,   !
18     !       no 4, pp 868-884, 2009.                                        !
19     !                                                                      !
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
21           SUBROUTINE RXNS_GS_DES1
22     
23           use constant, only: Pi
24     ! Flag: The fluid and discrete solids are explicitly coupled.
25           use discretelement, only: DES_EXPLICITLY_COUPLED
26           use discretelement, only: DTSOLID
27           use discretelement, only: PIJK
28           use discretelement, only: MAX_PIP
29     
30           use particle_filter, only: DES_INTERP_ON
31     
32           use particle_filter, only: FILTER_CELL
33           use particle_filter, only: FILTER_WEIGHT
34           use particle_filter, only: FILTER_SIZE
35     
36           use des_rxns, only: DES_R_gp, DES_R_gc
37           use des_rxns, only: DES_R_PHASE, DES_SUM_R_g
38           use des_rxns, only: DES_HOR_G
39     
40           use physprop, only: NMAX
41     
42           use functions, only: FLUID_AT
43           use functions, only: IS_NORMAL
44     
45           use param1, only: ZERO
46           use des_rxns, only: DES_R_s
47           use des_thermo, only: RXNS_Qs
48     
49           use param1, only: DIMENSION_LM
50           use param, only: DIMENSION_M
51     
52           IMPLICIT NONE
53     
54           INTEGER :: IJK, LC, NP
55     
56     
57     ! Local gas phase values.
58           DOUBLE PRECISION :: lRgp(NMAX(0)) ! Rate of species production
59           DOUBLE PRECISION :: lRgc(NMAX(0)) ! Rate of species consumption
60           DOUBLE PRECISION :: lHoRg         ! Heat of reaction
61           DOUBLE PRECISION :: lSUMRg        ! lSUMRg
62     
63     ! Interphase mass transfer
64           DOUBLE PRECISION :: lRPhase(DIMENSION_LM+DIMENSION_M-1)
65     
66           DOUBLE PRECISION :: WEIGHT
67     
68           DO NP=1,MAX_PIP
69              IF(.NOT.IS_NORMAL(NP)) CYCLE
70     
71     ! Avoid convection calculations in cells without fluid (cut-cell)
72              IF(.NOT.FLUID_AT(PIJK(NP,4))) THEN
73                 DES_R_s(NP,:) = ZERO
74                 RXNS_Qs(NP) = ZERO
75     
76     ! No additional calculations are needed for explictly coupled
77              ELSEIF(.NOT.DES_EXPLICITLY_COUPLED) THEN
78     
79     ! Calculate the heat transfer coefficient.
80                 CALL CALC_RRATES_DES(NP, lRgp, lRgc, lRPhase, lHoRg, lSUMRg)
81     
82     ! Integrate over solids time step and store in global array. This
83     ! needs updated when interpolation is reintroduced into thermo code.
84     !---------------------------------------------------------------------//
85     ! Store the gas phase source terms.
86                 IF(DES_INTERP_ON) THEN
87                    DO LC=1,FILTER_SIZE
88                       IJK = FILTER_CELL(LC,NP)
89                       WEIGHT = DTSOLID*FILTER_WEIGHT(LC,NP)
90     
91                       DES_R_gp(IJK,:) = DES_R_gp(IJK,:)+lRgp*WEIGHT
92                       DES_R_gc(IJK,:) = DES_R_gc(IJK,:)+lRgc*WEIGHT
93                       DES_R_PHASE(IJK,:) = DES_R_PHASE(IJK,:)+lRPhase*WEIGHT
94                       DES_HOR_G(IJK) = DES_HOR_G(IJK) + lHoRg*WEIGHT
95                       DES_SUM_R_g(IJK) = DES_SUM_R_g(IJK) + lSUMRg*WEIGHT
96                    ENDDO
97                 ELSE
98                    IJK = PIJK(NP,4)
99                    WEIGHT = DTSOLID
100     
101                    DES_R_gp(IJK,:) = DES_R_gp(IJK,:) + lRgp*WEIGHT
102                    DES_R_gc(IJK,:) = DES_R_gc(IJK,:) + lRgc*WEIGHT
103                    DES_R_PHASE(IJK,:) = DES_R_PHASE(IJK,:) + lRPhase*WEIGHT
104                    DES_HOR_G(IJK) = DES_HOR_G(IJK) + lHoRg*WEIGHT
105                    DES_SUM_R_g(IJK) = DES_SUM_R_g(IJK) + lSUMRg*WEIGHT
106                 ENDIF
107              ENDIF
108           ENDDO
109     
110     ! Note that MPI sync is managed at the end of des_time_march for
111     ! non-explicitly coupled cases that use interpolation.
112     
113           RETURN
114           END SUBROUTINE RXNS_GS_DES1
115     
116     
117     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
118     !                                                                      !
119     !  Subroutine: RXNS_GS_GAS1                                            !
120     !  Author: J.Musser                                   Date: 21-NOV-14  !
121     !                                                                      !
122     !  Purpose: This routine is called from the CONTINUUM. It calculates   !
123     !  the scalar cell center drag force acting on the fluid using         !
124     !  interpolated values for the gas velocity and volume fraction. The   !
125     !  The resulting sources are interpolated back to the fluid grid.      !
126     !                                                                      !
127     !  NOTE: The loop over particles includes ghost particles so that MPI  !
128     !  communications are needed to distribute overlapping force between   !
129     !  neighboring grid cells. This is possible because only cells "owned" !
130     !  by the current process will have non-zero weights.                  !
131     !                                                                      !
132     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
133           SUBROUTINE RXNS_GS_GAS1
134     
135     ! Size of particle array on this process.
136           use discretelement, only: MAX_PIP
137     ! Flag to use interpolation
138           use particle_filter, only: DES_INTERP_ON
139     ! Interpolation cells and weights
140           use particle_filter, only: FILTER_CELL,FILTER_WEIGHT,FILTER_SIZE
141     ! IJK of fluid cell containing particles center
142           use discretelement, only: PIJK
143     ! Gas phase mass, species, and energy equation sources
144           use des_rxns, only: DES_R_gp, DES_R_gc, DES_SUM_R_g
145           use des_rxns, only: DES_R_PHASE, DES_HOR_g
146     ! Number of species for each phase
147           use physprop, only: NMAX
148     ! Funtion for identifying fluid cells and normal particles.
149           use functions, only: FLUID_AT
150           use functions, only: IS_NORMAL
151     ! MPI function for collecting interpolated data from ghost cells.
152           use sendrecvnode, only: DES_COLLECT_gDATA
153     ! MPI wrapper for halo exchange.
154           use sendrecv, only: SEND_RECV
155     ! Fluid time step size.
156           use run, only: DT
157     ! Flag to use stiff chemistry solver
158           use stiff_chem, only: stiff_chemistry
159     ! Routine to sync stiff solver results
160           use stiff_chem, only: FINALIZE_STIFF_SOLVER
161     
162     ! Global Parameters:
163     !---------------------------------------------------------------------//
164     ! Double precision values.
165           use param1, only: ZERO, ONE
166           use param1, only: DIMENSION_LM
167           use param, only: DIMENSION_M
168     
169           IMPLICIT NONE
170     
171     ! Loop counters: Particle, fluid cell, neighbor cells
172           INTEGER :: NP, IJK, LC
173     ! Loop bound for filter
174           DOUBLE PRECISION :: GAMMAxTp
175     ! Local gas phase values.
176           DOUBLE PRECISION :: lRgp(NMAX(0)) ! Rate of species production
177           DOUBLE PRECISION :: lRgc(NMAX(0)) ! Rate of species consumption
178           DOUBLE PRECISION :: lHoRg         ! Heat of reaction
179           DOUBLE PRECISION :: lSUMRg        ! lSUMRg
180     
181     ! Interphase mass transfer
182           DOUBLE PRECISION :: lRPhase(DIMENSION_LM+DIMENSION_M-1)
183     
184           DOUBLE PRECISION :: WEIGHT
185     
186     ! Initialize fluid cell values.
187           DES_R_gp = ZERO
188           DES_R_gc = ZERO
189           DES_R_PHASE = ZERO
190           DES_HOR_G = ZERO
191           DES_SUM_R_g = ZERO
192     
193     ! Calculate the gas phase forces acting on each particle.
194     
195           DO NP=1,MAX_PIP
196     
197     ! Only calculate chemical reactions for normal particles that are inside
198     ! fluid cells. Ex: Skip ghost particles or particles that are in a cut-
199     ! cell dead space.
200              IF(.NOT.IS_NORMAL(NP)) CYCLE
201              IF(.NOT.FLUID_AT(PIJK(NP,4))) CYCLE
202     
203              IJK = PIJK(NP,4)
204     
205     ! Calculate the rates of species formation/consumption.
206              CALL CALC_RRATES_DES(NP, lRgp, lRgc, lRPhase, lHoRg, lSUMRg)
207     
208     ! Directly integrate fluid cell reaction sources
209              IF(STIFF_CHEMISTRY) THEN
210                 CALL DES_STIFF_CHEM(IJK, DT, lRgp, lRgc, lHoRg, lSUMRg)
211     
212     ! Store the gas phase source terms.
213              ELSEIF(DES_INTERP_ON) THEN
214                 DO LC=1,FILTER_SIZE
215                    IJK = FILTER_CELL(LC,NP)
216                    WEIGHT = FILTER_WEIGHT(LC,NP)
217     
218                    DES_R_gp(IJK,:) = DES_R_gp(IJK,:)+lRgp*WEIGHT
219                    DES_R_gc(IJK,:) = DES_R_gc(IJK,:)+lRgc*WEIGHT
220                    DES_R_PHASE(IJK,:) = DES_R_PHASE(IJK,:)+lRPhase*WEIGHT
221                    DES_HOR_G(IJK) = DES_HOR_G(IJK) + lHoRg*WEIGHT
222                    DES_SUM_R_g(IJK) = DES_SUM_R_g(IJK) + lSUMRg*WEIGHT
223                 ENDDO
224              ELSE
225                 DES_R_gp(IJK,:) = DES_R_gp(IJK,:) + lRgp
226                 DES_R_gc(IJK,:) = DES_R_gc(IJK,:) + lRgc
227                 DES_R_PHASE(IJK,:) = DES_R_PHASE(IJK,:) + lRPhase
228                 DES_HOR_G(IJK) = DES_HOR_G(IJK) + lHoRg
229                 DES_SUM_R_g(IJK) = DES_SUM_R_g(IJK) + lSUMRg
230              ENDIF
231           ENDDO
232     
233     ! Make the necessary send/recv calls for field vars
234           IF(STIFF_CHEMISTRY) THEN
235              CALL FINALIZE_STIFF_SOLVER
236           ELSE
237     ! Add in data stored in ghost cells from interpolation. This call must
238     ! preceed the SEND_RECV to avoid overwriting ghost cell data.
239              IF(DES_INTERP_ON) THEN
240                 CALL DES_COLLECT_gDATA(DES_R_gp)
241                 CALL DES_COLLECT_gDATA(DES_R_gc)
242                 CALL DES_COLLECT_gDATA(DES_R_PHASE)
243                 CALL DES_COLLECT_gDATA(DES_HOR_g)
244                 CALL DES_COLLECT_gDATA(DES_SUM_R_g)
245              ENDIF
246     
247     ! Update the species mass sources in ghost layers.
248              CALL SEND_RECV(DES_R_gp, 2)
249              CALL SEND_RECV(DES_R_gc, 2)
250              CALL SEND_RECV(DES_R_PHASE, 2)
251              CALL SEND_RECV(DES_HOR_g, 2)
252              CALL SEND_RECV(DES_SUM_R_g, 2)
253           ENDIF
254     
255           RETURN
256           END SUBROUTINE RXNS_GS_GAS1
257     
258     
259     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
260     !                                                                      !
261     !  Subroutine: DES_STIFF_CHEM                                          !
262     !  Author: J.Musser                                   Date:  7-FEB-16  !
263     !                                                                      !
264     !  Purpose: This routine updates a fluid cell by directly integrating  !
265     !  the source terms due to reaction. Note that this routine is called  !
266     !  for each particle to prevent over consumption of reactants.         !
267     !                                                                      !
268     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
269           SUBROUTINE DES_STIFF_CHEM(pIJK, pDT, pRgp, pRgc, pHoRg, pSUMRg)
270     
271           USE constant, only: GAS_CONST
272           USE fldvar, only: EP_g, RO_g, ROP_g, T_g, X_g, P_g
273           USE geometry, only: VOL
274           USE param1, only: ZERO, SMALL_NUMBER, ONE
275           USE physprop, only: C_pg, MW_g, MW_MIX_g
276           use toleranc, only: ZERO_X_gs
277           use physprop, only: NMAX
278     
279           use compar, only: myPE
280           IMPLICIT NONE
281     
282     ! Passed variables
283     !---------------------------------------------------------------------//
284           INTEGER, INTENT(IN) :: pIJK   ! fluid cell index
285           DOUBLE PRECISION, INTENT(IN) :: pDT           ! fluid time step
286           DOUBLE PRECISION, INTENT(IN) :: pRgp(NMAX(0)) ! Rate of production
287           DOUBLE PRECISION, INTENT(IN) :: pRgc(NMAX(0)) ! Rate of consumption
288           DOUBLE PRECISION, INTENT(IN) :: pHoRg         ! Heat of reaction
289           DOUBLE PRECISION, INTENT(IN) :: pSUMRg        ! Net mass change
290     
291     ! Local variables
292     !---------------------------------------------------------------------//
293           DOUBLE PRECISION :: lRg(NMAX(0)), sumXg
294           DOUBLE PRECISION :: lDToVOL
295           INTEGER :: N
296     !......................................................................!
297     
298     ! Time step divided by cell volume.
299           lDToVOL = pDT/VOL(pIJK)
300     ! Net change in species mass
301           lRg = pRgp - pRgc
302     
303     ! Gas phase density: ROP_g
304           ROP_G(pIJK) = ROP_G(pIJK) + lDToVOL*pSUMRg
305     
306     ! Temperature: T_g
307           T_g(pIJK) = T_g(pIJK) - lDToVOL*(pHORg/(ROP_g(pIJK)*C_pg(pIJK)))
308     
309     ! Species mass fractions: X_g
310           DO N=1,NMAX(0)
311              X_g(pIJK,N) = X_g(pIJK,N) + lDToVOL * &
312                 (lRg(N) - X_g(pIJK,N)*pSUMRg)/(ROP_g(pIJK))
313           ENDDO
314     
315     ! Cleanup over/under shoots.
316           X_g(pIJK,:) = min(max(X_g(pIJK,:), ZERO), ONE)
317           sumXg = sum(X_g(pIJK,:))
318           X_g(pIJK,:) = X_g(pIJK,:)/sumXg
319     
320     ! Gas phase bulk density is updated within the stiff solver (lVar(1)).
321     ! Now that the gas phase volume fraction is updated, the gas phase
322     ! density can be backed out. RO_g * EP_g = ROP_g
323           IF(EP_g(pIJK) > SMALL_NUMBER) THEN
324              RO_g(pIJK) = ROP_g(pIJK) / EP_g(pIJK)
325           ELSE
326     ! This case shouldn't happen, however HUGE is used to aid in tracking
327     ! errors should this somehow become and issue.
328              RO_g(pIJK) = huge(0.0)
329           ENDIF
330     
331     ! Calculate the mixture molecular weight.
332           MW_MIX_G(pIJK) = sum(X_G(pIJK,1:NMAX(0))/MW_g(1:NMAX(0)))
333           MW_MIX_G(pIJK) = ONE/MW_MIX_G(pIJK)
334     
335     ! Calculate the gas phase pressure.
336           P_G(pIJK) = (RO_G(pIJK)*GAS_CONST*T_G(pIJK))/MW_MIX_G(pIJK)
337     
338           RETURN
339           END SUBROUTINE DES_STIFF_CHEM
340