File: RELATIVE:/../../../mfix.git/model/des/drag_gp_des.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DES_DRAG_GP                                             C
4     !  Purpose: Calculate the gas-particle drag coefficient using          C
5     !           the gas velocity interpolated to the particle position     C
6     !           and the particle velocity.                                 C
7     !           Invoked from des_drag_gs and calc_des_drag_gs              C
8     !                                                                      C
9     !  Comments: The BVK drag model and all drag models with the           C
10     !            polydisperse correction factor (i.e., suffix _PCF)        C
11     !            require an average particle diameter. This has been       C
12     !            loosely defined for discrete particles based on their     C
13     !            solids phase                                              C
14     !                                                                      C
15     !  Variables referenced:                                               C
16     !  Variables modified:                                                 C
17     !  Local variables:                                                    C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20           SUBROUTINE DES_DRAG_GP(NP, PARTICLE_VEL, FLUID_VEL, EPg)
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE compar
26           USE constant
27           USE discretelement
28           USE drag
29           USE fldvar
30           USE functions
31           USE geometry
32           USE indices
33           USE machine, only: start_log, end_log
34           USE param
35           USE param1
36           USE physprop
37           USE run
38           USE sendrecv
39           USE ur_facs
40     
41           IMPLICIT NONE
42     !-----------------------------------------------
43     ! Dummy arguments
44     !-----------------------------------------------
45     ! particle number id.
46           INTEGER , INTENT(IN) :: NP
47     ! particle velocity
48           DOUBLE PRECISION, INTENT(IN) :: PARTICLE_VEL(3)
49     ! fluid velocity interpolated to particle position
50           DOUBLE PRECISION, INTENT(IN) :: FLUID_VEL(3)
51     ! Gas phase volume fraction.
52           DOUBLE PREcISION, INTENT(IN) :: EPg
53     
54     !-----------------------------------------------
55     ! Local variables
56     !-----------------------------------------------
57     ! indices, associated with current particle
58           INTEGER :: IJK
59     ! solids phase index, associated with current particle
60           INTEGER :: M
61     ! Slip velocity and its magnitude
62           DOUBLE PRECISION :: VSLP(3), VREL
63     ! gas laminar viscosity redefined here to set viscosity at pressure
64     ! boundaries
65           DOUBLE PRECISION :: Mu
66     ! drag coefficient
67           DOUBLE PRECISION :: DgA
68     ! current value of F_gs (i.e., without underrelaxation)
69           DOUBLE PRECISION F_gstmp
70     ! indices of solids phases (continuous, discrete)
71           INTEGER :: lM
72     ! correction factors for implementing polydisperse drag model
73     ! proposed by van der Hoef et al. (2005)
74           DOUBLE PRECISION :: F_cor, tSUM
75     ! average particle diameter in polydisperse systems
76           DOUBLE PRECISION :: DPA
77     ! diameter ratio in polydisperse systems
78           DOUBLE PRECISION :: Y_i
79     ! total solids volume fraction
80           DOUBLE PRECISION :: PHIS
81     ! aliases for void fraction, gas density, gas bulk density,
82     ! solids volume fraction, particle diameter, particle density
83           DOUBLE PRECISION :: ROg, ROPg, DPM
84     !-----------------------------------------------
85     
86     ! values based on current particle
87           IJK = PIJK(NP,4)
88     ! solids phase index of current particle
89           M = PIJK(NP,5)
90     ! Gas material and bulk densities
91           ROg = RO_G(IJK)
92           ROPg = RO_G(IJK) * EPg
93     ! Laminar viscosity.
94           Mu = MU_G(IJK)
95     ! Slip velocity and its magnitude
96           VSLP = FLUID_VEL - PARTICLE_VEL
97           VREL = SQRT(dot_product(VSLP, VSLP))
98     ! assign variables for short dummy arguments
99           DPM = 2.0d0*DES_RADIUS(NP)
100     ! Total solids volume fraction.
101           PHIS = ONE - EPg
102     
103     ! determine the drag coefficient
104           SELECT CASE(DRAG_TYPE_ENUM)
105     
106           CASE (SYAM_OBRIEN)
107              CALL DRAG_SYAM_OBRIEN(DgA,EPG,Mu,ROg,VREL,DPM)
108     
109           CASE (GIDASPOW)
110              CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
111     
112           CASE (GIDASPOW_BLEND)
113              CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPM)
114     
115           CASE (WEN_YU)
116              CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPM)
117     
118           CASE (KOCH_HILL)
119              CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPM,phis)
120     
121           CASE (USER_DRAG)
122              CALL DRAG_USR(IJK,NP,DgA,EPg,Mu,ROg,VREL,DPM,DES_RO_S(M), &
123                 FLUID_VEL(1), FLUID_VEL(2), FLUID_VEL(3))
124     
125           CASE DEFAULT
126     
127     ! calculate the average particle diameter and particle ratio
128              tSUM = ZERO
129              DO lM = 1,DES_MMAX
130                 IF(PHIS > ZERO) THEN
131                    tSUM = tSUM + DES_ROP_S(IJK,lM) / &
132                       (PHIS*DES_RO_S(lM)*DES_D_p0(lM))
133                  ELSE
134                    tSUM = tSUM + ONE/DES_D_p0(lM)
135                  ENDIF
136              ENDDO
137              IF(DES_CONTINUUM_HYBRID) THEN
138                 DO lM = 1,SMAX
139                    IF (phis .GT. ZERO) THEN
140                       tSUM = tSUM + EP_S(IJK,lM)/(PHIS*D_P(IJK,lM))
141                    ELSE
142                       tSUM = tSUM + ONE/D_P(IJK,lM)
143                    ENDIF
144                 ENDDO
145              ENDIF
146     
147              DPA = ONE / tSUM
148              Y_i = DPM * tSUM
149     
150              SELECT CASE(DRAG_TYPE_ENUM)
151     
152     
153              CASE (GIDASPOW_PCF)
154                 CALL DRAG_GIDASPOW(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
155              CASE (GIDASPOW_BLEND_PCF)
156                 CALL DRAG_GIDASPOW_BLEND(DgA,EPg,Mu,ROg,ROPg,VREL,DPA)
157              CASE (WEN_YU_PCF)
158                 CALL DRAG_WEN_YU(DgA,EPg,Mu,ROPg,VREL,DPA)
159              CASE (KOCH_HILL_PCF)
160                 CALL DRAG_KOCH_HILL(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
161              CASE (BVK)
162                 CALL DRAG_BVK(DgA,EPg,Mu,ROPg,VREL,DPM,DPA,phis)
163     
164              CASE DEFAULT
165                 CALL START_LOG
166                 IF(DMP_LOG) WRITE (*, '(A,A)') &
167                    'Unknown DRAG_TYPE: ', DRAG_TYPE
168                 WRITE (UNIT_LOG, '(A,A)') 'Unknown DRAG_TYPE: ', DRAG_TYPE
169                 CALL END_LOG
170                 CALL mfix_exit(myPE)
171              END SELECT   ! end selection of drag_type
172     
173     ! Modify drag coefficient to account for possible corrections and for
174     ! differences between Model B and Model A; see erratum Beetstra (2007)
175              IF(Model_B) THEN
176                 IF (M == 1) THEN
177                    F_cor = (EPg*Y_i + PHIS*Y_i**2)
178                 ELSE
179                    F_cor = (EPg*Y_i + PHIS*Y_i**2 + &
180                       0.064d0*EPg*Y_i**3)
181                 ENDIF
182              ELSE
183                 F_cor = Y_i
184              ENDIF
185              DgA = ONE/(Y_i*Y_i) * DgA * F_cor
186     
187           END SELECT   ! end selection of drag_type
188     
189     
190     
191     ! Calculate the drag coefficient (Model B coeff = Model A coeff/EP_g)
192           IF(MODEL_B) THEN
193              F_gstmp = DgA * PVOL(NP)/EP_G(IJK)
194           ELSE
195              F_gstmp = DgA * PVOL(NP)
196           ENDIF
197     
198     ! Determine drag force coefficient accounting for any under relaxation
199     ! f_gp() =  single particle drag excluding vector(v_g - v_p)
200           F_gp(NP) = (ONE - UR_F_gs) * F_gp(NP) + UR_F_gs * F_gstmp
201     
202           RETURN
203           END SUBROUTINE DES_DRAG_GP
204