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