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