File: N:\mfix\model\des\calc_grad_des.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_GRAD_DES(PHI, DEL_PHI)
15     
16     ! Flag for cut-cells
17           use cutcell, only: CARTESIAN_GRID
18     
19           use param, only: DIMENSION_3
20     
21           IMPLICIT NONE
22     
23     ! Dummy Arguments:
24     !---------------------------------------------------------------------//
25           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
26           DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,DIMENSION_3)
27     !......................................................................!
28     
29           IF(CARTESIAN_GRID) THEN
30              CALL CALC_GRAD_DES_CG(PHI, DEL_PHI)
31           ELSE
32              CALL CALC_GRAD_DES_STD(PHI, DEL_PHI)
33           ENDIF
34     
35           RETURN
36           END SUBROUTINE CALC_GRAD_DES
37     
38     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
39     !                                                                      !
40     !  Subroutine: CALC_GRAD_DES_STD                                       !
41     !                                                                      !
42     !  Purpose: Calculate cell centered gradient (DEL_PHI) of scalar PHI.  !
43     !           This routine calculates the average scalar value at the    !
44     !           cell faces to calculate the graident across the cell.      !
45     !                                                                      !
46     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
47           SUBROUTINE CALC_GRAD_DES_STD(PHI, DEL_PHI)
48     
49     ! Modules
50     !-----------------------------------------------
51     
52           use geometry, only: DO_K
53           use geometry, only: oDX, oDY, oDZ
54           USE geometry, only: IMIN1, JMIN1, KMIN1
55           USE geometry, only: IMAX1, JMAX1, KMAX1
56           USE indices, only: I_OF, J_OF, K_OF
57           USE compar, only: IJKSTART3, IJKEND3
58     
59           use functions, only: FLUID_AT
60           use functions, only: IM_OF, JM_OF, KM_OF
61           use functions, only: IP_OF, JP_OF, KP_OF
62           use fun_avg, only: AVG_X, AVG_Y, AVG_Z
63     
64           USE param, only: DIMENSION_3
65           USE param1, only: ZERO
66     
67           IMPLICIT NONE
68     
69     ! Dummy Arguments:
70     !---------------------------------------------------------------------//
71           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
72           DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,DIMENSION_3)
73     
74     ! Local variables
75     !---------------------------------------------------------------------//
76     ! general i, j, k indices
77           INTEGER :: I, J, K, IJK
78           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
79     !......................................................................!
80     
81     
82           DO IJK = IJKSTART3, IJKEND3
83              DEL_PHI(:,IJK) = ZERO
84              IF(.NOT.FLUID_AT(IJK)) CYCLE
85     
86              I = I_OF(IJK)
87              IMJK = IM_OF(IJK)
88              IPJK = IP_OF(IJK)
89     
90              IF(IMIN1 == IMAX1) THEN
91              ELSEIF((I>IMIN1).AND.(I<IMAX1)) THEN
92                 DEL_PHI(1,IJK) = oDX(I)*(AVG_X(PHI(IJK),PHI(IPJK),I) -     &
93                    AVG_X(PHI(IMJK),PHI(IJK),I-1))
94              ELSEIF(I == IMIN1) THEN
95                 DEL_PHI(1,IJK) = 2.0d0*oDX(I) *                            &
96                    (AVG_X(PHI(IJK),PHI(IPJK),I) -  PHI(IJK))
97              ELSEIF(I == IMAX1) THEN
98                 DEL_PHI(1,IJK) = 2.0d0*oDX(I) *                            &
99                    (PHI(IJK) - AVG_X(PHI(IMJK), PHI(IJK), I-1))
100              ELSE
101                 DEL_PHI(1,IJK) = ZERO
102              ENDIF
103     
104     
105              J = J_OF(IJK)
106              IJMK = JM_OF(IJK)
107              IJPK = JP_OF(IJK)
108     
109              IF(JMIN1 == JMAX1) THEN
110              ELSEIF((J>JMIN1) .AND. (J<JMAX1)) THEN
111                 DEL_PHI(2,IJK) = oDY(J)*(AVG_Y(PHI(IJK),PHI(IJPK),J) -     &
112                    AVG_Y(PHI(IJMK),PHI(IJK),J-1))
113              ELSEIF(J == JMIN1) THEN
114                 DEL_PHI(2,IJK) = 2.0d0*oDY(J) *                            &
115                    (AVG_Y(PHI(IJK),PHI(IJPK),J) - PHI(IJK))
116              ELSEIF(J == JMAX1) THEN
117                 DEL_PHI(2,IJK) = 2.0d0*oDY(J) *                            &
118                    (PHI(IJK)- AVG_Y(PHI(IJMK),PHI(IJK),J-1))
119              ELSE
120                 DEL_PHI(2,IJK) = ZERO
121              ENDIF
122     
123              IF(DO_K) THEN
124     
125                 K = K_OF(IJK)
126                 IJKM = KM_OF(IJK)
127                 IJKP = KP_OF(IJK)
128     
129                 IF(KMIN1 == KMAX1) THEN
130                 ELSEIF((K>KMIN1) .AND. (K<KMAX1)) THEN
131                    DEL_PHI(3,IJK) = oDZ(K)*(AVG_Z(PHI(IJK),PHI(IJKP),K) -  &
132                       AVG_Z(PHI(IJKM),PHI(IJK),K-1))
133                 ELSEIF(K == KMIN1) THEN
134                    DEL_PHI(3,IJK) = 2.0d0*oDZ(K) *                         &
135                       (AVG_Z(PHI(IJK),PHI(IJKP),K) - PHI(IJK))
136                 ELSEIF(K == KMAX1) THEN
137                    DEL_PHI(3,IJK) = 2.0d0*oDZ(K) *                         &
138                       (PHI(IJK) - AVG_Z(PHI(IJKM),PHI(IJK),K-1))
139                 ENDIF
140              ENDIF
141           ENDDO
142     
143           RETURN
144           END SUBROUTINE CALC_GRAD_DES_STD
145     
146     
147     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
148     !                                                                      !
149     !  Subroutine: CALC_GRAD_DES_CG                                        !
150     !                                                                      !
151     !  Purpose: Calculate cell centered gradient (DEL_PHI) of scalar PHI.  !
152     !           This routine calculates the graidents at the cell face     !
153     !           and averages the face values to the cell center.           !
154     !                                                                      !
155     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
156           SUBROUTINE CALC_GRAD_DES_CG(PHI, DEL_PHI)
157     
158     ! Modules
159     !---------------------------------------------------------------------//
160           use geometry, only: DO_K
161           use geometry, only: DX, DY, DZ
162           use indices, only: I_OF, J_OF, K_OF
163           use compar, only: IJKSTART3, IJKEND3
164     
165           use functions, only: FLUID_AT
166           use functions, only: IM_OF, JM_OF, KM_OF
167           use functions, only: IP_OF, JP_OF, KP_OF
168     
169           use param, only: DIMENSION_3
170           use param1, only: ZERO
171     
172           IMPLICIT NONE
173     
174     ! Dummy Arguments:
175     !---------------------------------------------------------------------//
176           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
177           DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,DIMENSION_3)
178     
179     ! Local variables
180     !---------------------------------------------------------------------//
181     ! General i, j, k indices
182           INTEGER :: I, J, K, IJK
183           INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
184     ! Counter
185           DOUBLE PRECISION :: dLC
186     !......................................................................!
187     
188     
189           DO IJK = IJKSTART3, IJKEND3
190     
191              I = I_OF(IJK)
192              J = J_OF(IJK)
193              K = K_OF(IJK)
194              DEL_PHI(:,IJK) = ZERO
195     
196              IF(.NOT.FLUID_AT(IJK)) CYCLE
197     
198              dLC = 0.0d0
199              IJKE = IP_OF(IJK)
200              IJKW = IM_OF(IJK)
201     
202              IF(FLUID_AT(IJKE)) THEN
203                 DEL_PHI(1,IJK) = DEL_PHI(1,IJK) +                          &
204                    2.0d0*(PHI(IJKE) - PHI(IJK))/(DX(I) + DX(I_OF(IJKE)))
205                 dLC = dLC + 1.0d0
206              ENDIF
207              IF(FLUID_AT(IJKW)) THEN
208                 DEL_PHI(1,IJK) = DEL_PHI(1,IJK) +                          &
209                    2.0d0*(PHI(IJK) - PHI(IJKW))/(DX(I) + DX(I_OF(IJKW)))
210                 dLC = dLC + 1.0d0
211              ENDIF
212              DEL_PHI(1,IJK) = DEL_PHI(1,IJK)/max(1.0d0,dLC)
213     
214     
215              dLC = 0.0d0
216              IJKN = JP_OF(IJK)
217              IJKS = JM_OF(IJK)
218     
219              IF(FLUID_AT(IJKN)) THEN
220                 DEL_PHI(2,IJK) = DEL_PHI(2,IJK) +                          &
221                    2.0d0*(PHI(IJKN) - PHI(IJK))/(DY(J) + DY(J_OF(IJKN)))
222                 dLC = dLC + 1.0d0
223              ENDIF
224     
225              IF(FLUID_AT(IJKS)) THEN
226                 DEL_PHI(2,IJK) = DEL_PHI(2,IJK) +                          &
227                    2.d0*(PHI(IJK) - PHI(IJKS))/(DY(J) + DY(J_OF(IJKS)))
228                 dLC = dLC + 1.0d0
229              ENDIF
230              DEL_PHI(2,IJK) = DEL_PHI(2,IJK)/max(1.0d0, dLC)
231     
232              IF(DO_K) THEN
233     
234                 dLC = 0.0d0
235                 IJKT = KP_OF(IJK)
236                 IJKB = KM_OF(IJK)
237     
238                 IF(FLUID_AT(IJKT)) THEN
239                    DEL_PHI(3,IJK) = DEL_PHI(3,IJK) +                       &
240                       2.0d0*(PHI(IJKT) - PHI(IJK))/(DZ(K) + DZ(K_OF(IJKT)))
241                    dLC = dLC + 1.0d0
242                 ENDIF
243                 IF(FLUID_AT(IJKB)) THEN
244                    DEL_PHI(3,IJK) = DEL_PHI(3,IJK) +                       &
245                       2.0d0*(PHI(IJK) - PHI(IJKB))/(DZ(K) + DZ(K_OF(IJKB)))
246                    dLC = dLC + 1.0d0
247                 ENDIF
248                 DEL_PHI(3,IJK) = DEL_PHI(3,IJK)/max(1.0d0, dLC)
249              ENDIF
250           ENDDO
251     
252           RETURN
253           END SUBROUTINE CALC_GRAD_DES_CG
254