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