File: /nfs/home/0/users/jenkins/mfix.git/model/kintheory_energy_dissipation_ss.f
1
2
3
4
5
6
7
8
9
10
11
12 SUBROUTINE CALC_IA_ENERGY_DISSIPATION_SS(M)
13
14
15
16
17 USE compar
18 USE constant
19 USE fldvar
20 USE functions
21 USE geometry
22 USE indices
23 USE kintheory
24 USE param
25 USE param1
26 USE physprop
27 USE rdf
28 USE run
29 USE toleranc
30 IMPLICIT NONE
31
32
33
34
35 INTEGER, INTENT(IN) :: M
36
37
38
39
40 INTEGER :: IJK, I, J, K
41
42
43 INTEGER :: L
44
45
46
47 INTEGER :: LM
48
49
50 DOUBLE PRECISION :: ED_common_term
51 DOUBLE PRECISION :: EDvel_sL, EDvel_sM
52 DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, &
53 D_PL, DPSUMo2
54 DOUBLE PRECISION :: Ap_lm, Dp_lm, R1p_lm, R10p_lm, R3p_lm, &
55 R4p_lm, R5p_lm, Bp_lm
56
57
58 DO IJK = ijkstart3, ijkend3
59 I = I_OF(IJK)
60 J = J_OF(IJK)
61 K = K_OF(IJK)
62
63 IF ( FLUID_AT(IJK) ) THEN
64
65 D_PM = D_P(IJK,M)
66 M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
67 NU_PM = ROP_S(IJK,M)/M_PM
68
69 DO L = 1, MMAX
70 LM = FUNLM(L,M)
71 D_PL = D_P(IJK,L)
72 M_PL = (PI/6.d0)*(D_PL**3)*RO_S(IJK,L)
73
74 MPSUM = M_PM + M_PL
75 DPSUMo2 = (D_PM+D_PL)/2.d0
76 NU_PL = ROP_S(IJK,L)/M_PL
77
78 ED_common_term = (3.d0/4.d0)*(DPSUMo2*DPSUMo2)*&
79 (1.d0+C_E)*G_0(IJK,M,L)*NU_PM*NU_PL*(M_PM*&
80 M_PL/MPSUM)*((M_PM*M_PL)**1.5)
81
82 IF (M .eq. L) THEN
83 Ap_lm = MPSUM/(2.d0)
84 Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
85 R1p_lm = 1.d0/( (Ap_lm**1.5)*(Dp_lm**3) )
86 R3p_lm = 1.d0/( (Ap_lm**1.5)*(Dp_lm**3.5) )
87
88
89
90
91
92 (IJK,LM) = ZERO
93
94
95 (IJK,M,L) = -ED_common_term* (1.d0-C_E)*&
96 (M_PL/MPSUM)*(DSQRT(PI)/6.d0)*R1p_lm*&
97 (Theta_m(IJK,M)**1.5)
98
99
100
101
102 = ED_common_term*((1.d0-C_E)*(M_PL/MPSUM)*&
103 (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm)
104 EDvel_sL_ip(IJK,M,L) = EDvel_sL*(Theta_m(IJK,L))
105
106
107
108
109
110
111
112
113
114 (IJK,M,L) = EDvel_sL_ip(IJK,M,L)
115
116 ELSE
117
118 Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/&
119 2.d0
120 Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
121 Theta_m(IJK,M) ))/(2.d0*MPSUM)
122 Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+M_PL*&
123 Theta_m(IJK,L) ))/(2.d0*MPSUM*MPSUM)
124
125 R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
126 ((9.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**4))+&
127 ((30.d0*Bp_lm**4)/(2.d0*Ap_lm**3.5 * Dp_lm**5))
128
129 R3p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3.5)))+&
130 ((21.d0*Bp_lm*Bp_lm)/(2.d0 * Ap_lm**2.5 * Dp_lm**4.5))+&
131 ((315.d0*Bp_lm**4)/(8.d0 * Ap_lm**3.5 *Dp_lm**5.5))
132
133 R4p_lm = (3.d0/( Ap_lm**2.5 * Dp_lm**3.5))+&
134 ((35.d0*Bp_lm*Bp_lm)/(2.d0 * Ap_lm**3.5 * Dp_lm**4.5))+&
135 ((441.d0*Bp_lm**4)/(8.d0 * Ap_lm**4.5 * Dp_lm**5.5))
136
137 R5p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3))+ &
138 ((5.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
139 ((14.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5))
140
141 R10p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**2.5))+&
142 ((25.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**3.5 * Dp_lm**3.5))+&
143 ((1225.d0*Bp_lm**4)/(24.d0* Ap_lm**4.5 * Dp_lm**4.5))
144
145
146
147 (IJK,LM) = ED_common_term*DSQRT(PI)*(M_PM*M_PL/&
148 (2.d0*MPSUM))*R5p_lm*( (Theta_M(IJK,M)*Theta_M(IJK,L))**3 )
149
150
151 (IJK,M,L) = -ED_common_term* (1.d0-C_E)*(M_PL/MPSUM)*&
152 (DSQRT(PI)/6.d0)*R1p_lm*( (Theta_m(IJK,M)*Theta_m(IJK,L))**3 )
153
154
155
156 = ED_common_term*( ((3.d0*DPSUMo2*PI/40.d0)*M_PL*&
157 R10p_lm)+( (DPSUMo2*PI/4.d0)*(M_PM*M_PL/MPSUM)*Bp_lm*&
158 R4p_lm)-( (1.d0-C_E)*(M_PL/MPSUM)*(DPSUMo2*PI/16.d0)*&
159 M_PL*Bp_lm*R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*&
160 (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm))
161
162 EDvel_sL_ip(IJK,M,L) = EDvel_sL*( Theta_m(IJK,M)**3.5 *&
163 Theta_m(IJK,L)**2.5 )
164
165
166
167 = ED_common_term*( (-(3.d0*DPSUMo2*PI/40.d0)*M_PM*&
168 R10p_lm)+( (DPSUMo2*PI/4.d0)*(M_PM*M_PL/MPSUM)*Bp_lm*&
169 R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*(DPSUMo2*PI/16.d0)*&
170 M_PM*Bp_lm*R4p_lm)+( (1.d0-C_E)*(M_PL/MPSUM)*&
171 (DPSUMo2*PI/48.d0)*(M_PM*M_PL/MPSUM)*R3p_lm))
172
173 EDvel_sM_ip(IJK,M,L) = EDvel_sM*( Theta_m(IJK,M)**2.5 *&
174 Theta_m(IJK,L)**3.5 )
175
176 ENDIF
177
178 ENDDO
179
180 ENDIF
181 ENDDO
182
183 RETURN
184 END SUBROUTINE CALC_IA_ENERGY_DISSIPATION_SS
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199 SUBROUTINE CALC_GD_99_ENERGY_DISSIPATION_SS(M, IER)
200
201
202
203
204 USE compar
205 USE constant
206 USE fldvar
207 USE fun_avg
208 USE functions
209 USE geometry
210 USE indices
211 USE param
212 USE param1
213 USE physprop
214 USE rdf
215 USE run
216 USE toleranc
217 use kintheory
218 IMPLICIT NONE
219
220
221
222
223 INTEGER, INTENT(INOUT) :: IER
224
225 INTEGER, INTENT(IN) :: M
226
227
228
229
230 INTEGER :: IJK, I, J, K
231
232 DOUBLE PRECISION :: press_star, c_star, zeta0_star, nu_gamma_star, &
233 lambda_num, cd_num, zeta1
234 DOUBLE PRECISION :: D_PM, EP_SM
235 DOUBLE PRECISION :: nu0, Chi
236
237
238 DO IJK = ijkstart3, ijkend3
239 I = I_OF(IJK)
240 J = J_OF(IJK)
241 K = K_OF(IJK)
242
243 IF ( FLUID_AT(IJK) ) THEN
244
245
246
247
248 = G_0(IJK,M,M)
249 EP_SM = EP_s(IJK,M)
250 D_PM = D_P(IJK,M)
251
252
253 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
254
255 press_star = 1.d0 + 2.d0*(1.d0+C_E)*EP_SM*Chi
256
257 c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
258 / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
259
260 zeta0_star = (5.d0/12.d0)*Chi*(1.d0 - C_E*C_E) &
261 * (1.d0 + (3.d0/32.d0)*c_star)
262
263 nu_gamma_star = ((1.d0+C_E)/48.d0)*Chi*(128.d0-96.d0*C_E + &
264 15.d0*C_E*C_E-15.d0*C_E*C_E*C_E+ (c_star/64.d0) * (15.d0* &
265 C_E*C_E*C_E-15.d0*C_E*C_E+498.d0*C_E-434.d0))
266
267 lambda_num = (3.d0/8.d0)*( (1.d0-C_E)*(5.d0*C_E*C_E+4.d0*C_E-1.d0)+ &
268 (c_star/12.d0)*(159.d0*C_E+3.d0*C_E*C_E-19.d0*C_E- &
269 15.d0*C_E*C_E*C_E) ) * (1.d0+C_E)
270
271
272
273 = ( (4.d0/15.d0)*lambda_num*EP_SM*Chi + &
274 (press_star-1.d0)*(2.d0/3.d0-C_E)*c_star ) / &
275 ( 0.5d0*zeta0_star+nu_gamma_star + (5.d0*c_star/64.d0) * &
276 (1.d0+(3.d0*c_star/64.d0))*Chi*(1.d0-C_E*C_E))
277
278
279
280
281 = -(1.d0-C_E)*(press_star-1.d0) + (5.d0/32.d0) * &
282 (1.d0-C_E*C_E)*(1.d0+(3.d0*c_star/64.d0))*Chi*cd_num
283
284
285
286 (IJK,M,M) = (3.d0/2.d0)*ROP_s(IJK,M)*zeta1
287
288
289
290 (IJK,M,M) = (3.d0/2.d0)*ROP_s(IJK,M)*nu0*zeta0_star
291
292 ENDIF
293 ENDDO
294
295 RETURN
296 END SUBROUTINE CALC_GD_99_ENERGY_DISSIPATION_SS
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311 SUBROUTINE CALC_GTSH_ENERGY_DISSIPATION_SS(M)
312
313
314
315
316 USE compar
317 USE constant
318 USE fldvar
319 USE fun_avg
320 USE functions
321 USE geometry
322 USE indices
323 USE param
324 USE param1
325 USE physprop
326 USE rdf
327 USE run
328 USE toleranc
329 use kintheory
330 IMPLICIT NONE
331
332
333
334
335 INTEGER, INTENT(IN) :: M
336
337
338
339
340 INTEGER :: IJK, I, J, K, IMJK, IJMK, IJKM
341
342 DOUBLE PRECISION :: D_PM, EP_SM, V_p, N_p, M_p
343 DOUBLE PRECISION :: nu0, Chi
344 DOUBLE PRECISION :: UGC, VGC, WGC, USCM, VSCM, WSCM, VREL
345 DOUBLE PRECISION :: Re_m, Re_T
346 DOUBLE PRECISION :: zeta_star, mu2_0, mu4_0, mu4_1
347 DOUBLE PRECISION :: omega, nu_j, rho_10, rho_11
348
349
350
351
352 DOUBLE PRECISION S_star
353 DOUBLE PRECISION G_gtsh
354
355
356 DO IJK = ijkstart3, ijkend3
357 I = I_OF(IJK)
358 J = J_OF(IJK)
359 K = K_OF(IJK)
360 IMJK = IM_OF(IJK)
361 IJMK = JM_OF(IJK)
362 IJKM = KM_OF(IJK)
363
364 IF ( FLUID_AT(IJK) ) THEN
365
366
367 = G_0(IJK,M,M)
368 EP_SM = EP_s(IJK,M)
369 D_PM = D_P(IJK,M)
370 V_p = pi*D_PM**3/6d0
371 n_p = EP_SM/V_p
372 M_p = V_p * ro_s(ijk,m)
373
374 nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
375
376
377
378 = AVG_X_E(U_G(IMJK),U_G(IJK),I)
379 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
380 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
381 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
382 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
383 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
384 VREL = DSQRT((UGC - USCM)**2 + (VGC - VSCM)**2 + (WGC - WSCM)**2)
385
386 Re_m = D_PM*VREL*ROP_g(ijk)/Mu_g(ijk)
387 = ro_g(ijk)*D_PM*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
388
389
390
391 (ijk) = one/6d0*D_PM*VREL**2*(3d0*pi*mu_g(ijk)*&
392 D_PM/M_p)**2 / dsqrt(pi*theta_m(ijk,m)) * &
393 S_star(EP_SM, Chi)
394
395
396 = dsqrt(2d0*pi) * Chi * (one-C_E**2)
397
398
399 = (4.5d0+C_E**2) * mu2_0
400
401
402
403
404 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + &
405 2d0*dsqrt(2d0*pi)* Chi*(one+C_E)
406
407 A2_gtsh(ijk) = zero
408 if(EP_SM> small_number) then
409
410 = 4.5d0*dsqrt(2d0*Pi)*&
411 (ro_g(ijk)/ro_s(ijk,m))**2*Re_m**2 * &
412 S_star(EP_SM,Chi) / (EP_SM*(one-EP_SM)**2 * Re_T**4)
413
414
415
416
417 (ijk) = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
418 (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
419 endif
420
421
422
423
424
425
426 (ijk,M,M) = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
427 (one+0.1875d0*A2_gtsh(ijk))*n_p*D_PM**2*dsqrt(theta_m(ijk,m))
428
429
430
431 = (one+C_E)*nu0*((one-C_E**2)*(5d0*C_E-one) - &
432 A2_gtsh(ijk)/ 6d0 * (15d0*C_E**3-3d0*C_E**2+81d0*C_E-61d0))
433 nu_j = (one+C_E)/192d0*Chi*nu0* &
434 (241d0-177d0*C_E+30d0*C_E**2-30d0*C_E**3)
435
436
437 = 2d0*Chi*EP_SM*(C_E**2-one)
438 rho_11 = 25d0/1024d0*EP_SM*Chi**2*(one-C_E**2)* &
439 (one+3d0/128d0*A2_gtsh(ijk)) * (omega/10d0 - &
440 (one+C_E)*nu0*(one/3d0-C_E)*A2_gtsh(ijk)/2d0) / &
441 (nu_j + G_gtsh(EP_SM, Chi, IJK, M)/m_p + 1.5d0* &
442 xsi_gtsh(ijk)/theta_m(ijk,m) -1.5d0*EDT_s_ip(ijk,M,M))
443
444
445 (IJK,M,M) = rho_10 + rho_11
446
447 ENDIF
448 ENDDO
449
450 RETURN
451 END SUBROUTINE CALC_GTSH_ENERGY_DISSIPATION_SS
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470 DOUBLE PRECISION FUNCTION G_gtsh (EP_SM, Chi, IJK, M)
471
472
473
474
475 USE param
476 USE param1
477 USE physprop
478 USE fldvar
479 USE constant
480 IMPLICIT NONE
481
482
483
484
485 DOUBLE PRECISION, INTENT(IN) :: EP_SM
486
487 DOUBLE PRECISION, INTENT(IN) :: Chi
488
489 INTEGER, INTENT(IN) :: IJK
490
491
492 INTEGER, INTENT(IN) :: M
493
494
495
496 DOUBLE PRECISION :: Re_T
497 DOUBLE PRECISION :: Rdiss, RdissP
498
499
500
501 DOUBLE PRECISION :: K_phi
502
503
504 if(EP_SM <= 0.1d0) then
505 RdissP = one+3d0*dsqrt(EP_SM/2d0)
506 else
507 RdissP = &
508 one + 3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*dlog(EP_SM) + &
509 11.26d0*EP_SM*(one-5.1d0*EP_SM+16.57d0*EP_SM**2-21.77d0* &
510 EP_SM**3) - EP_SM*Chi*dlog(epM)
511 endif
512
513 Re_T = ro_g(ijk)*d_p(ijk,m)*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
514
515 Rdiss = RdissP + Re_T * K_phi(EP_SM)
516
517 G_gtsh = 3d0*Pi*mu_g(ijk)*d_p(ijk,m)*Rdiss
518
519 RETURN
520 END FUNCTION G_gtsh
521
522
523
524
525
526
527
528
529 DOUBLE PRECISION FUNCTION K_phi (phi)
530 IMPLICIT NONE
531
532
533
534
535 DOUBLE PRECISION, INTENT(IN) :: phi
536
537
538 = (0.096d0 + 0.142d0*phi**0.212d0) / (1d0-phi)**4.454d0
539 K_phi = 0.0d0
540 RETURN
541 END FUNCTION K_phi
542
543
544
545
546
547
548
549
550 DOUBLE PRECISION FUNCTION R_d (phi)
551 IMPLICIT NONE
552
553
554
555
556 DOUBLE PRECISION, INTENT(IN) :: phi
557
558
559 = 1.0d0
560 if((phi > 1d-15) .and. (phi <= 0.4d0)) then
561 R_d = (1d0+3d0*dsqrt(phi/2d0)+135d0/64d0*phi*dlog(phi)+17.14d0*phi) / &
562 (1d0+0.681d0*phi-8.48*phi**2+8.16d0*phi**3)
563 elseif(phi > 0.4d0) then
564 R_d = 10d0*phi/(1d0-phi)**3 + 0.7d0
565 endif
566 RETURN
567 END FUNCTION R_d
568
569
570
571
572
573
574
575
576 DOUBLE PRECISION FUNCTION S_star (phi, Chi)
577 IMPLICIT NONE
578
579
580
581
582 DOUBLE PRECISION, INTENT(IN) :: phi
583
584 DOUBLE PRECISION, INTENT(IN) :: Chi
585
586
587
588 DOUBLE PRECISION R_d
589
590
591 = 1.0d0
592 if(phi >= 0.1d0) &
593 S_star = R_d(phi)**2/(Chi*(1d0+3.5d0*dsqrt(phi)+5.9*phi))
594 RETURN
595 END FUNCTION S_Star
596