3490
Ind. Eng. Chem. Res. 1998, 37, 3490-3496
A Theoretical Correction of the Ouchiyama and Tanaka Formula for Predicting Average Porosity of Packed Beds Consisting of Nonuniform Spheres M. Song, Karl T. Chuang,* and K. Nandakumar Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6
To predict the overall average porosity of a packing of unequal spheres, Ouchiyama and Tanaka (Ind. Eng. Chem. Fundam. 1981, 20, 66.) established a simple packing model and derived an analytical formula. In this model two uniformity assumptions were introduced, which lead to systematic errors in the model, making the prediction of the average porosity higher or the prediction of the packing density lower than experimental data. This paper proposes a theoretical correction for this model. For binary packings of spheres, the correction parameter is fixed. The prediction by our correction is in excellent agreement with experimental data for a large range of j0, the porosity of a packing of equal spheres, and k, the ratio of the smaller to the larger diameter of spheres. Introduction Establishing theoretical and numerical models to predict the characteristics of a random packing made of unequal spheres is an important research subject in particulate systems (Dodds, 1980; Powell, 1980; Ouchiyama and Tanaka, 1980, 1981, 1984, 1986, 1988, 1989; Gyenis et al., 1982; Suzuki and Oshima, 1985; Suzuki et al., 1986; Yu and Standish, 1987, 1988, 1991; Nolan and Kavanagh, 1992, 1993, 1994). Although a lot of numerical models are available in the literature, Ouchiyama and Tanaka’s (1981) model is the sole analytical model from which analytical formulas are derived to express the relative coordination number (Ouchiyama and Tanaka, 1980) and the overall average porosity (Ouchiyama and Tanaka, 1981) in a packing mixture of unequal spheres. Since they published their first two papers (Ouchiyama and Tanaka, 1980, 1981), many investigators, including themselves, have applied or extended this model to predict the average porosity, or its equivalent packing density under different packing conditions in different engineering areas (Ouchiyama and Tanaka, 1984, 1986, 1988, 1989; Cross et al., 1985; Gupta and Seshadri, 1986; Poslinski et al., 1988; Fiske et al., 1994; Mayadunne et al., 1996). Results from the model of Ouchiyama and Tanaka (1981) are in reasonable agreement with experimental data. Nevertheless, usually the prediction of the average porosity is higher, or the prediction of the packing density is lower than experimental data. It seems that there are some systematic errors in the model. The objective of this paper is to correct these systematic errors so that the prediction using this model can better match experimental data. * To whom correspondence should be addressed. Phone: 403-492-4676. Fax: 403-492-2881. E-mail: KarlT.Chuang@ Ualberta.ca.
Figure 1. Ouchiyama and Tanaka’s packing model (1981).
Ouchiyama and Tanaka’s Model and Its Systematic Errors Figure 1 schematically shows the packing model proposed by Ouchiyama and Tanaka (1981). The first assumption introduced by them is that an arbitrary sphere of diameter D in a packing is in direct contact with its neighbors having the average diameter D h . We call this assumption the first uniformity assumption. Consider the hypothetical sphere of diameter D + D h, which encloses the sphere of diameter D in its core. Define the hypothetical spherical shell whose inner and outer diameters are D ∼ D h and D + D h , respectively, where
D∼D h )
{
0 D-D h
DeD h D>D h
(1)
The hypothetical sphere shares the space of the abovedefined shell in common with the other hypothetical spheres around the contacting neighbors. To get the analytical formulas, Ouchiyama and Tanaka (1981)
S0888-5885(98)00036-0 CCC: $15.00 © 1998 American Chemical Society Published on Web 07/17/1998
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3491
where
y+ )
x+ ) xD h 2 - y+2
(4)
(D + D h )2 + 4D h 2 - D2 4(D + D h)
(5)
Simplifying eqs 2 and 3, we have
∆VA )
Figure 2. Contact of the central sphere with a neighbor.
introduced the second assumption that the part of the volume shared by n hypothetical spheres is equally allocated to each sphere, independent of the size of these spheres. We call this assumption the second uniformity assumption. The two assumptions cause model systematic errors. The uniformization makes the overall average porosity higher, or the packing density lower as the packing bed of equal spheres has a larger average porosity or a smaller packing density than that of unequal spheres. The larger the nonuniformity, the smaller the average. That is why the prediction of the average porosity using this model is higher. In accordance with the above analysis, this work presents a theoretical correction for the Ouchiyama and Tanaka (1981) model.
πD h3 (20D - 7D h) 96(D + D h)
(6)
πD h3 (8D + 5D h) 96(D + D h)
(7)
∆VC )
According to our consideration, the part of ∆Vn allocated to the central sphere should be
g(∆VA) ∆VnA ) ∆Vn g(∆VA) + (n - 1)g(∆VC)
(8)
The total volume VC(D) of the space allocated to the specified sphere of diameter D can be expressed as
g(∆VA) (D ∼ D h )3 + + ∆V1 + ∆V2 6 g(∆VA) + g(∆VC) g(∆VA) + ... (9) ∆V3 g(∆VA) + 2g(∆VC)
VC(D) ≡ π Derivation of Correction Formulas Examine the hypothetical spherical shell around the central sphere. Let the volume of the space in common occupied by n hypothetical spheres be ∆Vn. In the Ouchiyama and Tanaka (1981) model, ∆Vn/n is assumed to be equally allocated to each sphere. Actually, the larger portion of ∆Vn should be allocated to a larger sphere. Figure 2 shows the contact of the central sphere, named sphere A, with a neighbor, named sphere C. We assume that the hypothetical spheres around sphere A and sphere C are marked as sphere B and sphere E, respectively. ∆Vn is shared by the central sphere (sphere A) and (n - 1) contacting spheres (like sphere C). We assume that the portion of ∆Vn allocated to a sphere is proportional to some monotonically increasing function of ∆V, the volume of the part of the sharing sphere, cut by the hypothetical sphere around its contacting sphere. We define this function as g(∆V). Let the volume of the part of sphere A, cut by sphere E, be ∆VA, and the volume of the part of sphere C, cut by sphere B, be ∆VC. The following formulas can be readily derived:
{
)}
2 D3 D+D h 3 ∆VA ) π D h3 + - y+3 + y+ 3 8 2 1 π(D + D h )x+2 (2) 2
(
∆VC )
h 1 3 D πD h3 1 12 8D+D h
(
)
(3)
This can be written as
(D ∼ D h )3
VC(D) ≡ π
∞
+
6
∑ ∆Vng(∆V n)1
A)
g(∆VA) + (n - 1)g(∆VC) (10)
Define n j so that ∞
g(∆VA) ) ∆Vn g(∆VA) + (n - 1)g(∆VC) n)1 g(∆VA) Vm(D) (11) g(∆VA) + (n j - 1)g(∆VC)
∑
where Vm(D) is the total volume of the hypothetical spherical shell, that is ∞
Vm(D) )
∑ ∆Vn n)1
π h )3 - (D ∼ D ) [(D + D h )3] 6
(12)
3492 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998
n j is the average value of n and can be expressed by
1 ∞
∑ n)1 n j)1+
[
-1
∆Vn(D)
]
(13)
The total volume of the packing bed should be
(14)
where f(D) is the frequency distribution function of spheres based on number and N is the total number of spheres in the packing. Therefore, the average porosity is
∫0
∫0 6
[
∞
π (D ∼ D h )3 + 6
Vm(D) g(∆VC)
1 + (n j - 1)
g(∆VA)
]
(15)
Nf(D) dD
∫0
[
∞
∫0∞D3f(D) dD
(D ∼ D h )3 +
]
(D + D h )3 - (D ∼ D h )3 f(D) dD g(∆VC) 1 + (n j - 1) g(∆VA) (16)
Let the local porosity within the hypothetical shell around the central sphere be m(D). Clearly,
π 3 h )3] + C(D)∆VC [D - (D ∼ D 6 (17) m(D) ) 1 Vm(D) where
D+D h 2 C(D) ) 16(1 - A) 2D
(
)
(18)
is the coordination number of the central sphere and A is the surface porosity. A can be evaluated from the volume porosity j0 of a packing of equal spheres (Ouchiyama and Tanaka, 1981):
16j0 - 1 A ) 13
(21)
∫
∞
4(7 - 8j0) 3 D h D h (D + D h )2 1 13 8D+D h × g(∆VC) 1 + (n j - 1) g(∆VC)
(
[D3 - (D ∼ D h ) 3] +
0
f(D) dD )
∫ [D ∞
3
0
)
- (D ∼ D h )3]f(D) dD (22)
From eq 22, we can solve n j . Substituting n j into eq 16, we can obtain the average porosity j for a packing characterized by the frequency distribution function f(D) of sphere sizes based on number. j is solely determined by j0 and f(D). Calculations and Results
or
j ) 1 -
From the volume balance,
Substituting eqs 17-20 into eq 21, we obtain
∫0∞VC(D)Nf(D) dD
D3Nf(D) dD
g(∆VA)
∫0∞{π6(D ∼ Dh )3 + [1 - m(D)]∆VS}Nf(D) dD ) ∫0∞π6D3Nf(D) dD
g(∆VA)
∞π
(20)
g(∆VC)
1 + (n j - 1)
g(∆VC)
j ) 1 -
Vm(D)
∆VS )
g(∆VC)
Vm(D) 1 + (n - 1) g(∆VA)
VB )
Because the spherical shell around the central sphere is shared on average by (n j - 1) spherical shells, the volume allocated to the central shell should be
(19)
In eq 22, n j cannot be expressed in an explicit form. In our calculation n j is solved with the help of an iteration method. We design an iteration process leading to a fast convergence as follows. Assume n1 to be the value of n j from the model of Ouchiyama and Tanaka (1981), that is
n1 ) 1 +
∫0∞(D + Dh )2(1 - 83 D +Dh Dh )f(D) dD ∫0∞[D3 - (D ∼ Dh )3]f(D) dD
4(7 - 8j0) D h 13
(23)
Approximating eq 22, we have
[∫ { ∞
0
4(7 - 8j0) D h (D + D h )2 1 13 ∞g(∆VC) 1 - (n j - 1) 0 f(D) dD ) g(∆VA)
(
[D3 - (D ∼ D h )3] +
)}
]/[
h 3 D f(D) dD 8D+D h
∫
∫0∞[D3 - (D ∼ Dh )3]f(D) dD
]
(24)
Assume n2 to be the solution of n j from eq 22, that is,
n2 ) 1 +
n1 - 1 ∞g(∆VC) f(D) dD 0 g(∆V ) A
∫
(25)
Define F(n j ) as the difference of both sides of eq 22. Letting n j ) n1 and n j ) n2, calculate F(n j 1) and F(n j 2). If j 2) is less than δ, a specified small quantity F(n j 1) or F(n
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3493
Figure 3. Comparison of the average porosity predicted by our model with Ouchiyama and Tanaka’s (1981) model and experiments for j0 ) 0.4 and various k values.
(e.g., 10-4), n j ) n1 or n j ) n2 can be considered as the solution of n j . If not, calculate
n3 ) n j 1 - (n j2 - n j 1)
F(n j 1) F(n j 2) - F(n j 1)
Figure 4. Comparison of the average porosity predicted by our model with Ouchiyama and Tanaka’s (1981) model and experiments for j0 ) 0.5 and various k values.
(26)
j 2 ) n3. If F(n j 2) < δ, n j ) n3. Let n j 1 ) n2 and n Otherwise, continue until F(n j 2) < δ is satisfied. Taking δ ) 10-4, only three to five iterations are needed. Since the frequency distribution function fw(D) of sphere sizes based on weight (or volume) is often known, the frequency distribution function f(D) based on number can be calculated by
fw(D) f(D) )
D3 ∞fw(D)
∫0
D3
(27) dD
To ensure that the local porosity of both sides at the contacting point is equal, we take
g(∆V) ) ∆V
(28)
in equations (16) and (22). For the binary packing system studied by Furnas (1928, 1929), the average porosity is a function of j0, the porosity of a packing of equal size spheres, and k, the ratio of the smaller to the larger diameter of spheres. Taking j0 ) 0.4, 0.5, and 0.6, the calculation results are shown for k ) 0.5, 0.4, 0.3, 0.2, and 0.1 in Figures 3a-e, 4a-e, and 5a-e, respectively. The experimental data from Furnas (1928, 1929) and the results from the Ouchiyama and Tanaka formula (1981) are plotted for comparison. Obviously, the prediction of the average porosity is improved to a certain extent by our correction formulas (16) and (22). This shows that our view on the systematic errors caused by the assumptions in the Ouchiyama and Tanaka model is correct. However, our correction is not good enough for the theoretical results to be in satisfactory agreement with experimental data because we just correct the systematic error caused by
Figure 5. Comparison of the average porosity predicted by our model with Ouchiyama and Tanaka’s (1981) model and experiments for j0 ) 0.6 and various k values.
the second uniformity assumption. The first one also creates a systematic error leading to a higher average porosity or a lower packing density. Try the exponential function
g(∆V) ) ∆Vm
(29)
The model of Ouchiyama and Tanaka can be viewed as if m ) 0 is taken. m ) 1 is not sufficient to correct the error well. Therefore, we let m be larger than 1 so as to compensate for the model systematic error caused by the first uniformity assumption. Through numerical calculations, we get
m)
10k 210j0-5
(30)
3494 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 Table 1. Comparison of the Prediction Error between Our Correction and Ouchiyama and Tanaka’s Model
a
case
Ouchiyama and Tanaka’s model
our correction with m ) 10k/210j0-5
error ratio of Ouchiyama and Tanaka’s model to ours
j0 ) 0.4, k ) 0.5 j0 ) 0.4, k ) 0.4 j0 ) 0.4, k ) 0.3 j0 ) 0.4, k ) 0.2 j0 ) 0.4, k ) 0.1
0.0198 0.0416 0.0409 0.0202 0.00789
0.004 74 0.012 4 0.007 89 0.007 93 0.003 32
4.19 3.36 5.18 2.55 2.38
j0 ) 0.4, k ) 0.5 j0 ) 0.4, k ) 0.4 j0 ) 0.4, k ) 0.3 j0 ) 0.4, k ) 0.2 j0 ) 0.4, k ) 0.1
0.00738 0.0307 0.0302 0.0187 0.0114
0.017 0 0.004 82 0.003 88 0.001 27 0.005 15
0.434a 6.36 7.78 14.8 2.20
j0 ) 0.4, k ) 0.5 j0 ) 0.4, k ) 0.4 j0 ) 0.4, k ) 0.3 j0 ) 0.4, k ) 0.2 j0 ) 0.4, k ) 0.1
0.00963 0.0163 0.0191 0.0119 0.00745
0.006 79 0.003 63 0.005 59 0.005 73 0.004 87
1.42 4.50 3.42 2.08 1.53
This is the only case in which our prediction is worse than Ouchiyama and Tanaka’s.
Using eq 30, the calculated results are shown in Figures 3-5. Except for the case of j0 ) 0.5 and k ) 0.5, the agreement between our correction and the experimental data is fairly good. To quantify the improvement of prediction using our model, we define the relative average error as
x∑ P
E)
(p,exp - p,pre)2/P
p)1
j0
(31)
where P is the number of experimental data and p is obtained at wp ) 0.5(p - 1), p ) 1, 2, ..., P. The calculated results are listed in Table 1. The last column shows the ratio of error in Ouchiyama and Tanaka’s (1981) model to that in our model. Values greater than unity indicate an improvement in the prediction. It is clear that our model is superior in all except one case. For the cases of k ) 0.0 at all j0 results from the Ouchiyama and Tanaka (1981) formula, our corrections with m ) 1 and m expressed by eq 30 are all in excellent agreement with the experimental data. In the calculations, k ) 0.001 is set for k ) 0. Further Discussions Our correction for the Ouchiyama and Tanaka formula gives an extremely close agreement with Furnas’ experimental data for all cases except one (j0 ) 0.5, k ) 0.5). We suspect that there might be an error in this one set of experimental data. This belief is based on the expectation that changes in porosity should conform to a pattern of uniform changes in value when k and j0 are changed continuously. Plotting the experimental porosity versus k at different weight fractions of larger spheres, we get Figure 6a-c for j0 ) 0.4, 0.5, and 0.6, respectively. From these data we compute the local slope to check the trend using the slope defined as
R ) tan-1
i - i-1 ki - k-1
(32)
The results are shown in Figure 7. From this figure we can see a big jump in the slope in the range of k ) 0.4 to k ) 0.5 for j0 ) 0.5 because of
Figure 6. Experimental porosity vs k. (a) j0 ) 0.4, (b) j0 ) 0.5, and (c) j0 ) 0.6.
possible error in the experimental data for the case of j0 ) 0.5 k ) 0.5. Conclusions On the basis of our analysis of the systematic errors caused by the two uniformity assumptions introduced by Ouchiyama and Tanaka in establishing their analytical model, this investigation proposes a theoretical correction to predict the overall average porosity or the packing density of a packing bed of unequal spheres. For the binary packing system an index of the correction function is obtained. Our correction makes the predic-
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3495
Figure 7. Slope at different j0 for various w.
tion results match the experimental data for a large range of j0 (