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