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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CALC_GS_DES1                                            !
4     !  Author: J.Musser                                   Date: 21-NOV-14  !
5     !                                                                      !
6     !  Purpose: This routine is called from the DISCRETE side to calculate !
7     !  the gas-based forces acting on each particle using interpolated     !
8     !  values for gas velocity, gas pressure, and gas volume fraction.     !
9     !                                                                      !
10     !  Notes:                                                              !
11     !                                                                      !
12     !   F_gp is obtained from des_drag_gp subroutine is given as:          !
13     !    F_GP = beta*VOL_P/EP_s where VOL_P is the particle volume.        !
14     !                                                                      !
15     !  The drag force on each particle is equal to:                        !
16     !    D_FORCE = beta*VOL_P/EP_s*(Ug - Us) = F_GP *(Ug - Us)             !
17     !                                                                      !
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
19           SUBROUTINE DRAG_GS_DES1
20     
21     ! Modules
22     !---------------------------------------------------------------------//
23     ! Gas phase volume fraction
24           use fldvar, only: EP_G
25     ! Gas phase velocities
26           use fldvar, only: U_G, V_G, W_G
27     ! Function to deterine if a cell contains fluid.
28           use functions, only: FLUID_AT
29           use functions, only: is_normal
30     ! Flag for 3D simulatoins.
31           use geometry, only: DO_K
32     ! Size of particle arrays on this processor.
33           use discretelement, only: MAX_PIP
34     ! Flag: The fluid and discrete solids are explicitly coupled.
35           use discretelement, only: DES_EXPLICITLY_COUPLED
36     ! IJK of fluid cell containing particles center
37           use discretelement, only: PIJK
38     ! Drag force on each particle
39           use discretelement, only: F_GP
40     ! Particle velocity
41           use discretelement, only: DES_VEL_NEW
42     ! Total forces acting on particle
43           use discretelement, only: FC
44     ! Gas pressure force by fluid cell
45           use discretelement, only: P_FORCE
46     ! Particle volume.
47           use discretelement, only: PVOL
48     ! Particle drag force
49           use discretelement, only: DRAG_FC
50     ! Flag for MPPIC runs.
51           use mfix_pic, only: MPPIC
52     ! Flag to use implicit drag for MPPIC
53           use mfix_pic, only: MPPIC_PDRAG_IMPLICIT
54     ! Double precision values.
55           use param1, only: ZERO
56           use param, only: dimension_3
57     ! Flag to use interpolation
58           use particle_filter, only: DES_INTERP_ON
59     ! Interpolation cells and weights
60           use particle_filter, only: FILTER_CELL, FILTER_WEIGHT
61     ! Model B momentum equation
62           use run, only: MODEL_B
63           use param, only: DIMENSION_3
64           use compar
65           IMPLICIT NONE
66     
67     ! Local variables
68     !---------------------------------------------------------------------//
69     !     Cell-center gas velocities.
70           DOUBLE PRECISION :: UGC(DIMENSION_3)
71           DOUBLE PRECISION :: VGC(DIMENSION_3)
72           DOUBLE PRECISION :: WGC(DIMENSION_3)
73     ! Loop counters: Particle, fluid cell, neighbor cells
74           INTEGER :: NP, IJK, LC
75     ! Interpolation weight
76           DOUBLE PRECISION :: WEIGHT
77     ! Interpolated gas phase quantities and particle velocity.
78           DOUBLE PRECISION :: lEPg, VELFP(3), lPF(3), vel_p(3)
79     ! Drag force acting on each particle.
80           DOUBLE PRECISION :: D_FORCE(3)
81     ! Flag for Model A momentum equation
82           LOGICAL :: MODEL_A
83     ! Loop bound for filter
84           INTEGER :: LP_BND
85     !......................................................................!
86     
87     ! Set flag for Model A momentum equation.
88           MODEL_A = .NOT.MODEL_B
89     ! Loop bounds for interpolation.
90           LP_BND = merge(27,9,DO_K)
91     
92     ! Calculate the cell center gas velocities.
93           CALL CALC_CELL_CENTER_GAS_VEL(U_G, V_G, W_G, UGC, VGC, WGC)
94     
95     ! Calculate the gas phase forces acting on each particle.
96     
97     !$omp parallel default(none) private(np,lepg,velfp,ijk,weight,lpf,d_force,vel_p)    &
98     !$omp          shared(max_pip,des_interp_on,lp_bnd,filter_cell,filter_weight, &
99     !$omp          ep_g,pijk,des_vel_new,f_gp,mppic,ugc,vgc,wgc,p_force,    &
100     !$omp          des_explicitly_coupled,drag_fc,mppic_pdrag_implicit,fc,model_a,pvol)
101     !$omp do
102           DO NP=1,MAX_PIP
103              IF(.NOT.IS_NORMAL(NP)) CYCLE
104     
105              lEPG = ZERO
106              VELFP = ZERO
107              lPF = ZERO
108     
109     ! Calculate the gas volume fraction, velocity, and pressure force at
110     ! the particle's position.
111              IF(DES_INTERP_ON) THEN
112                 DO LC=1,LP_BND
113                    IJK = FILTER_CELL(LC,NP)
114                    WEIGHT = FILTER_WEIGHT(LC,NP)
115     ! Gas phase volume fraction.
116                    lEPG = lEPG + EP_G(IJK)*WEIGHT
117     ! Gas phase velocity.
118                    VELFP(1) = VELFP(1) + UGC(IJK)*WEIGHT
119                    VELFP(2) = VELFP(2) + VGC(IJK)*WEIGHT
120                    VELFP(3) = VELFP(3) + WGC(IJK)*WEIGHT
121     ! Gas pressure force.
122                    lPF = lPF + P_FORCE(:,IJK)*WEIGHT
123                 ENDDO
124              ELSE
125                 IJK = PIJK(NP,4)
126                 lEPG = EP_G(IJK)
127                 VELFP(1) = UGC(IJK)
128                 VELFP(2) = VGC(IJK)
129                 VELFP(3) = WGC(IJK)
130                 lPF = P_FORCE(:,IJK)
131              ENDIF
132     
133     ! Avoid drag calculations in cells without fluid (cut-cell)
134              IF(.NOT.FLUID_AT(PIJK(NP,4))) THEN
135                 DRAG_FC(NP,:) = 0.0d0
136                 F_GP(NP) = 0.0d0
137     
138     ! For explicit coupling, use the drag coefficient calculated for the
139     ! gas phase drag calculations.
140              ELSEIF(DES_EXPLICITLY_COUPLED) THEN
141                 DRAG_FC(NP,:) = F_GP(NP)*(VELFP - DES_VEL_NEW(NP,:))
142              ELSE
143     
144     ! Calculate the drag coefficient.
145                 vel_p = DES_VEL_NEW(NP,:)
146                 CALL DES_DRAG_GP(NP, vel_p, VELFP, lEPg)
147     
148     ! Calculate the gas-solids drag force on the particle
149                 IF(MPPIC .AND. MPPIC_PDRAG_IMPLICIT) THEN
150     ! implicit treatment of the drag term for mppic
151                    D_FORCE = F_GP(NP)*VELFP
152                 ELSE
153     ! default case
154                    D_FORCE = F_GP(NP)*(VELFP - DES_VEL_NEW(NP,:))
155                 ENDIF
156     
157     ! Update the contact forces (FC) on the particle to include gas
158     ! pressure and gas-solids drag
159                 FC(NP,:) = FC(NP,:) + D_FORCE(:)
160                 IF(MODEL_A) FC(NP,:) = FC(NP,:) + lPF*PVOL(NP)
161     
162              ENDIF
163     
164           ENDDO
165     !$omp end parallel
166     
167           RETURN
168           END SUBROUTINE DRAG_GS_DES1
169     
170     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
171     !                                                                      !
172     !  Subroutine: DRAG_GS_GAS1                                            !
173     !  Author: J.Musser                                   Date: 21-NOV-14  !
174     !                                                                      !
175     !                                                                      !
176     !  Purpose: This routine is called from the CONTINUUM. It calculates   !
177     !  the scalar cell center drag force acting on the fluid using         !
178     !  interpolated values for the gas velocity and volume fraction. The   !
179     !  The resulting sources are interpolated back to the fluid grid.      !
180     !                                                                      !
181     !  NOTE: The loop over particles includes ghost particles so that MPI  !
182     !  communications are needed to distribute overlapping force between   !
183     !  neighboring grid cells. This is possible because only cells "owned" !
184     !  by the current process will have non-zero weights.                  !
185     !                                                                      !
186     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
187           SUBROUTINE DRAG_GS_GAS1
188     
189     ! Flag: The fluid and discrete solids are explicitly coupled.
190           use discretelement, only: DES_EXPLICITLY_COUPLED
191     ! Size of particle array on this process.
192           use discretelement, only: MAX_PIP
193     ! IJK of fluid cell containing particles center
194           use discretelement, only: PIJK
195     ! Drag force on each particle
196           use discretelement, only: F_GP
197     ! Particle velocity
198           use discretelement, only: DES_VEL_NEW
199     ! Contribution to gas momentum equation due to drag
200           use discretelement, only: DRAG_BM
201     ! Scalar cell center total drag force
202           use discretelement, only: F_GDS
203     ! Gas phase volume fraction
204           use fldvar, only: EP_G
205     ! Gas phase velocities
206           use fldvar, only: U_G, V_G, W_G
207           use fldvar, only: U_GO, V_GO, W_GO
208           use functions, only: FLUID_AT
209           use functions, only: IS_NORMAL
210     ! Volume of scalar cell.
211           use geometry, only: VOL
212     ! Flag for 3D simulatoins.
213           use geometry, only: DO_K
214     ! Flag for MPPIC runs
215           use mfix_pic, only: MPPIC
216     ! Statical weight of each MPPIC parcel
217           use mfix_pic, only: DES_STAT_WT
218     ! Double precision values.
219           use param1, only: ZERO, ONE
220     ! Flag to use interpolation
221           use particle_filter, only: DES_INTERP_ON
222     ! Interpolation cells and weights
223           use particle_filter, only: FILTER_CELL, FILTER_WEIGHT
224     ! Array size for fluid variables
225           use param, only: DIMENSION_3
226     ! MPI function for collecting interpolated data from ghost cells.
227           use sendrecvnode, only: DES_COLLECT_gDATA
228     ! MPI wrapper for halo exchange.
229           use sendrecv, only: SEND_RECV
230     
231           IMPLICIT NONE
232     
233     ! Local variables
234     !---------------------------------------------------------------------//
235     !     Cell-center gas velocities.
236           DOUBLE PRECISION :: UGC(DIMENSION_3)
237           DOUBLE PRECISION :: VGC(DIMENSION_3)
238           DOUBLE PRECISION :: WGC(DIMENSION_3)
239     ! Loop counters: Particle, fluid cell, neighbor cells
240           INTEGER :: NP, IJK, LC
241     ! Interpolation weight
242           DOUBLE PRECISION :: WEIGHT
243     ! Interpolated gas phase quantities and particle velocity
244           DOUBLE PRECISION :: lEPg, VELFP(3), vel_p(3)
245     ! Loop bound for filter
246           INTEGER :: LP_BND
247     ! Drag force (intermediate calculation)
248           DOUBLE PRECISION :: lFORCE
249     ! Drag sources for fluid (intermediate calculation)
250           DOUBLE PRECISION :: lDRAG_BM(3)
251     !......................................................................!
252     
253     ! Initialize fluid cell values.
254           F_GDS = ZERO
255           DRAG_BM = ZERO
256     
257     ! Loop bounds for interpolation.
258           LP_BND = merge(27,9,DO_K)
259     
260     ! Calculate the cell center gas velocities.
261           CALL CALC_CELL_CENTER_GAS_VEL(U_G, V_G, W_G, UGC, VGC, WGC)
262     
263     ! Calculate the gas phase forces acting on each particle.
264     
265     !$omp parallel default(none) private(np,lepg,velfp,ijk,weight,ldrag_bm,lforce,vel_p) &
266     !$omp          shared(max_pip,des_interp_on,lp_bnd,filter_cell,filter_weight, &
267     !$omp          ep_g,pijk,des_vel_new,f_gp,vol,des_stat_wt,mppic,drag_bm,f_gds,ugc,vgc,wgc)
268     !$omp do
269           DO NP=1,MAX_PIP
270     
271              IF(.NOT.IS_NORMAL(NP)) CYCLE
272              IF(.NOT.FLUID_AT(PIJK(NP,4))) CYCLE
273     
274              lEPG = ZERO
275              VELFP = ZERO
276     
277     ! Calculate the gas volume fraction, velocity, and at the
278     ! particle's position.
279              IF(DES_INTERP_ON) THEN
280                 DO LC=1,LP_BND
281                    IJK = FILTER_CELL(LC,NP)
282                    WEIGHT = FILTER_WEIGHT(LC,NP)
283     ! Gas phase volume fraction.
284                    lEPG = lEPG + EP_G(IJK)*WEIGHT
285     ! Gas phase velocity.
286                    VELFP(1) = VELFP(1) + UGC(IJK)*WEIGHT
287                    VELFP(2) = VELFP(2) + VGC(IJK)*WEIGHT
288                    VELFP(3) = VELFP(3) + WGC(IJK)*WEIGHT
289                 ENDDO
290              ELSE
291                 IJK = PIJK(NP,4)
292                 lEPG = EP_G(IJK)
293                 VELFP(1) = UGC(IJK)
294                 VELFP(2) = VGC(IJK)
295                 VELFP(3) = WGC(IJK)
296              ENDIF
297     
298     ! This avoids FP exceptions for some ghost particles.
299              IF(lEPg == ZERO) lEPG = EP_g(PIJK(NP,4))
300     
301     ! Calculate drag coefficient
302              vel_p = DES_VEL_NEW(NP,:)
303              CALL DES_DRAG_GP(NP, vel_p, VELFP, lEPg)
304     
305              lFORCE = F_GP(NP)
306              IF(MPPIC) lFORCE = lFORCE*DES_STAT_WT(NP)
307     
308              lDRAG_BM = lFORCE*DES_VEL_NEW(NP,:)
309     
310              IF(DES_INTERP_ON) THEN
311                 DO LC=1,LP_BND
312                    IJK = FILTER_CELL(LC,NP)
313                    WEIGHT = FILTER_WEIGHT(LC,NP)/VOL(IJK)
314     
315     !$omp atomic
316                    DRAG_BM(IJK,1) = DRAG_BM(IJK,1) + lDRAG_BM(1)*WEIGHT
317     !$omp atomic
318                    DRAG_BM(IJK,2) = DRAG_BM(IJK,2) + lDRAG_BM(2)*WEIGHT
319     !$omp atomic
320                    DRAG_BM(IJK,3) = DRAG_BM(IJK,3) + lDRAG_BM(3)*WEIGHT
321     !$omp atomic
322                    F_GDS(IJK) = F_GDS(IJK) + lFORCE*WEIGHT
323                 ENDDO
324              ELSE
325                 IJK = PIJK(NP,4)
326                 WEIGHT = ONE/VOL(IJK)
327     
328     !$omp atomic
329                 DRAG_BM(IJK,1) = DRAG_BM(IJK,1) + lDRAG_BM(1)*WEIGHT
330     !$omp atomic
331                 DRAG_BM(IJK,2) = DRAG_BM(IJK,2) + lDRAG_BM(2)*WEIGHT
332     !$omp atomic
333                 DRAG_BM(IJK,3) = DRAG_BM(IJK,3) + lDRAG_BM(3)*WEIGHT
334     
335     !$omp atomic
336                 F_GDS(IJK) = F_GDS(IJK) + lFORCE*WEIGHT
337              ENDIF
338     
339           ENDDO
340     !$omp end parallel
341     
342     
343     ! Add in data stored in ghost cells from interpolation. This call must
344     ! preceed the SEND_RECV to avoid overwriting ghost cell data.
345           IF(DES_INTERP_ON) THEN
346              CALL DES_COLLECT_gDATA(F_GDS)
347              CALL DES_COLLECT_gDATA(DRAG_BM)
348           ENDIF
349     
350     ! Update the drag force and sources in ghost layers.
351           CALL SEND_RECV(F_GDS, 2)
352           CALL SEND_RECV(DRAG_BM, 2)
353     
354           RETURN
355           END SUBROUTINE DRAG_GS_GAS1
356     
357     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
358     !                                                                      !
359     !  Subroutine: CALC_CELL_CENTER_GAS_VEL                                !
360     !  Author: J.Musser                                   Date: 07-NOV-14  !
361     !                                                                      !
362     !  Purpose: Calculate the scalar cell center gas velocity. This code   !
363     !  is common to the DEM and GAS calls for non-interpolated drag        !
364     !  routines.                                                           !
365     !                                                                      !
366     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
367           SUBROUTINE CALC_CELL_CENTER_GAS_VEL(lUg, lVg, lWg, UGC, VGC, WGC)
368     
369     ! Modules
370     !---------------------------------------------------------------------//
371     ! Fluid grid loop bounds.
372           use compar, only: IJKStart3, IJKEnd3
373     ! Flags and correction factors for cut momentum cells.
374           use cutcell, only: CUT_U_TREATMENT_AT, THETA_UE, THETA_UE_BAR
375           use cutcell, only: CUT_V_TREATMENT_AT, THETA_VN, THETA_VN_BAR
376           use cutcell, only: CUT_W_TREATMENT_AT, THETA_WT, THETA_WT_BAR
377     ! Functions to average momentum to scalar cell center.
378           use fun_avg, only: AVG_X_E, AVG_Y_N, AVG_Z_T
379     ! Functions to lookup adjacent cells by index.
380           use functions, only: IM_OF, JM_OF, KM_OF
381     ! Function to deterine if a cell contains fluid.
382           use functions, only: FLUID_AT
383     ! Flag for 3D simulatoins.
384           use geometry, only: DO_K
385           use indices, only: I_OF
386     ! Double precision parameters.
387           use param, only: DIMENSION_3
388           use param1, only: ZERO
389     
390           IMPLICIT NONE
391     
392     ! Dummy arguments
393     !---------------------------------------------------------------------//
394           DOUBLE PRECISION, INTENT(IN) :: lUg(DIMENSION_3)
395           DOUBLE PRECISION, INTENT(IN) :: lVg(DIMENSION_3)
396           DOUBLE PRECISION, INTENT(IN) :: lWg(DIMENSION_3)
397           DOUBLE PRECISION, INTENT(OUT) :: UGC(DIMENSION_3)
398           DOUBLE PRECISION, INTENT(OUT) :: VGC(DIMENSION_3)
399           DOUBLE PRECISION, INTENT(OUT) :: WGC(DIMENSION_3)
400     ! Local variables:
401     !---------------------------------------------------------------------//
402     ! Indices of adjacent cells
403           INTEGER :: IJK, IMJK, IJMK, IJKM
404     !......................................................................!
405     
406     
407     ! Calculate the cell center gas velocity components.
408           DO IJK=IJKSTART3, IJKEND3
409              IF(FLUID_AT(IJK)) THEN
410                 IMJK = IM_OF(IJK)
411                 IF(CUT_U_TREATMENT_AT(IMJK)) THEN
412                    UGC(IJK) = (THETA_UE_BAR(IMJK)*lUG(IMJK) +              &
413                       THETA_UE(IMJK)*lUg(IJK))
414                 ELSE
415                    UGC(IJK) = AVG_X_E(lUG(IMJK),lUG(IJK),I_OF(IJK))
416                 ENDIF
417     
418                 IJMK = JM_OF(IJK)
419                 IF(CUT_V_TREATMENT_AT(IJMK)) THEN
420                    VGC(IJK) = (THETA_VN_BAR(IJMK)*lVG(IJMK) +              &
421                       THETA_VN(IJMK)*lVg(IJK))
422                 ELSE
423                    VGC(IJK) = AVG_Y_N(lVg(IJMK),lVg(IJK))
424                 ENDIF
425     
426                 IF(DO_K) THEN
427                    IJKM = KM_OF(IJK)
428                    IF(CUT_W_TREATMENT_AT(IJKM)) THEN
429                       WGC(IJK) = (THETA_WT_BAR(IJKM)*lWg(IJKM) +           &
430                          THETA_WT(IJKM)* lWg(IJK))
431                    ELSE
432                       WGC(IJK) = AVG_Z_T(lWg(IJKM),lWg(IJK))
433                    ENDIF
434                 ELSE
435                    WGC(IJK) = ZERO
436                 ENDIF
437              ELSE
438                 UGC(IJK) = ZERO
439                 VGC(IJK) = ZERO
440                 WGC(IJK) = ZERO
441              ENDIF
442           ENDDO
443     
444           RETURN
445           END SUBROUTINE CALC_CELL_CENTER_GAS_VEL
446