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