File: /nfs/home/0/users/jenkins/mfix.git/model/drag_gs.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21 SUBROUTINE DRAG_GS(M, IER)
22
23
24
25
26 USE param
27 USE param1
28 USE parallel
29 USE fldvar
30 USE geometry
31 USE indices
32 USE physprop
33 USE run
34 USE constant
35 USE compar
36 USE drag
37 USE sendrecv
38 USE discretelement
39 USE ur_facs
40 USE funits
41 USE mms
42 USE fun_avg
43 USE functions
44 IMPLICIT NONE
45
46
47
48
49 INTEGER, INTENT(IN) :: M
50
51 INTEGER, INTENT(INOUT) :: IER
52
53
54
55
56
57
58
59
60 INTEGER :: I, IJK, IMJK, IJMK, IJKM
61
62 DOUBLE PRECISION :: USCM, VSCM, WSCM
63
64 DOUBLE PRECISION :: USCM_HYS, VSCM_HYS, WSCM_HYS
65
66 DOUBLE PRECISION :: UGC, VGC, WGC
67
68 DOUBLE PRECISION :: VREL
69
70
71 DOUBLE PRECISION :: Mu
72
73 DOUBLE PRECISION :: DgA
74
75 DOUBLE PRECISION F_gstmp
76
77 INTEGER :: CM, DM, L
78
79
80 INTEGER :: MAXM
81
82
83 DOUBLE PRECISION :: DP_loc(2*DIM_M)
84
85
86 DOUBLE PRECISION :: EPs_loc(2*DIM_M)
87
88
89 DOUBLE PRECISION :: ROs_loc(2*DIM_M)
90
91
92 DOUBLE PRECISION :: F_cor, tmp_sum, tmp_fac
93
94 DOUBLE PRECISION :: DPA
95
96 DOUBLE PRECISION :: Y_i
97
98 DOUBLE PRECISION :: phis
99
100
101
102 DOUBLE PRECISION :: EPG, ROg, ROPg, EP_SM, DPM, ROs
103
104
105
106
107
108
109
110
111
112
113 DO IJK = ijkstart3, ijkend3
114
115 IF (FLUIDorP_FLOW_AT(IJK)) THEN
116
117 I = I_OF(IJK)
118 IMJK = IM_OF(IJK)
119 IJMK = JM_OF(IJK)
120 IJKM = KM_OF(IJK)
121
122 MAXM = SMAX
123 DO CM = 1,SMAX
124 DP_loc(CM) = D_p(IJK,CM)
125 EPs_loc(CM) = EP_S(IJK,CM)
126 ROs_loc(CM) = RO_S(IJK,CM)
127 ENDDO
128
129
130
131 IF(DES_CONTINUUM_HYBRID) THEN
132 MAXM = SMAX + DES_MMAX
133 DO DM = 1,DES_MMAX
134 L = SMAX + DM
135 DP_loc(L) = DES_D_p0(DM)
136 EPs_loc(L) = DES_ROP_S(IJK,DM)/DES_RO_S(DM)
137 ROs_loc(L) = DES_RO_S(DM)
138 ENDDO
139 ENDIF
140
141
142 = AVG_X_E(U_G(IMJK),U_G(IJK),I)
143 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
144 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
145
146 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
147 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
148 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
149
150
151 = SQRT((UGC-USCM)**2 + (VGC-VSCM)**2 + (WGC-WSCM)**2)
152
153
154
155
156
157 IF (P_FLOW_AT(IJK)) THEN
158 IF( FLUID_AT(EAST_OF(IJK)) ) THEN
159 Mu = MU_G(EAST_OF(IJK))
160 ELSEIF ( FLUID_AT(WEST_OF(IJK)) ) THEN
161 Mu = MU_G(WEST_OF(IJK))
162 ELSEIF ( FLUID_AT(NORTH_OF(IJK)) ) THEN
163 Mu = MU_G(NORTH_OF(IJK))
164 ELSEIF ( FLUID_AT(SOUTH_OF(IJK)) ) THEN
165 Mu = MU_G(SOUTH_OF(IJK))
166 ELSEIF ( FLUID_AT(TOP_OF(IJK)) ) THEN
167 Mu = MU_G(TOP_OF(IJK))
168 ELSEIF ( FLUID_AT(BOTTOM_OF(IJK)) ) THEN
169 Mu = MU_G(BOTTOM_OF(IJK))
170 ENDIF
171 ELSE
172 Mu = MU_G(IJK)
173 ENDIF
174
175
176 = ZERO
177 DO L = 1, MAXM
178
179 = phis + EPs_loc(L)
180 ENDDO
181
182
183 = ZERO
184 tmp_sum = ZERO
185 tmp_fac = ZERO
186 DO L = 1, MAXM
187 IF (phis .GT. ZERO) THEN
188 tmp_fac = EPs_loc(L)/phis
189 tmp_sum = tmp_sum + tmp_fac/DP_loc(L)
190 ELSE
191 tmp_sum = tmp_sum + ONE/DP_loc(L)
192 ENDIF
193 ENDDO
194 DPA = ONE / tmp_sum
195 Y_i = DP_loc(M)/DPA
196
197
198 = EP_G(IJK)
199 ROg = RO_G(IJK)
200 ROPg = ROP_G(IJK)
201 EP_SM = EPs_loc(M)
202 DPM = DP_loc(M)
203 ROs = ROs_loc(M)
204
205
206
207 IF (EP_SM <= ZERO) THEN
208 DgA = ZERO
209 ELSEIF (EPg == ZERO) THEN
210
211
212
213 = ZERO
214
215 ELSEIF (USE_MMS) THEN
216 DgA = ZERO
217 ELSE
218 SELECT CASE(DRAG_TYPE_ENUM)
219
220 CASE (SYAM_OBRIEN)
221 CALL DRAG_SYAM_OBRIEN(DgA,EPG,Mu,ROg,VREL,DPM)
222
223 CASE (GIDASPOW)
224 CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
225
226 CASE (GIDASPOW_PCF)
227 CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
228
229 CASE (GIDASPOW_BLEND)
230 CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
231
232 CASE (GIDASPOW_BLEND_PCF)
233 CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
234
235 CASE (WEN_YU)
236 CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPM)
237
238 CASE (WEN_YU_PCF)
239 CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPA)
240
241 CASE (KOCH_HILL)
242 CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPM,phis)
243
244 CASE (KOCH_HILL_PCF)
245 CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
246
247 CASE (BVK)
248 CALL DRAG_BVK(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
249
250 CASE (USER_DRAG)
251 CALL DRAG_USR(IJK, M, DgA, EPg, Mu, ROg, VREL, DPM, &
252 ROs, UGC, VGC, WGC)
253
254 CASE (HYS)
255
256 = ZERO
257 VSCM_HYS = ZERO
258 WSCM_HYS = ZERO
259 IF(phis > ZERO) THEN
260 DO L = 1, MAXM
261 USCM_HYS = USCM_HYS + EPs_loc(L)*(UGC - &
262 AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I))
263 VSCM_HYS = VSCM_HYS + EPs_loc(L)*(VGC - &
264 AVG_Y_N(V_S(IJMK,L),V_S(IJK,L)))
265 WSCM_HYS = WSCM_HYS + EPs_loc(L)*(WGC - &
266 AVG_Z_T(W_S(IJKM,L),W_S(IJK,L)))
267 ENDDO
268 USCM_HYS = USCM_HYS/phis
269 VSCM_HYS = VSCM_HYS/phis
270 WSCM_HYS = WSCM_HYS/phis
271 ENDIF
272
273 = SQRT(USCM_HYS**2 +VSCM_HYS**2 +WSCM_HYS**2)
274
275 CALL DRAG_HYS(DgA,EPg,Mu,ROPg,VREL,&
276 DP_loc(:),DPA,Y_i,EPs_loc(:),phis,M,MAXM,IJK)
277
278 CASE DEFAULT
279 CALL START_LOG
280 IF(.NOT.DMP_LOG) call open_pe_log(ier)
281 IF(DMP_LOG) WRITE (*, '(A,A)') &
282 'Unknown DRAG_TYPE: ', DRAG_TYPE
283 WRITE (UNIT_LOG, '(A,A)') 'Unknown DRAG_TYPE: ', DRAG_TYPE
284 CALL END_LOG
285 CALL mfix_exit(myPE)
286 END SELECT
287 ENDIF
288
289
290 IF (SUBGRID_TYPE_ENUM .ne. UNDEFINED_SUBGRID_TYPE) THEN
291 IF (SUBGRID_TYPE_ENUM .EQ. IGCI) THEN
292 CALL SUBGRID_DRAG_IGCI(DgA,EPg,Mu,ROg,DPM,ROs,IJK)
293 ELSEIF (SUBGRID_TYPE_ENUM .EQ. MILIOLI) THEN
294 CALL SUBGRID_DRAG_MILIOLI(DgA,EPg,Mu,ROg,VREL,&
295 DPM,ROs,IJK)
296 ENDIF
297 ENDIF
298
299
300
301
302 IF(DRAG_TYPE_ENUM == HYS) THEN
303
304 IF(Model_B)THEN
305 F_gstmp = DgA/EPg
306 ELSE
307 F_gstmp = DgA
308 ENDIF
309 ELSE
310 IF(DRAG_TYPE_ENUM == GIDASPOW_PCF .OR. &
311 DRAG_TYPE_ENUM == GIDASPOW_BLEND_PCF .OR. &
312 DRAG_TYPE_ENUM == WEN_YU_PCF .OR. &
313 DRAG_TYPE_ENUM == KOCH_HILL_PCF .OR. &
314 DRAG_TYPE_ENUM == BVK) THEN
315
316
317
318
319 IF(Model_B) THEN
320 IF (M == 1) THEN
321 F_cor = (EPg*Y_i + phis*Y_i**2)
322 ELSE
323 F_cor = (EPg*Y_i + phis*Y_i**2 + &
324 0.064d0*EPg*Y_i**3)
325 ENDIF
326 ELSE
327 F_cor = Y_i
328 ENDIF
329 DgA = ONE/(Y_i*Y_i) * DgA * F_cor
330 ENDIF
331
332
333 IF(Model_B)THEN
334 F_gstmp = DgA * EP_SM/EPg
335 ELSE
336 F_gstmp = DgA * EP_SM
337 ENDIF
338 ENDIF
339
340
341 (IJK,M) = (ONE - UR_F_gs)*F_GS(IJK, M) + &
342 UR_F_gs*F_gstmp
343
344 IF(KT_TYPE_ENUM == GHD_2007) THEN
345 IF(M==1) THEN
346 F_gs(IJK,MMAX) = F_gs(IJK,M)
347 ELSE
348 F_gs(IJK,MMAX) = F_gs(IJK,MMAX) + F_gs(IJK,M)
349 ENDIF
350 ENDIF
351
352 ELSE
353
354 (IJK,M) = ZERO
355 IF(KT_TYPE_ENUM == GHD_2007) F_gs(IJK, MMAX) = ZERO
356
357 ENDIF
358
359 ENDDO
360
361 RETURN
362 END SUBROUTINE DRAG_GS
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378 DOUBLE PRECISION FUNCTION C_DSXRE_DV(RE)
379 USE param
380 USE param1
381 IMPLICIT NONE
382 DOUBLE PRECISION, INTENT(IN) :: RE
383
384 = (0.63D0*SQRT(RE) + 4.8D0)**2
385 RETURN
386 END FUNCTION C_DSXRE_DV
387
388
389
390
391 DOUBLE PRECISION FUNCTION C_DSXRE_TL(RE)
392 USE param
393 USE param1
394 IMPLICIT NONE
395 DOUBLE PRECISION, INTENT(IN) :: RE
396
397 = 24.D0*(1.D0 + 0.173D0*RE**0.657D0) + &
398 0.413D0*RE**2.09D0/(RE**1.09D0 + 16300.D0)
399 RETURN
400 END FUNCTION C_DSXRE_TL
401
402
403
404
405
406 DOUBLE PRECISION FUNCTION C_DS_SN(RE)
407 USE param
408 USE param1
409 IMPLICIT NONE
410 DOUBLE PRECISION, INTENT(IN) :: RE
411
412 = 24.D0*(1.D0 + 0.15D0*RE**0.687D0)/(RE+SMALL_NUMBER)
413 RETURN
414 END FUNCTION C_DS_SN
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430 SUBROUTINE DRAG_SYAM_OBRIEN(lDgA,EPg,Mug,ROg,VREL,&
431 DPM)
432
433
434
435
436 USE param
437 USE param1
438 USE constant, only : drag_c1, drag_d1, c
439 IMPLICIT NONE
440
441
442
443
444 DOUBLE PRECISION, INTENT(OUT) :: ldGA
445
446 DOUBLE PRECISION, INTENT(IN) :: EPg
447
448 DOUBLE PRECISION, INTENT(IN) :: Mug
449
450 DOUBLE PRECISION, INTENT(IN) :: ROg
451
452 DOUBLE PRECISION, INTENT(IN) :: VREL
453
454 DOUBLE PRECISION, INTENT(IN) :: DPM
455
456
457
458
459
460
461
462
463 DOUBLE PRECISION, PARAMETER :: A2 = 0.005D0
464 DOUBLE PRECISION, PARAMETER :: A3 = 90.0D0
465 DOUBLE PRECISION, PARAMETER :: RE_C = 5.D0
466 DOUBLE PRECISION, PARAMETER :: EP_C = 0.92D0
467
468
469
470
471 DOUBLE PRECISION :: A, B
472
473
474 DOUBLE PRECISION :: V_rm
475
476 DOUBLE PRECISION :: RE
477
478
479
480
481 DOUBLE PRECISION, EXTERNAL :: C_DSXRE_DV
482
483
484 IF(Mug > ZERO) THEN
485 RE = DPM*VREL*ROg/Mug
486 ELSE
487 RE = LARGE_NUMBER
488 ENDIF
489
490
491 = EPg**4.14D0
492 IF (EPg <= 0.85D0) THEN
493 B = drag_c1*EPg**1.28D0
494 ELSE
495 B = EPg**drag_d1
496 ENDIF
497
498 V_RM=HALF*(A-0.06D0*RE+&
499 SQRT((3.6D-3)*RE*RE+0.12D0*RE*(2.D0*B-A)+A*A) )
500
501
502
503
504
505
506
507
508 = 0.75D0*Mug*EPg*C_DSXRE_DV(RE/V_RM)/(V_RM*DPM*DPM)
509
510 IF (RE == ZERO) lDgA = ZERO
511
512 RETURN
513 END SUBROUTINE DRAG_SYAM_OBRIEN
514
515
516
517
518
519
520
521
522
523
524
525
526 SUBROUTINE DRAG_GIDASPOW(lDgA,EPg,Mug,ROg,ROPg,VREL, DPM)
527
528
529
530
531 USE param
532 USE param1
533 IMPLICIT NONE
534
535
536
537
538 DOUBLE PRECISION, INTENT(OUT) :: lDgA
539
540 DOUBLE PRECISION, INTENT(IN) :: EPg
541
542 DOUBLE PRECISION, INTENT(IN) :: Mug
543
544 DOUBLE PRECISION, INTENT(IN) :: ROg
545
546 DOUBLE PRECISION, INTENT(IN) :: ROPg
547
548 DOUBLE PRECISION, INTENT(IN) :: VREL
549
550
551 DOUBLE PRECISION, INTENT(IN) :: DPM
552
553
554 DOUBLE PRECISION, EXTERNAL :: C_DS_SN
555
556
557
558
559
560 DOUBLE PRECISION :: RE
561
562 DOUBLE PRECISION :: C_d
563
564
565
566
567 = merge(DPM*VREL*ROPg/Mug, LARGE_NUMBER, MUg > ZERO)
568
569
570 IF(EPg <= 0.8D0) THEN
571 lDgA = 150D0*(ONE-EPg)*Mug / (EPg*DPM**2) + &
572 1.75D0*ROg*VREL/DPM
573 ELSE
574
575 IF(RE <= 1000D0)THEN
576
577 = C_DS_SN(RE)
578 ELSE
579 C_d = 0.44D0
580 ENDIF
581 lDgA = 0.75D0*C_d*VREL*ROPg*EPg**(-2.65D0) / DPM
582 ENDIF
583
584 IF (RE == ZERO) lDgA = ZERO
585
586 RETURN
587 END SUBROUTINE DRAG_GIDASPOW
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607 SUBROUTINE DRAG_GIDASPOW_BLEND(lDgA,EPg,Mug,ROg,ROPg,VREL,&
608 DPM)
609
610
611
612
613 USE param
614 USE param1
615 USE constant, only : PI
616 IMPLICIT NONE
617
618
619
620
621 DOUBLE PRECISION, INTENT(OUT) :: lDgA
622
623 DOUBLE PRECISION, INTENT(IN) :: EPg
624
625 DOUBLE PRECISION, INTENT(IN) :: Mug
626
627 DOUBLE PRECISION, INTENT(IN) :: ROg
628
629 DOUBLE PRECISION, INTENT(IN) :: ROPg
630
631 DOUBLE PRECISION, INTENT(IN) :: VREL
632
633
634 DOUBLE PRECISION, INTENT(IN) :: DPM
635
636
637
638
639 DOUBLE PRECISION :: RE
640
641 DOUBLE PRECISION :: C_d
642
643 DOUBLE PRECISION :: Ergun
644 DOUBLE PRECISION :: WenYu
645 DOUBLE PRECISION :: PHI_gs
646
647
648 IF(Mug > ZERO) THEN
649
650 = DPM*VREL*ROPg/Mug
651 ELSE
652 RE = LARGE_NUMBER
653 ENDIF
654
655
656 = 150D0*(ONE-EPg)*Mug / (EPg*DPM**2) + &
657 1.75D0*ROg*VREL/DPM
658
659 IF(RE <= 1000D0)THEN
660 C_d = (24.D0/(RE+SMALL_NUMBER)) * &
661 (ONE + 0.15D0 * RE**0.687D0)
662 ELSE
663 C_d = 0.44D0
664 ENDIF
665 WenYu = 0.75D0*C_d*VREL*ROPg*EPg**(-2.65D0) / DPM
666
667
668 = ATAN(150.D0*1.75D0*(EPg - 0.8D0))/PI + 0.5D0
669
670
671 = (1.D0-PHI_gs)*Ergun + PHI_gs*WenYu
672 IF (RE == ZERO) lDgA = ZERO
673
674 RETURN
675 END SUBROUTINE DRAG_GIDASPOW_BLEND
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690 SUBROUTINE DRAG_WEN_YU(lDgA,EPg,Mug,ROPg,VREL,DPM)
691
692
693
694
695 USE param
696 USE param1
697 IMPLICIT NONE
698
699
700
701
702 DOUBLE PRECISION, INTENT(OUT) :: lDgA
703
704 DOUBLE PRECISION, INTENT(IN) :: EPg
705
706 DOUBLE PRECISION, INTENT(IN) :: Mug
707
708 DOUBLE PRECISION, INTENT(IN) :: ROPg
709
710 DOUBLE PRECISION, INTENT(IN) :: VREL
711
712
713 DOUBLE PRECISION, INTENT(IN) :: DPM
714
715
716
717
718 DOUBLE PRECISION :: RE
719
720 DOUBLE PRECISION :: C_d
721
722
723 IF(Mug > ZERO) THEN
724
725 = DPM*VREL*ROPg/Mug
726 ELSE
727 RE = LARGE_NUMBER
728 ENDIF
729
730 IF(RE <= 1000.0D0)THEN
731 C_d = (24.D0/(RE+SMALL_NUMBER)) * (ONE + 0.15D0*RE**0.687D0)
732 ELSE
733 C_d = 0.44D0
734 ENDIF
735
736 lDgA = 0.75D0 * C_d * VREL * ROPg * EPg**(-2.65D0) / DPM
737 IF (RE == ZERO) lDgA = ZERO
738
739 RETURN
740 END SUBROUTINE DRAG_WEN_YU
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766 SUBROUTINE SUBGRID_DRAG_IGCI(lDgA,EPg,Mug,ROg,DPM,ROs,IJK)
767
768
769
770
771 USE param
772 USE param1
773 USE run, only : filter_size_ratio, SUBGRID_WALL
774 USE constant, only : GRAVITY
775 USE geometry, only : VOL,AXY,DO_K
776 IMPLICIT NONE
777
778
779
780
781 DOUBLE PRECISION, INTENT(INOUT) :: lDgA
782
783 DOUBLE PRECISION, INTENT(IN) :: EPg
784
785 DOUBLE PRECISION, INTENT(IN) :: Mug
786
787 DOUBLE PRECISION, INTENT(IN) :: ROg
788
789 DOUBLE PRECISION, INTENT(IN) :: DPM
790
791 DOUBLE PRECISION, INTENT(IN) :: ROs
792
793 INTEGER, INTENT(IN) :: IJK
794
795
796
797
798 DOUBLE PRECISION :: F_Subgrid
799
800
801 DOUBLE PRECISION :: F_SubGridWall
802
803 DOUBLE PRECISION :: vt
804
805 DOUBLE PRECISION :: filtersize
806
807 DOUBLE PRECISION :: Inv_Froude
808
809 DOUBLE PRECISION :: EPs
810
811 DOUBLE PRECISION :: GG_phip, h_phip, h_phip2, c_function,&
812 f_filter
813
814
815
816 = ONE
817 F_SubgridWall = ONE
818
819
820 = GRAVITY*DPM*DPM*(ROs - ROg) / (18.0d0*Mug)
821
822 IF(DO_K) THEN
823 filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
824 ELSE
825 filtersize = filter_size_ratio * DSQRT(AXY(IJK))
826 ENDIF
827
828
829 IF(ABS(vt) > SMALL_NUMBER) THEN
830 Inv_Froude = filtersize * GRAVITY / vt**2
831 ELSE
832 Inv_Froude = LARGE_NUMBER
833 ENDIF
834
835
836 = ONE - EPg
837
838 IF (EPs .LT. 0.0012d0) THEN
839 h_phip = 2.7d0*(EPs**0.234)
840 ELSEIF (EPs .LT. 0.014d0) THEN
841 h_phip = -0.019d0/(EPs**0.455) + 0.963d0
842 ELSEIF (EPs .LT. 0.25d0) THEN
843 h_phip = 0.868d0*EXP((-0.38*EPs)) - &
844 0.176d0*EXP((-119.2*EPs))
845 ELSEIF (EPs .LT. 0.455d0) THEN
846 h_phip = -4.59d-5*EXP((19.75*EPs)) + &
847 0.852d0*EXP((-0.268*EPs))
848 ELSEIF (EPs .LE. 0.59d0) THEN
849 h_phip = (EPs - 0.59d0) * (-1501.d0*(EPs**3) + &
850 2203.d0*(EPs**2) - 1054.d0*EPs + 162.d0)
851 ELSE
852 h_phip=ZERO
853 ENDIF
854
855 IF (EPs .LT. 0.18d0) THEN
856 GG_phip = (EPs**0.24)*(1.48d0 + EXP(-18.0*EPs))
857 ELSE
858 GG_phip = ONE
859 ENDIF
860
861
862 = (Inv_Froude**1.6) / ((Inv_Froude**1.6)+0.4d0)
863 h_phip2=h_phip*GG_phip
864 c_function=-h_phip2*f_filter
865 F_Subgrid =(ONE + c_function)
866
867 IF (SUBGRID_WALL) THEN
868 CALL SUBGRID_DRAG_WALL(F_SubgridWall,vt,IJK)
869 ENDIF
870
871 lDgA = F_SubgridWall*F_Subgrid * lDgA
872
873 RETURN
874 END SUBROUTINE SUBGRID_DRAG_IGCI
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901 SUBROUTINE SUBGRID_DRAG_MILIOLI(lDgA,EPg,Mug,ROg,VREL,DPM,ROs,&
902 IJK)
903
904
905
906
907 USE param
908 USE param1
909 USE run, only : filter_size_ratio, SUBGRID_WALL
910 USE constant, only : GRAVITY
911 USE geometry, only : VOL,AXY,DO_K
912 IMPLICIT NONE
913
914
915
916
917 DOUBLE PRECISION, INTENT(INOUT) :: lDgA
918
919 DOUBLE PRECISION, INTENT(IN) :: EPg
920
921 DOUBLE PRECISION, INTENT(IN) :: Mug
922
923 DOUBLE PRECISION, INTENT(IN) :: ROg
924
925 DOUBLE PRECISION, INTENT(IN) :: VREL
926
927 DOUBLE PRECISION, INTENT(IN) :: DPM
928
929 DOUBLE PRECISION, INTENT(IN) :: ROs
930
931 INTEGER, INTENT(IN) :: IJK
932
933
934
935
936 DOUBLE PRECISION :: F_Subgrid
937
938
939 DOUBLE PRECISION :: F_SubGridWall
940
941 DOUBLE PRECISION :: vt
942
943 DOUBLE PRECISION :: filtersize
944
945 DOUBLE PRECISION :: Inv_Froude
946
947 DOUBLE PRECISION :: vslip
948
949 DOUBLE PRECISION :: EPs
950
951 DOUBLE PRECISION :: h1, henv, hlin
952
953
954
955 = ONE
956 F_SubgridWall = ONE
957
958
959 = GRAVITY*DPM*DPM*(ROs - ROg) / (18.0d0*Mug)
960
961 IF(DO_K) THEN
962 filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
963 ELSE
964 filtersize = filter_size_ratio * DSQRT(AXY(IJK))
965 ENDIF
966
967 IF(ABS(vt) > SMALL_NUMBER) THEN
968 Inv_Froude = filtersize * GRAVITY / vt**2
969 ELSE
970 Inv_Froude = LARGE_NUMBER
971 ENDIF
972
973 = ONE - EPg
974
975 = VREL / vt
976
977 IF (Inv_Froude .LE. 1.028d0) THEN
978 h1 = (1.076d0 + 0.12d0*Vslip - (0.02d0/(Vslip+0.01d0)))*EPs + &
979 (0.084d0 + 0.09d0*Vslip - (0.01d0/(0.1d0*Vslip+0.01d0)))
980 IF (EPs .LE. 0.53d0) THEN
981 henv = (6.8d0*(ONE+EPs)*(EPs**0.3)) / &
982 (10.d0*(EPs**1.5) + 5.d0)
983 ELSEIF (EPs .GT. 0.53d0 .AND. EPs .LE. 0.65d0) THEN
984 henv = (2.23d0*((0.65d0-EPs)**(0.45))) / &
985 ((ONE/EPs)-ONE)
986 ELSEIF (EPs .GT. 0.65d0) THEN
987 henv=ZERO
988 ENDIF
989
990 ELSEIF (Inv_Froude .GT. 1.028d0 .AND. &
991 Inv_Froude .LE. 2.056d0) THEN
992 h1 = (1.268d0 - (0.2d0*Vslip) + (0.14d0/(Vslip+0.01d0)))*EPs + &
993 (0.385d0 + 0.09d0*Vslip - (0.05d0/(0.2d0*Vslip+0.01d0)))
994 IF (EPs .LE. 0.53d0) THEN
995 henv = (8.6d0*(ONE+EPs)*(EPs**0.2)) / (10.d0*EPs + 6.3d0)
996 ELSEIF (EPs .GT. 0.53d0 .AND. EPs .LE. 0.65d0) THEN
997 henv = (0.423d0*((0.65d0-EPs)**0.3)) / (ONE-(EPs**0.4))
998 ELSEIF (EPs .GT. 0.65d0) THEN
999 henv=ZERO
1000 ENDIF
1001
1002 ELSEIF (Inv_Froude .GT. 2.056d0 .AND. &
1003 Inv_Froude .LE. 4.112d0) THEN
1004 h1 = ((0.018d0*Vslip + 0.1d0)/(0.14d0*Vslip + 0.01d0))*EPs + &
1005 (0.9454d0 - (0.09d0/(0.2d0*Vslip + 0.01d0)))
1006 IF (EPs .LE. 0.5d0) THEN
1007 henv = (7.9d0*(ONE+EPs)*(EPs**0.2)) / &
1008 (10.d0*(EPs**0.9) + 5.d0)
1009 ELSEIF (EPs .GT. 0.5d0 .AND. EPs .LE. 0.63d0) THEN
1010 henv = (0.705d0*((0.63d0-EPs)**0.3)) / (ONE-(EPs**0.7))
1011 ELSEIF (EPs .GT. 0.63d0) THEN
1012 henv=ZERO
1013 ENDIF
1014
1015 ELSEIF (Inv_Froude .GT. 4.112d0 .AND. &
1016 Inv_Froude .LE. 8.224d0) THEN
1017 h1 = ((0.05d0*Vslip+0.3d0)/(0.4d0*Vslip+0.06d0))*EPs + &
1018 (0.9466d0 - (0.05d0/(0.11d0*Vslip+0.01d0)))
1019 IF (EPs .LE. 0.45d0) THEN
1020 henv = (7.9d0*(ONE+EPs)*(EPs**0.2)) / &
1021 ((10.d0*(EPs**0.6)) + 3.6d0)
1022 ELSEIF (EPs .GT. 0.45d0 .AND. EPs .LE. 0.57d0) THEN
1023 henv = (0.78d0*((0.57d0-EPs)**0.2)) / (ONE-(EPs**0.9))
1024 ELSEIF (EPs .GT. 0.57d0) THEN
1025 henv=ZERO
1026 ENDIF
1027
1028 ELSEIF (Inv_Froude .GT. 8.224d0 .AND. &
1029 Inv_Froude .LE. 12.336d0) THEN
1030 h1 = ((1.3d0*Vslip+2.2d0)/(5.2d0*Vslip+0.07d0))*EPs + &
1031 (0.9363d0-(0.11d0/(0.3d0*Vslip+0.01d0)))
1032 IF (EPs .LE. 0.35d0) THEN
1033 henv = (7.6d0*(ONE+EPs)*(EPs**0.2)) / &
1034 ((10.d0*(EPs**0.6)) + 3.3d0)
1035 ELSEIF (EPs .GT. 0.35d0 .AND. EPs .LE. 0.55d0) THEN
1036 henv = (0.81d0*((0.55d0-EPs)**0.3)) / (ONE-(EPs**0.7))
1037 ELSEIF (EPs .GT. 0.55d0) THEN
1038 henv=ZERO
1039 ENDIF
1040
1041 ELSEIF (Inv_Froude .GT. 12.336d0 .AND. &
1042 Inv_Froude .LE. 16.448d0) THEN
1043 h1 = ((2.6d0*Vslip+4.d0)/(10.d0*Vslip+0.08d0))*EPs + &
1044 (0.926d0-(0.17d0/(0.5d0*Vslip+0.01d0)))
1045 IF (EPs .LE. 0.25d0) THEN
1046 henv = (8.4d0*(ONE+EPs)*(EPs**0.2)) / &
1047 ((10.d0*(EPs**0.5)) + 3.3d0)
1048 ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1049 henv = (1.01d0*((0.52d0-EPs)**0.03))/(ONE-(EPs**0.9))
1050 ELSEIF (EPs .GT. 0.52d0) THEN
1051 henv=ZERO
1052 ENDIF
1053
1054 ELSEIF (Inv_Froude .GT. 16.448d0 .AND. &
1055 Inv_Froude .LE. 20.56d0) THEN
1056 h1 = ((2.5d0*Vslip+4.d0)/(10.d0*Vslip+0.08d0))*EPs + &
1057 (0.9261d0-(0.17d0/(0.5d0*Vslip+0.01d0)))
1058 IF (EPs .LE. 0.25d0) THEN
1059 henv = (8.4d0*(ONE+EPs)*(EPs**0.2)) / &
1060 ((10.d0*(EPs**0.5)) + 3.3d0)
1061 ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1062 henv = (1.065d0*((0.52d0-EPs)**0.3))/(ONE-EPs)
1063 ELSEIF (EPs .GT.0.52d0) THEN
1064 henv=ZERO
1065 ENDIF
1066
1067 ELSEIF (Inv_Froude .GT. 20.56d0) THEN
1068 h1 = ((1.6d0*Vslip+4.d0)/(7.9d0*Vslip+0.08d0))*EPs + &
1069 (0.9394d0 - (0.22d0/(0.6d0*Vslip+0.01d0)))
1070 IF (EPs .LE. 0.25d0) THEN
1071 henv = (9.d0*(ONE+EPs)*(EPs**0.15)) / &
1072 (10.d0*(EPs**0.45) + 4.2d0)
1073 ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1074 henv = (0.91d0*((0.52d0-EPs)**0.4))/(ONE-(EPs**0.6))
1075 ELSEIF (EPs .GT. 0.52d0) THEN
1076 henv=ZERO
1077 ENDIF
1078 ENDIF
1079
1080 IF (h1 .GT. ZERO) THEN
1081 hlin=h1
1082 ELSE
1083 hlin=ZERO
1084 ENDIF
1085
1086 IF (Inv_Froude .LT. 1.028d0) THEN
1087
1088
1089 = ONE
1090 ELSE
1091
1092 = ONE - MIN(henv,hlin)
1093 ENDIF
1094
1095
1096
1097
1098 IF (SUBGRID_WALL) THEN
1099 CALL SUBGRID_DRAG_WALL(F_SubgridWall,vt,IJK)
1100 ENDIF
1101
1102 lDgA = F_SubgridWall*F_Subgrid * lDgA
1103
1104
1105 RETURN
1106 END SUBROUTINE SUBGRID_DRAG_MILIOLI
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132 SUBROUTINE SUBGRID_DRAG_WALL(lSubgridWall,vt,IJK)
1133
1134
1135
1136
1137 USE param
1138 USE param1
1139 USE constant, only : GRAVITY
1140 USE cutcell, only : DWALL
1141 IMPLICIT NONE
1142
1143
1144
1145
1146
1147 DOUBLE PRECISION, INTENT(OUT) :: lSubGridWall
1148
1149 DOUBLE PRECISION, INTENT(IN) :: vt
1150
1151 INTEGER, INTENT(IN) :: IJK
1152
1153
1154
1155
1156 DOUBLE PRECISION, PARAMETER :: a22=6.0d0, b22=0.295d0
1157
1158
1159
1160
1161 DOUBLE PRECISION :: x_d
1162
1163
1164 = ONE
1165
1166
1167 = DWALL(IJK) * GRAVITY / vt**2
1168
1169
1170
1171 = ONE / ( ONE + a22 * (EXP(-b22*x_d)) )
1172
1173 RETURN
1174 END SUBROUTINE SUBGRID_DRAG_WALL
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198 SUBROUTINE DRAG_KOCH_HILL(lDgA,EPg,Mug,ROPg,VREL,&
1199 DPM,DPA,PHIS)
1200
1201
1202
1203
1204 USE param
1205 USE param1
1206 IMPLICIT NONE
1207
1208
1209
1210
1211 DOUBLE PRECISION, INTENT(OUT) :: lDgA
1212
1213 DOUBLE PRECISION, INTENT(IN) :: EPg
1214
1215 DOUBLE PRECISION, INTENT(IN) :: Mug
1216
1217 DOUBLE PRECISION, INTENT(IN) :: ROPg
1218
1219 DOUBLE PRECISION, INTENT(IN) :: VREL
1220
1221 DOUBLE PRECISION, INTENT(IN) :: DPM
1222
1223 DOUBLE PRECISION, INTENT(IN) :: DPA
1224
1225 DOUBLE PRECISION, INTENT(IN) :: PHIS
1226
1227
1228
1229
1230 DOUBLE PRECISION :: RE
1231
1232 DOUBLE PRECISION :: Re_Trans_1, Re_Trans_2
1233
1234 DOUBLE PRECISION :: F_STOKES
1235
1236 DOUBLE PRECISION :: F_0
1237
1238 DOUBLE PRECISION :: F_1
1239
1240 DOUBLE PRECISION :: F_2
1241
1242 DOUBLE PRECISION :: F_3
1243
1244 DOUBLE PRECISION :: F
1245
1246 DOUBLE PRECISION :: w
1247
1248
1249
1250 IF(Mug > ZERO) THEN
1251
1252 = 0.5D0*DPA*VREL*ROPG/Mug
1253 ELSE
1254 RE = LARGE_NUMBER
1255 ENDIF
1256
1257 F_STOKES = 18.D0*Mug*EPg*EPg/DPM**2
1258 = EXP(-10.0D0*(0.4D0-phis)/phis)
1259
1260 IF(phis > 0.01D0 .AND. phis < 0.4D0) THEN
1261 F_0 = (1.0D0-w) * (1.0D0 + 3.0D0*dsqrt(phis/2.0D0) + &
1262 135.0D0/64.0D0*phis*LOG(phis) + 17.14D0*phis) / &
1263 (1.0D0 + 0.681D0*phis - 8.48D0*phis*phis + &
1264 8.16D0*phis**3) + w*10.0D0*phis/(1.0D0-phis)**3
1265 ELSEIF(phis >= 0.4D0) THEN
1266 F_0 = 10.0D0*phis/(1.0D0-phis)**3
1267 ENDIF
1268
1269 IF(phis > 0.01D0 .AND. phis <= 0.1D0) THEN
1270 F_1 = dsqrt(2.0D0/phis) / 40.0D0
1271 ELSE IF(phis > 0.1D0) THEN
1272 F_1 = 0.11D0 + 5.1D-04 * exp(11.6D0*phis)
1273 ENDIF
1274
1275 IF(phis < 0.4D0) THEN
1276 F_2 = (1.0D0-w) * (1.0D0 + 3.0D0*dsqrt(phis/2.0D0) + &
1277 135.0D0/64.0D0*phis*LOG(phis) + 17.89D0*phis) / &
1278 (1.0D0 + 0.681D0*phis - 11.03D0*phis*phis + &
1279 15.41D0*phis**3)+ w*10.0D0*phis/(1.0D0-phis)**3
1280 ELSE
1281 F_2 = 10.0D0*phis/(1.0D0-phis)**3
1282 ENDIF
1283
1284 IF(phis < 0.0953D0) THEN
1285 F_3 = 0.9351D0*phis + 0.03667D0
1286 ELSE
1287 F_3 = 0.0673D0 + 0.212D0*phis +0.0232D0/(1.0-phis)**5
1288 ENDIF
1289
1290 Re_Trans_1 = (F_2 - 1.0D0)/(3.0D0/8.0D0 - F_3)
1291 Re_Trans_2 = (F_3 + dsqrt(F_3*F_3 - 4.0D0*F_1 &
1292 *(F_0-F_2))) / (2.0D0*F_1)
1293
1294 IF(phis <= 0.01D0 .AND. RE <= Re_Trans_1) THEN
1295 F = 1.0D0 + 3.0D0/8.0D0*RE
1296 ELSEIF(phis > 0.01D0 .AND. RE <= Re_Trans_2) THEN
1297 F = F_0 + F_1*RE*RE
1298 ELSEIF(phis <= 0.01D0 .AND. RE > Re_Trans_1 .OR. &
1299 phis > 0.01D0 .AND. RE > Re_Trans_2) THEN
1300 F = F_2 + F_3*RE
1301 ELSE
1302 F = zero
1303 ENDIF
1304
1305 lDgA = F * F_STOKES
1306 IF (RE == ZERO) lDgA = ZERO
1307
1308 RETURN
1309 END SUBROUTINE DRAG_KOCH_HILL
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324 SUBROUTINE DRAG_BVK(lDgA,EPg,Mug,ROPg,VREL,&
1325 DPM,DPA,PHIS)
1326
1327
1328
1329
1330 USE param
1331 USE param1
1332 IMPLICIT NONE
1333
1334
1335
1336
1337 DOUBLE PRECISION, INTENT(OUT) :: lDgA
1338
1339 DOUBLE PRECISION, INTENT(IN) :: EPg
1340
1341 DOUBLE PRECISION, INTENT(IN) :: Mug
1342
1343 DOUBLE PRECISION, INTENT(IN) :: ROPg
1344
1345 DOUBLE PRECISION, INTENT(IN) :: VREL
1346
1347 DOUBLE PRECISION, INTENT(IN) :: DPM
1348
1349 DOUBLE PRECISION, INTENT(IN) :: DPA
1350
1351 DOUBLE PRECISION, INTENT(IN) :: PHIS
1352
1353
1354
1355
1356 DOUBLE PRECISION :: RE
1357
1358 DOUBLE PRECISION :: F_STOKES
1359
1360 DOUBLE PRECISION :: F
1361
1362
1363 IF(Mug > ZERO) THEN
1364
1365 = DPA*VREL*ROPg/Mug
1366 ELSE
1367 RE = LARGE_NUMBER
1368 ENDIF
1369
1370
1371
1372 = 18D0*Mug*EPg/DPM**2
1373
1374 = 10d0*phis/EPg**2 + EPg**2*(ONE+1.5d0*DSQRT(phis))
1375 F = F + 0.413d0*RE/(24.d0*EPg**2) * &
1376 (ONE/EPg + 3d0*EPg*phis + 8.4d0/RE**0.343) / &
1377 (ONE+10.d0**(3d0*phis)/RE**(0.5+2.d0*phis))
1378
1379 lDgA = F*F_STOKES
1380 IF (RE == ZERO) lDgA = ZERO
1381
1382 RETURN
1383 END SUBROUTINE DRAG_BVK
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398 SUBROUTINE DRAG_HYS(lDgA,EPg,Mug,ROPg,VREL,&
1399 DPM,DPA,Y_I,EP_sM,PHIS,M,MAXM,IJK)
1400
1401
1402
1403
1404 USE param
1405 USE param1
1406 USE drag, only : beta_ij
1407 USE run, only : LAM_HYS
1408 IMPLICIT NONE
1409
1410
1411
1412
1413 DOUBLE PRECISION, INTENT(OUT) :: lDgA
1414
1415 DOUBLE PRECISION, INTENT(IN) :: EPg
1416
1417 DOUBLE PRECISION, INTENT(IN) :: Mug
1418
1419 DOUBLE PRECISION, INTENT(IN) :: ROPg
1420
1421 DOUBLE PRECISION, INTENT(IN) :: VREL
1422
1423 DOUBLE PRECISION :: DPM(2*DIM_M)
1424
1425 DOUBLE PRECISION, INTENT(IN) :: DPA
1426
1427 DOUBLE PRECISION, INTENT(IN) :: Y_i
1428
1429 DOUBLE PRECISION :: EP_sM(2*DIM_M)
1430
1431 DOUBLE PRECISION, INTENT(IN) :: PHIS
1432
1433 INTEGER, INTENT(IN) :: M, IJK
1434
1435 INTEGER, INTENT(IN) :: MAXM
1436
1437
1438
1439
1440 DOUBLE PRECISION :: DPmin
1441
1442 INTEGER :: L
1443
1444 DOUBLE PRECISION :: RE
1445
1446 DOUBLE PRECISION :: F_STOKES
1447
1448 DOUBLE PRECISION :: F
1449
1450 DOUBLE PRECISION :: a_YS
1451
1452 DOUBLE PRECISION :: alpha_YS
1453
1454 DOUBLE PRECISION :: beta_i_HYS
1455
1456 DOUBLE PRECISION :: beta_j_HYS
1457
1458 DOUBLE PRECISION :: FSTOKES_j
1459
1460 DOUBLE PRECISION :: Y_i_J
1461
1462 DOUBLE PRECISION :: F_D_BVK
1463
1464 DOUBLE PRECISION :: F_YS
1465
1466
1467 IF (Mug > ZERO) THEN
1468
1469 = DPA*VREL*ROPg/Mug
1470 ELSE
1471 RE = LARGE_NUMBER
1472 ENDIF
1473
1474
1475
1476 = 18D0*Mug*EPg*EP_SM(M)/DPM(M)**2
1477
1478
1479 = DPM(1)
1480 IF (MAXM > 1) THEN
1481 DO L=2,MAXM
1482 Dpmin = MIN(Dpmin,DPM(L))
1483 ENDDO
1484 ENDIF
1485
1486 a_YS = 1d0 - 2.66d0*phis + 9.096d0*phis**2 - 11.338d0*phis**3
1487
1488
1489
1490 = 1.313d0*LOG10(DPmin/lam_HYS) - 1.249d0
1491
1492
1493
1494 = ZERO
1495 F = 10.d0*phis/EPg**2 + EPg**2*(ONE+1.5d0*DSQRT(phis))
1496
1497 IF(RE > ZERO) F_D_BVK = F + 0.413d0*RE/(24.d0*EPg**2)*&
1498 (ONE/EPg + 3.d0*EPg*phis + 8.4d0/RE**0.343d0) / &
1499 (ONE+10**(3.d0*phis)/RE**(0.5d0+2.d0*phis))
1500
1501
1502 = 1d0/EPg + (F_D_BVK - 1d0/EPg)*&
1503 (a_YS*Y_i+(1d0-a_YS)*Y_i**2)
1504 F_YS = F_YS*F_STOKES
1505 beta_i_HYS = F_YS
1506
1507 DO L= 1,MAXM
1508 IF (L /= M) THEN
1509 Y_i_J = DPM(L)/DPA
1510 beta_j_HYS = 1.d0/EPg + (F_D_BVK - 1.d0/EPg) * &
1511 (a_YS*Y_i_J + (1d0-a_YS)*Y_i_J**2)
1512 FSTOKES_j = 18.D0*Mug*EP_sM(L)*EPg/&
1513 DPM(L)**2
1514
1515 beta_j_HYS = beta_j_HYS*FSTOKES_j
1516
1517
1518 (IJK,M,L) = ZERO
1519
1520
1521 IF (EP_sM(M) > ZERO .AND. EP_SM(L) > ZERO) &
1522 beta_ij(IJK,M,L) = (2.d0*alpha_YS*EP_sM(M)*EP_sM(L))/ &
1523 (EP_sM(M)/beta_i_HYS + EP_sM(L)/beta_j_HYS)
1524 F_YS = F_YS + beta_ij(IJK,M,L)
1525
1526 ENDIF
1527 ENDDO
1528
1529
1530 = F_YS
1531
1532 RETURN
1533 END SUBROUTINE DRAG_HYS
1534