File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/calc_external_forces.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21 SUBROUTINE CALC_EXTERNAL_FORCES()
22
23
24
25
26 USE param
27 USE param1
28 USE geometry
29 USE compar
30 USE fldvar
31 USE indices
32 USE ghdtheory
33 USE physprop
34 USE run
35 USE constant
36 USE drag
37 USE bc
38 use scales
39 use bodyforce
40 USE fun_avg
41 USE functions
42 IMPLICIT NONE
43
44
45
46
47 INTEGER IJK, I, J, K
48
49
50 INTEGER IJKE, IJKN, IJKT
51
52
53 INTEGER M ,IM
54
55
56 DOUBLE PRECISION NjC, NjE, NjN, NjT
57
58
59 DOUBLE PRECISION Mj, Vj
60
61
62 DOUBLE PRECISION dragFc, dragFe, dragFn, dragFt
63
64
65 DOUBLE PRECISION PGE, PGN, PGT, SDPx, SDPy, SDPz
66
67
68 DOUBLE PRECISION avgDragx, avgDragy, avgDragz
69
70
71
72
73
74 DO 200 IJK = ijkstart3, ijkend3
75 I = I_OF(IJK)
76 J = J_OF(IJK)
77 K = K_OF(IJK)
78
79
80 IF ( FLUID_AT(IJK) ) THEN
81
82 IJKE = EAST_OF(IJK)
83 IJKN = NORTH_OF(IJK)
84 IJKT = TOP_OF(IJK)
85
86
87 = P_G(IJKE)
88 PGN = P_G(IJKN)
89 PGT = P_G(IJKT)
90 IF (CYCLIC_X_PD) THEN
91 IF (CYCLIC_AT_E(IJK)) PGE = P_G(IJKE) - DELP_X
92 ENDIF
93 IF (CYCLIC_Y_PD) THEN
94 IF (CYCLIC_AT_N(IJK)) PGN = P_G(IJKN) - DELP_Y
95 ENDIF
96 IF (CYCLIC_Z_PD) THEN
97 IF (CYCLIC_AT_T(IJK)) PGT = P_G(IJKT) - DELP_Z
98 ENDIF
99 SDPx = -P_SCALE*(PGE - P_G(IJK)) * oDX_E(I)
100 SDPy = -P_SCALE*(PGN - P_G(IJK)) * oDY_N(J)
101 SDPz = -P_SCALE*(PGT - P_G(IJK)) * (oX_E(I)*oDZ_T(K))
102
103 DO M = 1, SMAX
104 Vj = (PI/6.d0)*D_P(IJK,M)**3
105 = Vj * RO_S(IJK,M)
106
107 = ROP_s(IJK,M) / Mj
108 NjE = ROP_S(IJKE,M) / Mj
109 NjN = ROP_S(IJKN,M) / Mj
110 NjT = ROP_S(IJKT,M) / Mj
111
112
113 = zero
114 dragFe = zero
115 dragFn = zero
116 dragFt = zero
117 if(NjC > zero) dragFc = F_GS(IJK ,M)/NjC
118 if(NjE > zero) dragFe = F_GS(IJKE,M)/NjE
119 if(NjN > zero) dragFn = F_GS(IJKN,M)/NjN
120 if(NjT > zero) dragFt = F_GS(IJKT,M)/NjT
121
122 dragFxflux(IJK,M) = AVG_X(dragFc,dragFe,I) * (U_g(IJK) - U_s(IJK,M))
123 dragFyflux(IJK,M) = AVG_Y(dragFc,dragFn,J) * (V_g(IJK) - V_s(IJK,M))
124 dragFzflux(IJK,M) = AVG_Z(dragFc,dragFt,K) * (W_g(IJK) - W_s(IJK,M))
125
126 dragFx(IJK,M) = AVG_X(dragFc,dragFe,I) * (U_g(IJK))
127 dragFy(IJK,M) = AVG_Y(dragFc,dragFn,J) * (V_g(IJK))
128 dragFz(IJK,M) = AVG_Z(dragFc,dragFt,K) * (W_g(IJK))
129
130 beta_cell_X(IJK,M)=AVG_X(dragFc,dragFe,I)
131 beta_cell_Y(IJK,M)=AVG_Y(dragFc,dragFn,J)
132 beta_cell_Z(IJK,M)=AVG_Z(dragFc,dragFt,K)
133
134 IF(DRAG_TYPE_ENUM == HYS)THEN
135 DO IM=1,SMAX
136 IF(IM /= M)THEN
137
138
139 = zero
140 dragFe = zero
141 dragFn = zero
142 dragFt = zero
143 if(NjC > zero) dragFc = beta_ij(IJK,M,IM) /NjC
144 if(NjE > zero) dragFe = beta_ij(IJKE,M,IM)/NjE
145 if(NjN > zero) dragFn = beta_ij(IJKN,M,IM)/NjN
146 if(NjT > zero) dragFt = beta_ij(IJKT,M,IM)/NjT
147
148 avgDragx = AVG_X(dragFc,dragFe,I)
149 dragFx(IJK,M) = dragFx(IJK,M) - avgDragx*(U_g(IJK))
150 dragFxflux(IJK,M) = dragFxflux(IJK,M) - avgDragx*(U_g(IJK) - U_s(IJK,IM))
151 beta_ij_cell_X(IJK,M,IM)= AVG_X(dragFc,dragFe,I)
152
153 avgDragy = AVG_Y(dragFc,dragFn,J)
154 dragFy(IJK,M) = dragFy(IJK,M) - avgDragy*(V_g(IJK))
155 dragFyflux(IJK,M) = dragFyflux(IJK,M) - avgDragy*(V_g(IJK) - V_s(IJK,IM))
156 beta_ij_cell_Y(IJK,M,IM)=AVG_Y(dragFc,dragFn,J)
157
158 avgDragz = AVG_Z(dragFc,dragFt,K)
159 dragFz(IJK,M) = dragFz(IJK,M) - avgDragz*(W_g(IJK))
160 dragFzflux(IJK,M) = dragFzflux(IJK,M) - avgDragz*(W_g(IJK) - W_s(IJK,IM))
161 beta_ij_cell_Z(IJK,M,IM)=AVG_Z(dragFc,dragFt,K)
162
163 ENDIF
164 ENDDO
165 ENDIF
166
167
168 FiXvel(IJK,M) = (Mj * BFX_S(IJK,M)+dragFx(IJK,M) +Vj*SDPx)
169 FiYvel(IJK,M) = (Mj * BFY_S(IJK,M)+dragFy(IJK,M) +Vj*SDPy)
170 FiZvel(IJK,M) = (Mj * BFZ_S(IJK,M)+dragFz(IJK,M) +Vj*SDPz)
171
172 FiX(IJK,M) = (Mj * BFX_S(IJK,M)+dragFxflux(IJK,M) +Vj*SDPx)
173 FiY(IJK,M) = (Mj * BFY_S(IJK,M)+dragFyflux(IJK,M) +Vj*SDPy)
174 FiZ(IJK,M) = (Mj * BFZ_S(IJK,M)+dragFzflux(IJK,M) +Vj*SDPz)
175
176 FiMinusDragX(IJK,M) = (Mj * BFX_S(IJK,M) + Vj*SDPx)
177 FiMinusDragY(IJK,M) = (Mj * BFY_S(IJK,M) + Vj*SDPy)
178 FiMinusDragZ(IJK,M) = (Mj * BFZ_S(IJK,M) + Vj*SDPz)
179 ENDDO
180 ENDIF
181 CONTINUE
182
183 RETURN
184 END SUBROUTINE CALC_EXTERNAL_FORCES
185