File: /nfs/home/0/users/jenkins/mfix.git/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     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !     Gera, D., Syamlal, M., O'Brien T.J. 2004. International Journal  C
12     !        of Multiphase Flow, 30, p419-428.                             C
13     !                                                                      C
14     !  Variables referenced:                                               C
15     !  Variables modified:                                                 C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20           SUBROUTINE DRAG_SS(L, M)
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE compar
26           USE constant
27           USE discretelement
28           USE drag
29           USE fldvar
30           USE fun_avg
31           USE functions
32           USE geometry
33           USE indices
34           USE parallel
35           USE param
36           USE param1
37           USE physprop
38           USE rdf
39           USE sendrecv
40           IMPLICIT NONE
41     !-----------------------------------------------
42     ! Dummy Arguments
43     !-----------------------------------------------
44     ! Index of solids phases
45           INTEGER, INTENT(IN) :: L, M
46     !-----------------------------------------------
47     ! Local Variables
48     !-----------------------------------------------
49     ! indices
50           INTEGER :: I, IJK, IMJK, IJMK, IJKM
51           INTEGER :: CM, DM
52     ! index for storing solids-solids drag coefficients in the upper
53     ! triangle of the matrix
54           INTEGER :: LM
55     ! cell center value of U_sm, U_sl, V_sm, V_sl, W_sm, W_sl
56           DOUBLE PRECISION :: USCM, USCL, VSCM, VSCL, WSCM, WSCL
57     ! relative velocity between solids phase m and l
58           DOUBLE PRECISION :: VREL
59     ! particle diameters of phase M and phase L
60           DOUBLE PRECISION :: D_pm, D_pl
61     ! particle densities of phase M and phase L
62           DOUBLE PRECISION :: RO_M, RO_L
63     ! radial distribution function between phases M and L
64           DOUBLE PRECISION :: G0_ML
65     ! void fraction and solids volume fraction
66           DOUBLE PRECISION :: EPg, EPs
67     ! sum over all phases of ratio volume fraction over particle diameter
68           DOUBLE PRECISION :: EPSoDP
69     ! solid-solid drag coefficient
70           DOUBLE PRECISION :: lDss
71     !-----------------------------------------------
72     
73           LM = FUNLM(L,M)
74     
75           DO IJK = ijkstart3, ijkend3
76     
77              IF (.NOT.WALL_AT(IJK)) THEN
78     ! Evaluate at all flow boundaries and fluid cells
79     ! This is unlike the fluid-solid drag coefficient, which is only
80     ! evluated in fluid cells and pressure inflow cells
81     
82                 I = I_OF(IJK)
83                 IMJK = IM_OF(IJK)
84                 IJMK = JM_OF(IJK)
85                 IJKM = KM_OF(IJK)
86     
87     ! calculating velocity components at i, j, k (cell center)
88                 USCL = AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I)
89                 VSCL = AVG_Y_N(V_S(IJMK,L),V_S(IJK,L))
90                 WSCL = AVG_Z_T(W_S(IJKM,L),W_S(IJK,L))
91     
92                 USCM = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
93                 VSCM = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M))
94                 WSCM = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M))
95     
96     ! magnitude of solids-solids relative velocity
97                 VREL = SQRT((USCL - USCM)**2 + (VSCL - VSCM)**2 + &
98                             (WSCL - WSCM)**2)
99     
100     ! setting aliases for easy reference
101                 D_PM = D_P(IJK,M)
102                 D_PL = D_P(IJK,L)
103                 RO_M = RO_S(IJK,M)
104                 RO_L = RO_S(IJK,L)
105     
106                 IF (DES_CONTINUUM_HYBRID) THEN
107     ! evaluating g0 - taken from G_0.f subroutine (lebowitz form)
108     ! this section is needed to account for all solids phases until g0 for
109     ! multiple solids types (i.e. discrete & continuum) can be addressed
110     ! more effectively.
111                    EPSoDP = ZERO
112                    DO CM = 1, MMAX
113                       EPS = EP_s(IJK, CM)
114                       EPSoDP = EPSoDP + EPS / D_p(IJK,CM)
115                    ENDDO
116                    DO DM = 1, DES_MMAX
117                       EPS = DES_ROP_S(IJK,DM)/DES_RO_S(DM)
118                       EPSoDP = EPSoDP + EPS / DES_D_p0(DM)
119                    ENDDO
120                    EPg = EP_g(IJK)
121                    G0_ML = ONE/EPg + 3.0d0*EPSoDP*D_pM*D_PL / &
122                       (EPg*EPg *(D_pM + D_pL))
123                 ELSE
124                    G0_ML = G_0(IJK,L,M)
125                 ENDIF
126     
127     ! determining the solids-solids 'drag coefficient'
128                 CALL DRAG_SS_SYAM(lDss,D_PM,D_PL,RO_M,RO_L,G0_ML,VREL)
129     
130                 F_SS(IJK,LM) = lDss*ROP_S(IJK,M)*ROP_S(IJK,L)
131     
132     ! Gera: accounting for particle-particle drag due to enduring contact in a
133     ! close-packed system. Note the check for mmax >= 2 below is unnecessary
134     ! since this routine will not be entered unless mmax >=2
135                 IF(CLOSE_PACKED(M) .AND. CLOSE_PACKED(L) .AND. &
136                   (MMAX >= 2))  F_SS(IJK,LM) = F_SS(IJK,LM) + &
137                    SEGREGATION_SLOPE_COEFFICIENT*P_star(IJK)
138     
139               ELSE   ! elseif (.not.wall_at(ijk))
140     
141                 F_SS(IJK,LM) = ZERO
142     
143              ENDIF   ! end if (.not.wall_at(ijk))
144     
145           ENDDO    ! end do (ijk=ijkstart3,ijkend3)
146     
147           RETURN
148           END SUBROUTINE DRAG_SS
149     
150     
151     
152     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
153     !                                                                      C
154     !  Subroutine: DRAG_SS_SYAM                                            C
155     !  Purpose: Calculate the solids-solids drag coefficient between a     C
156     !           continuous solids phase and discrete solids                C
157     !                                                                      C
158     !  Literature/Document References:                                     C
159     !     M. Syamlal. 1987. The particle-particle drag term in a           C
160     !        multiparticle model of fluidization. Technical Report.        C
161     !        DOE/MC/21353-2373. Office of Fossil Energy, Morgantown        C
162     !        Energy Technology Center, Morgantown, West Virginia.          C
163     !                                                                      C
164     !                                                                      C
165     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
166     
167           SUBROUTINE DRAG_SS_SYAM(lDss,D_PM,D_PL,RO_M,RO_L, G0_ML, VREL)
168     
169     !-----------------------------------------------
170     ! Modules
171     !-----------------------------------------------
172           USE param
173           USE param1
174           USE constant
175           use run
176           IMPLICIT NONE
177     !-----------------------------------------------
178     ! Dummy arguments
179     !-----------------------------------------------
180     ! drag coefficient
181           DOUBLE PRECISION, INTENT(OUT) :: ldss
182     ! particle diameter of solids phase M
183           DOUBLE PRECISION, INTENT(IN) :: D_PM
184     ! particle diameter of solids phase L
185           DOUBLE PRECISION, INTENT(IN) :: D_PL
186     ! particle density of solids phase M
187           DOUBLE PRECISION, INTENT(IN) :: RO_M
188     ! particle density of solids phase L
189           DOUBLE PRECISION, INTENT(IN) :: RO_L
190     ! radial distribution function between phases M and L
191           DOUBLE PRECISION, INTENT(IN) :: G0_ML
192     ! relative velocity between solids phase m and l
193           DOUBLE PRECISION, INTENT(IN) :: VREL
194     !-----------------------------------------------
195     ! Local variables
196     !-----------------------------------------------
197     ! Sum of particle diameters
198           DOUBLE PRECISION :: DPSUM
199     ! Intermediate calculation
200           DOUBLE PRECISION :: const
201     !-----------------------------------------------
202     ! External functions
203     !-----------------------------------------------
204     !-----------------------------------------------
205     
206           DPSUM = D_PL + D_PM
207     
208           const = 3.d0*(ONE + C_E)*(PI/2.d0 + C_F*PI*PI/8.d0)*&
209              DPSUM**2/(2.d0*PI*(RO_L*D_PL**3+RO_M*D_PM**3))
210     
211           ldss = const * G0_ML * VREL
212     
213           RETURN
214           END SUBROUTINE DRAG_SS_SYAM
215     
216     
217     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
218     !                                                                      !
219     !  Subroutine: DRAG_ss_IA                                              !
220     !  Purpose: Compute the coefficient of drag between solids phase m     !
221     !           and solids phase l using Iddir Arastoopour (2005) kinetic  !
222     !           theory model                                               !
223     !                                                                      !
224     !  Literature/Document References:                                     !
225     !    Iddir, Y.H., "Modeling of the multiphase mixture of particles     !
226     !       using the kinetic theory approach," PhD Thesis, Illinois       !
227     !       Institute of Technology, Chicago, Illinois, 2004               !
228     !    Iddir, Y.H., & H. Arastoopour, "Modeling of Multitype particle    !
229     !      flow using the kinetic theory approach," AIChE J., Vol 51,      !
230     !      no. 6, June 2005                                                !
231     !                                                                      !
232     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
233     
234           SUBROUTINE DRAG_SS_IA(L, M)
235     
236     !-----------------------------------------------
237     ! Modules
238     !-----------------------------------------------
239           USE compar
240           USE constant
241           USE drag
242           USE fldvar
243           USE functions
244           USE geometry
245           USE indices
246           USE kintheory
247           USE param1
248           USE physprop
249           USE rdf
250           USE sendrecv
251     
252           IMPLICIT NONE
253     !-----------------------------------------------
254     ! Dummy arguments
255     !-----------------------------------------------
256     ! Solids phase index
257           INTEGER, INTENT(IN) :: M, L
258     !-----------------------------------------------
259     ! Local variables
260     !-----------------------------------------------
261     ! Indices
262           INTEGER :: IJK
263     ! Index for storing solids-solids drag coefficients
264     ! in the upper triangle of the matrix
265           INTEGER :: LM
266     ! Particle diameters
267           DOUBLE PRECISION :: D_PM, D_PL
268     ! Sum of particle diameters
269           DOUBLE PRECISION :: DPSUM
270     !
271           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, DPSUMo2, NU_PL, NU_PM
272           DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R2p_lm, R3p_lm, &
273                               R4p_lm, R10p_lm, Bp_lm
274           DOUBLE PRECISION :: Fss_ip, Fnus_ip, FTsM_ip, FTsL_ip, &
275                               F_common_term
276     !-----------------------------------------------
277     
278           DO IJK = ijkstart3, ijkend3
279     
280               IF (.NOT.WALL_AT(IJK)) THEN
281     
282                    LM = FUNLM(L,M)
283     
284                    IF (M == L) THEN
285                         F_SS(IJK,LM) = ZERO
286                         Fnu_s_ip(IJK,M,L) = ZERO
287                         FT_sM_ip(IJK,M,L) = ZERO
288                         FT_sL_ip(IJK,M,L) = ZERO
289     
290                    ELSE
291                         D_PM = D_P(IJK,M)
292                         D_PL = D_P(IJK,L)
293                         DPSUM = D_PL + D_PM
294                         M_PM = (Pi/6.d0) * D_PM**3 *RO_S(IJK,M)
295                         M_PL = (Pi/6.d0) * D_PL**3 *RO_S(IJK,L)
296                         MPSUM = M_PM + M_PL
297                         DPSUMo2 = DPSUM/2.d0
298                         NU_PM = ROP_S(IJK,M)/M_PM
299                         NU_PL = ROP_S(IJK,L)/M_PL
300     
301                      IF(Theta_m(IJK,M) > ZERO .AND. Theta_m(IJK,L) > ZERO) THEN
302     
303                         Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/&
304                               2.d0
305                         Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-Theta_m(IJK,M) ))/&
306                              (2.d0*MPSUM)
307                         Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+M_PL*Theta_m(IJK,L) ))/&
308                              (2.d0*MPSUM*MPSUM)
309     
310                         R0p_lm = ( 1.d0/( Ap_lm**1.5 * Dp_lm**2.5 ) )+ &
311                                   ( (15.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**2.5 * Dp_lm**3.5 ) )+&
312                                   ( (175.d0*(Bp_lm**4))/( 8.d0*Ap_lm**3.5 * Dp_lm**4.5 ) )
313     
314                         R2p_lm = ( 1.d0/( 2.d0*Ap_lm**1.5 * Dp_lm*Dp_lm ) )+&
315                                   ( (3.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**3 ) )+&
316                                   ( (15.d0*Bp_lm**4)/( 2.d0*Ap_lm**3.5 * Dp_lm**4 ) )
317     
318                         R3p_lm = ( 1.d0/( (Ap_lm**1.5)*(Dp_lm**3.5) ) )+&
319                                   ( (21.d0*Bp_lm*Bp_lm)/( 2.d0 * Ap_lm**2.5 * Dp_lm**4.5 ) )+&
320                                   ( (315.d0*Bp_lm**4)/( 8.d0 * Ap_lm**3.5 *Dp_lm**5.5 ) )
321     
322                         R4p_lm = ( 3.d0/( Ap_lm**2.5 * Dp_lm**3.5 ) )+&
323                                   ( (35.d0*Bp_lm*Bp_lm)/( 2.d0 * Ap_lm**3.5 * Dp_lm**4.5 ) )+&
324                                   ( (441.d0*Bp_lm**4)/( 8.d0 * Ap_lm**4.5 * Dp_lm**5.5 ) )
325     
326                         R10p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**2.5 ) )+&
327                                   ( (25.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**3.5 * Dp_lm**3.5 ) )+&
328                                   ( (1225.d0*Bp_lm**4)/( 24.d0* Ap_lm**4.5 * Dp_lm**4.5 ) )
329     
330                         F_common_term = (DPSUMo2*DPSUMo2/4.d0)*(M_PM*M_PL/MPSUM)*&
331                              G_0(IJK,M,L)*(1.d0+C_E)*(M_PM*M_PL)**1.5
332     
333     ! Momentum source associated with relative velocity between solids
334     ! phase m and solid phase l
335                         Fss_ip = F_common_term*NU_PM*NU_PL*DSQRT(PI)*R2p_lm*&
336                                (Theta_m(IJK,M)*Theta_m(IJK,L))**2
337     
338     ! Momentum source associated with the difference in the gradients in
339     ! number density of solids phase m and all other solids phases
340                         Fnus_ip = F_common_term*(PI*DPSUMo2/12.d0)*R0p_lm*&
341                                (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
342     
343     ! Momentum source associated with the gradient in granular temperature
344     ! of solid phase M
345                         FTsM_ip = F_common_term*NU_PM*NU_PL*DPSUMo2*PI*&
346                                   (Theta_m(IJK,M)**1.5 * Theta_m(IJK,L)**2.5) *&
347                                   (  (-1.5d0/12.d0*R0p_lm)+&
348                                   Theta_m(IJK,L)/16.d0*(  (-M_PM*R10p_lm) - &
349                                   ((5.d0*M_PL*M_PL*M_PM/(192.d0*MPSUM*MPSUM))*R3p_lm)+&
350                                   ((5.d0*M_PM*M_PL)/(96.d0*MPSUM)*R4p_lm*Bp_lm)  )  )
351     
352     ! Momentum source associated with the gradient in granular temperature
353     ! of solid phase L ! no need to recompute (sof Aug 30 2006)
354                         FTsL_ip = F_common_term*NU_PM*NU_PL*DPSUMo2*PI*&
355                                  (Theta_m(IJK,L)**1.5 * Theta_m(IJK,M)**2.5) *&
356                                   (  (1.5d0/12.d0*R0p_lm)+&
357                                   Theta_m(IJK,M)/16.d0*(  (M_PL*R10p_lm)+&
358                                   (5.d0*M_PM*M_PM*M_PL/(192.d0*MPSUM*MPSUM)*R3p_lm)+&
359                                   (5.d0*M_PM*M_PL/(96.d0*MPSUM) *R4p_lm*Bp_lm)  )  )
360     
361                         F_SS(IJK,LM) = Fss_ip
362     
363                         Fnu_s_ip(IJK,M,L) = Fnus_ip
364     
365     ! WARNING: the following two terms have caused some convergence problems
366     ! earlier. Set them to ZERO for debugging in case of converegence
367     ! issues. (sof)
368                         FT_sM_ip(IJK,M,L) = FTsM_ip  ! ZERO
369     
370                         FT_sL_ip(IJK,M,L) = FTsL_ip  ! ZERO
371                      ELSE
372                        F_SS(IJK,LM) = ZERO
373                         Fnu_s_ip(IJK,M,L) = ZERO
374                         FT_sM_ip(IJK,M,L) = ZERO
375                         FT_sL_ip(IJK,M,L) = ZERO
376                      ENDIF
377                    ENDIF
378               ENDIF
379           ENDDO
380     
381           RETURN
382           END SUBROUTINE DRAG_SS_IA
383     
384