File: N:\mfix\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_NORMAL
45     
46     ! Global Parameters:
47     !---------------------------------------------------------------------//
48     ! Double precision values.
49           use param1, only: ZERO
50     
51           implicit none
52     
53     ! Loop counters: Particle, fluid cell, neighbor cells
54           INTEGER :: NP, IJK, LC
55     ! Interpolation weight
56           DOUBLE PRECISION :: WEIGHT
57     ! Interpolated gas phase quanties.
58           DOUBLE PRECISION :: lPF(3)
59     ! Loop bound for
60           INTEGER :: LP_BND
61     ! mean pressure gradient for the case of periodic boundaries
62           DOUBLE PRECISION :: cPG(3)
63     !......................................................................!
64     
65     ! Calculate the gas phase pressure gradient. (dP/dx)
66           CALL CALC_GRAD_DES(P_G, P_FORCE)
67     
68     
69     ! Add in cyclic BC pressure drop.
70           cPG(1) = merge(DELP_X/XLENGTH, ZERO, CYCLIC_X_PD)
71           cPG(2) = merge(DELP_Y/YLENGTH, ZERO, CYCLIC_Y_PD)
72           cPG(3) = merge(DELP_Z/ZLENGTH, ZERO, CYCLIC_Z_PD)
73     
74           DO IJK=IJKSTART3, IJKEND3
75              P_FORCE(:,IJK) = cPG - P_FORCE(:,IJK)
76           ENDDO
77     
78     
79           IF(DES_EXPLICITLY_COUPLED .AND. .NOT.MODEL_B) THEN
80     
81     ! Loop bounds for interpolation.
82              LP_BND = merge(27,9,DO_K)
83     
84     ! Calculate the gas phase forces acting on each particle.
85     
86     !$omp  parallel do default(none) &
87     !$omp  private(NP,lPF,lc,ijk,weight) &
88     !$omp  shared(MAX_PIP,DES_INTERP_ON,LP_BND,p_force,drag_fc, &
89     !$omp     filter_cell,filter_weight,pijk,pvol)
90              DO NP=1,MAX_PIP
91     
92                 IF(IS_NORMAL(NP)) THEN
93                    IF(.NOT.FLUID_AT(PIJK(NP,4))) CYCLE
94     
95                    IF(DES_INTERP_ON) THEN
96                       lPF = ZERO
97                       DO LC=1,LP_BND
98                          IJK = FILTER_CELL(LC,NP)
99                          WEIGHT = FILTER_WEIGHT(LC,NP)
100                          lPF = lPF + P_FORCE(:,IJK)*WEIGHT
101                       ENDDO
102                    ELSE
103                       lPF = P_FORCE(:,PIJK(NP,4))
104                    ENDIF
105     
106     ! Include gas pressure and gas-solids drag
107                    DRAG_FC(NP,:) = DRAG_FC(NP,:) + lPF*PVOL(NP)
108                 ENDIF
109     
110              ENDDO
111           ENDIF
112     
113           RETURN
114           END SUBROUTINE CALC_PG_GRAD
115