File: RELATIVE:/../../../mfix.git/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((I>IMIN1).AND.(I<IMAX1)) THEN
91                 DEL_PHI(1,IJK) = oDX(I)*(AVG_X(PHI(IJK),PHI(IPJK),I) -     &
92                    AVG_X(PHI(IMJK),PHI(IJK),I-1))
93              ELSEIF(I == IMIN1) THEN
94                 DEL_PHI(1,IJK) = 2.0d0*oDX(I) *                            &
95                    (AVG_X(PHI(IJK),PHI(IPJK),I) -  PHI(IJK))
96              ELSEIF(I == IMAX1) THEN
97                 DEL_PHI(1,IJK) = 2.0d0*oDX(I) *                            &
98                    (PHI(IJK) - AVG_X(PHI(IMJK), PHI(IJK), I-1))
99              ELSE
100                 DEL_PHI(1,IJK) = ZERO
101              ENDIF
102     
103     
104              J = J_OF(IJK)
105              IJMK = JM_OF(IJK)
106              IJPK = JP_OF(IJK)
107     
108              IF((J>JMIN1) .AND. (J<JMAX1)) THEN
109                 DEL_PHI(2,IJK) = oDY(J)*(AVG_Y(PHI(IJK),PHI(IJPK),J) -     &
110                    AVG_Y(PHI(IJMK),PHI(IJK),J-1))
111              ELSEIF(J == JMIN1) THEN
112                 DEL_PHI(2,IJK) = 2.0d0*oDY(J) *                            &
113                    (AVG_Y(PHI(IJK),PHI(IJPK),J) - PHI(IJK))
114              ELSEIF(J == JMAX1) THEN
115                 DEL_PHI(2,IJK) = 2.0d0*oDY(J) *                            &
116                    (PHI(IJK)- AVG_Y(PHI(IJMK),PHI(IJK),J-1))
117              ELSE
118                 DEL_PHI(2,IJK) = ZERO 
119              ENDIF
120     
121              IF(DO_K) THEN
122     
123                 K = K_OF(IJK)
124                 IJKM = KM_OF(IJK)
125                 IJKP = KP_OF(IJK)
126     
127                 IF((K>KMIN1) .AND. (K<KMAX1)) THEN
128                    DEL_PHI(3,IJK) = oDZ(K)*(AVG_Z(PHI(IJK),PHI(IJKP),K) -  &
129                       AVG_Z(PHI(IJKM),PHI(IJK),K-1))
130                 ELSEIF(K == KMIN1) THEN
131                    DEL_PHI(3,IJK) = 2.0d0*oDZ(K) *                         &
132                       (AVG_Z(PHI(IJK),PHI(IJKP),K) - PHI(IJK))
133                 ELSEIF(K == KMAX1) THEN
134                    DEL_PHI(3,IJK) = 2.0d0*oDZ(K) *                         &
135                       (PHI(IJK) - AVG_Z(PHI(IJKM),PHI(IJK),K-1))
136                 ENDIF
137              ENDIF
138           ENDDO
139     
140           RETURN
141           END SUBROUTINE CALC_GRAD_DES_STD
142     
143     
144     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
145     !                                                                      !
146     !  Subroutine: CALC_GRAD_DES_CG                                        !
147     !                                                                      !
148     !  Purpose: Calculate cell centered gradient (DEL_PHI) of scalar PHI.  !
149     !           This routine calculates the graidents at the cell face     !
150     !           and averages the face values to the cell center.           !
151     !                                                                      !
152     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
153           SUBROUTINE CALC_GRAD_DES_CG(PHI, DEL_PHI)
154     
155     ! Modules
156     !---------------------------------------------------------------------//
157           use geometry, only: DO_K
158           use geometry, only: DX, DY, DZ
159           use indices, only: I_OF, J_OF, K_OF
160           use compar, only: IJKSTART3, IJKEND3
161     
162           use functions, only: FLUID_AT
163           use functions, only: IM_OF, JM_OF, KM_OF
164           use functions, only: IP_OF, JP_OF, KP_OF
165     
166           use param, only: DIMENSION_3
167           use param1, only: ZERO
168     
169           IMPLICIT NONE
170     
171     ! Dummy Arguments:
172     !---------------------------------------------------------------------//
173           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
174           DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,DIMENSION_3)
175     
176     ! Local variables
177     !---------------------------------------------------------------------//
178     ! General i, j, k indices
179           INTEGER :: I, J, K, IJK
180           INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
181     ! Counter
182           DOUBLE PRECISION :: dLC
183     !......................................................................!
184     
185     
186           DO IJK = IJKSTART3, IJKEND3
187     
188              I = I_OF(IJK)
189              J = J_OF(IJK)
190              K = K_OF(IJK)
191              DEL_PHI(:,IJK) = ZERO
192     
193              IF(.NOT.FLUID_AT(IJK)) CYCLE
194     
195              dLC = 0.0d0
196              IJKE = IP_OF(IJK)
197              IJKW = IM_OF(IJK)
198     
199              IF(FLUID_AT(IJKE)) THEN
200                 DEL_PHI(1,IJK) = DEL_PHI(1,IJK) +                          &
201                    2.0d0*(PHI(IJKE) - PHI(IJK))/(DX(I) + DX(I_OF(IJKE)))
202                 dLC = dLC + 1.0d0
203              ENDIF
204              IF(FLUID_AT(IJKW)) THEN
205                 DEL_PHI(1,IJK) = DEL_PHI(1,IJK) +                          &
206                    2.0d0*(PHI(IJK) - PHI(IJKW))/(DX(I) + DX(I_OF(IJKW)))
207                 dLC = dLC + 1.0d0
208              ENDIF
209              DEL_PHI(1,IJK) = DEL_PHI(1,IJK)/max(1.0d0,dLC)
210     
211     
212              dLC = 0.0d0
213              IJKN = JP_OF(IJK)
214              IJKS = JM_OF(IJK)
215     
216              IF(FLUID_AT(IJKN)) THEN
217                 DEL_PHI(2,IJK) = DEL_PHI(2,IJK) +                          &
218                    2.0d0*(PHI(IJKN) - PHI(IJK))/(DY(J) + DY(J_OF(IJKN)))
219                 dLC = dLC + 1.0d0
220              ENDIF
221     
222              IF(FLUID_AT(IJKS)) THEN
223                 DEL_PHI(2,IJK) = DEL_PHI(2,IJK) +                          &
224                    2.d0*(PHI(IJK) - PHI(IJKS))/(DY(J) + DY(J_OF(IJKS)))
225                 dLC = dLC + 1.0d0
226              ENDIF
227              DEL_PHI(2,IJK) = DEL_PHI(2,IJK)/max(1.0d0, dLC)
228     
229              IF(DO_K) THEN
230     
231                 dLC = 0.0d0
232                 IJKT = KP_OF(IJK)
233                 IJKB = KM_OF(IJK)
234     
235                 IF(FLUID_AT(IJKT)) THEN
236                    DEL_PHI(3,IJK) = DEL_PHI(3,IJK) +                       &
237                       2.0d0*(PHI(IJKT) - PHI(IJK))/(DZ(K) + DZ(K_OF(IJKT)))
238                    dLC = dLC + 1.0d0
239                 ENDIF
240                 IF(FLUID_AT(IJKB)) THEN
241                    DEL_PHI(3,IJK) = DEL_PHI(3,IJK) +                       &
242                       2.0d0*(PHI(IJK) - PHI(IJKB))/(DZ(K) + DZ(K_OF(IJKB)))
243                    dLC = dLC + 1.0d0
244                 ENDIF
245                 DEL_PHI(3,IJK) = DEL_PHI(3,IJK)/max(1.0d0, dLC)
246              ENDIF
247           ENDDO
248     
249           RETURN
250           END SUBROUTINE CALC_GRAD_DES_CG
251