File: N:\mfix\model\qmomk\qmomk_gas_drag.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: QMOMK_GAS_DRAG(A_M, B_M, IER, UV, VV, WV)              C
4     !  Purpose: QMOMK - Accounting for the equal and opposite drag force   C
5     !           on gas due to particles by introducing the drag            C
6     !           as a source term. Face centered                            C
7     !                                                                      C
8     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
9     !  Reviewer: Alberto Passalacqua                      Date:            C
10     !            Re-used to develop QMOMK_GAS_DRAG from the                C
11     !            original DES_GAS_DRAG                                     C
12     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
13     
14           SUBROUTINE QMOMK_GAS_DRAG(A_M, B_M, IER, UV, VV, WV)
15     !-----------------------------------------------
16     !     M o d u l e s
17     !-----------------------------------------------
18     
19           USE compar, only: ijkstart3, ijkend3
20           USE discretelement, only: dimn
21           USE fun_avg, only: avg_x, avg_y, avg_z
22           USE functions, only: fluid_at, east_of, north_of, top_of
23           USE indices, only: i_of, j_of, k_of
24           USE geometry, only: vol_u, vol_v, vol_w
25           USE param
26           USE param1
27           USE physprop, only: mmax
28           USE qmom_kinetic_equation, only: qmomk_nn
29           USE qmom_kinetic_equation
30     
31           IMPLICIT NONE
32     !-----------------------------------------------
33     !     G l o b a l   P a r a m e t e r s
34     !-----------------------------------------------
35     !-----------------------------------------------
36     !     D u m m y   A r g u m e n t s
37     !-----------------------------------------------
38     !
39     !     Error index
40           INTEGER          IER
41     !
42     !     Indices
43           INTEGER          IJK
44     !
45     !     Phase index
46           INTEGER          M, UV, VV, WV
47     !     Averaging Factor
48           DOUBLE PRECISION :: AVG_FACTOR
49     !     Grid indices
50           INTEGER          I,J,K,IN
51     !
52           DOUBLE PRECISION USFCM, VSFCM, WSFCM
53           DOUBLE PRECISION A_M(DIMENSION_3, -3:3, 0:DIMENSION_M)
54           DOUBLE PRECISION B_M(DIMENSION_3, 0:DIMENSION_M)
55           DOUBLE PRECISION tmp_A, tmp_B
56     
57     !-----------------------------------------------
58           AVG_FACTOR = 0.25D0*(DIMN-2) + 0.5D0*(3-DIMN)
59           IF(UV.EQ.1) THEN
60              DO M = 1, MMAX
61                 DO IJK = IJKSTART3, IJKEND3
62                    IF(FLUID_AT(IJK)) THEN
63                       I = I_OF(IJK)
64                       J = J_OF(IJK)
65                       K = K_OF(IJK)
66     
67                       DO IN = 1, QMOMK_NN
68                         USFCM = AVG_X(QMOMK_U1(IN,IJK,M),QMOMK_U1(IN,EAST_OF(IJK),M),I_OF(IJK))
69                         tmp_A = -AVG_X(QMOMK_F_GS(IN,IJK,M),QMOMK_F_GS(IN,EAST_OF(IJK),M),I)*VOL_U(IJK)
70                         tmp_B = tmp_A*USFCM
71     
72                         A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A
73                         B_M(IJK,0) = B_M(IJK,0) + tmp_B
74                       END DO
75                    END IF
76                 END DO
77              END DO
78     
79           ELSE IF(VV.EQ.1) THEN
80     
81              DO M = 1, MMAX
82                 DO IJK = IJKSTART3, IJKEND3
83                    IF(FLUID_AT(IJK)) THEN
84                       I = I_OF(IJK)
85                       J = J_OF(IJK)
86                       K = K_OF(IJK)
87     
88                       DO IN = 1, QMOMK_NN
89                          VSFCM = AVG_Y(QMOMK_V1(IN,IJK,M), QMOMK_V1(IN,NORTH_OF(IJK),M),J_OF(IJK))
90                          tmp_A = -AVG_Y(QMOMK_F_GS(IN,IJK,M),QMOMK_F_GS(IN,NORTH_OF(IJK),M),J)*VOL_V(IJK)
91                          tmp_B = tmp_A*VSFCM
92     
93                          A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A
94                          B_M(IJK,0) = B_M(IJK,0) + tmp_B
95                       END DO
96                    END IF
97     
98                 END DO
99              END DO
100              close(900)
101     
102           ELSE IF(WV.EQ.1) THEN
103              DO M = 1, MMAX
104                 DO IJK = IJKSTART3, IJKEND3
105                    IF(FLUID_AT(IJK)) THEN
106                       I = I_OF(IJK)
107                       J = J_OF(IJK)
108                       K = K_OF(IJK)
109     
110                       DO IN = 1, QMOMK_NN
111                          WSFCM = AVG_Z(QMOMK_W1(IN,IJK,M),QMOMK_W1(IN,TOP_OF(IJK),M),K_OF(IJK))
112                          tmp_A = -AVG_Z(QMOMK_F_GS(IN,IJK,M),QMOMK_F_GS(IN,TOP_OF(IJK),M),K)*VOL_W(IJK)
113                          tmp_B = tmp_A*WSFCM
114     
115                          A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A
116                          B_M(IJK,0) = B_M(IJK,0) + tmp_B
117                       END DO
118                    END IF
119                 END DO
120              END DO
121     
122           END IF
123     
124           RETURN
125           END SUBROUTINE QMOMK_GAS_DRAG
126