File: RELATIVE:/../../../mfix.git/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           USE param
19           USE param1
20           USE parallel
21           USE matrix
22           USE scales
23           USE constant
24           USE physprop
25           USE fldvar
26           USE visc_g
27           USE rxns
28           USE run
29           USE toleranc
30           USE geometry
31           USE indices
32           USE is
33           USE tau_g
34           USE bc
35           USE compar
36           USE sendrecv
37           USE discretelement
38           USE qmom_kinetic_equation
39           USE drag
40           USE fun_avg
41           USE functions
42     
43           IMPLICIT NONE
44     !-----------------------------------------------
45     !     G l o b a l   P a r a m e t e r s
46     !-----------------------------------------------
47     !-----------------------------------------------
48     !     D u m m y   A r g u m e n t s
49     !-----------------------------------------------
50     !
51     !
52     !     Error index
53           INTEGER          IER
54     !
55     !     Indices
56           INTEGER          IJK
57     !
58     !     Phase index
59           INTEGER          M, UV, VV, WV
60     !     Averaging Factor
61           DOUBLE PRECISION :: AVG_FACTOR
62     !     Grid indices
63           INTEGER          I,J,K,IN
64     !
65           DOUBLE PRECISION USFCM, VSFCM, WSFCM
66           DOUBLE PRECISION A_M(DIMENSION_3, -3:3, 0:DIMENSION_M)
67           DOUBLE PRECISION B_M(DIMENSION_3, 0:DIMENSION_M)
68           DOUBLE PRECISION tmp_A, tmp_B
69     
70     !-----------------------------------------------
71           AVG_FACTOR = 0.25D0*(DIMN-2) + 0.5D0*(3-DIMN)
72           IF(UV.EQ.1) THEN
73              DO M = 1, MMAX
74                 DO IJK = IJKSTART3, IJKEND3
75                    IF(FLUID_AT(IJK)) THEN
76                       I = I_OF(IJK)
77                       J = J_OF(IJK)
78                       K = K_OF(IJK)
79     
80                       DO IN = 1, QMOMK_NN
81                         USFCM = AVG_X(QMOMK_U1(IN,IJK,M),QMOMK_U1(IN,EAST_OF(IJK),M),I_OF(IJK))
82                         tmp_A = -AVG_X(QMOMK_F_GS(IN,IJK,M),QMOMK_F_GS(IN,EAST_OF(IJK),M),I)*VOL_U(IJK)
83                         tmp_B = tmp_A*USFCM
84     
85                         A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A
86                         B_M(IJK,0) = B_M(IJK,0) + tmp_B
87                       END DO
88                    END IF
89                 END DO
90              END DO
91     
92     
93           ELSE IF(VV.EQ.1) THEN
94     
95              DO M = 1, MMAX
96                 DO IJK = IJKSTART3, IJKEND3
97                    IF(FLUID_AT(IJK)) THEN
98                       I = I_OF(IJK)
99                       J = J_OF(IJK)
100                       K = K_OF(IJK)
101     
102                       DO IN = 1, QMOMK_NN
103                          VSFCM = AVG_Y(QMOMK_V1(IN,IJK,M), QMOMK_V1(IN,NORTH_OF(IJK),M),J_OF(IJK))
104                          tmp_A = -AVG_Y(QMOMK_F_GS(IN,IJK,M),QMOMK_F_GS(IN,NORTH_OF(IJK),M),J)*VOL_V(IJK)
105                          tmp_B = tmp_A*VSFCM
106     
107                          A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A
108                          B_M(IJK,0) = B_M(IJK,0) + tmp_B
109                       END DO
110                    END IF
111     
112                 END DO
113              END DO
114     90       FORMAT(3(1x,i2),10(2x,E12.5))
115              close(900)
116     
117           ELSE IF(WV.EQ.1) THEN
118              DO M = 1, MMAX
119                 DO IJK = IJKSTART3, IJKEND3
120                    IF(FLUID_AT(IJK)) THEN
121                       I = I_OF(IJK)
122                       J = J_OF(IJK)
123                       K = K_OF(IJK)
124     
125                       DO IN = 1, QMOMK_NN
126                          WSFCM = AVG_Z(QMOMK_W1(IN,IJK,M),QMOMK_W1(IN,TOP_OF(IJK),M),K_OF(IJK))
127                          tmp_A = -AVG_Z(QMOMK_F_GS(IN,IJK,M),QMOMK_F_GS(IN,TOP_OF(IJK),M),K)*VOL_W(IJK)
128                          tmp_B = tmp_A*WSFCM
129     
130                          A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A
131                          B_M(IJK,0) = B_M(IJK,0) + tmp_B
132                       END DO
133                    END IF
134                 END DO
135              END DO
136     
137           END IF
138     
139           RETURN
140           END SUBROUTINE QMOMK_GAS_DRAG
141