File: /nfs/home/0/users/jenkins/mfix.git/model/calc_mu_s.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56 SUBROUTINE CALC_MU_s(M, IER)
57
58
59
60
61
62
63 USE run, only: granular_energy
64
65 USE run, only: kt_type_enum
66 USE run, only: lun_1984
67 USE run, only: simonin_1996
68 USE run, only: ahmadi_1995
69 USE run, only: gd_1999
70 USE run, only: gtsh_2012
71 USE run, only: ia_2005
72 USE run, only: ghd_2007
73 USE run, only: kt_type
74
75 USE run, only: subgrid_type_enum, milioli, igci, undefined_subgrid_type
76
77 USE run, only: friction, schaeffer
78
79 USE run, only: blending_stress
80
81 USE run, only: shear
82
83
84 USE visc_s, only: mu_s, mu_s_c, mu_s_v, mu_s_p, mu_s_f
85 USE visc_s, only: lambda_s, lambda_s_c, lambda_s_v, lambda_s_p, lambda_s_f
86
87 USE fldvar, only: P_s, p_s_c, p_s_v, p_s_p, p_s_f
88
89
90 USE physprop, only: mu_s0, k_s0
91
92 USE physprop, only: Kth_s
93
94
95 USE physprop, only: close_packed
96
97 USE physprop, only: mmax
98
99 USE qmom_kinetic_equation, only: qmomk
100
101 USE mms, only: use_mms
102
103
104 USE compar
105 USE geometry
106 USE indices
107 USE functions
108
109 IMPLICIT NONE
110
111
112
113
114 INTEGER, INTENT(IN) :: M
115
116 INTEGER, INTENT(INOUT) :: IER
117
118
119
120
121 INTEGER :: IJK
122
123 DOUBLE PRECISION :: BLEND
124
125
126
127 DOUBLE PRECISION, EXTERNAL :: BLEND_FUNCTION
128
129
130
131 IF (KT_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) RETURN
132
133 IF (SHEAR) CALL add_shear(M)
134
135
136 CALL INIT_MU_S(M, IER)
137
138 IF (SHEAR) call remove_shear(M)
139
140
141
142
143 IF( MU_s0 /= UNDEFINED) THEN
144
145 IF(USE_MMS) THEN
146 DO IJK = ijkstart3, ijkend3
147 MU_S(IJK,M) = MU_S0
148 LAMBDA_S(IJK,M) = -2.D0/3.D0*MU_S(IJK,M)
149 KTH_S(IJK,M) = K_S0(M)
150 ENDDO
151 ENDIF
152 RETURN
153 ENDIF
154
155
156
157
158 IF (.NOT. QMOMK) THEN
159
160
161 IF(.NOT.GRANULAR_ENERGY) THEN
162 IF(SUBGRID_TYPE_ENUM .ne. UNDEFINED_SUBGRID_TYPE) THEN
163 IF (SUBGRID_TYPE_ENUM .EQ. IGCI) THEN
164 CALL SUBGRID_STRESS_IGCI(M)
165 ELSEIF (SUBGRID_TYPE_ENUM .EQ. MILIOLI) THEN
166 CALL SUBGRID_STRESS_MILIOLI(M)
167 ENDIF
168 ELSE
169 call gt_algebraic(M)
170 ENDIF
171 ELSE
172 SELECT CASE(KT_TYPE_ENUM)
173 CASE (IA_2005)
174 CALL gt_pde_ia(M)
175 CASE (GD_1999)
176 CALL gt_pde_gd(M)
177 CASE (GTSH_2012)
178 CALL gt_pde_gtsh(M)
179 CASE (GHD_2007)
180 CALL TRANSPORT_COEFF_GHD(M)
181 CASE(LUN_1984)
182 CALL gt_pde_lun(M)
183 CASE(SIMONIN_1996)
184 CALL gt_pde_simonin(M)
185 CASE (AHMADI_1995)
186 CALL gt_pde_ahmadi(M)
187 CASE DEFAULT
188
189 WRITE (*, '(A)') 'CALC_MU_S'
190 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
191 call mfix_exit(myPE)
192 END SELECT
193 ENDIF
194 ENDIF
195
196
197
198
199
200
201 IF (SCHAEFFER .AND. CLOSE_PACKED(M)) call friction_schaeffer(M)
202
203 IF (FRICTION .AND. CLOSE_PACKED(M)) call friction_princeton(M)
204
205
206
207
208 IF(BLENDING_STRESS) THEN
209 DO IJK = ijkstart3, ijkend3
210 blend = blend_function(IJK)
211 Mu_s_c(IJK,M) = Mu_s_v(IJK)
212 Mu_s(IJK,M) = (ONE-blend)*Mu_s_p(IJK) &
213 + blend*Mu_s_v(IJK) + Mu_s_f(IJK)
214
215
216 (IJK,M)= Lambda_s_v(IJK)
217 LAMBDA_s(IJK,M) = (ONE-blend)*Lambda_s_p(IJK) + &
218 blend*Lambda_s_v(IJK) + Lambda_s_f(IJK)
219
220
221
222 (IJK,M) = P_s_v(IJK)
223 P_s(IJK,M) = (ONE-blend)*P_s_p(IJK) + &
224 blend*P_s_v(IJK) + P_s_f(IJK)
225
226 ENDDO
227 ELSE
228
229
230 (:,M) = Mu_s_v(:)
231 Mu_s(:,M) = Mu_s_p(:) + Mu_s_v(:) + Mu_s_f(:)
232
233
234 (:,M)= Lambda_s_v(:)
235 LAMBDA_s(:,M) = Lambda_s_p(:) + Lambda_s_v(:) + Lambda_s_f(:)
236
237
238
239
240
241
242 (:,M) = P_s_v(:)
243 P_s(:,M) = P_s_v(:) + P_s_f(:)
244 ENDIF
245
246 RETURN
247 END SUBROUTINE CALC_MU_s
248
249
250
251
252
253
254
255
256
257
258
259
260
261 Subroutine friction_schaeffer(M)
262
263
264
265
266 USE param
267 USE param1
268 USE geometry
269 USE compar
270 USE fldvar
271 USE indices
272 USE visc_s
273 USE physprop
274 USE run
275 USE constant
276 USE functions
277 IMPLICIT NONE
278
279
280
281
282 INTEGER, INTENT(IN) :: M
283
284
285
286
287 DOUBLE PRECISION :: MAX_MU_s
288 PARAMETER (MAX_MU_s = 1000.D0)
289
290
291
292
293 INTEGER :: IJK
294
295 INTEGER :: MM
296
297 DOUBLE PRECISION :: SUM_EPS_CP
298
299 DOUBLE PRECISION :: qxP_s
300
301
302 DO 200 IJK = ijkstart3, ijkend3
303 IF ( FLUID_AT(IJK) ) THEN
304
305 IF(EP_g(IJK) .LT. EP_g_blend_end(IJK)) THEN
306
307
308
309
310
311 =0.0
312 DO MM=1,SMAX
313 IF (CLOSE_PACKED(MM)) SUM_EPS_CP=SUM_EPS_CP+EP_S(IJK,MM)
314 END DO
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338 = SQRT( (4.D0 * Sin2_Phi) * I2_devD_s(IJK))
339 MU_s_p(IJK) = P_star(IJK) * Sin2_Phi&
340 / (qxP_s + SMALL_NUMBER) &
341 *(EP_S(IJK,M)/SUM_EPS_CP)
342
343 (IJK) = MIN(MU_s_p(IJK), to_SI*MAX_MU_s)
344
345 LAMBDA_s_p(IJK) = ZERO
346
347
348
349 IF(.NOT.GRANULAR_ENERGY) THETA_m(IJK, M) = ZERO
350 ENDIF
351 ENDIF
352
353
354 CONTINUE
355
356 RETURN
357 END SUBROUTINE FRICTION_SCHAEFFER
358
359
360
361
362
363
364
365 Subroutine Gt_algebraic (M)
366
367
368
369
370 USE compar
371 USE constant
372 USE fldvar
373 USE functions
374 USE geometry
375 USE indices
376 USE param
377 USE param1
378 USE physprop
379 USE rdf
380 USE run
381 USE trace
382 USE visc_s
383 IMPLICIT NONE
384
385
386
387
388 INTEGER, INTENT(IN) :: M
389
390
391
392
393 INTEGER :: IJK
394
395 DOUBLE PRECISION :: aq, bq, cq
396
397 DOUBLE PRECISION :: K_1m
398
399 DOUBLE PRECISION :: K_2m
400
401 DOUBLE PRECISION :: K_3m
402
403 DOUBLE PRECISION :: K_4m
404
405 DOUBLE PRECISION :: K_5m
406
407 DOUBLE PRECISION :: EP_sxSQRTHETA
408
409 DOUBLE PRECISION :: EP_s2xTHETA, temp_local
410
411
412
413
414
415
416 DO IJK = ijkstart3, ijkend3
417
418 IF ( FLUID_AT(IJK) ) THEN
419
420 IF(EP_g(IJK) .GE. EP_g_blend_start(IJK)) THEN
421
422
423 = 2.D0 * (ONE + C_e) * RO_S(IJK,M) * G_0(IJK, M, M)
424 K_3m = HALF * D_p(IJK,M) * RO_S(IJK,M) * (&
425 ( (SQRT_PI / (3.D0*(3.D0 - C_e))) *&
426 (HALF*(3d0*C_e+ONE) + 0.4D0*(ONE + C_e)*(3.D0*C_e - ONE)*&
427 EP_s(IJK,M)*G_0(IJK, M,M)) ) + 8.D0*EP_s(IJK,M)&
428 *G_0(IJK, M,M)*(ONE + C_e)/ (5.D0*SQRT_PI) )
429 K_2m = 4.D0 * D_p(IJK,M) * RO_S(IJK,M) * (ONE + C_e) *&
430 EP_s(IJK,M) * G_0(IJK, M,M) / (3.D0 * SQRT_PI) - 2.D0/3.D0 * K_3m
431 K_4m = 12.D0 * (ONE - C_e*C_e) *&
432 RO_S(IJK,M) * G_0(IJK, M,M) / (D_p(IJK,M) * SQRT_PI)
433 aq = K_4m*EP_s(IJK,M)
434 bq = K_1m*EP_s(IJK,M)*trD_s_C(IJK,M)
435 cq = -(K_2m*trD_s_C(IJK,M)*trD_s_C(IJK,M)&
436 + 2.D0*K_3m*trD_s2(IJK,M))
437
438
439 IF(V_ex .NE. ZERO) THEN
440 K_5m = 0.4 * (ONE + C_e) * G_0(IJK, M,M) * RO_S(IJK,M) *&
441 ( (V_ex * D_p(IJK,M)) / (ONE - EP_s(IJK,M) * V_ex) )**2
442 bq = bq + EP_s(IJK,M) * K_5m * (trM_s(IJK) + 2.D0 * trDM_s(IJK))
443 ELSE
444 K_5m = ZERO
445 ENDIF
446
447
448 = bq**2 - 4.D0 * aq * cq
449 EP_sxSQRTHETA = (-bq + SQRT(temp_local))&
450 / ( 2.D0 * K_4m )
451 EP_s2xTHETA = EP_sxSQRTHETA * EP_sxSQRTHETA
452
453 IF(EP_s(IJK,M) > SMALL_NUMBER)THEN
454
455 (IJK,M) = EP_s2xTHETA/(EP_s(IJK,M)*EP_s(IJK,M))
456 ELSE
457 THETA_m(IJK,M) = ZERO
458 ENDIF
459
460
461 (IJK) = K_1m * EP_s2xTHETA
462
463
464 (IJK) = K_2m * EP_sxSQRTHETA
465
466
467 (IJK) = K_3m * EP_sxSQRTHETA
468
469
470 (IJK, M) = -K_5m * EP_s2xTHETA
471 ENDIF
472
473 ENDIF
474 ENDDO
475
476
477 RETURN
478 END SUBROUTINE GT_ALGEBRAIC
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500 Subroutine gt_pde_lun (M)
501
502
503
504
505 USE param
506 USE param1, only: zero, one, small_number
507
508
509 USE fldvar, only: ro_g
510
511 USE fldvar, only: d_p, rop_s, ro_s
512
513 USE fldvar, only: theta_m
514
515 USE fldvar, only: p_s_v, ep_s
516
517 USE fldvar, only: ep_s
518
519 USE visc_s, only: mu_s_v, mu_b_v, lambda_s_v
520
521 USE physprop, only: smax
522
523 USE physprop, only: Kth_s, kphi_s
524
525 USE constant, only: switch
526 USE constant, only: alpha
527 USE constant, only: pi
528 USe constant, only: C_e, eta
529
530 USE toleranc, only: dil_ep_s
531
532 USE drag, only: f_gs, dga_s
533
534 USE compar
535 USE functions
536
537 USE geometry
538 USE indices
539 USE rdf
540 IMPLICIT NONE
541
542
543
544
545 INTEGER, INTENT(IN) :: M
546
547
548
549
550 INTEGER :: IJK, I, J, K
551
552 INTEGER :: MM
553
554 DOUBLE PRECISION :: Mu_star, Kth, Kth_star
555
556 DOUBLE PRECISION :: SUM_EpsGo
557
558
559 DO 200 IJK = ijkstart3, ijkend3
560 I = I_OF(IJK)
561 J = J_OF(IJK)
562 K = K_OF(IJK)
563
564 IF ( FLUID_AT(IJK) ) THEN
565
566
567
568
569
570
571
572 = ZERO
573 DO MM = 1, SMAX
574 SUM_EpsGo = SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
575 ENDDO
576
577
578 (IJK) = ROP_s(IJK,M)*(1.d0 + 4.d0 * Eta *&
579 SUM_EpsGo)*Theta_m(IJK,M)
580
581
582 (IJK) = (5.d0*DSQRT(Pi*Theta_m(IJK,M))*D_p(IJK,M)*&
583 RO_S(IJK,M))/96.d0
584 Mu_b_v(IJK) = (256.d0*Mu_s_v(IJK)*EP_s(IJK,M)*SUM_EpsGo)&
585 /(5.d0*Pi)
586
587
588 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
589 = Mu_s_v(IJK)
590 ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
591 Mu_star = ZERO
592 ELSEIF(EP_S(IJK,M) < DIL_EP_S) THEN
593 Mu_star = RO_S(IJK,M)*EP_s(IJK,M)* G_0(IJK,M,M)*&
594 Theta_m(IJK,M)* Mu_s_v(IJK)/ &
595 (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M) &
596 + 2.d0*DgA_S(IJK,M)/RO_S(IJK,M)* Mu_s_v(IJK))
597 ELSE
598 Mu_star = RO_S(IJK,M)*EP_s(IJK,M)* G_0(IJK,M,M)*&
599 Theta_m(IJK,M)*Mu_s_v(IJK)/ &
600 (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M)+ &
601 (2.d0*F_gs(IJK,M)*Mu_s_v(IJK)/&
602 (RO_S(IJK,M)*EP_s(IJK,M) )) )
603 ENDIF
604
605
606 (IJK) =&
607 ((2.d0+ALPHA)/3.d0)*((Mu_star/(Eta*(2.d0-Eta)*&
608 G_0(IJK,M,M)))*(ONE+1.6d0*Eta*SUM_EpsGo)&
609 *(ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
610 SUM_EpsGo)+(0.6d0*Mu_b_v(IJK)*Eta))
611
612
613 (IJK) = Eta*Mu_b_v(IJK) - (2d0*Mu_s_v(IJK))/3d0
614
615 Kth=75d0*RO_S(IJK,M)*D_p(IJK,M)*DSQRT(Pi*Theta_m(IJK,M))/&
616 (48d0*Eta*(41d0-33d0*Eta))
617
618 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
619 =Kth
620 ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
621 Kth_star = ZERO
622 ELSEIF(EP_S(IJK,M) < DIL_EP_S) THEN
623 Kth_star = RO_S(IJK,M)*EP_s(IJK,M)* G_0(IJK,M,M)*Theta_m(IJK,M)* Kth/ &
624 (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M) &
625 + 1.2d0*DgA_s(IJK,M)/RO_S(IJK,M)* Kth)
626 ELSE
627 Kth_star = RO_S(IJK,M)*EP_s(IJK,M)* G_0(IJK,M,M)*Theta_m(IJK,M)*Kth/ &
628 (RO_S(IJK,M)*SUM_EpsGo*Theta_m(IJK,M)+ &
629 (1.2d0*F_gs(IJK,M)*Kth/(RO_S(IJK,M)*EP_s(IJK,M))) )
630 ENDIF
631
632
633 (IJK,M) = Kth_star/G_0(IJK,M,M)*(&
634 ( ONE + (12d0/5.d0)*Eta*SUM_EpsGo )&
635 * ( ONE + (12d0/5.d0)*Eta*Eta*(4d0*Eta-3d0)* SUM_EpsGo )&
636 + (64d0/(25d0*Pi)) * (41d0-33d0*Eta) * (Eta*SUM_EpsGo)**2 )
637
638
639
640
641
642
643
644
645 (IJK,M) = ZERO
646
647
648
649
650
651
652
653 ENDIF
654 CONTINUE
655
656
657 RETURN
658 END SUBROUTINE GT_PDE_LUN
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676 Subroutine gt_pde_ahmadi (M)
677
678
679
680
681 USE param
682 USE param1
683 USE geometry
684 USE compar
685 USE fldvar
686 USE indices
687 USE visc_s
688 USE physprop
689 USE run
690 USE constant
691 USE toleranc
692 USE turb
693 USE drag
694 USE functions
695
696 USE rdf
697 IMPLICIT NONE
698
699
700
701
702 INTEGER, INTENT(IN) :: M
703
704
705
706
707 INTEGER :: IJK, I, J, K
708
709 INTEGER :: MM
710
711 DOUBLE PRECISION :: Tau_12_st
712 DOUBLE PRECISION :: Tmp_Ahmadi_Const
713
714 DOUBLE PRECISION :: SUM_EpsGo
715
716
717 DO 200 IJK = ijkstart3, ijkend3
718 I = I_OF(IJK)
719 J = J_OF(IJK)
720 K = K_OF(IJK)
721
722 IF ( FLUID_AT(IJK) ) THEN
723
724
725
726
727 IF(Ep_s(IJK,M) > DIL_EP_S .AND. F_GS(IJK,1) > small_number) THEN
728 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1)
729 ELSE
730 = RO_S(IJK,M)/DgA_S(IJK,1)
731 ENDIF
732
733
734
735
736
737
738
739
740 = ZERO
741 DO MM = 1, SMAX
742 SUM_EpsGo = SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
743 ENDDO
744
745 P_s_v(IJK) = ROP_s(IJK,M)*Theta_m(IJK,M) * &
746 ( (ONE + 4.d0*SUM_EpsGo ) + HALF*(ONE - C_e*C_e) )
747
748 IF(EP_s(IJK,M) < (ONE-EP_star_array(ijk))) THEN
749 Tmp_Ahmadi_Const = &
750 ONE/(ONE+ Tau_1(ijk)/Tau_12_st * &
751 (ONE-EP_s(IJK,M)/(ONE-EP_star_array(ijk)))**3)
752 ELSE
753 Tmp_Ahmadi_Const = ONE
754 ENDIF
755
756
757
758
759 (IJK) = Tmp_Ahmadi_Const &
760 *0.1045d0*(ONE/G_0(IJK,M,M)+3.2d0*EP_s(IJK,M)+12.1824d0* &
761 G_0(IJK,M,M)*EP_s(IJK,M)*EP_s(IJK,M))*D_p(IJK,M)*RO_S(IJK,M)* &
762 DSQRT(Theta_m(IJK,M))
763
764
765
766
767
768
769 (IJK) = 5.d0/3.d0* Tmp_Ahmadi_Const &
770 *0.1045d0*(12.1824d0*G_0(IJK,M,M)*EP_s(IJK,M)*EP_s(IJK,M)) &
771 *D_p(IJK,M)*RO_S(IJK,M)* DSQRT(Theta_m(IJK,M))
772
773
774 (IJK) = Eta*Mu_b_v(IJK) - (2d0*Mu_s_v(IJK))/3d0
775
776
777
778 (IJK,M) = 0.1306D0*RO_S(IJK,M)*D_p(IJK,M)*(ONE+C_e**2)* &
779 (ONE/G_0(IJK,M,M)+4.8D0*EP_s(IJK,M)+12.1184D0 &
780 *EP_s(IJK,M)*EP_s(IJK,M)*G_0(IJK,M,M) ) &
781 *DSQRT(Theta_m(IJK,M))
782
783
784 ENDIF
785 CONTINUE
786
787
788 RETURN
789 END SUBROUTINE GT_PDE_AHMADI
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811 Subroutine gt_pde_simonin (M)
812
813
814
815
816 USE compar
817 USE constant
818 USE drag
819 USE fldvar
820 USE functions
821 USE geometry
822 USE indices
823 USE param
824 USE param1
825 USE physprop
826 USE rdf
827 USE run
828 USE toleranc
829 USE turb
830 USE visc_s
831 IMPLICIT NONE
832
833
834
835
836 INTEGER, INTENT(IN) :: M
837
838
839
840
841 INTEGER :: IJK, I, J, K
842
843 INTEGER :: MM
844
845 DOUBLE PRECISION :: Tau_12_st, Tau_2_c, Tau_2, Zeta_r, C_Beta
846 DOUBLE PRECISION :: Sigma_c, Zeta_c, Omega_c, Zeta_c_2, X_21, Nu_t
847 DOUBLE PRECISION :: MU_2_T_Kin, Mu_2_Col, Kappa_kin, Kappa_Col
848
849 DOUBLE PRECISION :: SUM_EpsGo
850
851
852 DO 200 IJK = ijkstart3, ijkend3
853 I = I_OF(IJK)
854 J = J_OF(IJK)
855 K = K_OF(IJK)
856
857 IF ( FLUID_AT(IJK) ) THEN
858
859
860
861
862
863
864
865
866 IF(Ep_s(IJK,M) > DIL_EP_S .AND. F_GS(IJK,1) > small_number) THEN
867 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1)
868 ELSE
869 Tau_12_st = RO_S(IJK,M)/DgA_s(IJK,1)
870 ENDIF
871
872
873 = 3.0d0 * VREL_array(IJK,1)**2 / &
874 (2.0d0*K_Turb_G(IJK)+small_number)
875
876
877
878 = 1.8d0 - 1.35d0*Cos_Theta(IJK)**2
879
880
881 (ijk) = Tau_1(ijk)/sqrt(ONE+C_Beta*Zeta_r)
882
883
884 IF(Ep_s(IJK,M) > DIL_EP_S) THEN
885 Tau_2_c = D_p(IJK,M)/(6.d0*Ep_s(IJK,M)*G_0(IJK,M,M) &
886 *DSQRT(16.d0*(Theta_m(ijk,m)+Small_number)/PI))
887 ELSE
888 = LARGE_NUMBER
889 ENDIF
890
891
892 = (ONE+ C_e)*(3.d0-C_e)/5.d0
893
894 = (ONE+ C_e)*(49.d0-33.d0*C_e)/100.d0
895 Omega_c = 3.d0*(ONE+ C_e)**2 *(2.d0*C_e-ONE)/5.d0
896 Zeta_c_2= 2./5.*(ONE+ C_e)*(3.d0*C_e-ONE)
897
898
899
900 = ONE/(2./Tau_12_st+Sigma_c/Tau_2_c)
901
902 = Tau_12(ijk)/Tau_12_st
903
904
905
906 = Ep_s(IJK,M)*RO_S(IJK,M)/(EP_g(IJK)*RO_g(IJK))
907
908
909
910
911
912 (ijk) = Nu_t / (ONE+Nu_t*(ONE+X_21)) * &
913 (2.d+0 *K_Turb_G(IJK) + 3.d+0 *X_21*theta_m(ijk,m))
914
915
916 IF(K_12(ijk) > DSQRT(6.0D0*K_Turb_G(IJK)*theta_m(ijk,m))) THEN
917 K_12(ijk) = DSQRT(6.0D0*K_Turb_G(IJK)*theta_m(ijk,m))
918 ENDIF
919
920
921
922
923
924
925
926 = ZERO
927 DO MM = 1, SMAX
928 SUM_EpsGo = SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
929 ENDDO
930
931
932
933
934 (IJK) = ROP_s(IJK,M) * &
935 (ONE + 4.d0*Eta*SUM_EpsGo)*Theta_m(IJK,M)
936
937
938
939 = (2.0d0/3.0d0*K_12(ijk)*Nu_t + Theta_m(IJK,M) * &
940 (ONE+ zeta_c_2*EP_s(IJK,M)*G_0(IJK,M,M)))*Tau_2
941 Mu_2_Col = 8.d0/5.d0*EP_s(IJK,M)*G_0(IJK,M,M)*Eta* (MU_2_T_Kin+ &
942 D_p(IJK,M)*DSQRT(Theta_m(IJK,M)/PI))
943 Mu_b_v(IJK) = 5.d0/3.d0*EP_s(IJK,M)*RO_S(IJK,M)*Mu_2_Col
944 Mu_s_v(IJK) = EP_s(IJK,M)*RO_S(IJK,M)*(MU_2_T_Kin + Mu_2_Col)
945
946
947
948 (IJK) = Eta*Mu_b_v(IJK) - (2d0*Mu_s_v(IJK))/3d0
949
950
951
952 = (9.d0/10.d0*K_12(ijk)*Nu_t + 3.0D0/2.0D0 * &
953 Theta_m(IJK,M)*(ONE+ Omega_c*EP_s(IJK,M)*G_0(IJK,M,M)))/&
954 (9.d0/(5.d0*Tau_12_st) + zeta_c/Tau_2_c)
955 Kappa_Col = 18.d0/5.d0*EP_s(IJK,M)*G_0(IJK,M,M)*Eta* &
956 (Kappa_kin+ 5.d0/9.d0*D_p(IJK,M)*DSQRT(Theta_m(IJK,M)/PI))
957 Kth_s(IJK,M) = EP_s(IJK,M)*RO_S(IJK,M)*(Kappa_kin + Kappa_Col)
958
959
960 ENDIF
961 CONTINUE
962
963
964 RETURN
965 END SUBROUTINE GT_PDE_SIMONIN
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983 Subroutine gt_pde_gd (M)
984
985
986
987
988 USE compar
989 USE constant
990 USE drag
991 USE fldvar
992 USE functions
993 USE geometry
994 USE indices
995 USE param
996 USE param1
997 USE physprop
998 USE rdf
999 USE run
1000 USE toleranc
1001 USE visc_s
1002 use kintheory
1003 IMPLICIT NONE
1004
1005
1006
1007
1008 INTEGER, INTENT(IN) :: M
1009
1010
1011
1012
1013 INTEGER :: IJK, I, J, K
1014
1015 DOUBLE PRECISION :: Mu_star, Kth_star
1016
1017 DOUBLE PRECISION :: D_PM, M_PM, NU_PM, EP_SM, RO_SM, ROP_SM
1018
1019 DOUBLE PRECISION :: chi, dChiOdphi
1020
1021 DOUBLE PRECISION :: c_star, zeta0_star, nu_eta_star, &
1022 gamma_star, eta_k_star, eta_star, eta0, &
1023 kappa0, nu_kappa_star, kappa_k_star, &
1024 qmu_k_star, qmu_star, kappa_star, press_star
1025
1026
1027 DO 200 IJK = ijkstart3, ijkend3
1028 I = I_OF(IJK)
1029 J = J_OF(IJK)
1030 K = K_OF(IJK)
1031
1032 IF ( FLUID_AT(IJK) ) THEN
1033
1034
1035 = D_P(IJK,M)
1036 RO_SM = RO_S(IJK,M)
1037 ROP_SM = ROP_S(IJK,M)
1038 EP_SM = EP_S(IJK,M)
1039 M_PM = (PI/6.d0)*D_PM**3 * RO_SM
1040 NU_PM = ROP_SM/M_PM
1041 Chi = G_0(IJK,M,M)
1042 dChiOdphi = DG_0DNU(EP_SM)
1043
1044
1045
1046
1047
1048
1049 = 1.d0 + 2.d0*(1.d0+C_E)*EP_SM*Chi
1050
1051
1052 (IJK) = ROP_sM*Theta_m(IJK,M)*press_star
1053
1054
1055 = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
1056 / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
1057
1058 zeta0_star = (5.d0/12.d0)*Chi*(1.d0 - C_E*C_E) &
1059 * (1.d0 + (3.d0/32.d0)*c_star)
1060
1061 nu_eta_star = Chi*(1.d0 - 0.25d0*(1.d0-C_E)*(1.d0-C_E)) &
1062 * (1.d0-(c_star/64.d0))
1063
1064 gamma_star = (4.d0/5.d0)*(32.d0/PI)*EP_SM*EP_SM &
1065 * Chi*(1.d0+C_E) * (1.d0 - (c_star/32.d0))
1066
1067 eta_k_star = (1.d0 - (2.d0/5.d0)*(1.d0+C_E)*(1.d0-3.d0*C_E) &
1068 * EP_SM*Chi ) / (nu_eta_star - 0.5d0*zeta0_star)
1069
1070 eta_star = eta_k_star*(1.d0 + (4.d0/5.d0)*EP_SM*Chi &
1071 * (1.d0+C_E) ) + (3.d0/5.d0)*gamma_star
1072
1073 eta0 = 5.0d0*M_PM*DSQRT(Theta_m(IJK,M)/PI) / (16.d0*D_PM*D_PM)
1074
1075 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1076 Mu_star = eta0
1077 ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
1078 Mu_star = ZERO
1079 ELSEIF(EP_SM < DIL_EP_S) THEN
1080 Mu_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*eta0 / &
1081 ( RO_S(IJK,M)*EP_SM*Chi*Theta_m(IJK,M) + &
1082 2.d0*DgA_s(IJK,M)*eta0/RO_S(IJK,M) )
1083 ELSE
1084 Mu_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*eta0 / &
1085 ( RO_SM*EP_SM*Chi*Theta_m(IJK,M) + &
1086 (2.d0*F_gs(IJK,M)*eta0/(RO_SM*EP_SM)) )
1087 ENDIF
1088
1089
1090 (IJK) = Mu_star * eta_star
1091 Mu_b_v(IJK) = Mu_star * gamma_star
1092
1093
1094 (IJK) = Mu_b_v(IJK) - (2.d0/3.d0)*Mu_s_v(IJK)
1095
1096
1097
1098
1099 = (15.d0/4.d0)*eta0
1100
1101 nu_kappa_star = (Chi/3.d0)*(1.d0+C_E) * ( 1.d0 + &
1102 (33.d0/16.d0)*(1.d0-C_E) + ((19.d0-3.d0*C_E)/1024.d0)*&
1103 c_star)
1104
1105
1106 = (2.d0/3.d0)*(1.d0 +0.5d0*(1.d0+press_star)*&
1107 c_star + (3.d0/5.d0)*EP_SM*Chi*(1.d0+C_E)*(1.d0+C_E) * &
1108 (2.d0*C_E - 1.d0 + ( 0.5d0*(1.d0+C_E) - 5.d0/&
1109 (3*(1.d0+C_E))) * c_star ) ) / (nu_kappa_star - &
1110 2.d0*zeta0_star)
1111
1112 kappa_star = kappa_k_star * (1.d0 + (6.d0/5.d0)*EP_SM* &
1113 Chi*(1.d0+C_E) ) + (256.d0/25.d0)*(EP_SM* &
1114 EP_SM/PI)*Chi*(1.d0+C_E)*(1.d0+(7.d0/32.d0)* &
1115 c_star)
1116
1117 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1118 Kth_star= kappa0
1119 ELSEIF(Theta_m(IJK,M) .LT. SMALL_NUMBER)THEN
1120 Kth_star = ZERO
1121 ELSEIF(EP_SM < DIL_EP_S) THEN
1122 Kth_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*kappa0/ &
1123 (RO_SM*EP_SM*Chi*Theta_m(IJK,M) + 1.2d0*DgA_s(IJK,M)* &
1124 kappa0/RO_SM)
1125 ELSE
1126 Kth_star = RO_SM*EP_SM*Chi*Theta_m(IJK,M)*kappa0/ &
1127 (RO_SM*EP_SM*Chi*Theta_m(IJK,M)+ &
1128 (1.2d0*F_gs(IJK,M)*kappa0/(RO_SM*EP_SM)) )
1129 ENDIF
1130
1131
1132 (IJK,M) = Kth_star * kappa_star
1133
1134
1135
1136 = 2.d0*( (1.d0+EP_SM*dChiOdphi)* &
1137 zeta0_star*kappa_k_star + ( (press_star/3.d0) + &
1138 (2.d0/3.d0)* EP_SM*(1.d0+C_E)*(Chi+EP_SM* dChiOdphi) )*&
1139 c_star - (4.d0/5.d0)*EP_SM*Chi* (1.d0+(EP_SM/2.d0)*&
1140 dChiOdphi)* (1.d0+C_E) * ( C_E*(1.d0-C_E)+0.25d0*&
1141 ((4.d0/3.d0)+C_E* (1.d0-C_E))*c_star ) ) / &
1142 (2.d0*nu_kappa_star-3.d0*zeta0_star)
1143
1144 qmu_star = qmu_k_star*(1.d0+(6.d0/5.d0)*EP_SM*Chi*&
1145 (1.d0+C_E) )
1146
1147 IF (EP_SM .LT. SMALL_NUMBER) THEN
1148 Kphi_s(IJK,M) = ZERO
1149 ELSE
1150 Kphi_s(IJK,M) = (Theta_m(IJK,M)*Kth_star/NU_PM)*qmu_star
1151 ENDIF
1152
1153 ENDIF
1154 CONTINUE
1155
1156 RETURN
1157 END SUBROUTINE GT_PDE_GD
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178 Subroutine gt_pde_gtsh (M)
1179
1180
1181
1182
1183 USE compar
1184 USE constant
1185 USE drag
1186 USE fldvar
1187 USE functions
1188 USE geometry
1189 USE indices
1190 USE param
1191 USE param1
1192 USE physprop
1193 USE rdf
1194 USE run
1195 USE toleranc
1196 USE visc_s
1197 use kintheory
1198 IMPLICIT NONE
1199
1200
1201
1202
1203 INTEGER, INTENT(IN) :: M
1204
1205
1206
1207
1208 INTEGER :: IJK, I, J, K
1209
1210 DOUBLE PRECISION :: D_PM, M_PM, NU_PM, EP_SM
1211
1212 DOUBLE PRECISION :: chi, dChiOdphi, eta0
1213
1214 DOUBLE PRECISION :: nu0, nuN, etaK
1215 DOUBLE PRECISION :: dZeta_dT, dGama_dT, NuK, Kth0, KthK
1216 DOUBLE PRECISION :: Rdissdphi, Kphidphi, Re_T, dGamadn, dRdphi
1217 DOUBLE PRECISION :: denom
1218 DOUBLE PRECISION :: dSdphi, R_dphi, Tau_st, dPsidn, MuK
1219
1220 DOUBLE PRECISION :: VREL
1221
1222
1223
1224
1225 DOUBLE PRECISION, EXTERNAL :: G_gtsh
1226 DOUBLE PRECISION, EXTERNAL :: S_star
1227 DOUBLE PRECISION, EXTERNAL :: K_phi
1228 DOUBLE PRECISION, EXTERNAL :: R_d
1229
1230
1231 DO 200 IJK = ijkstart3, ijkend3
1232 I = I_OF(IJK)
1233 J = J_OF(IJK)
1234 K = K_OF(IJK)
1235
1236 IF ( FLUID_AT(IJK) ) THEN
1237
1238
1239 = D_P(IJK,M)
1240 EP_SM = EP_S(IJK,M)
1241 M_PM = (PI/6.d0)*D_PM**3 * RO_S(IJK,M)
1242 NU_PM = ROP_S(IJK,M)/M_PM
1243 Chi = G_0(IJK,M,M)
1244 dChiOdphi = DG_0DNU(EP_SM)
1245 VREL = VREL_array(IJK,M)
1246
1247
1248
1249
1250
1251
1252 (IJK) = ROP_s(IJK,M)*Theta_m(IJK,M)* &
1253 (one+2d0*(one+C_E)*Chi*EP_SM)
1254
1255
1256
1257 = 0.3125d0/(dsqrt(pi)*D_PM**2)*M_pm*&
1258 dsqrt(theta_m(ijk,m))
1259
1260 nu0 = (96.d0/5.d0)*(EP_SM/D_PM)*DSQRT(Theta_m(IJK,M)/PI)
1261
1262
1263 = 0.25d0*nu0*Chi*(3d0-C_E)*(one+C_E) * &
1264 (one+0.4375d0*A2_gtsh(ijk))
1265
1266 = rop_s(ijk,m)*theta_m(ijk,m) / (nuN-0.5d0*( &
1267 EDT_s_ip(ijk,M,M)-xsi_gtsh(ijk)/theta_m(ijk,m) - &
1268 2d0*G_gtsh(EP_SM, Chi, IJK, M)/M_PM)) * (one -0.4d0 * &
1269 (one+C_E)*(one-3d0*C_E)*EP_SM*Chi)
1270
1271
1272 (IJK) = 25.6d0/pi * EP_SM**2 * Chi *(one+C_E) * &
1273 (one - A2_gtsh(ijk)/16d0)*eta0
1274
1275
1276 (IJK) = etaK*(one+0.8d0*EP_SM*Chi*(one+C_E)) + &
1277 0.6d0*Mu_b_v(IJK)
1278
1279
1280 (IJK) = Mu_b_v(IJK) - (2.d0/3.d0)*Mu_s_v(IJK)
1281
1282
1283
1284
1285
1286
1287
1288
1289 = -0.5d0*xsi_gtsh(ijk)/(M_pm*theta_m(ijk,m))
1290
1291 dGama_dT = 3d0*pi*D_PM**2*RO_g(ijk)*K_phi(EP_SM)/ &
1292 (2d0*M_pm*dsqrt(theta_m(ijk,m)))
1293
1294
1295 = nu0*(one+C_E)/3d0*Chi*( one+2.0625d0*(one-C_E)+ &
1296 ((947d0-579*C_E)/256d0*A2_gtsh(ijk)) )
1297
1298
1299 = 3.75d0*eta0/M_pm
1300
1301
1302
1303 = zero
1304 IF(EP_SM > SMALL_NUMBER) KthK = 2d0/3d0*Kth0*nu0 / (NuK - &
1305 2d0*EDT_s_ip(ijk,M,M) - 2d0*theta_m(ijk,m)*dGama_dT) * &
1306 (one+2d0*A2_gtsh(ijk)+0.6d0*EP_SM*Chi* &
1307 (one+C_E)**2*(2*C_E-one+A2_gtsh(ijk)*(one+C_E)))
1308
1309
1310 (IJK,M) = KthK*(one+1.2d0*EP_SM*Chi*(one+C_E)) + &
1311 (10.24d0/pi* EP_SM**2*Chi*(one+C_E)*(one+0.4375d0* &
1312 A2_gtsh(ijk))*Kth0)
1313
1314
1315
1316 (IJK,M) = M_pm * Kth_s(IJK,M)
1317
1318
1319
1320
1321
1322
1323
1324 = ZERO
1325 IF(EP_SM > SMALL_NUMBER) Rdissdphi = &
1326 1.5d0*dsqrt(EP_SM/2d0)+135d0/64d0*EP_SM*(dlog(EP_SM)+one) +&
1327 11.26d0*EP_SM*(one-10.2*EP_SM+49.71d0*EP_SM**2-87.08d0* &
1328 EP_SM**3) - EP_SM*dlog(epM)*(Chi+EP_SM*dChiOdphi)
1329
1330 = EP_SM*(0.212d0*0.142d0/(EP_SM**0.788d0*&
1331 (one-EP_SM)**4.454d0) + 4.454d0*K_phi(EP_SM)/(one-EP_SM))
1332 Kphidphi = zero
1333
1334 = ro_g(ijk)*D_PM*dsqrt(theta_m(ijk,m)) / mu_g(ijk)
1335
1336
1337 = 3d0*pi*D_pm*Mu_g(ijk)*(Rdissdphi+Re_T*Kphidphi)
1338
1339
1340
1341
1342
1343
1344
1345 = zero
1346 IF((EP_SM > SMALL_NUMBER) .AND. (EP_SM <= 0.4d0)) THEN
1347 denom = one+0.681d0*EP_SM-8.48d0*EP_SM**2+8.16d0*EP_SM**3
1348 dRdphi = (1.5d0*dsqrt(EP_SM/2d0)+135d0/64d0*EP_SM*&
1349 (dlog(EP_SM)+one)+ 17.14d0*EP_SM)/denom - EP_SM*&
1350 (one+3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*&
1351 dlog(EP_SM)+17.14*EP_SM)/denom**2 * &
1352 (0.681d0-16.96d0*EP_SM+24.48d0*EP_SM**2)
1353 ELSEIF(EP_SM > 0.4d0) THEN
1354 dRdphi = 10d0*(one+2d0*EP_SM)/(one-EP_SM)**4
1355 ENDIF
1356
1357
1358 = zero
1359 IF(EP_SM >= 0.1d0) THEN
1360 R_dphi = R_d(EP_SM)
1361 denom = one+3.5d0*dsqrt(EP_SM)+5.9d0*EP_SM
1362 dSdphi = 2d0*R_dphi*dRdphi/(Chi*denom) - &
1363 EP_SM*R_dphi**2 * (dChiOdphi/(Chi**2*denom) + &
1364 (1.75d0/dsqrt(EP_SM)+5.9d0)/(Chi*denom**2))
1365 ENDIF
1366
1367
1368 = M_pm/(3d0*pi*mu_g(ijk)*D_pm)
1369
1370
1371 = dsqrt(pi)*D_pm**4*VREL**2 / &
1372 (36d0*Tau_st**2*dsqrt(theta_m(ijk,m))) * dSdphi
1373
1374
1375 = ZERO
1376 IF(EP_SM > SMALL_NUMBER) Muk = Kth0*Nu0*M_pm*&
1377 theta_m(ijk,m)/NU_PM / (NuK-1.5d0*(EDT_s_ip(ijk,M,M)-&
1378 xsi_gtsh(ijk)/theta_m(ijk,m))) * ( KthK/(Kth0*Nu0)*&
1379 (2d0/M_pm*dGamadn-ro_s(IJK,m)/(M_pm* theta_m(ijk,m))*&
1380 dPsidn + EDT_s_ip(ijk,M,M)*(one+EP_SM/Chi*dChiOdphi)) +&
1381 2d0/3d0*A2_gtsh(ijk) + 0.8d0*EP_SM*Chi* (one+C_E)*&
1382 (one+0.5d0*EP_SM/Chi*dChiOdphi)*(C_E*(C_E-one)+ &
1383 A2_gtsh(ijk)/6d0*(16d0-3d0*C_E+3d0*C_E**2)))
1384
1385
1386 (IJK,M) = Muk*(one+1.2d0*EP_SM*Chi*(one+C_E))
1387
1388 ENDIF
1389 CONTINUE
1390
1391 RETURN
1392 END SUBROUTINE GT_PDE_GTSH
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414 Subroutine gt_pde_ia (M)
1415
1416
1417
1418
1419 USE compar
1420 USE constant
1421 USE drag
1422 USE fldvar
1423 USE functions
1424 USE geometry
1425 USE indices
1426 USE kintheory
1427 USE param
1428 USE param1
1429 USE physprop
1430 USE rdf
1431 USE run
1432 USE toleranc
1433 USE turb
1434 USE ur_facs
1435 USE visc_s
1436 IMPLICIT NONE
1437
1438
1439
1440
1441 INTEGER, INTENT(IN) :: M
1442
1443
1444
1445
1446 INTEGER :: IJK, I, J, K
1447
1448 INTEGER :: L
1449
1450 DOUBLE PRECISION :: Mu_star, Mu_s_dil, Kth_star, K_s_dil
1451
1452 DOUBLE PRECISION :: P_s_sum, P_s_MM, P_s_LM
1453 DOUBLE PRECISION :: MU_common_term, K_common_term
1454 DOUBLE PRECISION :: Mu_sM_sum, MU_s_MM, MU_s_LM, MU_sM_LM, MU_sL_LM
1455 DOUBLE PRECISION :: XI_sM_sum, XI_s_v
1456 DOUBLE PRECISION :: M_PM, M_PL, MPSUM, NU_PL, NU_PM, D_PM, D_PL, DPSUMo2
1457 DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R1p_lm, R8p_lm, R9p_lm, Bp_lm,&
1458 R5p_lm, R6p_lm, R7p_lm
1459 DOUBLE PRECISION :: K_s_sum, K_s_MM, K_s_LM
1460
1461 DOUBLE PRECISION :: SUM_EpsGo
1462
1463 DOUBLE PRECISION :: Kth_sL_iptmp
1464
1465
1466 DO 200 IJK = ijkstart3, ijkend3
1467 I = I_OF(IJK)
1468 J = J_OF(IJK)
1469 K = K_OF(IJK)
1470
1471 IF ( FLUID_AT(IJK) ) THEN
1472
1473
1474
1475 IF(SWITCH_IA) THEN
1476 SUM_EpsGo = ZERO
1477 DO L = 1, SMAX
1478 SUM_EpsGo = SUM_EpsGo+EP_s(IJK,L)*G_0(IJK,M,L)
1479 ENDDO
1480 ELSE
1481 SUM_EpsGo = EP_s(IJK,M)*G_0(IJK,M,M)
1482 ENDIF
1483
1484 P_s_sum = ZERO
1485 Mu_sM_sum = ZERO
1486 XI_sM_sum = ZERO
1487
1488 D_PM = D_P(IJK,M)
1489 M_PM = (PI/6.d0)*D_PM**3 * RO_S(IJK,M)
1490 NU_PM = ROP_S(IJK,M)/M_PM
1491
1492 P_s_MM = NU_PM*Theta_m(IJK,M)
1493
1494 MU_s_dil = (5.d0/96.d0)*D_PM* RO_S(IJK,M)*&
1495 DSQRT(PI*Theta_m(IJK,M)/M_PM)
1496
1497 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1498 Mu_star = MU_s_dil
1499 ELSEIF(Theta_m(IJK,M)/M_PM < SMALL_NUMBER)THEN
1500 Mu_star = ZERO
1501 ELSEIF(EP_S(IJK,M) <= DIL_EP_s) THEN
1502 Mu_star = MU_s_dil*EP_s(IJK,M)*G_0(IJK,M,M)/ &
1503 (SUM_EpsGo + 2.0d0*DgA_s(IJK,M)*MU_s_dil &
1504 / (RO_S(IJK,M)**2 *(Theta_m(IJK,M)/M_PM)))
1505 ELSE
1506 Mu_star = MU_s_dil*EP_S(IJK,M)*G_0(IJK,M,M)/ &
1507 (SUM_EpsGo + 2.0d0*F_gs(IJK,M)*MU_s_dil &
1508 / (RO_S(IJK,M)**2 *EP_s(IJK,M)*(Theta_m(IJK,M)/M_PM)))
1509 ENDIF
1510
1511 MU_s_MM = (Mu_star/G_0(IJK,M,M))*&
1512 (1.d0+(4.d0/5.d0)*(1.d0+C_E)*SUM_EpsGo)**2
1513
1514 DO L = 1, SMAX
1515 D_PL = D_P(IJK,L)
1516 M_PL = (PI/6.d0)*D_PL**3 * RO_S(IJK,L)
1517 MPSUM = M_PM + M_PL
1518 DPSUMo2 = (D_PM+D_PL)/2.d0
1519 NU_PL = ROP_S(IJK,L)/M_PL
1520
1521 IF ( L .eq. M) THEN
1522 Ap_lm = MPSUM/(2.d0)
1523 Dp_lm = M_PL*M_PM/(2.d0*MPSUM)
1524 R0p_lm = ONE/( Ap_lm**1.5 * Dp_lm**2.5 )
1525 R1p_lm = ONE/( Ap_lm**1.5 * Dp_lm**3 )
1526
1527 P_s_LM = PI*(DPSUMo2**3 / 48.d0)*G_0(IJK,M,L)*&
1528 (M_PM*M_PL/MPSUM)* (M_PM*M_PL)**1.5 *&
1529 NU_PM*NU_PL*(1.d0+C_E)*R0p_lm*Theta_m(IJK,M)
1530
1531 MU_s_LM = DSQRT(PI)*( DPSUMo2**4 / 240d0 )*&
1532 G_0(IJK,M,L)*(M_PL*M_PM/MPSUM)**2 *&
1533 (M_PL*M_PM)**1.5 * NU_PM*NU_PL*&
1534 (1.d0+C_E) * R1p_lm * DSQRT(Theta_m(IJK,M))
1535
1536
1537 (IJK,M,L) = (MU_s_MM + MU_s_LM)
1538
1539
1540 (IJK,M,L) = MU_s_LM
1541
1542
1543
1544 (IJK,M,L) = (5.d0/3.d0)*MU_s_LM
1545
1546
1547
1548 (IJK,M,L) = (5.d0/3.d0)*MU_s_LM
1549
1550 ELSE
1551 Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*&
1552 Theta_m(IJK,M))/2.d0
1553 Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
1554 Theta_m(IJK,M) ))/(2.d0*MPSUM)
1555 Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+&
1556 M_PL*Theta_m(IJK,L) ))/&
1557 (2.d0*MPSUM*MPSUM)
1558 R0p_lm = (1.d0/(Ap_lm**1.5 * Dp_lm**2.5))+ &
1559 ((15.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**2.5 *&
1560 Dp_lm**3.5))+&
1561 ((175.d0*(Bp_lm**4))/(8.d0*Ap_lm**3.5 * &
1562 Dp_lm**4.5))
1563 R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
1564 ((9.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**4))+&
1565 ((30.d0*Bp_lm**4) /( 2.d0*Ap_lm**3.5 * &
1566 Dp_lm**5))
1567
1568 P_s_LM = PI*(DPSUMo2**3 / 48.d0)*G_0(IJK,M,L)*&
1569 (M_PM*M_PL/MPSUM)* (M_PM*M_PL)**1.5 *&
1570 NU_PM*NU_PL*(1.d0+C_E)*R0p_lm* &
1571 (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
1572
1573 MU_common_term = DSQRT(PI)*( DPSUMo2**4 / 240d0 )*&
1574 G_0(IJK,M,L)*(M_PL*M_PM/MPSUM)**2 *&
1575 (M_PL*M_PM)**1.5 * NU_PM*NU_PL*&
1576 (1.d0+C_E) * R1p_lm
1577 MU_sM_LM = MU_common_term * Theta_m(IJK,M)**2 *&
1578 Theta_m(IJK,L)**3
1579 MU_sL_LM = MU_common_term * Theta_m(IJK,L)**2 *&
1580 Theta_m(IJK,M)**3
1581
1582
1583
1584 (IJK,M,L) = MU_sM_LM
1585
1586
1587
1588 (IJK,M,L) = MU_sL_LM
1589
1590
1591
1592 (IJK,M,L) = (5.d0/3.d0)*MU_sM_LM
1593
1594
1595
1596 (IJK,M,L) = (5.d0/3.d0)*MU_sL_LM
1597 ENDIF
1598
1599 P_s_sum = P_s_sum + P_s_LM
1600 MU_sM_sum = MU_sM_sum + MU_sM_ip(IJK,M,L)
1601 XI_sM_sum = XI_sM_sum + XI_sM_ip(IJK,M,L)
1602 ENDDO
1603
1604
1605
1606 (IJK) = P_s_sum + P_S_MM
1607
1608
1609
1610 (IJK) = MU_sM_sum + MU_sL_ip(IJK,M,M)
1611 XI_s_v = XI_sM_sum + XI_sL_ip(IJK,M,M)
1612
1613
1614 (IJK) = -(2.d0/3.d0)*Mu_s_v(IJK) + XI_s_v
1615
1616
1617
1618 = ZERO
1619
1620 K_s_dil = (75.d0/384.d0)*D_PM* RO_S(IJK,M)*&
1621 DSQRT(PI*Theta_m(IJK,M)/M_PM)
1622
1623 IF(SWITCH == ZERO .OR. RO_G(IJK) == ZERO) THEN
1624 Kth_star = K_s_dil
1625 ELSEIF(Theta_m(IJK,M)/M_PM < SMALL_NUMBER)THEN
1626 Kth_star = ZERO
1627
1628 ELSEIF(EP_S(IJK,M) <= DIL_EP_s) THEN
1629 Kth_star = K_s_dil*EP_s(IJK,M)*G_0(IJK,M,M)/ &
1630 (SUM_EpsGo+ 1.2d0*DgA_s(IJK,M)*K_s_dil &
1631 / (RO_S(IJK,M)**2 *(Theta_m(IJK,M)/M_PM)))
1632 ELSE
1633 Kth_star = K_s_dil*EP_S(IJK,M)*G_0(IJK,M,M)/ &
1634 (SUM_EpsGo+ 1.2d0*F_gs(IJK,M)*K_s_dil &
1635 / (RO_S(IJK,M)**2 *EP_s(IJK,M)*(Theta_m(IJK,M)/M_PM)))
1636 ENDIF
1637
1638
1639 = (Kth_star/(M_PM*G_0(IJK,M,M)))*&
1640 (1.d0+(3.d0/5.d0)*(1.d0+C_E)*(1.d0+C_E)*SUM_EpsGo)**2
1641
1642 DO L = 1, SMAX
1643 D_PL = D_P(IJK,L)
1644 M_PL = (PI/6.d0)*D_PL**3 *RO_S(IJK,L)
1645 MPSUM = M_PM + M_PL
1646 DPSUMo2 = (D_PM+D_PL)/2.d0
1647 NU_PL = ROP_S(IJK,L)/M_PL
1648
1649 IF ( L .eq. M) THEN
1650
1651
1652
1653
1654
1655 (IJK,M,L) = ZERO
1656
1657
1658
1659
1660
1661
1662
1663 (IJK,M,L) = ZERO
1664
1665
1666
1667 (IJK,M,L) = ZERO
1668
1669
1670
1671 = K_s_sum + K_s_MM
1672
1673 ELSE
1674 Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/2.d0
1675 Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-&
1676 Theta_m(IJK,M) ))/(2.d0*MPSUM)
1677 Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+&
1678 M_PL*Theta_m(IJK,L) ))/(2.d0*MPSUM*MPSUM)
1679
1680 R0p_lm = (1.d0/(Ap_lm**1.5 * Dp_lm**2.5))+&
1681 ((15.d0*Bp_lm*Bp_lm)/(2.d0* Ap_lm**2.5 * Dp_lm**3.5))+&
1682 ((175.d0*(Bp_lm**4))/(8.d0*Ap_lm**3.5 * Dp_lm**4.5))
1683
1684 R1p_lm = (1.d0/((Ap_lm**1.5)*(Dp_lm**3)))+ &
1685 ((9.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**4))+&
1686 ((30.d0*Bp_lm**4)/(2.d0*Ap_lm**3.5 * Dp_lm**5))
1687
1688 R5p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3 ) )+ &
1689 ((5.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
1690 ((14.d0*Bp_lm**4)/(Ap_lm**4.5 * Dp_lm**5))
1691
1692 R6p_lm = (1.d0/(Ap_lm**3.5 * Dp_lm**3))+ &
1693 ((7.d0*Bp_lm*Bp_lm)/(Ap_lm**4.5 * Dp_lm**4))+&
1694 ((126.d0*Bp_lm**4)/(5.d0*Ap_lm**5.5 * Dp_lm**5))
1695
1696 R7p_lm = (3.d0/(2.d0*Ap_lm**2.5 * Dp_lm**4))+ &
1697 ((10.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**5))+&
1698 ((35.d0*Bp_lm**4)/(Ap_lm**4.5 * Dp_lm**6))
1699
1700 R8p_lm = (1.d0/(2.d0*Ap_lm**1.5 * Dp_lm**4))+ &
1701 ((6.d0*Bp_lm*Bp_lm)/(Ap_lm**2.5 * Dp_lm**5))+&
1702 ((25.d0*Bp_lm**4)/(Ap_lm**3.5 * Dp_lm**6))
1703
1704 R9p_lm = (1.d0/(Ap_lm**2.5 * Dp_lm**3))+ &
1705 ((15.d0*Bp_lm*Bp_lm)/(Ap_lm**3.5 * Dp_lm**4))+&
1706 ((70.d0*Bp_lm**4)/(Ap_lm**4.5 * Dp_lm**5))
1707
1708 K_common_term = DPSUMo2**3 * M_PL*M_PM/(2.d0*MPSUM)*&
1709 (1.d0+C_E)*G_0(IJK,M,L) * (M_PM*M_PL)**1.5
1710
1711
1712
1713 = - K_common_term*NU_PM*NU_PL*(&
1714 ((DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)+&
1715 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1716 (3.d0/2.d0)*R1p_lm)-(&
1717 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM/8.d0)*Bp_lm*R6p_lm)+&
1718 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1719 8.d0)*M_PM*R9p_lm)+&
1720 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1721 M_PL*Bp_lm*R7p_lm)+&
1722 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1723 2.d0)*(M_PM/MPSUM)**2 * M_PL*R8p_lm)+&
1724 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1725 R9p_lm)-&
1726 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1727 (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*Theta_M(IJK,L) )*&
1728 (Theta_M(IJK,M)**2 * Theta_M(IJK,L)**3)
1729
1730
1731
1732
1733 = K_common_term*NU_PM*NU_PL*(&
1734 (-(DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)-&
1735 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1736 (3.d0/2.d0)*R1p_lm)+(&
1737 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL/8.d0)*Bp_lm*R6p_lm)+&
1738 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1739 8.d0)*M_PL*R9p_lm)+&
1740 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1741 M_PM*Bp_lm*R7p_lm)+&
1742 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1743 2.d0)*(M_PM/MPSUM)**2 * M_PM*R8p_lm)+&
1744 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1745 R9p_lm)+&
1746 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1747 (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*Theta_M(IJK,M) )*&
1748 (Theta_M(IJK,L)**2 * Theta_M(IJK,M)**3)
1749
1750 Kth_sL_ip(IJK,M,L) = (ONE-UR_Kth_sml)*Kth_sL_ip(IJK,M,L) +&
1751 UR_Kth_sml * Kth_sl_iptmp
1752
1753
1754
1755 (IJK,M,L) = K_common_term*NU_PM*NU_PL*&
1756 (M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(3.d0*PI/10.d0)*&
1757 R0p_lm * (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
1758
1759
1760
1761 (IJK,M,L) = K_common_term*(&
1762 ((DPSUMo2*DSQRT(PI)/16.d0)*Bp_lm*R5p_lm)+&
1763 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1764 R1p_lm) )* (Theta_m(IJK,M)*Theta_m(IJK,L))**3
1765
1766
1767 (IJK,M,L) = NU_PL * Knu_sL_ip(IJK,M,L)
1768 Knu_sL_ip(IJK,M,L) = NU_PM * Knu_sL_ip(IJK,M,L)
1769 K_s_sum = K_s_sum + K_s_LM
1770 ENDIF
1771 ENDDO
1772
1773
1774 (IJK,M) = K_s_sum
1775
1776 ENDIF
1777 CONTINUE
1778
1779 RETURN
1780 END SUBROUTINE GT_PDE_IA
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801 Subroutine Friction_princeton(M)
1802
1803
1804
1805
1806 USE compar
1807 USE constant
1808 USE fldvar
1809 USE functions
1810 USE geometry
1811 USE indices
1812 USE param
1813 USE param1
1814 USE physprop
1815 USE rdf
1816 USE run
1817 USE trace
1818 USE visc_s
1819 IMPLICIT NONE
1820
1821
1822
1823
1824 INTEGER, INTENT(IN) :: M
1825
1826
1827
1828
1829 INTEGER :: IJK
1830
1831 INTEGER :: MM
1832
1833 DOUBLE PRECISION :: Chi, Pc, Mu_zeta,PfoPc, N_Pff
1834 DOUBLE PRECISION :: ZETA
1835
1836 DOUBLE PRECISION :: SUM_EPS_CP
1837
1838 DOUBLE PRECISION :: dpc_dphi, dp_avg
1839
1840
1841 DO 200 IJK = ijkstart3, ijkend3
1842
1843 IF ( FLUID_AT(IJK) ) THEN
1844
1845
1846
1847
1848
1849 IF (EP_g(IJK) .LT. (ONE-eps_f_min)) THEN
1850
1851
1852 = ZERO
1853 dp_avg = ZERO
1854 DO MM=1,SMAX
1855 dp_avg = dp_avg + D_p(IJK,MM)
1856 IF (CLOSE_PACKED(MM)) SUM_EPS_CP=SUM_EPS_CP+EP_S(IJK,MM)
1857 END DO
1858 dp_avg = dp_avg/DBLE(SMAX)
1859
1860
1861 IF (SAVAGE.EQ.1) THEN
1862 = &
1863 ((2d0+ALPHA)/3d0)*((Mu_s_v(IJK)/(Eta*(2d0-Eta)*&
1864 G_0(IJK,M,M)))*(1d0+1.6d0*Eta*EP_s(IJK,M)*&
1865 G_0(IJK,M,M))*(1d0+1.6d0*Eta*(3d0*Eta-2d0)*&
1866 EP_s(IJK,M)*G_0(IJK,M,M))+(0.6d0*Mu_b_v(IJK)*Eta))
1867 ZETA = &
1868 ((48d0*Eta*(1d0-Eta)*RO_S(IJK,M)*EP_s(IJK,M)*&
1869 EP_s(IJK,M)*G_0(IJK,M,M)*&
1870 (Theta_m(IJK,M)**1.5d0))/&
1871 (SQRT_Pi*D_p(IJK,M)*2d0*Mu_zeta))**0.5d0
1872
1873 ELSEIF (SAVAGE.EQ.0) THEN
1874 = (SMALL_NUMBER +&
1875 trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
1876 trD_s_C(IJK,M))/3.d0))**0.5d0
1877
1878 ELSE
1879 IF(KT_TYPE_ENUM == GHD_2007) THEN
1880 ZETA = ((Theta_m(IJK,M)/dp_avg**2) +&
1881 (trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
1882 trD_s_C(IJK,M))/3.d0)))**0.5d0
1883 ELSE
1884 ZETA = ((Theta_m(IJK,M)/D_p(IJK,M)**2) +&
1885 (trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
1886 trD_s_C(IJK,M))/3.d0)))**0.5d0
1887 ENDIF
1888 ENDIF
1889
1890
1891 IF ((ONE-EP_G(IJK)) .GT. ((ONE-ep_star_array(ijk))-delta)) THEN
1892
1893 = (to_SI*Fr)*((delta**5)*(2d0*(ONE-&
1894 ep_star_array(IJK)-delta) - 2d0*eps_f_min)+&
1895 ((ONE-ep_star_array(ijk)-delta)-eps_f_min)*&
1896 (5*delta**4))/(delta**10)
1897
1898 Pc = (to_SI*Fr)*(((ONE-ep_star_array(IJK)-delta) -&
1899 EPS_f_min)**N_Pc)/(delta**D_Pc)
1900 Pc = Pc + dpc_dphi*((ONE-EP_G(IJK))+delta-(ONE-&
1901 ep_star_array(IJK)))
1902
1903
1904 ELSE
1905 Pc = (to_SI*Fr)*(((ONE-EP_G(IJK)) - EPS_f_min)**N_Pc) / &
1906 (((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK)) +&
1907 SMALL_NUMBER)**D_Pc)
1908 ENDIF
1909
1910 IF (trD_s_C(IJK,M) .GE. ZERO) THEN
1911 N_Pff = DSQRT(3d0)/(2d0*Sin_Phi)
1912 ELSE
1913 N_Pff = N_Pf
1914 ENDIF
1915
1916 IF ((trD_s_C(IJK,M)/(ZETA*N_Pff*DSQRT(2d0)&
1917 *Sin_Phi)) .GT. 1d0) THEN
1918 P_s_f(IJK) =ZERO
1919 PfoPc = ZERO
1920 ELSEIF(trD_s_C(IJK,M) == ZERO) THEN
1921 P_s_f(IJK) = Pc
1922 PfoPc = ONE
1923 ELSE
1924 P_s_f(IJK) = Pc*(1d0 - (trD_s_C(IJK,M)/(ZETA&
1925 *N_Pff*DSQRT(2d0)*Sin_Phi)))**(N_Pff-1d0)
1926
1927 PfoPc = (1d0 - (trD_s_C(IJK,M)/(ZETA&
1928 *N_Pff*DSQRT(2d0)*Sin_Phi)))**(N_Pff-1d0)
1929 ENDIF
1930
1931 Chi = DSQRT(2d0)*P_s_f(IJK)*Sin_Phi*(N_Pff - (N_Pff-1d0)*&
1932 (PfoPc)**(1d0/(N_Pff-1d0)))
1933
1934 IF (Chi < ZERO) THEN
1935 P_s_f(IJK) = Pc*((N_Pff/(N_Pff-1d0))**(N_Pff-1d0))
1936 Chi = ZERO
1937 ENDIF
1938
1939 Mu_s_f(IJK) = Chi/(2d0*ZETA)
1940 Lambda_s_f(IJK) = -2d0*Mu_s_f(IJK)/3d0
1941
1942
1943
1944 IF(KT_TYPE_ENUM /= GHD_2007) THEN
1945 P_s_f(IJK) = P_s_f(IJK) * (EP_S(IJK,M)/SUM_EPS_CP)
1946 Mu_s_f(IJK) = Mu_s_f(IJK) * (EP_S(IJK,M)/SUM_EPS_CP)
1947 Lambda_s_f(IJK) = Lambda_s_f(IJK) * (EP_S(IJK,M)/SUM_EPS_CP)
1948 ENDIF
1949
1950 ENDIF
1951 ENDIF
1952 CONTINUE
1953
1954
1955 RETURN
1956 END SUBROUTINE FRICTION_PRINCETON
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984 Subroutine subgrid_stress_igci(M)
1985
1986
1987
1988
1989 USE param
1990 USE param1
1991 USE geometry
1992 USE compar
1993 USE fldvar
1994 USE indices
1995 USE visc_s
1996 USE physprop
1997 USE run
1998 USE constant
1999 USE fun_avg
2000 USE functions
2001 IMPLICIT NONE
2002
2003
2004
2005
2006 INTEGER, INTENT(IN) :: M
2007
2008
2009
2010
2011 INTEGER :: IJK, I, J, K
2012
2013 DOUBLE PRECISION :: Mu_sub,ps_sub
2014 DOUBLE PRECISION :: pressurefac
2015 DOUBLE PRECISION :: viscosityfac,ps_kinetic,ps_total, extra_visfactor
2016 DOUBLE PRECISION :: mu_kinetic,mu_total
2017
2018 DOUBLE PRECISION :: wfactor_Ps, wfactor_mus
2019
2020 DOUBLE PRECISION :: Inv_Froude
2021
2022 DOUBLE PRECISION :: vt
2023
2024 DOUBLE PRECISION :: filtersize
2025
2026
2027 DO IJK = ijkstart3, ijkend3
2028
2029 IF ( FLUID_AT(IJK) ) THEN
2030 I = I_OF(IJK)
2031 J = J_OF(IJK)
2032 K = K_OF(IJK)
2033
2034
2035 = ONE
2036 = ONE
2037
2038
2039 = GRAVITY*D_p0(M)*D_p0(M)*(RO_S(IJK,M) - RO_g(IJK)) / &
2040 (18.0d0*MU_G(IJK))
2041
2042
2043 IF(DO_K) THEN
2044 filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
2045 ELSE
2046 filtersize = filter_size_ratio * DSQRT(AXY(IJK))
2047 ENDIF
2048
2049
2050
2051 = filtersize * GRAVITY / vt**2
2052
2053
2054 = 0.48d0*(Inv_Froude**0.86)*&
2055 (ONE-EXP(-Inv_Froude/1.4))
2056 viscosityfac = 0.37d0*(Inv_Froude**1.22)
2057 Extra_visfactor = ONE/(0.28d0*(Inv_Froude**0.43)+ONE)
2058
2059 IF (EP_s(IJK,M) .LE. 0.0131) THEN
2060 Ps_kinetic = -10.4d0*(EP_s(IJK,m)**2)+0.31d0*EP_s(IJK,m)
2061 ELSEIF (EP_s(IJK,M) .LE. 0.290) THEN
2062 Ps_kinetic = -0.185d0*(EP_s(IJK,m)**3)+&
2063 0.066d0*(EP_s(IJK,m)**2)-0.000183d0*EP_s(IJK,m)+&
2064 0.00232d0
2065 ELSEIF (EP_s(IJK,M) .LE. 0.595) THEN
2066 Ps_kinetic = -0.00978d0*EP_s(IJK,m)+0.00615d0
2067 ELSE
2068 Ps_kinetic = -6.62d0*(EP_s(IJK,m)**3)+&
2069 49.5d0*(EP_s(IJK,m)**2)-50.3d0*EP_s(IJK,m)+13.8d0
2070 ENDIF
2071
2072 Ps_sub = pressurefac*(EP_s(IJK,M)-0.59d0)*&
2073 (-1.69d0*EP_s(IJK,M)-4.61d0*(EP_s(IJK,M)**2)+&
2074 11.d0*(EP_s(IJK,M)**3))
2075
2076 IF (Ps_sub .GE. ZERO) THEN
2077 Ps_total=Ps_kinetic+Ps_sub
2078 ELSE
2079 Ps_total=Ps_kinetic
2080 ENDIF
2081
2082 IF (EP_s(IJK,M) .LE. 0.02) THEN
2083 Mu_kinetic = 1720.d0*(EP_s(IJK,m)**4)-&
2084 215.d0*(EP_s(IJK,m)**3) + 9.81d0*(EP_s(IJK,m)**2)-&
2085 0.207d0*EP_s(IJK,m)+0.00254d0
2086 ELSEIF (EP_s(IJK,M) .LE. 0.2) THEN
2087 Mu_kinetic = 2.72d0*(EP_s(IJK,m)**4)-&
2088 1.55d0*(EP_s(IJK,m)**3)+0.329d0*(EP_s(IJK,m)**2)-&
2089 0.0296d0*EP_s(IJK,m)+0.00136d0
2090 ELSEIF (EP_s(IJK,M) .LE. 0.6095) THEN
2091 Mu_kinetic = -0.0128d0*(EP_s(IJK,m)**3)+&
2092 0.0107d0*(EP_s(IJK,m)**2)-0.0005d0*EP_s(IJK,m)+&
2093 0.000335d0
2094 ELSE
2095 Mu_kinetic = 23.6d0*(EP_s(IJK,m)**2)-&
2096 28.0d0*EP_s(IJK,m)+8.30d0
2097 ENDIF
2098
2099 Mu_sub = Extra_visfactor*viscosityfac*&
2100 (EP_s(IJK,M)-0.59d0)*(-1.22d0*EP_s(IJK,M)-&
2101 0.7d0*(EP_s(IJK,M)**2)-2.d0*(EP_s(IJK,M)**3))
2102
2103 IF (Mu_sub .GE. ZERO) THEN
2104 Mu_total = Mu_kinetic+Mu_sub
2105 ELSE
2106 Mu_total = Mu_kinetic
2107 ENDIF
2108
2109 IF (SUBGRID_WALL) THEN
2110 CALL SUBGRID_STRESS_WALL(wfactor_Ps,wfactor_Mus,vt,IJK)
2111 ENDIF
2112
2113
2114 (IJK) = Ps_total * wfactor_Ps * (vt**2) * &
2115 RO_S(IJK,M)
2116
2117
2118 (IJK) = Mu_total * wfactor_mus * (vt**3) * &
2119 RO_S(IJK,M)/GRAVITY
2120
2121
2122
2123 IF (P_s_v(IJK) .LE. SMALL_NUMBER) P_s_v(IJK) = SMALL_NUMBER
2124 IF (Mu_s_v(IJK) .LE. SMALL_NUMBER) Mu_s_v(IJK)= SMALL_NUMBER
2125
2126
2127 (IJK) = (-2.0d0/3.0d0)*Mu_s_v(IJK)
2128
2129
2130 (IJK, M) = ZERO
2131
2132 ENDIF
2133
2134 ENDDO
2135
2136 RETURN
2137 END SUBROUTINE SUBGRID_STRESS_IGCI
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164 Subroutine subgrid_stress_MILIOLI(M)
2165
2166
2167
2168
2169 USE param
2170 USE param1
2171 USE geometry
2172 USE compar
2173 USE fldvar
2174 USE indices
2175 USE visc_s
2176 USE physprop
2177 USE run
2178 USE constant
2179 USE fun_avg
2180 USE functions
2181 IMPLICIT NONE
2182
2183
2184
2185
2186 INTEGER, INTENT(IN) :: M
2187
2188
2189
2190
2191 INTEGER :: IJK, I, J, K
2192
2193 DOUBLE PRECISION :: cvisc_pot,cvisc_num,cvisc_den,Cvisc
2194 DOUBLE PRECISION :: cpress_pot,cpress_num,cpress_den,Cpress
2195
2196 DOUBLE PRECISION :: wfactor_Ps, wfactor_mus
2197
2198 DOUBLE PRECISION :: Inv_Froude
2199
2200 DOUBLE PRECISION :: vt
2201
2202 DOUBLE PRECISION :: filtersize
2203
2204
2205 DO IJK = ijkstart3, ijkend3
2206
2207 IF ( FLUID_AT(IJK) ) THEN
2208 I = I_OF(IJK)
2209 J = J_OF(IJK)
2210 K = K_OF(IJK)
2211
2212
2213 = ONE
2214 = ONE
2215
2216
2217 = GRAVITY*D_p0(M)*D_p0(M)*(RO_S(IJK,M) - RO_g(IJK)) / &
2218 (18.0d0*MU_G(IJK))
2219
2220
2221 IF(DO_K) THEN
2222 filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
2223 ELSE
2224 filtersize = filter_size_ratio * DSQRT(AXY(IJK))
2225 ENDIF
2226
2227
2228 = filtersize * GRAVITY / vt**2
2229
2230
2231 = (0.59d0-(ONE-EP_g(IJK)))
2232 cvisc_num = (0.7d0*(ONE-EP_g(IJK))*cvisc_pot)
2233 cvisc_den = (0.8d0+(17.d0*cvisc_pot*cvisc_pot*cvisc_pot))
2234 IF ((ONE-EP_g(IJK)) .GE. ZERO .AND. &
2235 (ONE-EP_g(IJK)) .LE. 0.59) THEN
2236 cvisc=(cvisc_num/cvisc_den)
2237 ELSE
2238 cvisc=ZERO
2239 ENDIF
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253 = (0.59d0-(ONE-EP_g(IJK)))
2254 cpress_num = (0.4d0*(ONE-EP_g(IJK))*cpress_pot)
2255 cpress_den = (0.5d0+(13.d0*cpress_pot*cpress_pot*cpress_pot))
2256 IF ((ONE-EP_g(IJK)) .GE. ZERO .AND. &
2257 (ONE-EP_g(IJK)) .LE. 0.59) THEN
2258 cpress = (cpress_num/cpress_den)
2259 ELSE
2260 cpress = ZERO
2261 ENDIF
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274 IF (SUBGRID_WALL) THEN
2275 CALL SUBGRID_STRESS_WALL(wfactor_Ps,wfactor_Mus,vt,IJK)
2276 ENDIF
2277
2278
2279 (IJK) = RO_S(IJK,M) * Inv_froude**(2/7) * &
2280 filtersize**2 * DSQRT( I2_devD_s(IJK) )**2 * &
2281 cpress * wfactor_Ps
2282
2283
2284 (IJK) = RO_S(IJK,M) * filtersize**2 * &
2285 DSQRT( I2_devD_s(IJK) ) * cvisc * wfactor_mus
2286
2287
2288
2289 IF (P_s_v(IJK) .LE. SMALL_NUMBER) P_s_v(IJK) = SMALL_NUMBER
2290 IF (Mu_s_v(IJK) .LE. SMALL_NUMBER) Mu_s_v(IJK)= SMALL_NUMBER
2291
2292
2293 (IJK) = (-2.0d0/3.0d0)*Mu_s_v(IJK)
2294
2295
2296 (IJK, M) = ZERO
2297
2298 ENDIF
2299
2300
2301 ENDDO
2302
2303 RETURN
2304 END SUBROUTINE SUBGRID_STRESS_MILIOLI
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330 Subroutine subgrid_stress_wall(lfactor_ps, lfactor_mus, vt, &
2331 IJK)
2332
2333
2334
2335
2336 USE param1, only: ONE
2337 USE constant, only : GRAVITY
2338 USE cutcell, only : DWALL
2339 IMPLICIT NONE
2340
2341
2342
2343
2344 DOUBLE PRECISION, INTENT(OUT) :: lfactor_ps
2345
2346 DOUBLE PRECISION, INTENT(OUT) :: lfactor_mus
2347
2348 DOUBLE PRECISION, INTENT(IN) :: vt
2349
2350 INTEGER, INTENT(IN) :: IJK
2351
2352
2353
2354
2355 DOUBLE PRECISION, PARAMETER :: aps=9.14d0, bps=0.345d0,&
2356 amus=5.69d0, bmus=0.228d0
2357
2358
2359
2360
2361 DOUBLE PRECISION :: x_d
2362
2363
2364
2365 = ONE
2366 lfactor_mus = ONE
2367
2368
2369 = DWALL(IJK) * GRAVITY / vt**2
2370
2371 = ONE / ( ONE + aps * (EXP(-bps*x_d)) )
2372
2373 = ONE / ( ONE + amus * (EXP(-bmus*x_d)) )
2374
2375 RETURN
2376 END SUBROUTINE SUBGRID_STRESS_WALL
2377
2378
2379
2380
2381
2382
2383
2384 Subroutine add_shear(M)
2385
2386
2387
2388
2389
2390 USE fldvar, only: v_s
2391
2392 USE vshear, only: vsh
2393
2394 USE compar
2395 USE geometry
2396 USE indices
2397 USE functions
2398 IMPLICIT NONE
2399
2400
2401
2402
2403 INTEGER, INTENT(IN) :: M
2404
2405
2406
2407
2408 INTEGER :: IJK
2409
2410
2411
2412 DO IJK= ijkstart3, ijkend3
2413 IF (FLUID_AT(IJK)) THEN
2414 V_s(ijk,m)=V_s(IJK,m)+VSH(IJK)
2415 ENDIF
2416 ENDDO
2417
2418 RETURN
2419 END SUBROUTINE ADD_SHEAR
2420
2421
2422
2423
2424
2425
2426 Subroutine remove_shear(M)
2427
2428
2429
2430
2431
2432 USE fldvar, only: v_s
2433
2434 USE vshear, only: vsh
2435
2436 USE compar
2437 USE geometry
2438 USE indices
2439 USE functions
2440 IMPLICIT NONE
2441
2442
2443
2444
2445 INTEGER, INTENT(IN) :: M
2446
2447
2448
2449
2450 INTEGER :: IJK
2451
2452
2453
2454 DO IJK= ijkstart3, ijkend3
2455 IF (FLUID_AT(IJK)) THEN
2456 V_s(IJK,m)=V_s(IJK,m)-VSH(IJK)
2457 ENDIF
2458 ENDDO
2459
2460 RETURN
2461 END SUBROUTINE REMOVE_SHEAR
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474 Subroutine init_mu_s (M, IER)
2475
2476
2477
2478
2479
2480 USE param1, only: zero, one, half, undefined, small_number
2481
2482
2483 USE run, only: granular_energy
2484
2485 USE run, only: kt_type_enum
2486 USE run, only: simonin_1996, ahmadi_1995
2487 USE run, only: ia_2005, ghd_2007
2488
2489 USE run, only: shear, V_sh
2490
2491
2492 Use vshear, only: VSH
2493
2494 use drag, only: dgA_s
2495
2496
2497 USE visc_s, only: mu_s, mu_s_v, mu_s_p, mu_s_f, mu_b_v
2498 USE visc_s, only: lambda_s, lambda_s_v, lambda_s_p, lambda_s_f
2499
2500 USE visc_s, only: ep_star_array
2501
2502 USE visc_s, only: VREL_array
2503
2504 USE visc_s, only: I2_devD_s
2505
2506 USE visc_s, only: trm_s, trdM_s, alpha_s
2507
2508
2509 USE trace, only: trD_s_c
2510
2511 USE trace, only: trD_s2
2512
2513 USE trace, only: trD_s2_ip
2514
2515
2516 USE fldvar, only: P_s, p_s_v, p_s_p, p_s_f
2517
2518 USE fldvar, only: u_s, v_s, w_s
2519
2520 USE fldvar, only: u_g, v_g, w_g
2521
2522 USE fldvar, only: ep_g, rop_g, ro_g
2523
2524
2525 USE fldvar, only: rop_s, ro_s, d_p
2526
2527 USE fldvar, only: e_turb_g, k_turb_g
2528
2529 USE fldvar, only: ep_s
2530
2531
2532 USE constant, only: v_ex
2533
2534
2535 USE constant, only: switch
2536
2537
2538 USE physprop, only: smax
2539
2540 USE physprop, only: mu_s0
2541
2542 USE physprop, only: mu_g
2543
2544 USE is, only: any_is_defined
2545
2546
2547 USE turb, only: cos_theta, tau_1
2548
2549
2550 USE toleranc, only: zero_ep_s
2551
2552
2553 USE cutcell, only: cut_cell_at
2554
2555
2556 USE compar
2557 USE fun_avg
2558 USE functions
2559 USE geometry
2560 USE indices
2561 USE rdf
2562
2563 IMPLICIT NONE
2564
2565
2566
2567
2568 INTEGER, INTENT(IN) :: M
2569
2570 INTEGER, INTENT(INOUT) :: IER
2571
2572
2573
2574
2575 DOUBLE PRECISION, PARAMETER :: C_mu = 9.0D-02
2576
2577
2578
2579
2580 DOUBLE PRECISION :: D_s(3,3), D_sl(3,3)
2581
2582 DOUBLE PRECISION :: U_s_N, Usl_N
2583
2584 DOUBLE PRECISION :: U_s_S, Usl_S
2585
2586 DOUBLE PRECISION :: U_s_T, Usl_T
2587
2588 DOUBLE PRECISION :: U_s_B, Usl_B
2589
2590
2591 DOUBLE PRECISION :: U_s_C, Usl_C
2592
2593 DOUBLE PRECISION :: V_s_E, Vsl_E
2594
2595 DOUBLE PRECISION :: V_s_W, Vsl_W
2596
2597 DOUBLE PRECISION :: V_s_T, Vsl_T
2598
2599 DOUBLE PRECISION :: V_s_B, Vsl_B
2600
2601 DOUBLE PRECISION :: W_s_E, Wsl_E
2602
2603 DOUBLE PRECISION :: W_s_W, Wsl_W
2604
2605 DOUBLE PRECISION :: W_s_N, Wsl_N
2606
2607 DOUBLE PRECISION :: W_s_S, Wsl_S
2608
2609
2610 DOUBLE PRECISION :: W_s_C, Wsl_C
2611
2612 DOUBLE PRECISION :: USCM, VSCM, WSCM, &
2613 UGC, VGC, WGC,&
2614 SqrtVs, SqrtVgMinusVs
2615
2616 INTEGER :: I1, I2
2617
2618 DOUBLE PRECISION :: DEP_soDX
2619
2620 DOUBLE PRECISION :: DEP_soDY
2621
2622 DOUBLE PRECISION :: DEP_soXDZ
2623
2624 DOUBLE PRECISION :: M_s(3,3)
2625
2626 INTEGER :: I, J, K, IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP,&
2627 IJKW, IJKE, IJKS, IJKN, IJKB, IJKT,&
2628 IM, JM, KM
2629 INTEGER :: IMJPK, IMJMK, IMJKP, IMJKM, IPJKM, IPJMK, IJMKP,&
2630 IJMKM, IJPKM
2631
2632 INTEGER :: L
2633
2634 DOUBLE PRECISION :: SRT
2635
2636 DOUBLE PRECISION :: C_d, Re
2637
2638
2639 IF(MU_s0 == UNDEFINED) THEN
2640
2641 (:,M) = ZERO
2642 LAMBDA_s(:,M) = ZERO
2643 P_s(:,M) = ZERO
2644
2645 (:) = ZERO
2646 Mu_b_v(:) = ZERO
2647 LAMBDA_s_v(:) = ZERO
2648 P_s_v(:) = ZERO
2649
2650 (:) = ZERO
2651 LAMBDA_s_p(:) = ZERO
2652 (:) = ZERO
2653
2654 (:) = ZERO
2655 LAMBDA_s_f(:) = ZERO
2656 P_s_f(:) = ZERO
2657
2658 (:,M) = ZERO
2659 ENDIF
2660
2661
2662 (:,M) = ZERO
2663 vrel_array(:,M) = ZERO
2664
2665 IF (SHEAR) SRT=(2d0*V_sh/XLENGTH)
2666
2667
2668
2669
2670
2671
2672
2673
2674
2675
2676
2677
2678
2679 DO IJK = ijkstart3, ijkend3
2680
2681 IF ( FLUID_AT(IJK) ) THEN
2682 I = I_OF(IJK)
2683 J = J_OF(IJK)
2684 K = K_OF(IJK)
2685 IM = Im1(I)
2686 JM = Jm1(J)
2687 KM = Km1(K)
2688 IJKW = WEST_OF(IJK)
2689 IJKE = EAST_OF(IJK)
2690 IJKS = SOUTH_OF(IJK)
2691 IJKN = NORTH_OF(IJK)
2692 IJKB = BOTTOM_OF(IJK)
2693 IJKT = TOP_OF(IJK)
2694 IMJK = IM_OF(IJK)
2695 IPJK = IP_OF(IJK)
2696 IJMK = JM_OF(IJK)
2697 IJPK = JP_OF(IJK)
2698 IJKM = KM_OF(IJK)
2699 IJKP = KP_OF(IJK)
2700 IMJPK = IM_OF(IJPK)
2701 IMJMK = IM_OF(IJMK)
2702 IMJKP = IM_OF(IJKP)
2703 IMJKM = IM_OF(IJKM)
2704 IPJKM = IP_OF(IJKM)
2705 IPJMK = IP_OF(IJMK)
2706 IJMKP = JM_OF(IJKP)
2707 IJMKM = JM_OF(IJKM)
2708 IJPKM = JP_OF(IJKM)
2709
2710 U_s_N = AVG_Y( &
2711 (U_s(IMJK, M), U_s(IJK, M), I),&
2712 AVG_X_E(U_s(IMJPK, M), U_s(IJPK, M), I), J)
2713 U_s_S = AVG_Y( &
2714 (U_s(IMJMK, M), U_s(IJMK, M), I),&
2715 AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I), JM)
2716 U_s_T = AVG_Z( &
2717 (U_s(IMJK, M), U_s(IJK, M), I),&
2718 AVG_X_E(U_s(IMJKP, M), U_s(IJKP, M), I), K)
2719 U_s_B = AVG_Z( &
2720 (U_s(IMJKM, M), U_s(IJKM, M), I),&
2721 AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I), KM)
2722
2723 IF (SHEAR) THEN
2724 V_s_E = AVG_X( &
2725 (V_s(IJMK, M), V_s(IJK, M)),&
2726 AVG_Y_N((V_s(IPJMK, M)-VSH(IPJMK)+VSH(IJMK)&
2727 +SRT*ONE/oDX_E(I)),&
2728 (V_s(IPJK, M)-VSH(IPJK)+VSH(IJK)&
2729 +SRT*ONE/oDX_E(I))), I)
2730 V_s_W = AVG_X( &
2731 ((V_s(IMJMK, M)-VSH(IMJMK)+VSH(IJMK)&
2732 -SRT*ONE/oDX_E(IM1(I))),&
2733 (V_s(IMJK, M)-VSH(IMJK)+VSH(IJK)&
2734 -SRT*ONE/oDX_E(IM1(I)))),&
2735 AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)), IM)
2736 ELSE
2737 V_s_E = AVG_X( &
2738 (V_s(IJMK, M), V_s(IJK, M)),&
2739 AVG_Y_N(V_s(IPJMK, M), V_s(IPJK, M)), I )
2740 V_s_W = AVG_X( &
2741 (V_s(IMJMK, M), V_s(IMJK, M)),&
2742 AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)), IM )
2743 ENDIF
2744
2745 V_s_T = AVG_Z( &
2746 (V_s(IJMK, M), V_s(IJK, M)),&
2747 AVG_Y_N(V_s(IJMKP, M), V_s(IJKP, M)), K )
2748 V_s_B = AVG_Z( &
2749 (V_s(IJMKM, M), V_s(IJKM, M)),&
2750 AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)), KM )
2751 W_s_N = AVG_Y( &
2752 (W_s(IJKM, M), W_s(IJK, M)),&
2753 AVG_Z_T(W_s(IJPKM, M), W_s(IJPK, M)), J )
2754 W_s_S = AVG_Y( &
2755 (W_s(IJMKM, M), W_s(IJMK, M)),&
2756 AVG_Z_T(W_s(IJKM, M), W_s(IJK, M)), JM )
2757 W_s_E = AVG_X( &
2758 (W_s(IJKM, M), W_s(IJK, M)),&
2759 AVG_Z_T(W_s(IPJKM, M), W_s(IPJK, M)), I)
2760 W_s_W = AVG_X( &
2761 (W_s(IMJKM, M), W_s(IMJK, M)),&
2762 AVG_Z_T(W_s(IJKM, M), W_s(IJK, M)), IM )
2763
2764 IF(CYLINDRICAL) THEN
2765 U_s_C = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
2766 = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
2767 ELSE
2768 U_s_C = ZERO
2769 W_s_C = ZERO
2770 ENDIF
2771
2772
2773 IF(ANY_IS_DEFINED) THEN
2774 IF(IS_AT_N(IJK) .AND. .NOT.WALL_AT(IJPK)) &
2775 U_s_N = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
2776 IF(IS_AT_N(IJMK) .AND. .NOT.WALL_AT(IJMK)) &
2777 U_s_S = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
2778 IF(IS_AT_T(IJK) .AND. .NOT.WALL_AT(IJKP)) &
2779 U_s_T = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
2780 IF(IS_AT_T(IJKM) .AND. .NOT.WALL_AT(IJKM)) &
2781 U_s_B = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
2782 IF(IS_AT_E(IJK) .AND. .NOT.WALL_AT(IPJK)) &
2783 V_s_E = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
2784 IF(IS_AT_E(IMJK) .AND. .NOT.WALL_AT(IMJK)) &
2785 V_s_W = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
2786 IF(IS_AT_T(IJK) .AND. .NOT.WALL_AT(IJKP)) &
2787 V_s_T = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
2788 IF(IS_AT_T(IJKM) .AND. .NOT.WALL_AT(IJKM)) &
2789 V_s_B = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
2790 IF(IS_AT_N(IJK) .AND. .NOT.WALL_AT(IJPK)) &
2791 W_s_N = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
2792 IF(IS_AT_N(IJMK) .AND. .NOT.WALL_AT(IJMK)) &
2793 W_s_S = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
2794 IF(IS_AT_E(IJK) .AND. .NOT.WALL_AT(IPJK)) &
2795 W_s_E = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
2796 IF(IS_AT_E(IMJK) .AND. .NOT.WALL_AT(IMJK)) &
2797 W_s_W = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
2798 ENDIF
2799
2800
2801
2802 (1,1) = ( U_s(IJK,M) - U_s(IMJK,M) ) * oDX(I)
2803 D_s(1,2) = HALF * ( (U_s_N - U_s_S) * oDY(J) +&
2804 (V_s_E - V_s_W) * oDX(I) )
2805 D_s(1,3) = HALF * ( (W_s_E - W_s_W) * oDX(I) +&
2806 (U_s_T - U_s_B) * (oX(I)*oDZ(K)) - W_s_C * oX(I) )
2807 D_s(2,1) = D_s(1,2)
2808 D_s(2,2) = ( V_s(IJK,M) - V_s(IJMK,M) ) * oDY(J)
2809 D_s(2,3) = HALF * ( (V_s_T - V_s_B) * (oX(I)*oDZ(K)) +&
2810 (W_s_N - W_s_S) * oDY(J) )
2811 D_s(3,1) = D_s(1,3)
2812 D_s(3,2) = D_s(2,3)
2813 D_s(3,3) = ( W_s(IJK,M) - W_s(IJKM,M) ) * (oX(I)*oDZ(K)) +&
2814 U_s_C * oX(I)
2815
2816 IF(CUT_CELL_AT(IJK)) CALL CG_CALC_VEL_S_GRAD(IJK,M,D_s)
2817
2818
2819 (IJK,M) = D_s(1,1) + D_s(2,2) + D_s(3,3)
2820
2821
2822 (IJK,M) = ZERO
2823 DO I1 = 1,3
2824 DO I2 = 1,3
2825 trD_s2(IJK,M) = trD_s2(IJK,M) + D_s(I1,I2)*D_s(I1,I2)
2826 ENDDO
2827 ENDDO
2828
2829
2830 IF (trD_s2(IJK,M) == zero) trD_s_C(IJK,M) = zero
2831
2832
2833
2834 (IJK) = ( (D_s(1,1)-D_s(2,2))**2&
2835 +(D_s(2,2)-D_s(3,3))**2&
2836 +(D_s(3,3)-D_s(1,1))**2 )/6.&
2837 + D_s(1,2)**2 + D_s(2,3)**2 + D_s(3,1)**2
2838
2839
2840
2841
2842 = AVG_X_E(U_G(IMJK),U_G(IJK),I)
2843 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
2844 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
2845
2846 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
2847 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
2848 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
2849
2850
2851 (IJK,M) = SQRT((UGC - USCM)**2 + &
2852 (VGC - VSCM)**2 + (WGC - WSCM)**2)
2853
2854
2855
2856
2857
2858
2859 IF (KT_TYPE_ENUM /= GHD_2007) THEN
2860 RE = D_p(IJK,M)*VREL_array(IJK,M)*ROP_G(IJK)/&
2861 (MU_G(IJK) + SMALL_NUMBER)
2862 IF(RE .LE. 1000.d0)THEN
2863 C_d = (24.d0/(Re+SMALL_NUMBER)) * &
2864 (ONE + 0.15d0 * Re**0.687D0)
2865 ELSE
2866 C_d = 0.44d0
2867 ENDIF
2868 DgA_s(IJK,M) = 0.75d0 * C_d * VREL_array(IJK,M) * &
2869 ROP_g(IJK) / D_p(IJK,M)
2870
2871
2872 IF(VREL_array(IJK,M) == ZERO) DgA_S(IJK,M) = LARGE_NUMBER
2873 ENDIF
2874
2875
2876
2877 IF((KT_TYPE_ENUM == SIMONIN_1996 .OR. &
2878 KT_TYPE_ENUM == AHMADI_1995) .AND. M == 1) THEN
2879
2880 (ijk) = 3.d0/2.d0*C_MU*K_Turb_G(IJK)/&
2881 (E_Turb_G(IJK)+small_number)
2882
2883
2884
2885
2886 IF(KT_TYPE_ENUM == SIMONIN_1996) THEN
2887 SqrtVs = SQRT(USCM**2+VSCM**2+WSCM**2)
2888 SqrtVgMinusVs = VREL_array(IJK,M)
2889 IF(SqrtVs > Small_Number .AND. SqrtVgMinusVs > &
2890 Small_Number .AND. EP_S(IJK,M) > ZERO_EP_S) THEN
2891 Cos_Theta(IJK) = ((UGC-USCM)*USCM+&
2892 (VGC-VSCM)*VSCM+&
2893 (WGC-WSCM)*WSCM)/ (SqrtVgMinusVs * SqrtVs)
2894 ELSE
2895 Cos_Theta(IJK) = ZERO
2896 ENDIF
2897 ENDIF
2898 ENDIF
2899
2900
2901
2902
2903
2904 IF (KT_TYPE_ENUM == IA_2005) THEN
2905 DO L = 1,SMAX
2906 IF (L .NE. M) THEN
2907 IF (L > M) THEN
2908 = AVG_Y(&
2909 (U_s(IMJK, L), U_s(IJK, L), I),&
2910 AVG_X_E(U_s(IMJPK, L), U_s(IJPK, L), I), J)
2911 Usl_S = AVG_Y(&
2912 (U_s(IMJMK, L), U_s(IJMK, L), I),&
2913 AVG_X_E(U_s(IMJK, L), U_s(IJK, L), I), JM)
2914 Usl_T = AVG_Z(&
2915 (U_s(IMJK, L), U_s(IJK, L), I),&
2916 AVG_X_E(U_s(IMJKP, L), U_s(IJKP, L), I), K)
2917 Usl_B = AVG_Z(&
2918 (U_s(IMJKM, L), U_s(IJKM, L), I),&
2919 AVG_X_E(U_s(IMJK, L), U_s(IJK, L), I), KM)
2920
2921 IF (SHEAR) THEN
2922 Vsl_E = AVG_X(&
2923 (V_s(IJMK, L), V_s(IJK, L)),&
2924 AVG_Y_N((V_s(IPJMK, L)-VSH(IPJMK)+&
2925 VSH(IJMK)+SRT*ONE/oDX_E(I)),&
2926 (V_s(IPJK, L)-VSH(IPJK)+VSH(IJK)&
2927 +SRT*ONE/oDX_E(I))), I)
2928 Vsl_W = AVG_X(&
2929 ((V_s(IMJMK, L)-VSH(IMJMK)+&
2930 VSH(IJMK)-SRT*ONE/oDX_E(IM1(I))),&
2931 (V_s(IMJK, L)-VSH(IMJK)+VSH(IJK)&
2932 -SRT*ONE/oDX_E(IM1(I)))),&
2933 AVG_Y_N(V_s(IJMK, L), V_s(IJK, L)), IM)
2934 ELSE
2935 Vsl_E = AVG_X(&
2936 (V_s(IJMK, L), V_s(IJK, L)),&
2937 AVG_Y_N(V_s(IPJMK, L), V_s(IPJK, L)), I )
2938 Vsl_W = AVG_X(&
2939 (V_s(IMJMK, L), V_s(IMJK, L)),&
2940 AVG_Y_N(V_s(IJMK, L), V_s(IJK, L)), IM )
2941 ENDIF
2942
2943 Vsl_T = AVG_Z(&
2944 (V_s(IJMK, L), V_s(IJK, L)),&
2945 AVG_Y_N(V_s(IJMKP, L), V_s(IJKP, L)), K )
2946 Vsl_B = AVG_Z(&
2947 (V_s(IJMKM, L), V_s(IJKM, L)),&
2948 AVG_Y_N(V_s(IJMK, L), V_s(IJK, L)), KM )
2949 Wsl_N = AVG_Y(&
2950 (W_s(IJKM, L), W_s(IJK, L)),&
2951 AVG_Z_T(W_s(IJPKM, L), W_s(IJPK, L)), J )
2952 Wsl_S = AVG_Y(&
2953 (W_s(IJMKM, L), W_s(IJMK, L)),&
2954 AVG_Z_T(W_s(IJKM, L), W_s(IJK, L)), JM )
2955 Wsl_E = AVG_X(&
2956 (W_s(IJKM, L), W_s(IJK, L)),&
2957 AVG_Z_T(W_s(IPJKM, L), W_s(IPJK, L)), I)
2958 Wsl_W = AVG_X(&
2959 (W_s(IMJKM, L), W_s(IMJK, L)),&
2960 AVG_Z_T(W_s(IJKM, L), W_s(IJK, L)), IM )
2961
2962 IF(CYLINDRICAL) THEN
2963 Usl_C = AVG_X_E(U_s(IMJK, L), U_s(IJK, L), I)
2964 = AVG_Z_T(W_s(IJKM, L), W_s(IJK, L))
2965 ELSE
2966 Usl_C = ZERO
2967 Wsl_C = ZERO
2968 ENDIF
2969
2970
2971 IF(ANY_IS_DEFINED) THEN
2972 IF(IS_AT_N(IJK) .AND. .NOT.WALL_AT(IJPK)) &
2973 Usl_N = AVG_X_E(U_s(IMJK,L),U_s(IJK,L), I)
2974 IF(IS_AT_N(IJMK) .AND. .NOT.WALL_AT(IJMK)) &
2975 Usl_S = AVG_X_E(U_s(IMJK,L),U_s(IJK,L), I)
2976 IF(IS_AT_T(IJK) .AND. .NOT.WALL_AT(IJKP)) &
2977 Usl_T = AVG_X_E(U_s(IMJK,L),U_s(IJK,L), I)
2978 IF(IS_AT_T(IJKM) .AND. .NOT.WALL_AT(IJKM)) &
2979 Usl_B = AVG_X_E(U_s(IMJK,L),U_s(IJK,L), I)
2980 IF(IS_AT_E(IJK) .AND. .NOT.WALL_AT(IPJK)) &
2981 Vsl_E = AVG_Y_N(V_s(IJMK,L),V_s(IJK,L))
2982 IF(IS_AT_E(IMJK) .AND. .NOT.WALL_AT(IMJK)) &
2983 Vsl_W = AVG_Y_N(V_s(IJMK,L),V_s(IJK,L))
2984 IF(IS_AT_T(IJK) .AND. .NOT.WALL_AT(IJKP)) &
2985 Vsl_T = AVG_Y_N(V_s(IJMK,L),V_s(IJK,L))
2986 IF(IS_AT_T(IJKM) .AND. .NOT.WALL_AT(IJKM)) &
2987 Vsl_B = AVG_Y_N(V_s(IJMK,L),V_s(IJK,L))
2988 IF(IS_AT_N(IJK) .AND. .NOT.WALL_AT(IJPK)) &
2989 Wsl_N = AVG_Z_T(W_s(IJKM,L),W_s(IJK,L))
2990 IF(IS_AT_N(IJMK) .AND. .NOT.WALL_AT(IJMK)) &
2991 Wsl_S = AVG_Z_T(W_s(IJKM,L),W_s(IJK,L))
2992 IF(IS_AT_E(IJK) .AND. .NOT.WALL_AT(IPJK)) &
2993 Wsl_E = AVG_Z_T(W_s(IJKM,L),W_s(IJK,L))
2994 IF(IS_AT_E(IMJK) .AND. .NOT.WALL_AT(IMJK)) &
2995 Wsl_W = AVG_Z_T(W_s(IJKM,L),W_s(IJK,L))
2996 ENDIF
2997
2998
2999
3000 (1,1) = ( U_s(IJK,L) - U_s(IMJK,L) ) * oDX(I)
3001 D_sl(1,2) = HALF * ( (Usl_N - Usl_S) * oDY(J) +&
3002 (Vsl_E - Vsl_W) * oDX(I) )
3003 D_sl(1,3) = HALF * ( (Wsl_E - Wsl_W) * oDX(I) +&
3004 (Usl_T - Usl_B) * (oX(I)*oDZ(K)) - Wsl_C * oX(I) )
3005 D_sl(2,1) = D_sl(1,2)
3006 D_sl(2,2) = ( V_s(IJK,L) - V_s(IJMK,L) ) * oDY(J)
3007 D_sl(2,3) = HALF * ( (Vsl_T - Vsl_B) * (oX(I)*oDZ(K)) +&
3008 (Wsl_N - Wsl_S) * oDY(J) )
3009 D_sl(3,1) = D_sl(1,3)
3010 D_sl(3,2) = D_sl(2,3)
3011 D_sl(3,3) = ( W_s(IJK,L) - W_s(IJKM,L) ) * (oX(I)*oDZ(K)) +&
3012 Usl_C * oX(I)
3013
3014
3015
3016 (IJK,M,L) = ZERO
3017 DO I1 = 1,3
3018 DO I2 = 1,3
3019 trD_s2_ip(IJK,M,L) = trD_s2_ip(IJK,M,L)+&
3020 D_sl(I1,I2)*D_s(I1,I2)
3021 ENDDO
3022 ENDDO
3023
3024 ELSE
3025 (IJK,M,L) = trD_s2_ip(IJK,L,M)
3026 ENDIF
3027 ELSE
3028 (IJK,M,M) = trD_s2(IJK,M)
3029 ENDIF
3030 ENDDO
3031 ENDIF
3032
3033
3034
3035
3036
3037
3038
3039 IF(.NOT.GRANULAR_ENERGY) THEN
3040 IF(EP_g(IJK) .GE. EP_star_array(IJK)) THEN
3041 IF(V_ex .NE. ZERO) THEN
3042 DEP_soDX = ( EP_s(IJKE, M) - EP_s(IJK, M) ) * oDX_E(I)&
3043 * ( ONE / ( (oDX_E(IM)/oDX_E(I)) + ONE ) ) +&
3044 ( EP_s(IJK, M) - EP_s(IJKW, M) ) * oDX_E(IM)&
3045 * ( ONE / ( (oDX_E(I)/oDX_E(IM)) + ONE ) )
3046 DEP_soDY = ( EP_s(IJKN, M) - EP_s(IJK, M) ) * oDY_N(J)&
3047 * ( ONE / ( (oDY_N(JM)/oDY_N(J)) + ONE ) ) +&
3048 ( EP_s(IJK, M) - EP_s(IJKS, M) ) * oDY_N(JM)&
3049 * ( ONE / ( (oDY_N(J)/oDY_N(JM)) + ONE ) )
3050 DEP_soXDZ = (( EP_s(IJKT, M) - EP_s(IJK, M) )&
3051 * oX(I)*oDZ_T(K)&
3052 * ( ONE / ( (oDZ_T(KM)/oDZ_T(K)) + ONE ) ) +&
3053 ( EP_s(IJK, M) - EP_s(IJKB, M) ) * oX(I)*oDZ_T(KM)&
3054 * ( ONE / ( (oDZ_T(K)/oDZ_T(KM)) + ONE ) ) ) /&
3055 X(I)
3056 M_s(1,1) = DEP_soDX * DEP_soDX
3057 M_s(1,2) = DEP_soDX * DEP_soDY
3058 M_s(1,3) = DEP_soDX * DEP_soXDZ
3059 M_s(2,1) = DEP_soDX * DEP_soDY
3060 M_s(2,2) = DEP_soDY * DEP_soDY
3061 M_s(2,3) = DEP_soDY * DEP_soXDZ
3062 M_s(3,1) = DEP_soDX * DEP_soXDZ
3063 M_s(3,2) = DEP_soDY * DEP_soXDZ
3064 M_s(3,3) = DEP_soXDZ * DEP_soXDZ
3065 trM_s(IJK) = M_s(1,1) + M_s(2,2) + M_s(3,3)
3066 trDM_s(IJK) = ZERO
3067 DO I1 = 1,3
3068 DO I2 = 1,3
3069 trDM_s(IJK) = trDM_s(IJK) + D_s(I1,I2)*M_s(I1,I2)
3070 ENDDO
3071 ENDDO
3072 ENDIF
3073 ENDIF
3074 ENDIF
3075
3076
3077 ENDIF
3078 ENDDO
3079
3080
3081 RETURN
3082 END SUBROUTINE INIT_MU_S
3083
3084