File: /nfs/home/0/users/jenkins/mfix.git/model/drag_ss.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 SUBROUTINE DRAG_SS(L, M)
21
22
23
24
25 USE compar
26 USE constant
27 USE discretelement
28 USE drag
29 USE fldvar
30 USE fun_avg
31 USE functions
32 USE geometry
33 USE indices
34 USE parallel
35 USE param
36 USE param1
37 USE physprop
38 USE rdf
39 USE sendrecv
40 IMPLICIT NONE
41
42
43
44
45 INTEGER, INTENT(IN) :: L, M
46
47
48
49
50 INTEGER :: I, IJK, IMJK, IJMK, IJKM
51 INTEGER :: CM, DM
52
53
54 INTEGER :: LM
55
56 DOUBLE PRECISION :: USCM, USCL, VSCM, VSCL, WSCM, WSCL
57
58 DOUBLE PRECISION :: VREL
59
60 DOUBLE PRECISION :: D_pm, D_pl
61
62 DOUBLE PRECISION :: RO_M, RO_L
63
64 DOUBLE PRECISION :: G0_ML
65
66 DOUBLE PRECISION :: EPg, EPs
67
68 DOUBLE PRECISION :: EPSoDP
69
70 DOUBLE PRECISION :: lDss
71
72
73 = FUNLM(L,M)
74
75 DO IJK = ijkstart3, ijkend3
76
77 IF (.NOT.WALL_AT(IJK)) THEN
78
79
80
81
82 = I_OF(IJK)
83 IMJK = IM_OF(IJK)
84 IJMK = JM_OF(IJK)
85 IJKM = KM_OF(IJK)
86
87
88 = AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I)
89 VSCL = AVG_Y_N(V_S(IJMK,L),V_S(IJK,L))
90 WSCL = AVG_Z_T(W_S(IJKM,L),W_S(IJK,L))
91
92 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
93 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
94 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
95
96
97 = SQRT((USCL - USCM)**2 + (VSCL - VSCM)**2 + &
98 (WSCL - WSCM)**2)
99
100
101 = D_P(IJK,M)
102 D_PL = D_P(IJK,L)
103 RO_M = RO_S(IJK,M)
104 RO_L = RO_S(IJK,L)
105
106 IF (DES_CONTINUUM_HYBRID) THEN
107
108
109
110
111 = ZERO
112 DO CM = 1, MMAX
113 EPS = EP_s(IJK, CM)
114 EPSoDP = EPSoDP + EPS / D_p(IJK,CM)
115 ENDDO
116 DO DM = 1, DES_MMAX
117 EPS = DES_ROP_S(IJK,DM)/DES_RO_S(DM)
118 EPSoDP = EPSoDP + EPS / DES_D_p0(DM)
119 ENDDO
120 EPg = EP_g(IJK)
121 G0_ML = ONE/EPg + 3.0d0*EPSoDP*D_pM*D_PL / &
122 (EPg*EPg *(D_pM + D_pL))
123 ELSE
124 G0_ML = G_0(IJK,L,M)
125 ENDIF
126
127
128 CALL DRAG_SS_SYAM(lDss,D_PM,D_PL,RO_M,RO_L,G0_ML,VREL)
129
130 F_SS(IJK,LM) = lDss*ROP_S(IJK,M)*ROP_S(IJK,L)
131
132
133
134
135 IF(CLOSE_PACKED(M) .AND. CLOSE_PACKED(L) .AND. &
136 (MMAX >= 2)) F_SS(IJK,LM) = F_SS(IJK,LM) + &
137 SEGREGATION_SLOPE_COEFFICIENT*P_star(IJK)
138
139 ELSE
140
141 (IJK,LM) = ZERO
142
143 ENDIF
144
145 ENDDO
146
147 RETURN
148 END SUBROUTINE DRAG_SS
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167 SUBROUTINE DRAG_SS_SYAM(lDss,D_PM,D_PL,RO_M,RO_L, G0_ML, VREL)
168
169
170
171
172 USE param
173 USE param1
174 USE constant
175 use run
176 IMPLICIT NONE
177
178
179
180
181 DOUBLE PRECISION, INTENT(OUT) :: ldss
182
183 DOUBLE PRECISION, INTENT(IN) :: D_PM
184
185 DOUBLE PRECISION, INTENT(IN) :: D_PL
186
187 DOUBLE PRECISION, INTENT(IN) :: RO_M
188
189 DOUBLE PRECISION, INTENT(IN) :: RO_L
190
191 DOUBLE PRECISION, INTENT(IN) :: G0_ML
192
193 DOUBLE PRECISION, INTENT(IN) :: VREL
194
195
196
197
198 DOUBLE PRECISION :: DPSUM
199
200 DOUBLE PRECISION :: const
201
202
203
204
205
206 = D_PL + D_PM
207
208 const = 3.d0*(ONE + C_E)*(PI/2.d0 + C_F*PI*PI/8.d0)*&
209 DPSUM**2/(2.d0*PI*(RO_L*D_PL**3+RO_M*D_PM**3))
210
211 ldss = const * G0_ML * VREL
212
213 RETURN
214 END SUBROUTINE DRAG_SS_SYAM
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234 SUBROUTINE DRAG_SS_IA(L, M)
235
236
237
238
239 USE compar
240 USE constant
241 USE drag
242 USE fldvar
243 USE functions
244 USE geometry
245 USE indices
246 USE kintheory
247 USE param1
248 USE physprop
249 USE rdf
250 USE sendrecv
251
252 IMPLICIT NONE
253
254
255
256
257 INTEGER, INTENT(IN) :: M, L
258
259
260
261
262 INTEGER :: IJK
263
264
265 INTEGER :: LM
266
267 DOUBLE PRECISION :: D_PM, D_PL
268
269 DOUBLE PRECISION :: DPSUM
270
271 DOUBLE PRECISION :: M_PM, M_PL, MPSUM, DPSUMo2, NU_PL, NU_PM
272 DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R2p_lm, R3p_lm, &
273 R4p_lm, R10p_lm, Bp_lm
274 DOUBLE PRECISION :: Fss_ip, Fnus_ip, FTsM_ip, FTsL_ip, &
275 F_common_term
276
277
278 DO IJK = ijkstart3, ijkend3
279
280 IF (.NOT.WALL_AT(IJK)) THEN
281
282 LM = FUNLM(L,M)
283
284 IF (M == L) THEN
285 F_SS(IJK,LM) = ZERO
286 Fnu_s_ip(IJK,M,L) = ZERO
287 FT_sM_ip(IJK,M,L) = ZERO
288 FT_sL_ip(IJK,M,L) = ZERO
289
290 ELSE
291 D_PM = D_P(IJK,M)
292 D_PL = D_P(IJK,L)
293 DPSUM = D_PL + D_PM
294 M_PM = (Pi/6.d0) * D_PM**3 *RO_S(IJK,M)
295 M_PL = (Pi/6.d0) * D_PL**3 *RO_S(IJK,L)
296 MPSUM = M_PM + M_PL
297 DPSUMo2 = DPSUM/2.d0
298 NU_PM = ROP_S(IJK,M)/M_PM
299 NU_PL = ROP_S(IJK,L)/M_PL
300
301 IF(Theta_m(IJK,M) > ZERO .AND. Theta_m(IJK,L) > ZERO) THEN
302
303 Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/&
304 2.d0
305 Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-Theta_m(IJK,M) ))/&
306 (2.d0*MPSUM)
307 Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+M_PL*Theta_m(IJK,L) ))/&
308 (2.d0*MPSUM*MPSUM)
309
310 R0p_lm = ( 1.d0/( Ap_lm**1.5 * Dp_lm**2.5 ) )+ &
311 ( (15.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**2.5 * Dp_lm**3.5 ) )+&
312 ( (175.d0*(Bp_lm**4))/( 8.d0*Ap_lm**3.5 * Dp_lm**4.5 ) )
313
314 R2p_lm = ( 1.d0/( 2.d0*Ap_lm**1.5 * Dp_lm*Dp_lm ) )+&
315 ( (3.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**3 ) )+&
316 ( (15.d0*Bp_lm**4)/( 2.d0*Ap_lm**3.5 * Dp_lm**4 ) )
317
318 R3p_lm = ( 1.d0/( (Ap_lm**1.5)*(Dp_lm**3.5) ) )+&
319 ( (21.d0*Bp_lm*Bp_lm)/( 2.d0 * Ap_lm**2.5 * Dp_lm**4.5 ) )+&
320 ( (315.d0*Bp_lm**4)/( 8.d0 * Ap_lm**3.5 *Dp_lm**5.5 ) )
321
322 R4p_lm = ( 3.d0/( Ap_lm**2.5 * Dp_lm**3.5 ) )+&
323 ( (35.d0*Bp_lm*Bp_lm)/( 2.d0 * Ap_lm**3.5 * Dp_lm**4.5 ) )+&
324 ( (441.d0*Bp_lm**4)/( 8.d0 * Ap_lm**4.5 * Dp_lm**5.5 ) )
325
326 R10p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**2.5 ) )+&
327 ( (25.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**3.5 * Dp_lm**3.5 ) )+&
328 ( (1225.d0*Bp_lm**4)/( 24.d0* Ap_lm**4.5 * Dp_lm**4.5 ) )
329
330 F_common_term = (DPSUMo2*DPSUMo2/4.d0)*(M_PM*M_PL/MPSUM)*&
331 G_0(IJK,M,L)*(1.d0+C_E)*(M_PM*M_PL)**1.5
332
333
334
335 = F_common_term*NU_PM*NU_PL*DSQRT(PI)*R2p_lm*&
336 (Theta_m(IJK,M)*Theta_m(IJK,L))**2
337
338
339
340 = F_common_term*(PI*DPSUMo2/12.d0)*R0p_lm*&
341 (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
342
343
344
345 = F_common_term*NU_PM*NU_PL*DPSUMo2*PI*&
346 (Theta_m(IJK,M)**1.5 * Theta_m(IJK,L)**2.5) *&
347 ( (-1.5d0/12.d0*R0p_lm)+&
348 Theta_m(IJK,L)/16.d0*( (-M_PM*R10p_lm) - &
349 ((5.d0*M_PL*M_PL*M_PM/(192.d0*MPSUM*MPSUM))*R3p_lm)+&
350 ((5.d0*M_PM*M_PL)/(96.d0*MPSUM)*R4p_lm*Bp_lm) ) )
351
352
353
354 = F_common_term*NU_PM*NU_PL*DPSUMo2*PI*&
355 (Theta_m(IJK,L)**1.5 * Theta_m(IJK,M)**2.5) *&
356 ( (1.5d0/12.d0*R0p_lm)+&
357 Theta_m(IJK,M)/16.d0*( (M_PL*R10p_lm)+&
358 (5.d0*M_PM*M_PM*M_PL/(192.d0*MPSUM*MPSUM)*R3p_lm)+&
359 (5.d0*M_PM*M_PL/(96.d0*MPSUM) *R4p_lm*Bp_lm) ) )
360
361 F_SS(IJK,LM) = Fss_ip
362
363 Fnu_s_ip(IJK,M,L) = Fnus_ip
364
365
366
367
368 (IJK,M,L) = FTsM_ip
369
370 (IJK,M,L) = FTsL_ip
371 ELSE
372 F_SS(IJK,LM) = ZERO
373 Fnu_s_ip(IJK,M,L) = ZERO
374 FT_sM_ip(IJK,M,L) = ZERO
375 FT_sL_ip(IJK,M,L) = ZERO
376 ENDIF
377 ENDIF
378 ENDIF
379 ENDDO
380
381 RETURN
382 END SUBROUTINE DRAG_SS_IA
383
384