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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DRAG_ss                                                 C
4     !  Purpose: This module computes the coefficient of drag between       C
5     !  two solids phases (M and L).                                        C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 20-MAY-96  C
8     !                                                                      C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11           SUBROUTINE DRAG_SS(L, M)
12     
13     ! Modules
14     !---------------------------------------------------------------------//
15           use compar, only: mype
16           use exit, only: mfix_exit
17           use functions, only: funlm
18           use run, only: granular_energy, kt_type, kt_type_enum
19           use run, only: lun_1984, simonin_1996, ahmadi_1995
20           use run, only: ia_2005, gd_1999, ghd_2007, gtsh_2012
21           use usr_prop, only: usr_fss, calc_usr_prop
22           use usr_prop, only: solidssolids_drag
23           IMPLICIT NONE
24     
25     ! Dummy Arguments
26     !---------------------------------------------------------------------//
27     ! Index of solids phases
28           INTEGER, INTENT(IN) :: L, M
29     
30     ! Local Variables
31     !---------------------------------------------------------------------//
32     ! index for storing solids-solids drag coefficients in the upper
33     ! triangle of the matrix
34           INTEGER :: LM
35     
36     !---------------------------------------------------------------------//
37           LM = FUNLM(L,M)
38     
39           IF (USR_FSS(LM)) THEN
40              CALL CALC_USR_PROP(SolidsSolids_Drag,lL=L,lM=M)
41           ELSE
42              IF (GRANULAR_ENERGY) THEN
43                 SELECT CASE (KT_TYPE_ENUM)
44                    CASE(LUN_1984, SIMONIN_1996, AHMADI_1995)
45                       CALL DRAG_SS_SYAM(L, M)
46                    CASE (IA_2005)
47                       CALL DRAG_SS_IA(L, M)
48                    CASE (GD_1999, GTSH_2012)
49     ! strictly speaking gd and gtsh are monodisperse theories and so
50     ! do not have solids-solids drag
51                       RETURN
52     ! ghd theory is a polydisperse theory but is self contained and does
53     ! not invoke this routine
54                    CASE (GHD_2007)
55                       RETURN
56                 CASE DEFAULT
57     ! should never hit this
58                    WRITE (*, '(A)') 'DRAG_SS'
59                    WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
60                 call mfix_exit(myPE)
61                 END SELECT
62              ELSE
63                 CALL DRAG_SS_SYAM(L, M)
64              ENDIF
65           ENDIF
66     
67           RETURN
68           END SUBROUTINE DRAG_SS
69     
70     
71     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
72     !                                                                      C
73     !  Subroutine: DRAG_SS_SYAM                                            C
74     !  Purpose: Calculate the solids-solids drag coefficient between a     C
75     !           continuous solids phase and discrete solids                C
76     !                                                                      C
77     !  Literature/Document References:                                     C
78     !     M. Syamlal. 1987. The particle-particle drag term in a           C
79     !        multiparticle model of fluidization. Technical Report.        C
80     !        DOE/MC/21353-2373. Office of Fossil Energy, Morgantown        C
81     !        Energy Technology Center, Morgantown, West Virginia.          C
82     !     Gera, D., Syamlal, M., O'Brien T.J. 2004. International Journal  C
83     !        of Multiphase Flow, 30, p419-428.                             C
84     !                                                                      C
85     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
86           SUBROUTINE DRAG_SS_SYAM(L,M)
87     
88     ! Modules
89     !---------------------------------------------------------------------//
90           use compar, only: mype, ijkstart3, ijkend3
91           USE constant, only: segregation_slope_coefficient
92           USE drag, only: f_ss
93           USE fldvar, only: d_p, ro_s, rop_s, theta_m, p_star
94           USE fldvar, only: u_s, v_s, w_s
95           USE fun_avg, only: avg_x_e, avg_y_n, avg_z_t
96           use functions, only: wall_at, funlm
97           use functions, only: im_of, jm_of, km_of
98           USE indices, only: i_of
99           USE physprop, only: close_packed
100           USE rdf, only: g_0
101           IMPLICIT NONE
102     
103     ! Dummy arguments
104     !---------------------------------------------------------------------//
105     ! solids phase indices
106           INTEGER, INTENT(IN) :: M, L
107     
108     ! Local variables
109     !---------------------------------------------------------------------//
110     ! indices
111           INTEGER :: I, IJK, IMJK, IJMK, IJKM
112     ! index for storing solids-solids drag coefficients in the upper
113     ! triangle of the matrix
114           INTEGER :: LM
115     ! cell center value of U_sm, U_sl, V_sm, V_sl, W_sm, W_sl
116           DOUBLE PRECISION :: USCM, USCL, VSCM, VSCL, WSCM, WSCL
117     ! relative velocity between solids phase m and l
118           DOUBLE PRECISION :: VREL
119     ! particle diameters of phase M and phase L
120           DOUBLE PRECISION :: D_pm, D_pl
121     ! particle densities of phase M and phase L
122           DOUBLE PRECISION :: RO_M, RO_L
123     ! radial distribution function between phases M and L
124           DOUBLE PRECISION :: G0_ML
125     ! solid-solid drag coefficient
126           DOUBLE PRECISION :: lDss
127     
128     !---------------------------------------------------------------------//
129           LM = FUNLM(L,M)
130     
131           DO IJK = ijkstart3, ijkend3
132              IF (WALL_AT(IJK)) CYCLE
133     ! Evaluate at all flow boundaries and fluid cellls. This is unlike
134     ! calculation of the fluid-solid drag coefficient, which is only
135     ! evaluated in fluid cells and pressure inflow cells
136     
137              I = I_OF(IJK)
138              IMJK = IM_OF(IJK)
139              IJMK = JM_OF(IJK)
140              IJKM = KM_OF(IJK)
141     
142     ! calculating velocity components at i, j, k (cell center)
143              USCL = AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I)
144              VSCL = AVG_Y_N(V_S(IJMK,L),V_S(IJK,L))
145              WSCL = AVG_Z_T(W_S(IJKM,L),W_S(IJK,L))
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 solids-solids relative velocity
152              VREL = SQRT((USCL - USCM)**2 + (VSCL - VSCM)**2 + &
153                          (WSCL - WSCM)**2)
154     
155     ! setting aliases for easy reference
156              D_PM = D_P(IJK,M)
157              D_PL = D_P(IJK,L)
158              RO_M = RO_S(IJK,M)
159              RO_L = RO_S(IJK,L)
160              G0_ML = G_0(IJK,L,M)
161     
162     ! evaluating the solids-solids drag coefficient
163              CALL DRAG_SS_SYAM0(ldss, d_pm, d_pl, ro_m, ro_l, g0_ml, vrel)
164              F_SS(IJK,LM) = lDss*ROP_S(IJK,M)*ROP_S(IJK,L)
165     
166     ! Gera: accounting for particle-particle drag due to enduring contact
167     ! in a close-packed system.
168              IF(CLOSE_PACKED(M) .AND. CLOSE_PACKED(L)) &
169                 F_SS(IJK,LM) = F_SS(IJK,LM) + &
170                    SEGREGATION_SLOPE_COEFFICIENT*P_star(IJK)
171     
172           ENDDO    ! end do (ijk=ijkstart3,ijkend3)
173     
174           RETURN
175           END SUBROUTINE DRAG_SS_SYAM
176     
177     
178     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
179     !                                                                      C
180     !  Subroutine: DRAG_SS_SYAM0                                           C
181     !  Purpose: Created as a work-around for des/hybrid cases that need    C
182     !  to use this model while passing particle type information.          C
183     !                                                                      C
184     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
185           SUBROUTINE DRAG_SS_SYAM0(ldss, d_pm, d_pl, ro_m, ro_l, g0_ml, vrel)
186     
187     ! Modules 
188     !---------------------------------------------------------------------//
189           use param1, only: one
190           use constant, only: c_e, c_f, pi
191           IMPLICIT NONE
192     
193     ! Dummy arguments
194     !---------------------------------------------------------------------//
195     ! solid-solid drag coefficient
196           DOUBLE PRECISION, INTENT(OUT) :: lDss
197     ! particle diameters of phase M and phase L
198           DOUBLE PRECISION, INTENT(IN) :: D_pm, D_pl
199     ! particle densities of phase M and phase L
200           DOUBLE PRECISION, INTENT(IN) :: RO_M, RO_L
201     ! radial distribution function between phases M and L
202           DOUBLE PRECISION, INTENT(IN) :: G0_ML
203     ! relative velocity between solids phase m and l
204           DOUBLE PRECISION, INTENT(IN) :: VREL
205     
206     ! Local variables
207     !---------------------------------------------------------------------//
208     ! Sum of particle diameters
209           DOUBLE PRECISION :: DPSUM
210     ! Intermediate calculation
211           DOUBLE PRECISION :: const
212     !---------------------------------------------------------------------//
213     
214           DPSUM = D_PL + D_PM
215           const = 3.d0*(ONE + C_E)*(PI/2.d0 + C_F*PI*PI/8.d0)*&
216                   DPSUM**2/(2.d0*PI*(RO_L*D_PL**3+RO_M*D_PM**3))
217           ldss = const * G0_ML * VREL
218     
219           RETURN
220           END SUBROUTINE DRAG_SS_SYAM0
221     
222     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
223     !                                                                      !
224     !  Subroutine: DRAG_SS_IA                                              !
225     !  Purpose: Compute collisional momentum source terms betweem solids   !
226     !  phase M and solids phase L using Iddir Arastoopour (2005) kinetic   !
227     !  theory mdoel that is proportional to the relative velocity between  !
228     !  between the two phases (i.e. solids-solids drag).                   !
229     !                                                                      !
230     !  Literature/Document References:                                     !
231     !    Iddir, Y.H., "Modeling of the multiphase mixture of particles     !
232     !       using the kinetic theory approach," PhD Thesis, Illinois       !
233     !       Institute of Technology, Chicago, Illinois, 2004               !
234     !    Iddir, Y.H., & H. Arastoopour, "Modeling of Multitype particle    !
235     !      flow using the kinetic theory approach," AIChE J., Vol 51,      !
236     !      no. 6, June 2005                                                !
237     !     Gera, D., Syamlal, M., O'Brien T.J. 2004. International Journal  !
238     !        of Multiphase Flow, 30, p419-428.                             !
239     !                                                                      !
240     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
241           SUBROUTINE DRAG_SS_IA(L, M)
242     
243     ! Modules
244     !---------------------------------------------------------------------//
245           use compar, only: ijkstart3, ijkend3
246           USE constant, only: segregation_slope_coefficient
247           USE constant, only: c_e, pi
248           USE drag, only: f_ss
249           USE fldvar, only: d_p, ro_s, rop_s, theta_m, p_star
250           use functions, only: wall_at, funlm
251           USE param1, only: zero
252           USE physprop, only: close_packed
253           USE rdf, only: g_0
254           IMPLICIT NONE
255     
256     ! Dummy arguments
257     !---------------------------------------------------------------------//
258     ! solids phase indices
259           INTEGER, INTENT(IN) :: M, L
260     
261     ! Local variables
262     !---------------------------------------------------------------------//
263     ! indices
264           INTEGER :: IJK
265     ! index for storing solids-solids drag coefficients in the upper
266     ! triangle of the matrix
267           INTEGER :: LM
268     ! particle diameters of phase M and phase L
269           DOUBLE PRECISION :: D_pm, D_pl
270     ! particle densities of phase M and phase L
271           DOUBLE PRECISION :: RO_M, RO_L
272     ! radial distribution function between phases M and L
273           DOUBLE PRECISION :: G0_ML
274     ! granular temperature of phase M and phase L
275           DOUBLE PRECISION :: thetaM, thetaL
276     ! solid-solid drag coefficient
277           DOUBLE PRECISION :: lDss
278     ! Sum of particle diameters
279           DOUBLE PRECISION :: DPSUM
280     ! Intermediate calculation
281           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, DPSUMo2
282           DOUBLE PRECISION :: Ap_lm, Dp_lm, Bp_lm, R2p_lm
283           DOUBLE PRECISION :: F_common_term
284     !---------------------------------------------------------------------//
285           LM = FUNLM(L,M)
286     
287           DO IJK = ijkstart3, ijkend3
288              IF (WALL_AT(IJK)) CYCLE
289     ! Evaluate at all flow boundaries and fluid cellls. This is unlike
290     ! calculation of the fluid-solid drag coefficient, which is only
291     ! evaluated in fluid cells and pressure inflow cells
292     
293     ! setting aliases for easy reference
294              D_PM = D_P(IJK,M)
295              D_PL = D_P(IJK,L)
296              RO_M = RO_S(IJK,M)
297              RO_L = RO_S(IJK,L)
298              ThetaM = Theta_M(IJK,M)
299              ThetaL = Theta_M(IJK,L)
300              G0_ML = G_0(IJK,L,M)
301     
302              DPSUM = D_PL + D_PM
303              M_PM = (Pi/6.d0) * D_PM**3 *RO_M
304              M_PL = (Pi/6.d0) * D_PL**3 *RO_L
305              MPSUM = M_PM + M_PL
306              DPSUMo2 = DPSUM/2.d0
307              ldss = zero
308     
309     ! evaluating the solids-solids drag coefficient
310              IF(ThetaM > ZERO .AND. ThetaL > ZERO) THEN
311                 Ap_lm = (M_PM*ThetaL+M_PL*ThetaM)/2.d0
312                 Bp_lm = (M_PM*M_PL*(ThetaL-ThetaM ))/(2.d0*MPSUM)
313                 Dp_lm = (M_PL*M_PM*(M_PM*ThetaM+M_PL*ThetaL ))/&
314                     (2.d0*MPSUM*MPSUM)
315     
316                 R2p_lm = ( 1.d0/( 2.d0*Ap_lm**1.5 * Dp_lm*Dp_lm ) )+&
317                          ( (3.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**3 ) )+&
318                          ( (15.d0*Bp_lm**4)/( 2.d0*Ap_lm**3.5 * Dp_lm**4 ) )
319     
320                 F_common_term = (DPSUMo2*DPSUMo2/4.d0)*(M_PM*M_PL/MPSUM)*&
321                      G0_ML*(1.d0+C_E)*(M_PM*M_PL)**1.5
322     
323     ! Momentum source associated with relative velocity between solids
324     ! phase m and solid phase l
325     ! Factor of rop_sm and rop_sl included in calling routine
326                 ldss = F_common_term/(M_PM*M_PL)*DSQRT(PI)*R2p_lm*&
327                             (ThetaM*ThetaL)**2
328              ENDIF
329              F_SS(IJK,LM) = lDss*ROP_S(IJK,M)*ROP_S(IJK,L)
330     
331     ! Gera: accounting for particle-particle drag due to enduring contact
332     ! in a close-packed system.
333              IF(CLOSE_PACKED(M) .AND. CLOSE_PACKED(L)) &
334                 F_SS(IJK,LM) = F_SS(IJK,LM) + &
335                    SEGREGATION_SLOPE_COEFFICIENT*P_star(IJK)
336     
337           ENDDO    ! end do (ijk=ijkstart3,ijkend3)
338     
339           RETURN
340           END SUBROUTINE DRAG_SS_IA
341