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

1     !``````````````````````````````````````````````````````````````````````!
2     !  Module: rdf                                                         !
3     !                                                                      !
4     !  Calculate radial distribution functions.                            !
5     !  Note that routines G_0AVG, G_0, and DG_0DNU need to be modified to  !
6     !  effect a change in the radial distribution function g_0. The old    !
7     !  routine G_0EP has been replaced with G_0CS, which used only for     !
8     !  Carnahan-Starling g_0.                                              !
9     !......................................................................!
10     
11     MODULE rdf
12     
13     CONTAINS
14     
15     !``````````````````````````````````````````````````````````````````````!
16     !  Function: G_0AVG                                                    !
17     !  Author: M. Syamlal                                 Date: 05-JAN-05  !
18     !                                                                      !
19     !......................................................................!
20           DOUBLE PRECISION FUNCTION G_0AVG (IJK1, IJK2, DIR, L, M1, M2)
21     
22           USE compar
23           USE constant
24           use discretelement, only: des_mmax
25           USE fldvar
26           USE functions
27           USE geometry
28           USE indices
29           USE param
30           USE param1
31           USE physprop
32           USE visc_s
33           USE run
34           USE toleranc
35     
36           IMPLICIT NONE
37     
38     ! Dummy Arguments:
39     !---------------------------------------------------------------------//
40     ! Cell indices for neighboring cells (passed variable)
41           INTEGER, INTENT(IN) :: IJK1, IJK2
42     ! Direction (X, Y, or Z) (passed variable)
43           CHARACTER, INTENT(IN) :: DIR
44     ! Direction index (i, j, or k)
45           INTEGER, INTENT(IN) :: L
46     ! Solids phase index-1
47           INTEGER, INTENT(IN) :: M1
48     ! Solids phase index-2
49           INTEGER, INTENT(IN) :: M2
50     
51     ! Local Variables:
52     !---------------------------------------------------------------------//
53     ! Solids phase index
54           INTEGER :: Mx
55     ! Average solids volume fraction
56           DOUBLE PRECISION :: EPS
57     ! Average void fraction
58           DOUBLE PRECISION :: EPg, EPg_STAR_AVG
59     ! Sum over m (EP_sm/D_pm)
60           DOUBLE PRECISION :: EPSoDP
61     ! Sum over EP_s
62           DOUBLE PRECISION :: SUM_EPS
63     ! Average number density of solids phase mm
64           DOUBLE PRECISION :: NU_MM
65     ! Volume of a single particle of solids phase mm
66           DOUBLE PRECISION :: VOLP
67     ! Quantity employed in rdf calculation
68           DOUBLE PRECISION :: XI
69     ! Average D_P for phase M1 and M2
70           DOUBLE PRECISION :: DP_AVG_M1, DP_AVG_M2, DP_AVG
71     
72           IF(IJK1 == IJK2)THEN
73              G_0AVG = G_0(IJK1, M1, M2)
74           ELSE
75              SELECT CASE(RDF_TYPE_ENUM)
76     
77     ! Lebowitz, J.L. (1964) The Physical Review, A133, 895-899
78     !---------------------------------------------------------------------//
79              CASE(LEBOWITZ)
80                 DP_AVG_M1 = HALF*(D_p(IJK1,M1)+D_p(IJK2,M1))
81                 DP_AVG_M2 = HALF*(D_p(IJK1,M2)+D_p(IJK2,M2))
82                 EPSoDP = ZERO
83     
84                 DO Mx = 1, DES_MMAX+MMAX
85                    EPS = AVG_XYZ(EP_s(IJK1, Mx), EP_s(IJK2, Mx), DIR, L)
86                    EPSoDP = EPSoDP + 2d0*EPS / (D_p(IJK1,Mx)+D_p(IJK2,Mx))
87                 ENDDO
88     
89                 EPg = AVG_XYZ(EP_g(IJK1), EP_g(IJK2), DIR, L)
90                 G_0AVG = ONE / EPg + 3.0d0*EPSoDP*DP_AVG_M1*DP_AVG_M2 /    &
91                    (EPg*EPg *(DP_AVG_M1 + DP_AVG_M2))
92     
93     
94     ! REF: ?
95     !---------------------------------------------------------------------//
96              CASE(MODIFIED_LEBOWITZ)
97                 DP_AVG_M1 = HALF*(D_p(IJK1,M1)+D_p(IJK2,M1))
98                 DP_AVG_M2 = HALF*(D_p(IJK1,M2)+D_p(IJK2,M2))
99                 EPSoDP = ZERO
100                 SUM_EPS = ZERO
101     
102                 DO Mx = 1, DES_MMAX+MMAX
103                    EPS = AVG_XYZ(EP_s(IJK1, Mx), EP_s(IJK2, Mx), DIR, L)
104                    DP_AVG = HALF*( D_p(IJK1,Mx)+D_p(IJK2,Mx) )
105                    EPSoDP = EPSoDP + EPS/DP_AVG
106                    SUM_EPS = SUM_EPS + EPS
107                 ENDDO
108     
109                 EPg_STAR_AVG = AVG_XYZ(EP_star_array(IJK1),                &
110                    EP_star_array(IJK2), DIR, L)
111     
112     
113     ! Prevent G_0 from becoming negative when the solids volume fraction
114     ! approachs a maximum packing.  This may occur during non-converged
115     ! iterations when the local solids volume fraction exceeds the maximum.
116                 IF(SUM_EPS >= (ONE-EPg_STAR_AVG))                          &
117                    SUM_EPS = SUM_EPS - DIL_EP_s
118     
119     
120                 G_0AVG = (ONE/(ONE-SUM_EPS/(ONE-EPg_STAR_AVG))) + 3.0d0 *  &
121                    ((DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2))*EPSoDP
122     
123     
124     
125     ! Mansoori, GA, Carnahan N.F., Starling, K.E. Leland, T.W. (1971).
126     ! The Journal of Chemical Physics, Vol. 54:1523-1525.
127     !---------------------------------------------------------------------//
128     ! Extended Carnahan & Starling see Garzo & Dufty (1999) for details.
129              CASE(MANSOORI)
130     
131                 DP_AVG_M1 = HALF*(D_p(IJK1,M1)+D_p(IJK2,M1))
132                 DP_AVG_M2 = HALF*(D_p(IJK1,M2)+D_p(IJK2,M2))
133                 SUM_EPS = ZERO
134                 XI = ZERO
135     
136                 DO Mx = 1, DES_MMAX+MMAX
137                    EPS = AVG_XYZ(EP_s(IJK1, Mx), EP_s(IJK2, Mx), DIR, L)
138                    SUM_EPS = SUM_EPS + EPS
139     
140                    DP_AVG = HALF*( D_p(IJK1,Mx)+D_p(IJK2,Mx) )
141                    VOLP = (PI/6.0D0)*DP_AVG**3
142     
143                    IF(DP_AVG > ZERO) THEN
144                       NU_MM = EPS/VOLP
145                       XI = XI + NU_MM*DP_AVG*DP_AVG
146                    ENDIF
147                 ENDDO
148     
149                 XI = (PI/6.0D0)*XI
150     
151                 G_0AVG = (ONE/(ONE-SUM_EPS)) + (3.0D0)*&
152                    ( (DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2) )*        &
153                    ( XI/((ONE-SUM_EPS)*(ONE-SUM_EPS)) ) + (2.0D0) *        &
154                    ( (DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2) ) *       &
155                    ( (DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2) ) *       &
156                    ( (XI*XI)/((ONE-SUM_EPS)*(ONE-SUM_EPS)*(ONE-SUM_EPS)))
157     
158     
159     ! van Wachem, B.G.M., Schouten, J.C., van den Bleek, C.M., Krishna, R.
160     ! and Sinclair, J. L. (2001). AIChE Journal 47:1035–1051.
161     !---------------------------------------------------------------------//
162              CASE(MODIFIED_MANSOORI)
163     
164                 DP_AVG_M1 = HALF*(D_p(IJK1,M1)+D_p(IJK2,M1))
165                 DP_AVG_M2 = HALF*(D_p(IJK1,M2)+D_p(IJK2,M2))
166                 SUM_EPS = ZERO
167                 XI = ZERO
168     
169                 DO Mx = 1, DES_MMAX+MMAX
170                    EPS = AVG_XYZ(EP_s(IJK1, Mx), EP_s(IJK2, Mx), DIR, L)
171                    SUM_EPS = SUM_EPS + EPS
172     
173                    DP_AVG = HALF*( D_p(IJK1,Mx)+D_p(IJK2,Mx) )
174                    VOLP = (PI/6.0D0)*DP_AVG**3
175     
176                    IF (DP_AVG > ZERO) THEN
177                       NU_MM = EPS/VOLP
178                       XI = XI + NU_MM*DP_AVG*DP_AVG
179                    ENDIF
180                 ENDDO
181     
182                 XI = (PI/6.0D0)*XI
183     
184                 EPg_STAR_AVG = AVG_XYZ(EP_star_array(IJK1),                &
185                    EP_star_array(IJK2), DIR, L)
186     
187                 IF(SUM_EPS >= (ONE-EPg_STAR_AVG))                          &
188                    SUM_EPS = SUM_EPS - DIL_EP_s
189     
190                 G_0AVG = (ONE/(ONE-SUM_EPS/(ONE-EPg_STAR_AVG))) + (3.0D0)* &
191                    ( (DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2) )*        &
192                    ( XI/((ONE-SUM_EPS/(ONE-EPg_STAR_AVG))*                 &
193                    (ONE-SUM_EPS/(ONE-EPg_STAR_AVG))) ) + (2.0D0) *         &
194                    ( (DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2) ) *       &
195                    ( (DP_AVG_M1*DP_AVG_M2)/(DP_AVG_M1+DP_AVG_M2) ) *       &
196                    ( (XI*XI)/((ONE-SUM_EPS/(ONE-EPg_STAR_AVG))*            &
197                    (ONE-SUM_EPS/(ONE-EPg_STAR_AVG))*                       &
198                    (ONE-SUM_EPS/(ONE-EPg_STAR_AVG))) )
199     
200     
201     ! Carnahan, N.F. and Starling K.E. (1969).
202     ! The Journal of Chemical Physics, Vol. 51(2):635-636.
203              CASE(CARNAHAN_STARLING)
204     !---------------------------------------------------------------------//
205     
206     ! Do not use when there are more than one granular phase.
207                 EPS = AVG_XYZ(EP_s(IJK1, M1), EP_s(IJK2, M1), DIR, L)
208                 G_0AVG = G_0CS(EPS)
209              END SELECT
210     
211           ENDIF
212     
213           RETURN
214           END FUNCTION G_0AVG
215     
216     
217     !``````````````````````````````````````````````````````````````````````!
218     !  Function: G_0 (IJK, M1, M2)                                         !
219     !  Author: M. Syamlal                                 Date: 16-MAR-92  !
220     !                                                                      !
221     !  Purpose: Calculate radial distribution function at contact for a    !
222     !           mixture of spheres of different diameters                  !
223     !                                                                      !
224     !  References:                                                         !
225     !                                                                      !
226     !    Lebowitz, J.L., "Exact solution of generalized Percus-Yevick      !
227     !      equation for a mixture of hard spheres," The Physical Review,   !
228     !      A133, 895-899 (1964).                                           !
229     !                                                                      !
230     !    Iddir, Y.H., "Modeling of the multiphase mixture of particles     !
231     !         using the kinetic theory approach," PhD Thesis, Illinois     !
232     !         Institute of Technology, Chicago, Illinois, 2004:            !
233     !         chapter 2, equations 2-49 through 2-52.                      !
234     !                                                                      !
235     !    Mansoori et al. (1971)                                            !
236     !       This RDF expression is equivalent to that cited by Jenkins and !
237     !          Mancini (1987) & Garzo and Dufty (1999).                    !
238     !......................................................................!
239           DOUBLE PRECISION FUNCTION G_0 (IJK, M1, M2)
240     
241           USE compar
242           USE constant
243           use discretelement, only: des_mmax
244           USE fldvar
245           USE functions
246           USE geometry
247           USE indices
248           USE param
249           USE param1
250           USE physprop
251           USE run
252           USE toleranc
253           USE visc_s
254     
255           IMPLICIT NONE
256     
257     ! Dummy Arguments:
258     !---------------------------------------------------------------------//
259     ! Solids phase index-1 (passed variable)
260           INTEGER, INTENT(IN) :: M1
261     ! Solids phase index-2 (passed variable)
262           INTEGER, INTENT(IN) :: M2
263     ! Fluid cell index
264           INTEGER, INTENT(IN) :: IJK
265     
266     ! Local Variables:
267     !---------------------------------------------------------------------//
268     ! Local solids phase index
269           INTEGER :: Mx, MM
270     ! Average solids volume fraction
271           DOUBLE PRECISION :: EPS
272     ! Average void fraction
273           DOUBLE PRECISION :: EPg
274     ! Sum over m (EP_sm/D_pm)
275           DOUBLE PRECISION :: EPSoDP
276     ! Sum over EP_s
277           DOUBLE PRECISION :: SUM_EPS
278     ! Number density of solids phase mm
279           DOUBLE PRECISION :: NU_MM
280     ! Volume of a single particle of solids phase mm
281           DOUBLE PRECISION :: VOLP
282     ! Quantity employed in rdf calculation
283           DOUBLE PRECISION :: XI
284     !---------------------------------------------------------------------//
285     
286           SUM_EPS = ZERO
287           EPg = EP_G(IJK)
288           DO MM = 1, DES_MMAX+MMAX
289               EPS = EP_s(IJK, MM)
290               SUM_EPS = SUM_EPS + EPS
291           END DO
292     
293           SELECT CASE(RDF_TYPE_ENUM)
294     
295     ! Lebowitz, J.L. (1964) The Physical Review, A133, 895-899
296     !---------------------------------------------------------------------//
297           CASE(LEBOWITZ)
298     
299              EPSoDP = ZERO
300              DO Mx = 1, DES_MMAX+MMAX
301                 EPS = EP_s(IJK, Mx)
302                 EPSoDP = EPSoDP + EPS / D_p(IJK,Mx)
303              ENDDO
304              EPg = EP_g(IJK)
305     
306              G_0 = ONE/EPg + 3.0d0 * EPSoDP * D_p(IJK,M1) * D_p(IJK,M2) /  &
307                 (EPg*EPg *(D_p(IJK,M1) + D_p(IJK,M2)))
308     
309     ! REF: ?
310     !---------------------------------------------------------------------//
311           CASE(MODIFIED_LEBOWITZ)
312     
313              EPSoDP = ZERO
314              SUM_EPS = ZERO
315     
316              DO MM = 1, DES_MMAX+MMAX
317                 EPS = EP_s(IJK, MM)
318                 EPSoDP = EPSoDP + (EPS/D_p(IJK,MM))
319                 SUM_EPS = SUM_EPS + EPS
320              END DO
321     
322     ! Prevent G_0 from becoming negative when the solids volume fraction
323     ! approachs a maximum packing.  This may occur during non-converged
324     ! iterations when the local solids volume fraction exceeds the maximum.
325              IF(SUM_EPS >= (ONE-EP_star_array(IJK)))                       &
326                 SUM_EPS = SUM_EPS - DIL_EP_s
327     
328              G_0 = (ONE/(ONE-SUM_EPS/(ONE-EP_star_array(IJK)) )) + 3.0d0 * &
329                 ((D_p(IJK,M1)*D_p(IJK,M2))/(D_p(IJK,M1)+D_p(IJK,M2)))*EPSoDP
330     
331     ! Mansoori, GA, Carnahan N.F., Starling, K.E. Leland, T.W. (1971).
332     ! The Journal of Chemical Physics, Vol. 54:1523-1525.
333     !---------------------------------------------------------------------//
334     ! Extended Carnahan & Starling see Garzo & Dufty (1999) for details.
335           CASE(MANSOORI)
336     
337              SUM_EPS = ZERO
338              XI = ZERO
339     
340              DO MM = 1, DES_MMAX+MMAX
341                 EPS = EP_s(IJK, MM)
342                 SUM_EPS = SUM_EPS + EPS
343                 VOLP = (PI/6.0D0)*D_P(IJK,MM)**3.0
344     
345                 IF (D_P(IJK,MM) > ZERO) THEN
346                    NU_MM = EPS/VOLP
347                    XI = XI + NU_MM*D_P(IJK,MM)*D_P(IJK,MM)
348                 ENDIF
349              ENDDO
350              XI = (PI/6.0D0)*XI
351     
352              G_0 = (ONE/(ONE-SUM_EPS)) + (3.0D0)*&
353                 ( (D_P(IJK,M1)*D_P(IJK,M2))/(D_P(IJK,M1)+D_P(IJK,M2)) )*   &
354                 ( XI/((ONE-SUM_EPS)*(ONE-SUM_EPS)) ) + (2.0D0) *           &
355                 ( (D_P(IJK,M1)*D_P(IJK,M2))/(D_P(IJK,M1)+D_P(IJK,M2)) ) *  &
356                 ( (D_P(IJK,M1)*D_P(IJK,M2))/(D_P(IJK,M1)+D_P(IJK,M2)) ) *  &
357                 ( (XI*XI)/((ONE-SUM_EPS)*(ONE-SUM_EPS)*(ONE-SUM_EPS)) )
358     
359     ! van Wachem, B.G.M., Schouten, J.C., van den Bleek, C.M., Krishna, R.
360     ! and Sinclair, J. L. (2001). AIChE Journal 47:1035–1051.
361     !---------------------------------------------------------------------//
362              CASE(MODIFIED_MANSOORI)
363     
364              SUM_EPS = ZERO
365              XI = ZERO
366     
367              DO MM = 1, DES_MMAX+MMAX
368                 EPS = EP_s(IJK, MM)
369                 SUM_EPS = SUM_EPS + EPS
370                 VOLP = (PI/6.0D0)*D_P(IJK,MM)**3.0
371     
372                 IF (D_P(IJK,MM) > ZERO) THEN
373                    NU_MM = EPS/VOLP
374                    XI = XI + NU_MM*D_P(IJK,MM)*D_P(IJK,MM)
375                 ENDIF
376              ENDDO
377              XI = (PI/6.0D0)*XI
378     
379              IF(SUM_EPS >= (ONE-EP_star_array(IJK)) )                      &
380                 SUM_EPS = SUM_EPS - DIL_EP_s
381     
382              G_0 = (ONE/(ONE-SUM_EPS/(ONE-EP_star_array(IJK)))) + (3.0D0)* &
383                 ( (D_P(IJK,M1)*D_P(IJK,M2))/(D_P(IJK,M1)+D_P(IJK,M2)) )*   &
384                 ( XI/((ONE-SUM_EPS/(ONE-EP_star_array(IJK)))*              &
385                 (ONE-SUM_EPS/(ONE-EP_star_array(IJK)))) ) + (2.0D0) *      &
386                 ( (D_P(IJK,M1)*D_P(IJK,M2))/(D_P(IJK,M1)+D_P(IJK,M2)) ) *  &
387                 ( (D_P(IJK,M1)*D_P(IJK,M2))/(D_P(IJK,M1)+D_P(IJK,M2)) ) *  &
388                 ( (XI*XI)/((ONE-SUM_EPS/(ONE-EP_star_array(IJK)))*         &
389                 (ONE-SUM_EPS/(ONE-EP_star_array(IJK)))*                    &
390                 (ONE-SUM_EPS/(ONE-EP_star_array(IJK)))) )
391     
392     ! Carnahan, N.F. and Starling K.E. (1969).
393     ! The Journal of Chemical Physics, Vol. 51(2):635-636.
394           CASE(CARNAHAN_STARLING)
395     !---------------------------------------------------------------------//
396              G_0 = G_0CS(EP_S(IJK,M1))
397     
398           END SELECT
399     
400           RETURN
401           END FUNCTION G_0
402     
403     
404     !``````````````````````````````````````````````````````````````````````!
405     !  Module name: DG_0DNU (EPs)                                          !
406     !  Author: K. Agrawal                                 Date: 16-FEB-98  !
407     !                                                                      !
408     !  Purpose: Calculate derivative of radial distribution function at    !
409     !           w.r.t granular volume fraction                             !
410     !......................................................................!
411           DOUBLE PRECISION FUNCTION DG_0DNU (EPS)
412     
413     ! Global Parameters:
414     !---------------------------------------------------------------------//
415           USE param1, only: ZERO, ONE
416           USE physprop, only: MMAX
417     
418           IMPLICIT NONE
419     
420     ! Dummy Arguments:
421     !---------------------------------------------------------------------//
422           DOUBLE PRECISION, INTENT(IN) :: EPS
423     
424           DG_0DNU = ZERO
425     
426     ! Carnahan-Starling derivative of (G0) wrt EP_s:
427     ! This value is only needed for monodisper simulatoins.
428           IF(MMAX == 1) DG_0DNU = (2.5D0-EPS)/(ONE - EPS)**4
429     
430           RETURN
431           END FUNCTION DG_0DNU
432     
433     !``````````````````````````````````````````````````````````````````````!
434     !  Module name: G_0CS(EPs)                                             !
435     !                                                                      !
436     !  Purpose: Carnahan-Starling radial distribution function.            !
437     !......................................................................!
438           DOUBLE PRECISION FUNCTION G_0CS (EPS)
439     
440     ! Global Parameters:
441     !---------------------------------------------------------------------//
442           USE param1, only: ONE
443     
444           IMPLICIT NONE
445     
446     ! Dummy Arguments:
447     !---------------------------------------------------------------------//
448     ! Solids volume fraction
449           DOUBLE PRECISION :: EPS
450     
451           G_0CS = (ONE-0.5D0*EPS)/(ONE - EPS)**3
452     
453           RETURN
454           END FUNCTION G_0CS
455     
456     !``````````````````````````````````````````````````````````````````````!
457     !  Function: AVG_XYZ                                                   !
458     !                                                                      !
459     !  Purpose: Calculate the arithmetic average in X, Y, or Z direction.  !
460     !......................................................................!
461           DOUBLE PRECISION FUNCTION AVG_XYZ (V1, V2, DIR, L)
462     
463           USE param
464           USE param1
465           USE geometry
466           USE indices
467           USE compar
468           USE functions
469           USE fun_avg
470     
471           IMPLICIT NONE
472     
473     ! Dummy Arguments:
474     !---------------------------------------------------------------------//
475     ! Direction (X, Y, or Z) (passed variable)
476           CHARACTER, INTENT(IN) :: DIR
477     ! Direction index (i, j, or k)
478           INTEGER, INTENT(IN) :: L
479     ! The variables to be averaged
480           DOUBLE PRECISION, INTENT(IN) :: V1, V2
481     
482     
483     ! Local Variables:
484     !---------------------------------------------------------------------//
485     ! Fluid cell index
486     
487           IF(DIR == 'X')THEN
488             AVG_XYZ = AVG_X(V1, V2, L)
489     
490           ELSEIF(DIR == 'Y')THEN
491             AVG_XYZ = AVG_Y(V1, V2, L)
492     
493           ELSEIF(DIR == 'Z')THEN
494             AVG_XYZ = AVG_Z(V1, V2, L)
495     
496           ELSE
497             CALL WRITE_ERROR('AVG_XYZ', 'Unkown direction', 1)
498           ENDIF
499     
500           RETURN
501           END FUNCTION AVG_XYZ
502     
503         END MODULE rdf
504