File: /nfs/home/0/users/jenkins/mfix.git/model/des/calc_pg_grad.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CALC_PG_GRAD                                            !
4     !  Purpose: Calculate cell centered pressure force exerted on the      !
5     !           particles in the cell by the gas/fluid phase               !
6     !           Note that P_force is evaluated as -dp/dx                   !
7     !                                                                      !
8     !  Notes: This pressure force only needs to be calculated once during  !
9     !         the DEM loop (at the beginning) since the gas/fluid phase    !
10     !         is essentially static at that point (i.e., gas field is not  !
11     !         updated during DEM loop                                      !
12     !                                                                      !
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
14           SUBROUTINE CALC_PG_GRAD
15     
16           use cutcell, only: CARTESIAN_GRID
17     
18     ! Model B momentum equation
19           use run, only: MODEL_B
20     
21           use discretelement, only: MAX_PIP, PIJK, DES_EXPLICITLY_COUPLED
22           use particle_filter, only: FILTER_CELL
23           use particle_filter, only: FILTER_WEIGHT
24           use particle_filter, only: DES_INTERP_ON
25     
26     ! Flags indicating the state of particle
27           use discretelement, only: PEA
28     ! Particle volume.
29           use discretelement, only: PVOL
30     ! Gas pressure force by fluid cell
31           use discretelement, only: P_FORCE
32     ! Particle drag force
33           use discretelement, only: DRAG_FC
34     ! Flag for 3D simulatoins.
35           use geometry, only: DO_K
36     
37     
38     ! Global Parameters:
39     !---------------------------------------------------------------------//
40     ! Double precision values.
41           use param1, only: ZERO
42     
43     
44           implicit none
45     
46     
47     ! Loop counters: Particle, fluid cell, neighbor cells
48           INTEGER :: NP, IJK, LC
49     ! Interpolation weight
50           DOUBLE PRECISION :: WEIGHT
51     ! Interpolated gas phase quanties.
52           DOUBLE PRECISION :: lPF(3)
53     ! Loop bound for
54           INTEGER :: LP_BND
55     
56     
57           IF(CARTESIAN_GRID) THEN
58              CALL CALC_PG_GRAD_CG
59           ELSE
60              CALL CALC_PG_GRAD_STD
61           ENDIF
62     
63     
64           IF(DES_EXPLICITLY_COUPLED .AND. .NOT.MODEL_B) THEN
65     
66     ! Loop bounds for interpolation.
67              LP_BND = merge(27,9,DO_K)
68     
69     
70     ! Calculate the gas phase forces acting on each particle.
71     
72     !$omp  parallel do default(none) &
73     !$omp              private(NP,lPF,lc,ijk,weight) &
74     !$omp              shared(MAX_PIP,PEA,DES_INTERP_ON,LP_BND,p_force,drag_fc,filter_cell,filter_weight,pijk,pvol)
75              DO NP=1,MAX_PIP
76                 IF(.NOT.PEA(NP,1)) CYCLE
77                 IF(any(PEA(NP,2:3))) CYCLE
78     
79                 IF(DES_INTERP_ON) THEN
80                    lPF = ZERO
81                    DO LC=1,LP_BND
82                       IJK = FILTER_CELL(LC,NP)
83                       WEIGHT = FILTER_WEIGHT(LC,NP)
84                       lPF = lPF + P_FORCE(:,IJK)*WEIGHT
85                    ENDDO
86                 ELSE
87                    lPF = P_FORCE(:,PIJK(NP,4))
88                 ENDIF
89     
90     ! Include gas pressure and gas-solids drag
91                 DRAG_FC(:,NP) = DRAG_FC(:,NP) + lPF*PVOL(NP)
92     
93              ENDDO
94           ENDIF
95     
96           RETURN
97           END SUBROUTINE CALC_PG_GRAD
98     
99     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
100     !                                                                      !
101     !  Subroutine: CALC_PG_GRAD_STD                                        !
102     !  Purpose: Calculate cell centered pressure force exerted on the      !
103     !           particles in the cell by the gas/fluid phase               !
104     !           Note that P_force is evaluated as -dp/dx                   !
105     !                                                                      !
106     !  Notes: This pressure force only needs to be calculated once during  !
107     !         the DEM loop (at the beginning) since the gas/fluid phase    !
108     !         is essentially static at that point (i.e., gas field is not  !
109     !         updated during DEM loop                                      !
110     !                                                                      !
111     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
112           SUBROUTINE CALC_PG_GRAD_STD
113     
114     !-----------------------------------------------
115     ! Modules
116     !-----------------------------------------------
117           USE param
118           USE param1
119           USE physprop
120           USE fldvar
121           USE run
122           USE geometry
123           USE indices
124           USE bc
125           USE compar
126           USE sendrecv
127           USE discretelement
128           USE cutcell
129           USE fun_avg
130           USE functions
131           implicit none
132     !-----------------------------------------------
133     ! Local variables
134     !-----------------------------------------------
135     ! general i, j, k indices
136           INTEGER :: I, J, K, IJK
137           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
138     ! temporary variables used to calculate pressure at scalar cell edge
139           DOUBLE PRECISION :: TEMP1, TEMP2
140     ! mean pressure gradient for the case of periodic boundaries
141           DOUBLE PRECISION :: MPG_CYCLIC(3)
142     !-----------------------------------------------
143     
144     
145           MPG_CYCLIC(1:3) = ZERO
146     
147           IF(CYCLIC_X_PD) MPG_CYCLIC(1) = DELP_X/XLENGTH
148           IF(CYCLIC_Y_PD) MPG_CYCLIC(2) = DELP_Y/YLENGTH
149           IF(CYCLIC_Z_PD.AND.DO_K) MPG_CYCLIC(3) = DELP_Z/ZLENGTH
150     
151           DO IJK = IJKSTART3, IJKEND3
152              P_FORCE(:,IJK) = ZERO
153              IF(.NOT.FLUID_AT(IJK)) CYCLE
154     
155              I = I_OF(IJK)
156              J = J_OF(IJK)
157              K = K_OF(IJK)
158              IMJK = IM_OF(IJK)
159              IPJK = IP_OF(IJK)
160              IJMK = JM_OF(IJK)
161              IJPK = JP_OF(IJK)
162              IJKM = KM_OF(IJK)
163              IJKP = KP_OF(IJK)
164     
165              IF(IMIN1.EQ.IMAX1) THEN
166                 P_FORCE(1,IJK) = MPG_CYCLIC(1)
167              ELSEIF(I.EQ.IMIN1) THEN
168                 TEMP2 = AVG_X(P_G(IJK), P_G(IPJK), I)
169                 TEMP1 = P_G(IJK)
170                 P_FORCE(1,IJK) = 2.d0*(TEMP1-TEMP2)/DX(I)  +  MPG_CYCLIC(1)
171              ELSEIF(I.EQ.IMAX1) THEN
172                 TEMP2 = AVG_X(P_G(IMJK), P_G(IJK), I-1)
173                 TEMP1 = P_G(IJK)
174                 P_FORCE(1,IJK) = 2.d0*(TEMP2 - TEMP1)/DX(I) +  MPG_CYCLIC(1)
175              ELSEIF((I.GT.IMIN1).AND.(I.LT.IMAX1)) THEN
176                 TEMP2 = AVG_X(P_G(IJK),  P_G(IPJK), I)
177                 TEMP1 = AVG_X(P_G(IMJK), P_G(IJK),  I-1)
178                 P_FORCE(1,IJK) = (TEMP1 - TEMP2)/DX(I) + MPG_CYCLIC(1)
179              ENDIF
180     
181              IF(JMIN1.EQ.JMAX1) THEN
182                 P_FORCE(2,IJK) = MPG_CYCLIC(2)
183              ELSEIF(J.EQ.JMIN1) THEN
184                 TEMP2 = AVG_Y(P_G(IJK), P_G(IJPK), J)
185                 TEMP1 = P_G(IJK)
186                 P_FORCE(2,IJK) = 2.d0*(TEMP1 - TEMP2)/DY(J)  + MPG_CYCLIC(2)
187              ELSEIF(J.EQ.JMAX1) THEN
188                 TEMP2 = AVG_Y(P_G(IJMK), P_G(IJK), J-1)
189                 TEMP1 = P_G(IJK)
190                 P_FORCE(2,IJK) = 2.d0*(TEMP2 - TEMP1)/DY(J) + MPG_CYCLIC(2)
191              ELSEIF((J.GT.JMIN1).AND.(J.LT.JMAX1)) THEN
192                 TEMP2 = AVG_Y(P_G(IJK),  P_G(IJPK), J)
193                 TEMP1 = AVG_Y(P_G(IJMK), P_G(IJK),  J-1)
194                 P_FORCE(2,IJK) = (TEMP1 - TEMP2)/DY(J) +  MPG_CYCLIC(2)
195              ENDIF
196     
197              IF(DO_K) THEN
198                 IF(KMIN1.EQ.KMAX1) THEN
199                    P_FORCE(3,IJK) = MPG_CYCLIC(3)
200                 ELSEIF(K.EQ.KMIN1) THEN
201                    TEMP2 = AVG_Z(P_G(IJK), P_G(IJKP), K)
202                    TEMP1 = P_G(IJK)
203                    P_FORCE(3,IJK) = 2.d0*(TEMP1 - TEMP2)/DZ(K) +  MPG_CYCLIC(3)
204                 ELSEIF(K.EQ.KMAX1) THEN
205                    TEMP2 = AVG_Z(P_G(IJKM), P_G(IJK), K-1)
206                    TEMP1 = P_G(IJK)
207                    P_FORCE(3,IJK) = 2.d0*(TEMP2 - TEMP1)/DZ(K) +  MPG_CYCLIC(3)
208                 ELSEIF((K.GT.KMIN1).AND.(K.LT.KMAX1)) THEN
209                    TEMP2 = AVG_Z(P_G(IJK),  P_G(IJKP), K)
210                    TEMP1 = AVG_Z(P_G(IJKM), P_G(IJK),  K-1)
211                    P_FORCE(3,IJK) = (TEMP1 - TEMP2)/DZ(K) +  MPG_CYCLIC(3)
212                 ENDIF
213              ENDIF
214           ENDDO         ! end do loop over ijk
215     
216           RETURN
217           END SUBROUTINE CALC_PG_GRAD_STD
218     
219     
220     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
221     !                                                                      !
222     !  Subroutine: CALC_PG_GRAD_CG                                         !
223     !  Purpose: Calculate cell centered pressure force exerted on the      !
224     !           particles in the cell by the gas/fluid phase               !
225     !           (cut-cell version)                                         !
226     !                                                                      !
227     !  Notes: This pressure force needs to be calculated once in a DEM     !
228     !         time step (at the beginning) since the gas/fluid phase is    !
229     !         not updated (is static) during the DEM portion of the        !
230     !         simulation.                                                  !
231     !                                                                      !
232     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
233           SUBROUTINE CALC_PG_GRAD_CG
234     
235     !-----------------------------------------------
236     ! Modules
237     !-----------------------------------------------
238           USE param
239           USE param1
240           USE parallel
241           USE physprop
242           USE fldvar
243           USE run
244           USE geometry
245           USE indices
246           USE bc
247           USE compar
248           USE sendrecv
249           USE discretelement
250           USE fun_avg
251           USE functions
252           implicit none
253     !-----------------------------------------------
254     ! Local variables
255     !-----------------------------------------------
256     ! general i, j, k indices
257           INTEGER :: I, J, K, IJK
258           INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
259     ! counters
260           INTEGER :: X_COUNT, Y_COUNT, Z_COUNT
261     ! mean pressure gradient for the case of periodic boundaries
262           DOUBLE PRECISION :: MPG_CYCLIC(3)
263     !-----------------------------------------------
264           MPG_CYCLIC(1:3) = ZERO
265     
266           IF(CYCLIC_X_PD) MPG_CYCLIC(1) = DELP_X/XLENGTH
267           IF(CYCLIC_Y_PD) MPG_CYCLIC(2) = DELP_Y/YLENGTH
268           IF(CYCLIC_Z_PD.AND.DO_K) MPG_CYCLIC(3) = DELP_Z/ZLENGTH
269     
270           DO IJK = IJKSTART3, IJKEND3
271              I = I_OF(IJK)
272              J = J_OF(IJK)
273              K = K_OF(IJK)
274              P_FORCE(:,IJK) = ZERO
275     
276              IF(.NOT.FLUID_AT(IJK).OR..NOT.IS_ON_myPE_owns(I, J, K)) CYCLE
277     
278              X_COUNT = 0
279              Y_COUNT = 0
280              Z_COUNT = 0
281     
282              IJKE = IP_OF(IJK)
283              IJKW = IM_OF(IJK)
284              IJKN = JP_OF(IJK)
285              IJKS = JM_OF(IJK)
286              IJKT = KP_OF(IJK)
287              IJKB = KM_OF(IJK)
288     
289              IF(FLUID_AT(IJKE)) THEN
290                 X_COUNT = X_COUNT + 1
291                 P_FORCE(1,IJK) = P_FORCE(1,IJK) +                          &
292                    2.d0*(P_G(IJKE) - P_G(IJK))/(DX(I) + DX(I_OF(IJKE)))
293              ENDIF
294              IF(FLUID_AT(IJKW)) THEN
295                 X_COUNT = X_COUNT + 1
296                 P_FORCE(1,IJK) = P_FORCE(1,IJK) +                          &
297                    2.d0*(P_G(IJK) - P_G(IJKW))/(DX(I) + DX(I_OF(IJKW)))
298              ENDIF
299              X_COUNT = MAX(1, X_COUNT) !to prevent division from zero
300     ! P_FORCE (by convention) is stored as -dp/dx. MPG_CYCLIC is already -dp/dx.
301     ! therefore, P_force is multiplied by "-" for consistency
302              P_FORCE(1,IJK) = MPG_CYCLIC(1) - P_FORCE(1,IJK)/REAL(X_COUNT)
303     
304              IF(FLUID_AT(IJKN)) THEN
305                 Y_COUNT = Y_COUNT + 1
306                 P_FORCE(2,IJK) = P_FORCE(2,IJK) +                          &
307                    2.d0*(P_G(IJKN) - P_G(IJK))/(DY(J) + DY(J_OF(IJKN)))
308              ENDIF
309     
310              IF(FLUID_AT(IJKS)) THEN
311                 Y_COUNT = Y_COUNT + 1
312                 P_FORCE(2,IJK) = P_FORCE(2,IJK) +                          &
313                    2.d0*(P_G(IJK) - P_G(IJKS))/(DY(J) + DY(J_OF(IJKS)))
314              ENDIF
315              Y_COUNT = MAX(1, Y_COUNT) !to prevent division from zero
316              P_FORCE(2,IJK) = MPG_CYCLIC(2) - P_FORCE(2,IJK)/REAL(Y_COUNT)
317     
318              IF(DO_K) THEN
319                 IF(FLUID_AT(IJKT)) THEN
320                    Z_COUNT = Z_COUNT + 1
321                    P_FORCE(3,IJK) = P_FORCE(3,IJK) +                       &
322                       2.d0*(P_G(IJKT) - P_G(IJK))/(DZ(K) + DZ(K_OF(IJKT)))
323                 ENDIF
324                 IF(FLUID_AT(IJKB)) THEN
325                    Z_COUNT = Z_COUNT + 1
326                    P_FORCE(3,IJK) = P_FORCE(3,IJK) +                       &
327                       2.d0*(P_G(IJK) - P_G(IJKB))/(DZ(K) + DZ(K_OF(IJKB)))
328                 ENDIF
329                 Z_COUNT = MAX(1, Z_COUNT) !to prevent division from zero
330                 P_FORCE(3,IJK) = MPG_CYCLIC(3)-P_FORCE(3,IJK)/REAL(Z_COUNT)
331              ENDIF
332     
333           ENDDO
334     
335           RETURN
336           END SUBROUTINE CALC_PG_GRAD_CG
337