File: N:\mfix\model\drag_gs.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DRAG_gs                                                 C
4     !  Purpose: Calculate the gas-solids drag coefficient                  C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 20-MAY-96  C
7     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !  Literature/Document References:                                     C
10     !                                                                      C
11     !  Comments:                                                           C
12     !    If changes are made to this file, especially in regard to the     C
13     !    structure of any drag subroutine call, then such changes must     C
14     !    also be made to the analgous routine des_drag_gp in the file      C
15     !    des/drag_fgs.f for consistency!                                   C
16     !                                                                      C
17     !                                                                      C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20     
21           SUBROUTINE DRAG_GS(M, IER)
22     
23     !-----------------------------------------------
24     ! Modules
25     !-----------------------------------------------
26           USE compar
27           USE constant
28           USE discretelement
29           USE drag
30           USE fldvar
31           USE fun_avg
32           USE functions
33           USE funits
34           USE geometry
35           USE indices
36           USE machine, only: start_log, end_log
37           USE mms
38           USE parallel
39           USE param
40           USE param1
41           USE physprop
42           USE run, only: syam_obrien, gidaspow, gidaspow_pcf, gidaspow_blend, gidaspow_blend_pcf, wen_yu, wen_yu_pcf, koch_hill, koch_hill_pcf, bvk, user_drag, hys, drag_type, drag_type_enum, ghd_2007, igci, kt_type_enum, milioli, model_b, subgrid_type_enum, undefined_subgrid_type
43           USE sendrecv
44           USE ur_facs
45           IMPLICIT NONE
46     !-----------------------------------------------
47     ! Dummy arguments
48     !-----------------------------------------------
49     ! Solids phase index
50           INTEGER, INTENT(IN) :: M
51     ! Error index
52           INTEGER, INTENT(INOUT) :: IER
53     !-----------------------------------------------
54     ! Local parameters
55     !-----------------------------------------------
56     
57     !-----------------------------------------------
58     ! Local variables
59     !-----------------------------------------------
60     ! indices
61           INTEGER :: I, IJK, IMJK, IJMK, IJKM
62     ! cell center value of U_sm, V_sm, W_sm
63           DOUBLE PRECISION :: USCM, VSCM, WSCM
64     ! cell center value of x, y and z-particle velocity
65           DOUBLE PRECISION :: USCM_HYS, VSCM_HYS, WSCM_HYS
66     ! cell center value of U_g, V_g, W_g
67           DOUBLE PRECISION :: UGC, VGC, WGC
68     ! magnitude of gas-solids relative velocity
69           DOUBLE PRECISION :: VREL
70     ! gas laminar viscosity redefined here to set viscosity at pressure
71     ! boundaries
72           DOUBLE PRECISION :: Mu
73     ! drag coefficient
74           DOUBLE PRECISION :: DgA
75     ! current value of F_gs (i.e., without underrelaxation)
76           DOUBLE PRECISION F_gstmp
77     ! indices of solids phases (continuous, discrete)
78           INTEGER :: CM, DM, L
79     ! temporary shift of total number of solids phases to account for both
80     ! discrete and continuous solids phases used for the hybrid mdoel
81           INTEGER :: MAXM
82     ! tmp local variable for the particle diameter of solids
83     ! phase M (continuous or discrete)
84           DOUBLE PRECISION :: DP_loc(DIM_M)
85     ! tmp local variable for the solids volume fraction of solids
86     ! phase M (continuous or discrete)
87           DOUBLE PRECISION :: EPs_loc(DIM_M)
88     ! tmp local variable for the particle density of solids
89     ! phase M (continuous or discrete)
90           DOUBLE PRECISION :: ROs_loc(DIM_M)
91     ! correction factors for implementing polydisperse drag model
92     ! proposed by van der Hoef et al. (2005)
93           DOUBLE PRECISION :: F_cor, tmp_sum, tmp_fac
94     ! average particle diameter in polydisperse systems
95           DOUBLE PRECISION :: DPA
96     ! diameter ratio in polydisperse systems
97           DOUBLE PRECISION :: Y_i
98     ! total solids volume fraction
99           DOUBLE PRECISION :: phis
100     ! temporary local variables to use for dummy arguments in subroutines
101     ! void fraction, gas density, gas bulk density, solids volume fraction
102     ! particle diameter, particle density
103           DOUBLE PRECISION :: EPG, ROg, ROPg, EP_SM, DPM, ROs
104     !-----------------------------------------------
105     
106     !!$omp  parallel do default(shared)                                   &
107     !!$omp  private( I,  IJK, IMJK, IJMK, IJKM, DM, MAXM, CM, L,          &
108     !!$omp           UGC, VGC, WGC, USCM, VSCM, WSCM, VREL, USCM_HYS,     &
109     !!$omp           VSCM_HYS, WSCM_HYS, tmp_sum, tmp_fac, Y_i, F_cor,    &
110     !!$omp           EP_SM, EPs_loc, ROs_loc, DP_loc, DPA, DPM, ROs,      &
111     !!$omp           phis, EPg, ROg, ROPg, Mu, DgA, F_gstmp)
112     
113     
114           DO IJK = ijkstart3, ijkend3
115     
116              IF (FLUIDorP_FLOW_AT(IJK)) THEN
117     
118                 I = I_OF(IJK)
119                 IMJK = IM_OF(IJK)
120                 IJMK = JM_OF(IJK)
121                 IJKM = KM_OF(IJK)
122     
123                 DO L = 1,DES_MMAX+MMAX
124                    IF(KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX) THEN
125     ! set this to avoid issues later in routine
126                       DP_loc(L) = one
127                       EPS_loc(L) = zero
128                       ROs_loc(L) = zero
129                    ENDIF
130                    DP_loc(L) = D_p(IJK,L)
131                    EPs_loc(L) = EP_S(IJK,L)
132                    ROs_loc(L) = RO_S(IJK,L)
133                 ENDDO
134     
135     ! Calculate velocity components at i, j, k
136                 UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
137                 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
138                 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
139     
140                 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
141                 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
142                 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
143     
144     ! magnitude of gas-solids relative velocity
145                 VREL = SQRT((UGC-USCM)**2 + (VGC-VSCM)**2 + (WGC-WSCM)**2)
146     
147     ! Laminar viscosity at a pressure boundary is given the value of the
148     ! fluid cell next to it. This applies just to the calculation of the
149     ! drag, in other routines the value of viscosity at a pressure boundary
150     ! always has a zero value.
151                 IF (P_FLOW_AT(IJK)) THEN
152                    IF( FLUID_AT(EAST_OF(IJK)) ) THEN
153                       Mu = MU_G(EAST_OF(IJK))
154                    ELSEIF ( FLUID_AT(WEST_OF(IJK)) ) THEN
155                       Mu = MU_G(WEST_OF(IJK))
156                    ELSEIF ( FLUID_AT(NORTH_OF(IJK)) ) THEN
157                       Mu = MU_G(NORTH_OF(IJK))
158                    ELSEIF ( FLUID_AT(SOUTH_OF(IJK)) ) THEN
159                       Mu = MU_G(SOUTH_OF(IJK))
160                    ELSEIF ( FLUID_AT(TOP_OF(IJK)) ) THEN
161                       Mu = MU_G(TOP_OF(IJK))
162                    ELSEIF ( FLUID_AT(BOTTOM_OF(IJK)) ) THEN
163                       Mu = MU_G(BOTTOM_OF(IJK))
164                    ENDIF
165                 ELSE
166                    Mu = MU_G(IJK)
167                 ENDIF
168     
169     ! calculate the total solids volume fraction
170                 phis = ZERO
171                 DO L = 1, DES_MMAX+MMAX
172     ! this is slightly /= one-ep_g due to round-off
173                    phis = phis + EPs_loc(L)
174                 ENDDO
175     
176     ! calculate the average particle diameter and particle ratio
177                 DPA = ZERO
178                 tmp_sum = ZERO
179                 tmp_fac = ZERO
180                 DO L = 1, DES_MMAX+MMAX
181                    IF (phis .GT. ZERO) THEN
182                       tmp_fac = EPs_loc(L)/phis
183                       tmp_sum = tmp_sum + tmp_fac/DP_loc(L)
184                     ELSE
185                       tmp_sum = tmp_sum + ONE/DP_loc(L) ! not important, but will avoid NaN's in empty cells
186                     ENDIF
187                 ENDDO
188                 DPA = ONE / tmp_sum
189                 Y_i = DP_loc(M)/DPA
190     
191     ! assign aliases for easy reference
192                 EPg = EP_G(IJK)
193                 ROg = RO_G(IJK)
194                 ROPg = ROP_G(IJK)
195                 EP_SM = EPs_loc(M)
196                 DPM = DP_loc(M)
197                 ROs = ROs_loc(M)
198     
199     
200     ! determine the drag coefficient
201                 IF (EP_SM <= ZERO) THEN
202                    DgA = ZERO
203                 ELSEIF (EPg == ZERO) THEN
204     ! this case will already be caught in most drag subroutines whenever
205     ! RE==0 (for correlations in which RE includes EPg). however, this will
206     ! prevent potential divisions by zero in some models by setting it now.
207                    DgA = ZERO
208     ! Force a ZERO drag coefficient for MMS cases.
209                 ELSEIF (USE_MMS) THEN
210                    DgA = ZERO
211                 ELSE
212                    SELECT CASE(DRAG_TYPE_ENUM)
213     
214                    CASE (SYAM_OBRIEN)
215                       CALL DRAG_SYAM_OBRIEN(DgA,EPG,Mu,ROg,VREL,DPM)
216     
217                    CASE (GIDASPOW)
218                       CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
219     
220                    CASE (GIDASPOW_PCF)
221                       CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
222     
223                    CASE (GIDASPOW_BLEND)
224                       CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
225     
226                    CASE (GIDASPOW_BLEND_PCF)
227                       CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
228     
229                    CASE (WEN_YU)
230                       CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPM)
231     
232                    CASE (WEN_YU_PCF)
233                       CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPA)
234     
235                    CASE (KOCH_HILL)
236                       CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPM,phis)
237     
238                    CASE (KOCH_HILL_PCF)
239                       CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
240     
241                    CASE (BVK)
242                       CALL DRAG_BVK(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
243     
244                    CASE (USER_DRAG)
245                       CALL DRAG_USR(IJK, M, DgA, EPg, Mu, ROg, VREL, DPM,  &
246                          ROs, UGC, VGC, WGC)
247     
248                    CASE (HYS)
249     ! only over the continuous two fluid phases
250                       MAXM = SMAX
251     ! calculate velocity components of each solids phase
252                       USCM_HYS = ZERO
253                       VSCM_HYS = ZERO
254                       WSCM_HYS = ZERO
255                       IF(phis > ZERO) THEN
256                          DO L = 1, MAXM
257                             USCM_HYS = USCM_HYS + EPs_loc(L)*(UGC - &
258                                AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I))
259                             VSCM_HYS = VSCM_HYS + EPs_loc(L)*(VGC - &
260                                AVG_Y_N(V_S(IJMK,L),V_S(IJK,L)))
261                             WSCM_HYS = WSCM_HYS + EPs_loc(L)*(WGC - &
262                                AVG_Z_T(W_S(IJKM,L),W_S(IJK,L)))
263                          ENDDO
264                          USCM_HYS = USCM_HYS/phis
265                          VSCM_HYS = VSCM_HYS/phis
266                          WSCM_HYS = WSCM_HYS/phis
267                       ENDIF
268     ! magnitude of gas-solids relative velocity
269                       VREL = SQRT(USCM_HYS**2 +VSCM_HYS**2 +WSCM_HYS**2)
270     
271                       CALL DRAG_HYS(DgA,EPg,Mu,ROPg,VREL,&
272                            DP_loc(:),DPA,Y_i,EPs_loc(:),phis,M,MAXM,IJK)
273     
274                    CASE DEFAULT
275                       CALL START_LOG
276                       IF(.NOT.DMP_LOG) call open_pe_log(ier)
277                       IF(DMP_LOG) WRITE (*, '(A,A)') &
278                          'Unknown DRAG_TYPE: ', DRAG_TYPE
279                       WRITE (UNIT_LOG, '(A,A)') 'Unknown DRAG_TYPE: ', DRAG_TYPE
280                       CALL END_LOG
281                       CALL mfix_exit(myPE)
282                    END SELECT   ! end selection of drag_type
283                 ENDIF   ! end if/elseif/else (ep_sm <= zero, ep_g==0)
284     
285     ! modify the drag coefficient to account for subgrid domain effects
286                 IF (SUBGRID_TYPE_ENUM .ne. UNDEFINED_SUBGRID_TYPE) THEN
287                    IF (SUBGRID_TYPE_ENUM .EQ. IGCI) THEN
288                       CALL SUBGRID_DRAG_IGCI(DgA,EPg,Mu,ROg,DPM,ROs,IJK)
289                    ELSEIF (SUBGRID_TYPE_ENUM .EQ. MILIOLI) THEN
290                       CALL SUBGRID_DRAG_MILIOLI(DgA,EPg,Mu,ROg,VREL,&
291                            DPM,ROs,IJK)
292                    ENDIF
293                 ENDIF
294     
295     
296     ! Modify drag coefficient to account for possible corrections and
297     ! for differences between Model B and Model A
298                 IF(DRAG_TYPE_ENUM == HYS) THEN
299     ! this drag model is handled differently than the others
300                    IF(Model_B)THEN
301                       F_gstmp = DgA/EPg
302                    ELSE
303                       F_gstmp = DgA
304                    ENDIF
305                 ELSE
306                    IF(DRAG_TYPE_ENUM == GIDASPOW_PCF .OR. &
307                       DRAG_TYPE_ENUM == GIDASPOW_BLEND_PCF .OR. &
308                       DRAG_TYPE_ENUM == WEN_YU_PCF .OR. &
309                       DRAG_TYPE_ENUM == KOCH_HILL_PCF .OR. &
310                       DRAG_TYPE_ENUM == BVK) THEN
311     ! see erratum by Beetstra et al. (2007) : the correction factor differs
312     ! for model A versus model B.
313     ! application of the correction factor for model A is found from
314     ! the correction factor for model B and neglects the Y_i**3 term
315                       IF(Model_B) THEN
316                          IF (M == 1) THEN
317                             F_cor = (EPg*Y_i + phis*Y_i**2)
318                          ELSE
319                             F_cor = (EPg*Y_i + phis*Y_i**2 + &
320                                0.064d0*EPg*Y_i**3)
321                          ENDIF
322                       ELSE
323                          F_cor = Y_i
324                       ENDIF
325                       DgA = ONE/(Y_i*Y_i) * DgA * F_cor
326                    ENDIF
327     
328     ! Calculate the drag coefficient (Model B coeff = Model A coeff/EP_g); all other models, eg., Wen_Yu
329                    IF(Model_B)THEN
330                       F_gstmp = DgA * EP_SM/EPg
331                    ELSE
332                       F_gstmp = DgA * EP_SM                     !3D buoyancy model
333                    ENDIF
334                 ENDIF   !end if/else trim(drag_type=='hys')
335     
336     ! Determine drag force coefficient accounting for any under relaxation
337                 F_GS(IJK,M) = (ONE - UR_F_gs)*F_GS(IJK, M) + &
338                    UR_F_gs*F_gstmp
339     
340                 IF(KT_TYPE_ENUM == GHD_2007) THEN
341                    IF(M==1) THEN
342                       F_gs(IJK,MMAX) = F_gs(IJK,M)
343                    ELSE
344                       F_gs(IJK,MMAX) = F_gs(IJK,MMAX) + F_gs(IJK,M)
345                    ENDIF
346                 ENDIF
347     
348              ELSE   ! .not.(fluidorp_flow_at(ijk)) branch
349     
350                 F_GS(IJK,M) = ZERO
351                 IF(KT_TYPE_ENUM == GHD_2007) F_gs(IJK, MMAX) = ZERO
352     
353              ENDIF   ! end if (fluidorp_flow_at(ijk))
354     
355           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
356     
357           RETURN
358           END SUBROUTINE DRAG_GS
359     
360     !-----------------------------------------------------------------<<<
361     
362     ! Turton and Levenspiel (1986)
363     !----------------------------------------------------------------->>>
364           DOUBLE PRECISION FUNCTION C_DSXRE_TL(RE)
365           USE param
366           USE param1
367           IMPLICIT NONE
368           DOUBLE PRECISION, INTENT(IN) :: RE ! Reynolds number
369     
370           C_DSXRE_TL = 24.D0*(1.D0 + 0.173D0*RE**0.657D0) + &
371              0.413D0*RE**2.09D0/(RE**1.09D0 + 16300.D0)
372           RETURN
373           END FUNCTION C_DSXRE_TL
374     !-----------------------------------------------------------------<<<
375     
376     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
377     !                                                                      C
378     !  Subroutine: DRAG_SYAM_OBRIEN                                        C
379     !  Purpose: Calculate the gas-solids drag coefficient                  C
380     !                                                                      C
381     !  Literature/Document References:                                     C
382     !     Syamlal M, O'Brien TJ (1988). International Journal of           C
383     !        Multiphase Flow 14: 473-481.                                  C
384     !                                                                      C
385     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
386     
387           SUBROUTINE DRAG_SYAM_OBRIEN(lDgA,EPg,Mug,ROg,VREL,&
388                      DPM)
389     
390     !-----------------------------------------------
391     ! Modules
392     !-----------------------------------------------
393           USE constant, only : drag_c1, drag_d1
394           USE drag
395           USE param
396           USE param1
397           IMPLICIT NONE
398     !-----------------------------------------------
399     ! Dummy arguments
400     !-----------------------------------------------
401     ! drag coefficient
402           DOUBLE PRECISION, INTENT(OUT) :: ldGA
403     ! gas volume fraction
404           DOUBLE PRECISION, INTENT(IN) :: EPg
405     ! gas laminar viscosity
406           DOUBLE PRECISION, INTENT(IN) :: Mug
407     ! gas density
408           DOUBLE PRECISION, INTENT(IN) :: ROg
409     ! Magnitude of gas-solids relative velocity
410           DOUBLE PRECISION, INTENT(IN) :: VREL
411     ! particle diameter of solids phase M
412           DOUBLE PRECISION, INTENT(IN) :: DPM
413     !-----------------------------------------------
414     ! Local parameters
415     !-----------------------------------------------
416     !     Parameters in the Cluster-effect model
417     !     PARAMETER (a1 = 250.)   !for G_s = 98 kg/m^2.s
418     !     PARAMETER (a1 = 1500.)  !for G_s = 147 kg/m^2.s
419     !     a1 depends upon solids flux.  It has been represented by C(1)
420     !     defined in the data file.
421     !     DOUBLE PRECISION, PARAMETER :: A2 = 0.005D0
422     !     DOUBLE PRECISION, PARAMETER :: A3 = 90.0D0
423     !     DOUBLE PRECISION, PARAMETER :: RE_C = 5.D0
424     !     DOUBLE PRECISION, PARAMETER :: EP_C = 0.92D0
425     !-----------------------------------------------
426     ! Local variables
427     !-----------------------------------------------
428     ! Variables which are function of EP_g
429           DOUBLE PRECISION :: A, BB
430     ! Ratio of settling velocity of a multiparticle system to
431     ! that of a single particle
432           DOUBLE PRECISION :: V_rm
433     ! Reynolds number
434           DOUBLE PRECISION :: RE
435     !-----------------------------------------------
436     
437           IF(Mug > ZERO) THEN
438              RE = DPM*VREL*ROg/Mug
439           ELSE
440              RE = LARGE_NUMBER
441           ENDIF
442     
443     ! Calculate V_rm
444           A = EPg**4.14D0
445           IF (EPg <= 0.85D0) THEN
446              BB = drag_c1*EPg**1.28D0
447           ELSE
448             BB = EPg**drag_d1
449           ENDIF
450     
451           V_RM=HALF*(A-0.06D0*RE+&
452                SQRT((3.6D-3)*RE*RE+0.12D0*RE*(2.D0*BB-A)+A*A) )
453     
454     !------------------Begin cluster correction --------------------------
455     ! uncomment the following four lines ...
456     !       V_RM = V_RM * (ONE + C(1)*&
457     !                      EXP(-A2*(RE-RE_C)**2 - A3*(EPg-EP_C)**2)* &
458     !                      RE*(1. - EPg))
459     !------------------End cluster correction ----------------------------
460     
461           lDgA = 0.75D0*Mug*EPg*C_DSXRE_DV(RE/V_RM)/(V_RM*DPM*DPM)
462     
463           IF (RE == ZERO) lDgA = ZERO
464     
465           RETURN
466     
467           END SUBROUTINE DRAG_SYAM_OBRIEN
468     
469     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
470     !                                                                      C
471     !  Subroutine: DRAG_GIDASPOW                                           C
472     !  Purpose: Calculate the gas-solids drag coefficient                  C
473     !                                                                      C
474     !  Literature/Document References:                                     C
475     !     Ding J, Gidaspow D (1990). AIChE Journal 36: 523-538.            C
476     !                                                                      C
477     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
478     
479           SUBROUTINE DRAG_GIDASPOW(lDgA,EPg,Mug,ROg,ROPg,VREL, DPM)
480     
481     !-----------------------------------------------
482     ! Modules
483     !-----------------------------------------------
484           USE drag
485           USE param
486           USE param1
487           IMPLICIT NONE
488     !-----------------------------------------------
489     ! Dummy arguments
490     !-----------------------------------------------
491     ! drag coefficient
492           DOUBLE PRECISION, INTENT(OUT) :: lDgA
493     ! gas volume fraction
494           DOUBLE PRECISION, INTENT(IN) :: EPg
495     ! gas laminar viscosity
496           DOUBLE PRECISION, INTENT(IN) :: Mug
497     ! gas density
498           DOUBLE PRECISION, INTENT(IN) :: ROg
499     ! gas density*EP_g
500           DOUBLE PRECISION, INTENT(IN) :: ROPg
501     ! Magnitude of gas-solids relative velocity
502           DOUBLE PRECISION, INTENT(IN) :: VREL
503     ! particle diameter of solids phase M or
504     ! average particle diameter if PCF
505           DOUBLE PRECISION, INTENT(IN) :: DPM
506     
507     !-----------------------------------------------
508     ! Local variables
509     !-----------------------------------------------
510     ! Reynolds number
511           DOUBLE PRECISION :: RE
512     ! Single sphere drag coefficient
513           DOUBLE PRECISION :: C_d
514     !-----------------------------------------------
515     
516     ! Note the presence of gas volume fraction in ROPG
517           RE = merge(DPM*VREL*ROPg/Mug, LARGE_NUMBER, MUg > ZERO)
518     
519     ! Dense phase
520           IF(EPg <= 0.8D0) THEN
521              lDgA = 150D0*(ONE-EPg)*Mug / (EPg*DPM**2) + &
522                     1.75D0*ROg*VREL/DPM
523           ELSE
524     ! Dilute phase - EP_g >= 0.8
525              IF(RE <= 1000D0)THEN
526     ! this could be replaced with the function C_DS_SN
527                 C_d = C_DS_SN(RE)
528              ELSE
529                 C_d = 0.44D0
530              ENDIF
531              lDgA = 0.75D0*C_d*VREL*ROPg*EPg**(-2.65D0) / DPM
532           ENDIF
533     
534           IF (RE == ZERO) lDgA = ZERO
535     
536           RETURN
537     
538           END SUBROUTINE DRAG_GIDASPOW
539     
540     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
541     !                                                                      C
542     !  Subroutine: DRAG_GIDASPOW_BLEND                                     C
543     !  Purpose: Calculate the gas-solids drag coefficient                  C
544     !                                                                      C
545     !  Author: Charles E.A. Finney                        Date: 23-Mar-06  C
546     !  Reviewer: Sreekanth Pannala                        Date: 24-Mar-06  C
547     !                                                                      C
548     !  Literature/Document References:                                     C
549     !     original source unknown:                                         C
550     !     Lathouwers D, Bellan J (2000). Proceedings of the 2000 U.S. DOE  C
551     !        Hydrogen Program Review NREL/CP-570-28890. Available from     C
552     !     http://www.eere.energy.gov/hydrogenandfuelcells/pdfs/28890k.pdf. C
553     !                                                                      C
554     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
555     
556           SUBROUTINE DRAG_GIDASPOW_BLEND(lDgA,EPg,Mug,ROg,ROPg,VREL,&
557                      DPM)
558     
559     !-----------------------------------------------
560     ! Modules
561     !-----------------------------------------------
562           USE param
563           USE param1
564           USE constant, only : PI
565           IMPLICIT NONE
566     !-----------------------------------------------
567     ! Dummy arguments
568     !-----------------------------------------------
569     ! drag coefficient
570           DOUBLE PRECISION, INTENT(OUT) :: lDgA
571     ! gas volume fraction
572           DOUBLE PRECISION, INTENT(IN) :: EPg
573     ! gas laminar viscosity
574           DOUBLE PRECISION, INTENT(IN) :: Mug
575     ! gas density
576           DOUBLE PRECISION, INTENT(IN) :: ROg
577     ! gas density*EP_g
578           DOUBLE PRECISION, INTENT(IN) :: ROPg
579     ! Magnitude of gas-solids relative velocity
580           DOUBLE PRECISION, INTENT(IN) :: VREL
581     ! particle diameter of solids phase M or
582     ! average particle diameter if PCF
583           DOUBLE PRECISION, INTENT(IN) :: DPM
584     !-----------------------------------------------
585     ! Local variables
586     !-----------------------------------------------
587     ! Reynolds number
588           DOUBLE PRECISION :: RE
589     ! Single sphere drag coefficient
590           DOUBLE PRECISION :: C_d
591     ! Gidaspow switch function variables
592           DOUBLE PRECISION :: Ergun
593           DOUBLE PRECISION :: WenYu
594           DOUBLE PRECISION :: PHI_gs
595     !-----------------------------------------------
596     
597           IF(Mug > ZERO) THEN
598     ! Note the presence of gas volume fraction in ROPG
599              RE = DPM*VREL*ROPg/Mug
600           ELSE
601              RE = LARGE_NUMBER
602           ENDIF
603     
604     ! Dense phase - EP_g <= 0.8
605           Ergun = 150D0*(ONE-EPg)*Mug / (EPg*DPM**2) + &
606                    1.75D0*ROg*VREL/DPM
607     ! Dilute phase - EP_g >= 0.8
608           IF(RE <= 1000D0)THEN
609              C_d = (24.D0/(RE+SMALL_NUMBER)) * &
610                (ONE + 0.15D0 * RE**0.687D0)
611           ELSE
612              C_d = 0.44D0
613           ENDIF
614           WenYu = 0.75D0*C_d*VREL*ROPg*EPg**(-2.65D0) / DPM
615     
616     ! Switch function
617           PHI_gs = ATAN(150.D0*1.75D0*(EPg - 0.8D0))/PI + 0.5D0
618     
619     ! Blend the models
620           lDgA = (1.D0-PHI_gs)*Ergun + PHI_gs*WenYu
621           IF (RE == ZERO) lDgA = ZERO
622     
623           RETURN
624           END SUBROUTINE DRAG_GIDASPOW_BLEND
625     
626     
627     
628     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
629     !                                                                      C
630     !  Subroutine: DRAG_WEN_YU                                             C
631     !  Purpose: Calculate the gas-solids drag coefficient                  C
632     !                                                                      C
633     !  Literature/Document References:                                     C
634     !     Wen CY, Yu YH (1966). Chemical Engineering Progress Symposium    C
635     !        Series 62: 100-111.                                           C
636     !                                                                      C
637     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
638     
639           SUBROUTINE DRAG_WEN_YU(lDgA,EPg,Mug,ROPg,VREL,DPM)
640     
641     !-----------------------------------------------
642     ! Modules
643     !-----------------------------------------------
644           USE param
645           USE param1
646           IMPLICIT NONE
647     !-----------------------------------------------
648     ! Dummy arguments
649     !-----------------------------------------------
650     ! drag coefficient
651           DOUBLE PRECISION, INTENT(OUT) :: lDgA
652     ! gas volume fraction
653           DOUBLE PRECISION, INTENT(IN) :: EPg
654     ! gas laminar viscosity
655           DOUBLE PRECISION, INTENT(IN) :: Mug
656     ! gas density*EP_g
657           DOUBLE PRECISION, INTENT(IN) :: ROPg
658     ! Magnitude of gas-solids relative velocity
659           DOUBLE PRECISION, INTENT(IN) :: VREL
660     ! particle diameter of solids phase M or
661     ! average particle diameter if PCF
662           DOUBLE PRECISION, INTENT(IN) :: DPM
663     !-----------------------------------------------
664     ! Local variables
665     !-----------------------------------------------
666     ! Reynolds number
667           DOUBLE PRECISION :: RE
668     ! Single sphere drag coefficient
669           DOUBLE PRECISION :: C_d
670     !-----------------------------------------------
671     
672           IF(Mug > ZERO) THEN
673     ! Note the presence of gas volume fraction in ROPG
674              RE = DPM*VREL*ROPg/Mug
675           ELSE
676              RE = LARGE_NUMBER
677           ENDIF
678     
679           IF(RE <= 1000.0D0)THEN
680              C_d = (24.D0/(RE+SMALL_NUMBER)) * (ONE + 0.15D0*RE**0.687D0)
681           ELSE
682              C_d = 0.44D0
683           ENDIF
684     
685           lDgA = 0.75D0 * C_d * VREL * ROPg * EPg**(-2.65D0) / DPM
686           IF (RE == ZERO) lDgA = ZERO
687     
688           RETURN
689           END SUBROUTINE DRAG_WEN_YU
690     
691     
692     
693     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
694     !                                                                      C
695     !  Subroutine: SUBGRID_DRAG_IGCI                                       C
696     !  Purpose: Calculate subgrid correction to the gas-solids drag        C
697     !           coefficient developed by Wen-Yu                            C
698     !                                                                      C
699     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
700     !                                                                      C
701     !  Revision: 1                                                         C
702     !  Purpose: Minor changes, e.g., fix inconsistenty with analogous      C
703     !     calls in des_drag_gp and with new variable density feature       C
704     !  Author: Janine Carney, June 2013                                    C
705     !                                                                      C
706     !  Literature/Document References:                                     C
707     !     Igci, Y., Pannala, S., Benyahia, S., & Sundaresan S.,            C
708     !        Validation studies on filtered model equations for gas-       C
709     !        particle flows in risers, Industrial & Engineering Chemistry  C
710     !        Research, 2012, 51(4), 2094-2103                              C
711     !                                                                      C
712     !                                                                      C
713     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
714     
715           SUBROUTINE SUBGRID_DRAG_IGCI(lDgA,EPg,Mug,ROg,DPM,ROs,IJK)
716     
717     !-----------------------------------------------
718     ! Modules
719     !-----------------------------------------------
720           USE param
721           USE param1
722           USE run, only : filter_size_ratio, SUBGRID_WALL
723           USE constant, only : GRAVITY
724           USE geometry, only : VOL,AXY,DO_K
725           IMPLICIT NONE
726     !-----------------------------------------------
727     ! Dummy arguments
728     !-----------------------------------------------
729     ! drag coefficient
730           DOUBLE PRECISION, INTENT(INOUT) :: lDgA
731     ! gas volume fraction
732           DOUBLE PRECISION, INTENT(IN) :: EPg
733     ! gas laminar viscosity
734           DOUBLE PRECISION, INTENT(IN) :: Mug
735     ! gas density
736           DOUBLE PRECISION, INTENT(IN) :: ROg
737     ! particle diameter of solids phase M
738           DOUBLE PRECISION, INTENT(IN) :: DPM
739     ! particle density of solids phase M
740           DOUBLE PRECISION, INTENT(IN) :: ROs
741     ! current cell index
742           INTEGER, INTENT(IN) :: IJK
743     !-----------------------------------------------
744     ! Local variables
745     !-----------------------------------------------
746     ! factor to correct the drag for subgrid domain effects
747           DOUBLE PRECISION :: F_Subgrid
748     ! factor to correct the drag for subgrid domain effects arising from
749     ! wall
750           DOUBLE PRECISION :: F_SubGridWall
751     ! particle terminal settling velocity from stokes' formulation
752           DOUBLE PRECISION :: vt
753     ! filter size which is a function of each grid cell volume
754           DOUBLE PRECISION :: filtersize
755     ! inverse Froude number, or dimensionless filter size
756           DOUBLE PRECISION :: Inv_Froude
757     ! total solids volume fraction
758           DOUBLE PRECISION :: EPs
759     ! Variables for Igci model
760           DOUBLE PRECISION :: GG_phip, h_phip, h_phip2, c_function,&
761                               f_filter
762     !-----------------------------------------------
763     
764     ! initialize
765           F_Subgrid = ONE
766           F_SubgridWall = ONE
767     
768     ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
769           vt = GRAVITY*DPM*DPM*(ROs - ROg) / (18.0d0*Mug)
770     ! filter size calculation for each specific gridcell volume
771           IF(DO_K) THEN
772              filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
773           ELSE
774              filtersize = filter_size_ratio * DSQRT(AXY(IJK))
775           ENDIF
776     
777     ! dimensionless inverse of Froude number
778           IF(ABS(vt) > SMALL_NUMBER) THEN
779              Inv_Froude =  filtersize * GRAVITY / vt**2
780           ELSE
781              Inv_Froude =  LARGE_NUMBER
782           ENDIF
783     
784     ! total solids volume fraction
785           EPs = ONE - EPg
786     
787           IF (EPs .LT. 0.0012d0) THEN
788              h_phip = 2.7d0*(EPs**0.234)
789           ELSEIF (EPs .LT. 0.014d0) THEN
790              h_phip = -0.019d0/(EPs**0.455) + 0.963d0
791           ELSEIF (EPs .LT. 0.25d0) THEN
792              h_phip = 0.868d0*EXP((-0.38*EPs)) - &
793                 0.176d0*EXP((-119.2*EPs))
794           ELSEIF (EPs .LT. 0.455d0) THEN
795              h_phip = -4.59d-5*EXP((19.75*EPs)) + &
796                 0.852d0*EXP((-0.268*EPs))
797           ELSEIF (EPs .LE. 0.59d0) THEN
798              h_phip = (EPs - 0.59d0) * (-1501.d0*(EPs**3) + &
799                 2203.d0*(EPs**2) - 1054.d0*EPs + 162.d0)
800           ELSE
801              h_phip=ZERO
802           ENDIF
803     
804           IF (EPs .LT. 0.18d0) THEN
805               GG_phip = (EPs**0.24)*(1.48d0 + EXP(-18.0*EPs))
806           ELSE
807               GG_phip = ONE
808           ENDIF
809     
810     ! a filter function needed in Igci Filtered/subgrid Model [dimensionless]
811           f_filter = (Inv_Froude**1.6) / ((Inv_Froude**1.6)+0.4d0)
812           h_phip2=h_phip*GG_phip
813           c_function=-h_phip2*f_filter
814           F_Subgrid =(ONE + c_function)
815     
816           IF (SUBGRID_WALL) THEN
817              CALL SUBGRID_DRAG_WALL(F_SubgridWall,vt,IJK)
818           ENDIF
819     
820           lDgA = F_SubgridWall*F_Subgrid * lDgA
821     
822           RETURN
823           END SUBROUTINE SUBGRID_DRAG_IGCI
824     
825     
826     
827     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
828     !                                                                      C
829     !  Subroutine: SUBGRID_DRAG_MILIOLI                                    C
830     !  Purpose: Calculate subgrid correction to the gas-solids drag        C
831     !           coefficient developed by Wen-Yu                            C
832     !                                                                      C
833     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
834     !                                                                      C
835     !  Revision: 1                                                         C
836     !  Purpose: Minor changes, e.g., fix inconsistenty with analogous      C
837     !     calls in des_drag_gp and with new variable density feature       C
838     !  Author: Janine Carney, June 2013                                    C
839     !                                                                      C
840     !  Literature/Document References:                                     C
841     !     Milioli, C. C., et al., Filtered two-fluid models of fluidized   C
842     !        gas-particle flows: new constitutive relations, AICHE J,      C
843     !        doi: 10.1002/aic.14130                                        C
844     !                                                                      C
845     !  Comments:                                                           C
846     !     Still needs to be reviewed for accuracy with source material     C
847     !                                                                      C
848     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
849     
850           SUBROUTINE SUBGRID_DRAG_MILIOLI(lDgA,EPg,Mug,ROg,VREL,DPM,ROs,&
851              IJK)
852     
853     !-----------------------------------------------
854     ! Modules
855     !-----------------------------------------------
856           USE param
857           USE param1
858           USE run, only : filter_size_ratio, SUBGRID_WALL
859           USE constant, only : GRAVITY
860           USE geometry, only : VOL,AXY,DO_K
861           IMPLICIT NONE
862     !-----------------------------------------------
863     ! Dummy arguments
864     !-----------------------------------------------
865     ! drag coefficient
866           DOUBLE PRECISION, INTENT(INOUT) :: lDgA
867     ! gas volume fraction
868           DOUBLE PRECISION, INTENT(IN) :: EPg
869     ! gas laminar viscosity
870           DOUBLE PRECISION, INTENT(IN) :: Mug
871     ! gas density
872           DOUBLE PRECISION, INTENT(IN) :: ROg
873     ! Magnitude of gas-solids relative velocity
874           DOUBLE PRECISION, INTENT(IN) :: VREL
875     ! particle diameter of solids phase M
876           DOUBLE PRECISION, INTENT(IN) :: DPM
877     ! particle density of solids phase M
878           DOUBLE PRECISION, INTENT(IN) :: ROs
879     ! current cell index
880           INTEGER, INTENT(IN) :: IJK
881     !-----------------------------------------------
882     ! Local variables
883     !-----------------------------------------------
884     ! factor to correct the drag for subgrid domain effects
885           DOUBLE PRECISION :: F_Subgrid
886     ! factor to correct the drag for subgrid domain effects arising from
887     ! wall
888           DOUBLE PRECISION :: F_SubGridWall
889     ! particle terminal settling velocity from stokes' formulation
890           DOUBLE PRECISION :: vt
891     ! filter size which is a function of each grid cell volume
892           DOUBLE PRECISION :: filtersize
893     ! inverse Froude number, or dimensionless filter size
894           DOUBLE PRECISION :: Inv_Froude
895     ! dimensionless slip velocity = VREL/vt
896           DOUBLE PRECISION :: vslip
897     ! total solids volume fraction
898           DOUBLE PRECISION :: EPs
899     ! Variables for Milioli model
900           DOUBLE PRECISION :: h1, henv, hlin
901     !-----------------------------------------------
902     
903     ! initialize
904           F_Subgrid = ONE
905           F_SubgridWall = ONE
906     
907     ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
908           vt = GRAVITY*DPM*DPM*(ROs - ROg) / (18.0d0*Mug)
909     ! filter size calculation for each specific gridcell volume
910           IF(DO_K) THEN
911              filtersize = filter_size_ratio * (VOL(IJK)**(ONE/3.0d0))
912           ELSE
913              filtersize = filter_size_ratio * DSQRT(AXY(IJK))
914           ENDIF
915     ! dimensionless inverse of Froude number
916           IF(ABS(vt) > SMALL_NUMBER) THEN
917              Inv_Froude =  filtersize * GRAVITY / vt**2
918           ELSE
919              Inv_Froude =  LARGE_NUMBER
920           ENDIF
921     ! total solids volume fractionn
922           EPs = ONE - EPg
923     ! dimensionless slip velocity between gas and solids phase M
924           Vslip = VREL / vt
925     
926           IF (Inv_Froude .LE. 1.028d0) THEN
927              h1 = (1.076d0 + 0.12d0*Vslip - (0.02d0/(Vslip+0.01d0)))*EPs + &
928                 (0.084d0 + 0.09d0*Vslip - (0.01d0/(0.1d0*Vslip+0.01d0)))
929              IF (EPs .LE. 0.53d0) THEN
930                 henv = (6.8d0*(ONE+EPs)*(EPs**0.3)) / &
931                    (10.d0*(EPs**1.5) + 5.d0)
932              ELSEIF (EPs .GT. 0.53d0 .AND. EPs .LE. 0.65d0) THEN
933                 henv = (2.23d0*((0.65d0-EPs)**(0.45))) / &
934                    ((ONE/EPs)-ONE)
935              ELSEIF (EPs .GT. 0.65d0) THEN
936                 henv=ZERO
937              ENDIF
938     
939           ELSEIF (Inv_Froude .GT. 1.028d0 .AND. &
940                   Inv_Froude .LE. 2.056d0) THEN
941              h1 = (1.268d0 - (0.2d0*Vslip) + (0.14d0/(Vslip+0.01d0)))*EPs + &
942                 (0.385d0 + 0.09d0*Vslip - (0.05d0/(0.2d0*Vslip+0.01d0)))
943              IF (EPs .LE. 0.53d0) THEN
944                 henv = (8.6d0*(ONE+EPs)*(EPs**0.2)) / (10.d0*EPs + 6.3d0)
945              ELSEIF (EPs .GT. 0.53d0 .AND. EPs .LE. 0.65d0) THEN
946                 henv = (0.423d0*((0.65d0-EPs)**0.3)) / (ONE-(EPs**0.4))
947              ELSEIF (EPs .GT. 0.65d0) THEN
948                 henv=ZERO
949              ENDIF
950     
951           ELSEIF (Inv_Froude .GT. 2.056d0 .AND. &
952                   Inv_Froude .LE. 4.112d0) THEN
953              h1 = ((0.018d0*Vslip + 0.1d0)/(0.14d0*Vslip + 0.01d0))*EPs + &
954                 (0.9454d0 - (0.09d0/(0.2d0*Vslip + 0.01d0)))
955              IF (EPs .LE. 0.5d0) THEN
956                 henv = (7.9d0*(ONE+EPs)*(EPs**0.2)) / &
957                    (10.d0*(EPs**0.9) + 5.d0)
958              ELSEIF (EPs .GT. 0.5d0 .AND. EPs .LE. 0.63d0) THEN
959                 henv = (0.705d0*((0.63d0-EPs)**0.3)) / (ONE-(EPs**0.7))
960              ELSEIF (EPs .GT. 0.63d0) THEN
961                 henv=ZERO
962              ENDIF
963     
964           ELSEIF (Inv_Froude .GT. 4.112d0 .AND. &
965                   Inv_Froude .LE. 8.224d0) THEN
966              h1 = ((0.05d0*Vslip+0.3d0)/(0.4d0*Vslip+0.06d0))*EPs + &
967                 (0.9466d0 - (0.05d0/(0.11d0*Vslip+0.01d0)))
968              IF (EPs .LE. 0.45d0) THEN
969                 henv = (7.9d0*(ONE+EPs)*(EPs**0.2)) / &
970                    ((10.d0*(EPs**0.6)) + 3.6d0)
971              ELSEIF (EPs .GT. 0.45d0 .AND. EPs .LE. 0.57d0) THEN
972                 henv = (0.78d0*((0.57d0-EPs)**0.2)) / (ONE-(EPs**0.9))
973              ELSEIF (EPs .GT. 0.57d0) THEN
974                 henv=ZERO
975              ENDIF
976     
977           ELSEIF (Inv_Froude .GT. 8.224d0 .AND. &
978                   Inv_Froude .LE. 12.336d0) THEN
979              h1 = ((1.3d0*Vslip+2.2d0)/(5.2d0*Vslip+0.07d0))*EPs + &
980                 (0.9363d0-(0.11d0/(0.3d0*Vslip+0.01d0)))
981              IF (EPs .LE. 0.35d0) THEN
982                 henv = (7.6d0*(ONE+EPs)*(EPs**0.2)) / &
983                    ((10.d0*(EPs**0.6)) + 3.3d0)
984              ELSEIF (EPs .GT. 0.35d0 .AND. EPs .LE. 0.55d0) THEN
985                 henv = (0.81d0*((0.55d0-EPs)**0.3)) / (ONE-(EPs**0.7))
986              ELSEIF (EPs .GT. 0.55d0) THEN
987                 henv=ZERO
988              ENDIF
989     
990           ELSEIF (Inv_Froude .GT. 12.336d0 .AND. &
991                   Inv_Froude .LE. 16.448d0) THEN
992              h1 = ((2.6d0*Vslip+4.d0)/(10.d0*Vslip+0.08d0))*EPs + &
993                 (0.926d0-(0.17d0/(0.5d0*Vslip+0.01d0)))
994              IF (EPs .LE. 0.25d0) THEN
995                 henv = (8.4d0*(ONE+EPs)*(EPs**0.2)) / &
996                    ((10.d0*(EPs**0.5)) + 3.3d0)
997              ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
998                 henv = (1.01d0*((0.52d0-EPs)**0.03))/(ONE-(EPs**0.9))
999              ELSEIF (EPs .GT. 0.52d0) THEN
1000                 henv=ZERO
1001              ENDIF
1002     
1003           ELSEIF (Inv_Froude .GT. 16.448d0 .AND. &
1004                   Inv_Froude .LE. 20.56d0) THEN
1005              h1 = ((2.5d0*Vslip+4.d0)/(10.d0*Vslip+0.08d0))*EPs + &
1006                 (0.9261d0-(0.17d0/(0.5d0*Vslip+0.01d0)))
1007              IF (EPs .LE. 0.25d0) THEN
1008                 henv = (8.4d0*(ONE+EPs)*(EPs**0.2)) / &
1009                    ((10.d0*(EPs**0.5)) + 3.3d0)
1010              ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1011                 henv = (1.065d0*((0.52d0-EPs)**0.3))/(ONE-EPs)
1012              ELSEIF (EPs .GT.0.52d0) THEN
1013                 henv=ZERO
1014              ENDIF
1015     
1016           ELSEIF (Inv_Froude .GT. 20.56d0) THEN
1017              h1 = ((1.6d0*Vslip+4.d0)/(7.9d0*Vslip+0.08d0))*EPs + &
1018                 (0.9394d0 - (0.22d0/(0.6d0*Vslip+0.01d0)))
1019              IF (EPs .LE. 0.25d0) THEN
1020                 henv = (9.d0*(ONE+EPs)*(EPs**0.15)) / &
1021                    (10.d0*(EPs**0.45) + 4.2d0)
1022              ELSEIF (EPs .GT. 0.25d0 .AND. EPs .LE. 0.52d0) THEN
1023                 henv = (0.91d0*((0.52d0-EPs)**0.4))/(ONE-(EPs**0.6))
1024              ELSEIF (EPs .GT. 0.52d0) THEN
1025                 henv=ZERO
1026              ENDIF
1027           ENDIF
1028     
1029           IF (h1 .GT. ZERO) THEN
1030              hlin=h1
1031           ELSE
1032              hlin=ZERO
1033           ENDIF
1034     
1035           IF (Inv_Froude .LT. 1.028d0) THEN
1036     ! for very small filtered size, the drag wont be changed:
1037     ! F_Subgrid = 1.0 - H where H = 0.0
1038              F_Subgrid = ONE
1039           ELSE
1040     ! MIN(henv,hlin) is H in Milioli paper, 2013
1041              F_Subgrid = ONE - MIN(henv,hlin)
1042           ENDIF
1043     
1044     ! Filtered drag = (1 - H)*Microscopic_drag; it is strange Milioli takes EPs
1045     !     F_Subgrid = EPs*(ONE-hmili)
1046     
1047           IF (SUBGRID_WALL) THEN
1048              CALL SUBGRID_DRAG_WALL(F_SubgridWall,vt,IJK)
1049           ENDIF
1050     
1051           lDgA = F_SubgridWall*F_Subgrid * lDgA
1052     
1053     
1054           RETURN
1055           END SUBROUTINE SUBGRID_DRAG_MILIOLI
1056     
1057     
1058     
1059     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1060     !                                                                      C
1061     !  Subroutine: SUBGRID_DRAG_WALL                                       C
1062     !  Purpose: Calculate subgrid correction arising from wall to the      C
1063     !     gas-solids drag coefficient                                      C
1064     !                                                                      C
1065     !  Author: Sebastien Dartevelle, LANL, May 2013                        C
1066     !                                                                      C
1067     !  Revision: 1                                                         C
1068     !  Author: Janine Carney, June 2013                                    C
1069     !                                                                      C
1070     !  Literature/Document References:                                     C
1071     !     Igci, Y., and Sundaresan, S., Verification of filtered two-      C
1072     !        fluid models for gas-particle flows in risers, AICHE J.,      C
1073     !        2011, 57 (10), 2691-2707.                                     C
1074     !                                                                      C
1075     !  Comments: Currently only valid for free-slip wall but no checks     C
1076     !     are made to ensure user has selected free-slip wall when this    C
1077     !     option is invoked                                                C
1078     !                                                                      C
1079     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1080     
1081           SUBROUTINE SUBGRID_DRAG_WALL(lSubgridWall,vt,IJK)
1082     
1083     !-----------------------------------------------
1084     ! Modules
1085     !-----------------------------------------------
1086           USE param
1087           USE param1
1088           USE constant, only : GRAVITY
1089           USE cutcell, only : DWALL
1090           IMPLICIT NONE
1091     !-----------------------------------------------
1092     ! Dummy arguments
1093     !-----------------------------------------------
1094     ! factor to correct the drag for subgrid domain effects arising from
1095     ! wall
1096           DOUBLE PRECISION, INTENT(OUT) :: lSubGridWall
1097     ! particle terminal settling velocity from stokes' formulation
1098           DOUBLE PRECISION, INTENT(IN) :: vt
1099     ! current cell index
1100           INTEGER, INTENT(IN) :: IJK
1101     !-----------------------------------------------
1102     ! Local parameters
1103     !-----------------------------------------------
1104     ! values are only correct for FREE-Slip walls
1105           DOUBLE PRECISION, PARAMETER :: a22=6.0d0, b22=0.295d0
1106     !-----------------------------------------------
1107     ! Local variables
1108     !-----------------------------------------------
1109     ! dimensionless distance to wall
1110           DOUBLE PRECISION :: x_d
1111     !-----------------------------------------------
1112     ! initialize
1113           lSubgridWall = ONE
1114     
1115     ! dimensionless distance to the Wall
1116           x_d = DWALL(IJK) * GRAVITY / vt**2
1117     
1118     ! decrease exponentionally away from the wall
1119     ! more complex model could be implemented with JJ wall model
1120           lSubgridWall = ONE / ( ONE + a22 * (EXP(-b22*x_d)) )
1121     
1122           RETURN
1123           END SUBROUTINE SUBGRID_DRAG_WALL
1124     
1125     
1126     
1127     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1128     !                                                                      C
1129     !  Subroutine: DRAG_KOCH_HILL                                          C
1130     !  Purpose: Calculate the gas-solids drag coefficient                  C
1131     !                                                                      C
1132     !  Author: Clay Sutton (Lehigh University)            Date: 14-Jul-04  C
1133     !                                                                      C
1134     !  Revision: 1                                                         C
1135     !  Author: Sofiane Benyahia                           Date: 21-Jan-05  C
1136     !                                                                      C
1137     !  Literature/Document References:                                     C
1138     !     Benyahia S, Syamlal M, O'Brien TJ (2006). Powder Technology      C
1139     !        162: 166-174.                                                 C
1140     !     Hill RJ, Koch DL, Ladd JC (2001). Journal of Fluid Mechanics     C
1141     !        448: 213-241.                                                 C
1142     !     Hill RJ, Koch DL, Ladd JC (2001). Journal of Fluid Mechanics     C
1143     !        448: 243-278.                                                 C
1144     !                                                                      C
1145     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1146     
1147           SUBROUTINE DRAG_KOCH_HILL(lDgA,EPg,Mug,ROPg,VREL,&
1148                      DPM,DPA,PHIS)
1149     
1150     !-----------------------------------------------
1151     ! Modules
1152     !-----------------------------------------------
1153           USE param
1154           USE param1
1155           IMPLICIT NONE
1156     !-----------------------------------------------
1157     ! Dummy arguments
1158     !-----------------------------------------------
1159     ! drag coefficient
1160           DOUBLE PRECISION, INTENT(OUT) :: lDgA
1161     ! gas volume fraction
1162           DOUBLE PRECISION, INTENT(IN) :: EPg
1163     ! gas laminar viscosity
1164           DOUBLE PRECISION, INTENT(IN) :: Mug
1165     ! gas density*EP_g
1166           DOUBLE PRECISION, INTENT(IN) :: ROPg
1167     ! Magnitude of gas-solids relative velocity
1168           DOUBLE PRECISION, INTENT(IN) :: VREL
1169     ! particle diameter of solids phase M
1170           DOUBLE PRECISION, INTENT(IN) :: DPM
1171     ! average particle diameter if pcf otherwise DPM again
1172           DOUBLE PRECISION, INTENT(IN) :: DPA
1173     ! total solids volume fraction of solids phases
1174           DOUBLE PRECISION, INTENT(IN) :: PHIS
1175     !-----------------------------------------------
1176     ! Local variables
1177     !-----------------------------------------------
1178     ! Reynolds number
1179           DOUBLE PRECISION :: RE
1180     ! transition Reynolds numbers
1181           DOUBLE PRECISION :: Re_Trans_1, Re_Trans_2
1182     ! Stokes Drag Force
1183           DOUBLE PRECISION :: F_STOKES
1184     ! zero Re function for low Reynolds number
1185           DOUBLE PRECISION :: F_0
1186     ! inertial function for low Reynolds number
1187           DOUBLE PRECISION :: F_1
1188     ! zero Re function for high Reynolds number
1189           DOUBLE PRECISION :: F_2
1190     ! inertial function for high Reynolds number
1191           DOUBLE PRECISION :: F_3
1192     ! dimensionless drag force F
1193           DOUBLE PRECISION :: F
1194     ! weighting factor to compute F_0 and F_2
1195           DOUBLE PRECISION :: ww
1196     !-----------------------------------------------
1197     
1198     
1199           IF(Mug > ZERO) THEN
1200     ! Note the presence of gas volume fraction in ROPG and factor of 1/2
1201              RE = 0.5D0*DPA*VREL*ROPG/Mug        ! if pcf DPA otherwise DPM
1202           ELSE
1203              RE = LARGE_NUMBER
1204           ENDIF
1205     
1206           F_STOKES = 18.D0*Mug*EPg*EPg/DPM**2    ! use DPM
1207           ww = EXP(-10.0D0*(0.4D0-phis)/phis)
1208     
1209           IF(phis > 0.01D0 .AND. phis < 0.4D0) THEN
1210              F_0 = (1.0D0-ww) * (1.0D0 + 3.0D0*dsqrt(phis/2.0D0) + &
1211                 135.0D0/64.0D0*phis*LOG(phis) + 17.14D0*phis) / &
1212                 (1.0D0 + 0.681D0*phis - 8.48D0*phis*phis + &
1213                 8.16D0*phis**3) + ww*10.0D0*phis/(1.0D0-phis)**3
1214           ELSEIF(phis >= 0.4D0) THEN
1215              F_0 = 10.0D0*phis/(1.0D0-phis)**3
1216           ENDIF
1217     
1218           IF(phis > 0.01D0 .AND. phis <= 0.1D0) THEN
1219             F_1 = dsqrt(2.0D0/phis) / 40.0D0
1220           ELSE IF(phis > 0.1D0) THEN
1221             F_1 = 0.11D0 + 5.1D-04 * exp(11.6D0*phis)
1222           ENDIF
1223     
1224           IF(phis < 0.4D0) THEN
1225             F_2 = (1.0D0-ww) * (1.0D0 + 3.0D0*dsqrt(phis/2.0D0) + &
1226                135.0D0/64.0D0*phis*LOG(phis) + 17.89D0*phis) / &
1227                (1.0D0 + 0.681D0*phis - 11.03D0*phis*phis + &
1228                15.41D0*phis**3)+ ww*10.0D0*phis/(1.0D0-phis)**3
1229           ELSE
1230              F_2 = 10.0D0*phis/(1.0D0-phis)**3
1231           ENDIF
1232     
1233           IF(phis < 0.0953D0) THEN
1234              F_3 = 0.9351D0*phis + 0.03667D0
1235           ELSE
1236              F_3 = 0.0673D0 + 0.212D0*phis +0.0232D0/(1.0-phis)**5
1237           ENDIF
1238     
1239           Re_Trans_1 = (F_2 - 1.0D0)/(3.0D0/8.0D0 - F_3)
1240           Re_Trans_2 = (F_3 + dsqrt(F_3*F_3 - 4.0D0*F_1 &
1241                *(F_0-F_2))) / (2.0D0*F_1)
1242     
1243           IF(phis <= 0.01D0 .AND. RE <= Re_Trans_1) THEN
1244              F = 1.0D0 + 3.0D0/8.0D0*RE
1245           ELSEIF(phis > 0.01D0 .AND. RE <= Re_Trans_2) THEN
1246              F = F_0 + F_1*RE*RE
1247           ELSEIF(phis <= 0.01D0 .AND. RE > Re_Trans_1 .OR.   &
1248              phis >  0.01D0 .AND. RE > Re_Trans_2) THEN
1249              F = F_2 + F_3*RE
1250           ELSE
1251              F = zero
1252           ENDIF
1253     
1254           lDgA = F * F_STOKES
1255           IF (RE == ZERO) lDgA = ZERO
1256     
1257           RETURN
1258           END SUBROUTINE DRAG_KOCH_HILL
1259     
1260     
1261     
1262     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1263     !                                                                      C
1264     !  Subroutine: DRAG_BVK                                                C
1265     !  Purpose: Calculate the gas-solids drag coefficient                  C
1266     !                                                                      C
1267     !  Literature/Document References:                                     C
1268     !     Beetstra, van der Hoef, Kuipers, Chem. Eng. Science 62           C
1269     !     (Jan 2007)                                                       C
1270     !                                                                      C
1271     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1272     
1273           SUBROUTINE DRAG_BVK(lDgA,EPg,Mug,ROPg,VREL,&
1274                      DPM,DPA,PHIS)
1275     
1276     !-----------------------------------------------
1277     ! Modules
1278     !-----------------------------------------------
1279           USE param
1280           USE param1
1281           IMPLICIT NONE
1282     !-----------------------------------------------
1283     ! Dummy arguments
1284     !-----------------------------------------------
1285     ! drag coefficient
1286           DOUBLE PRECISION, INTENT(OUT) :: lDgA
1287     ! gas volume fraction
1288           DOUBLE PRECISION, INTENT(IN) :: EPg
1289     ! gas laminar viscosity
1290           DOUBLE PRECISION, INTENT(IN) :: Mug
1291     ! gas density*EP_g
1292           DOUBLE PRECISION, INTENT(IN) :: ROPg
1293     ! magnitude of gas-solids relative velocity
1294           DOUBLE PRECISION, INTENT(IN) :: VREL
1295     ! particle diameter of solids phase M or
1296           DOUBLE PRECISION, INTENT(IN) :: DPM
1297     ! average particle diameter
1298           DOUBLE PRECISION, INTENT(IN) :: DPA
1299     ! total solids volume fraction of solids phases
1300           DOUBLE PRECISION, INTENT(IN) :: PHIS
1301     !-----------------------------------------------
1302     ! Local variables
1303     !-----------------------------------------------
1304     ! Reynolds number
1305           DOUBLE PRECISION :: RE
1306     ! Stokes Drag Force
1307           DOUBLE PRECISION :: F_STOKES
1308     ! dimensionless drag force F
1309           DOUBLE PRECISION :: F
1310     !-----------------------------------------------
1311     
1312           IF(Mug > ZERO) THEN
1313     ! Note the presence of gas volume fraction in ROPG
1314              RE = DPA*VREL*ROPg/Mug        ! use DPA
1315           ELSE
1316              RE = LARGE_NUMBER
1317           ENDIF
1318     
1319     ! eq(9) BVK J. fluid. Mech. 528, 2005
1320     ! (this F_Stokes is /= of Koch_Hill by a factor of ep_g)
1321           F_STOKES = 18D0*Mug*EPg/DPM**2   ! use DPM
1322     
1323           F = 10d0*phis/EPg**2 + EPg**2*(ONE+1.5d0*DSQRT(phis))
1324           F = F + 0.413d0*RE/(24.d0*EPg**2) * &
1325                  (ONE/EPg + 3d0*EPg*phis + 8.4d0/RE**0.343) / &
1326                  (ONE+10.d0**(3d0*phis)/RE**(0.5+2.d0*phis))
1327     
1328           lDgA = F*F_STOKES
1329           IF (RE == ZERO) lDgA = ZERO
1330     
1331           RETURN
1332           END SUBROUTINE DRAG_BVK
1333     
1334     
1335     
1336     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1337     !                                                                      C
1338     !  Subroutine: DRAG_HYS                                                C
1339     !  Purpose: Calculate the gas-solids drag coefficient                  C
1340     !                                                                      C
1341     !  Literature/Document References:                                     C
1342     !     Yin, X, Sundaresan, S. (2009). AIChE Journal 55: no 6, 1352-     C
1343     !     1368                                                             C
1344     !                                                                      C
1345     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1346     
1347           SUBROUTINE DRAG_HYS(lDgA,EPg,Mug,ROPg,VREL,&
1348                      DPM,DPA,Y_I,EP_sM,PHIS,M,MAXM,IJK)
1349     
1350     !-----------------------------------------------
1351     ! Modules
1352     !-----------------------------------------------
1353           USE param
1354           USE param1
1355           USE drag, only : beta_ij
1356           USE run, only : LAM_HYS
1357           IMPLICIT NONE
1358     !-----------------------------------------------
1359     ! Dummy arguments
1360     !-----------------------------------------------
1361     ! drag coefficient
1362           DOUBLE PRECISION, INTENT(OUT) :: lDgA
1363     ! gas volume fraction
1364           DOUBLE PRECISION, INTENT(IN) :: EPg
1365     ! gas laminar viscosity
1366           DOUBLE PRECISION, INTENT(IN) :: Mug
1367     ! gas density*EP_g
1368           DOUBLE PRECISION, INTENT(IN) :: ROPg
1369     ! magnitude of gas-solids relative velocity
1370           DOUBLE PRECISION, INTENT(IN) :: VREL
1371     ! local variable for the particle diameter
1372           DOUBLE PRECISION :: DPM(DIM_M)
1373     ! average particle diameter
1374           DOUBLE PRECISION, INTENT(IN) :: DPA
1375     ! diameter ratio in polydisperse systems
1376           DOUBLE PRECISION, INTENT(IN) :: Y_i
1377     ! local variable for the solids volume fraction
1378           DOUBLE PRECISION :: EP_sM(DIM_M)
1379     ! total solids volume fraction of solids phases
1380           DOUBLE PRECISION, INTENT(IN) :: PHIS
1381     ! current solids phase index and fluid cell index
1382           INTEGER, INTENT(IN) :: M, IJK
1383     ! maximum number of solids phases
1384           INTEGER, INTENT(IN) :: MAXM
1385     !-----------------------------------------------
1386     ! Local variables
1387     !-----------------------------------------------
1388     ! minimum particle diameter in mixture
1389           DOUBLE PRECISION :: DPmin
1390     ! Index for particles of other species
1391           INTEGER :: L
1392     ! Reynolds number
1393           DOUBLE PRECISION :: RE
1394     ! Stokes Drag Force
1395           DOUBLE PRECISION :: F_STOKES
1396     ! dimensionless drag force F
1397           DOUBLE PRECISION :: F
1398     ! Polydisperse correction factor for YS drag relation
1399           DOUBLE PRECISION :: a_YS
1400     ! Lubrication interaction prefactor in YS drag relation
1401           DOUBLE PRECISION :: alpha_YS
1402     ! Friction coefficient for a particle of type i (HYS drag relation)
1403           DOUBLE PRECISION :: beta_i_HYS
1404     ! Friction coefficient for a particle of type j (HYS drag relation)
1405           DOUBLE PRECISION :: beta_j_HYS
1406     ! Stokes drag of a particle of type j
1407           DOUBLE PRECISION :: FSTOKES_j
1408     ! Diameter ratio for particle of type j
1409           DOUBLE PRECISION :: Y_i_J
1410     ! Variable for Beetstra et. al. drag relation
1411           DOUBLE PRECISION :: F_D_BVK
1412     ! Variable for YS drag relation
1413           DOUBLE PRECISION :: F_YS
1414     !-----------------------------------------------
1415     
1416           IF (Mug > ZERO) THEN
1417     ! Note the presence of gas volume fraction in ROPG
1418              RE = DPA*VREL*ROPg/Mug   ! use DPA
1419           ELSE
1420              RE = LARGE_NUMBER
1421           ENDIF
1422     
1423     ! (this F_Stokes is /= of Koch_Hill by a factor of ep_g/ep_sm)
1424     ! (this F_Stokes is /= of BVK by a factor of 1/ep_sm)
1425           F_STOKES = 18D0*Mug*EPg*EP_SM(M)/DPM(M)**2   ! use DPM
1426     
1427     ! Find smallest diameter if number of particle types is greater than 1
1428           Dpmin= DPM(1)
1429           IF (MAXM > 1) THEN
1430              DO L=2,MAXM
1431                 Dpmin = MIN(Dpmin,DPM(L))
1432              ENDDO
1433           ENDIF
1434     
1435           a_YS = 1d0 - 2.66d0*phis + 9.096d0*phis**2 - 11.338d0*phis**3
1436     
1437     ! Calculate the prefactor of the off-diagonal friction coefficient
1438     ! Use default value of lamdba if there are no particle asparities
1439           alpha_YS = 1.313d0*LOG10(DPmin/lam_HYS) - 1.249d0
1440     
1441     
1442     ! Beetstra correction for monodisperse drag
1443           F_D_BVK = ZERO
1444           F = 10.d0*phis/EPg**2 + EPg**2*(ONE+1.5d0*DSQRT(phis))
1445     
1446           IF(RE > ZERO) F_D_BVK = F + 0.413d0*RE/(24.d0*EPg**2)*&
1447              (ONE/EPg + 3.d0*EPg*phis + 8.4d0/RE**0.343d0) / &
1448              (ONE+10**(3.d0*phis)/RE**(0.5d0+2.d0*phis))
1449     
1450     ! YS correction for polydisperse drag
1451           F_YS = 1d0/EPg + (F_D_BVK - 1d0/EPg)*&
1452                            (a_YS*Y_i+(1d0-a_YS)*Y_i**2)
1453           F_YS = F_YS*F_STOKES
1454           beta_i_HYS = F_YS
1455     
1456           DO L= 1,MAXM
1457              IF (L /= M) THEN
1458                 Y_i_J = DPM(L)/DPA
1459                 beta_j_HYS = 1.d0/EPg + (F_D_BVK - 1.d0/EPg) * &
1460                    (a_YS*Y_i_J + (1d0-a_YS)*Y_i_J**2)
1461                 FSTOKES_j = 18.D0*Mug*EP_sM(L)*EPg/&
1462                    DPM(L)**2
1463     
1464                 beta_j_HYS = beta_j_HYS*FSTOKES_j
1465     
1466     ! Calculate off-diagonal friction coefficient
1467                 beta_ij(IJK,M,L) = ZERO
1468     
1469     ! This if statement prevents NaN values from appearing for beta_ij
1470                 IF (EP_sM(M) > ZERO .AND. EP_SM(L) > ZERO) &
1471                    beta_ij(IJK,M,L) = (2.d0*alpha_YS*EP_sM(M)*EP_sM(L))/ &
1472                       (EP_sM(M)/beta_i_HYS + EP_sM(L)/beta_j_HYS)
1473                 F_YS = F_YS + beta_ij(IJK,M,L)
1474     
1475              ENDIF   ! end if (J/=M)
1476           ENDDO   ! end do (J=1,MAXM)
1477     
1478     
1479           lDgA = F_YS
1480     
1481           RETURN
1482           END SUBROUTINE DRAG_HYS
1483