File: RELATIVE:/../../../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     ! Model B momentum equation
17           use run, only: MODEL_B
18     ! Particle volume.
19           use discretelement, only: PVOL
20     ! Gas pressure force by fluid cell
21           use discretelement, only: P_FORCE
22     ! Particle drag force
23           use discretelement, only: DRAG_FC
24     ! Flag for 3D simulatoins.
25           use geometry, only: DO_K
26     ! Loop bounds for fluid grid
27           USE compar, only: IJKSTART3, IJKEND3
28     ! Flags for cyclic BC with pressure drop 
29           use geometry, only: CYCLIC_X_PD, CYCLIC_Y_PD, CYCLIC_Z_PD
30     ! Specified pressure drop
31           use bc, only: DELP_X, DELP_Y, DELP_Z
32     ! Domain length
33           use geometry, only: XLENGTH, YLENGTH, ZLENGTH
34     ! Gas phase pressure
35           use fldvar, only: P_G
36     
37           use discretelement, only: MAX_PIP, PIJK, DES_EXPLICITLY_COUPLED
38           use particle_filter, only: FILTER_CELL
39           use particle_filter, only: FILTER_WEIGHT
40           use particle_filter, only: DES_INTERP_ON
41     
42           use functions, only: FLUID_AT
43     
44           use functions, only: IS_NONEXISTENT
45           use functions, only: IS_ENTERING, IS_ENTERING_GHOST
46           use functions, only: IS_EXITING, IS_EXITING_GHOST
47     
48     ! Global Parameters:
49     !---------------------------------------------------------------------//
50     ! Double precision values.
51           use param1, only: ZERO
52     
53           implicit none
54     
55     ! Loop counters: Particle, fluid cell, neighbor cells
56           INTEGER :: NP, IJK, LC
57     ! Interpolation weight
58           DOUBLE PRECISION :: WEIGHT
59     ! Interpolated gas phase quanties.
60           DOUBLE PRECISION :: lPF(3)
61     ! Loop bound for
62           INTEGER :: LP_BND
63     ! mean pressure gradient for the case of periodic boundaries
64           DOUBLE PRECISION :: cPG(3)
65     !......................................................................!
66     
67     ! Calculate the gas phase pressure gradient. (dP/dx)
68           CALL CALC_GRAD_DES(P_G, P_FORCE)
69     
70     
71     ! Add in cyclic BC pressure drop.
72           cPG(1) = merge(DELP_X/XLENGTH, ZERO, CYCLIC_X_PD)
73           cPG(2) = merge(DELP_Y/YLENGTH, ZERO, CYCLIC_Y_PD)
74           cPG(3) = merge(DELP_Z/ZLENGTH, ZERO, CYCLIC_Z_PD)
75     
76           DO IJK=IJKSTART3, IJKEND3
77              P_FORCE(:,IJK) = cPG - P_FORCE(:,IJK)
78           ENDDO
79     
80     
81           IF(DES_EXPLICITLY_COUPLED .AND. .NOT.MODEL_B) THEN
82     
83     ! Loop bounds for interpolation.
84              LP_BND = merge(27,9,DO_K)
85     
86     ! Calculate the gas phase forces acting on each particle.
87     
88     !$omp  parallel do default(none) &
89     !$omp  private(NP,lPF,lc,ijk,weight) &
90     !$omp  shared(MAX_PIP,DES_INTERP_ON,LP_BND,p_force,drag_fc, &
91     !$omp     filter_cell,filter_weight,pijk,pvol)
92              DO NP=1,MAX_PIP
93     
94                 IF(IS_NONEXISTENT(NP) .or.                              &
95                    IS_ENTERING(NP) .or. IS_ENTERING_GHOST(NP) .or.      &
96                    IS_EXITING(NP)  .or. IS_EXITING_GHOST(NP)) CYCLE
97     
98                 IF(.NOT.FLUID_AT(PIJK(NP,4))) CYCLE
99     
100                 IF(DES_INTERP_ON) THEN
101                    lPF = ZERO
102                    DO LC=1,LP_BND
103                       IJK = FILTER_CELL(LC,NP)
104                       WEIGHT = FILTER_WEIGHT(LC,NP)
105                       lPF = lPF + P_FORCE(:,IJK)*WEIGHT
106                    ENDDO
107                 ELSE
108                    lPF = P_FORCE(:,PIJK(NP,4))
109                 ENDIF
110     
111     ! Include gas pressure and gas-solids drag
112                 DRAG_FC(:,NP) = DRAG_FC(:,NP) + lPF*PVOL(NP)
113     
114              ENDDO
115           ENDIF
116     
117           RETURN
118           END SUBROUTINE CALC_PG_GRAD
119