MFIX  2016-1
rdf_mod.f
Go to the documentation of this file.
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
visc_s
Definition: visc_s_mod.f:1
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
param1::one
double precision, parameter one
Definition: param1_mod.f:29
write_error
subroutine write_error(NAME, LINE, LMAX)
Definition: write_error.f:21
rdf::dg_0dnu
double precision function dg_0dnu(EPS)
Definition: rdf_mod.f:412
rdf::g_0
double precision function g_0(IJK, M1, M2)
Definition: rdf_mod.f:240
visc_s::ep_star_array
double precision, dimension(:), allocatable ep_star_array
Definition: visc_s_mod.f:54
fun_avg
Definition: fun_avg_mod.f:1
fldvar::d_p
double precision, dimension(:,:), allocatable d_p
Definition: fldvar_mod.f:57
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
rdf::avg_xyz
double precision function avg_xyz(V1, V2, DIR, L)
Definition: rdf_mod.f:462
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
param
Definition: param_mod.f:2
toleranc::dil_ep_s
double precision, parameter dil_ep_s
Definition: toleranc_mod.f:24
fldvar::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
rdf
Definition: rdf_mod.f:11
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
rdf::g_0avg
double precision function g_0avg(IJK1, IJK2, DIR, L, M1, M2)
Definition: rdf_mod.f:21
rdf::g_0cs
double precision function g_0cs(EPS)
Definition: rdf_mod.f:439
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27