Critical Choice of Variables in the Modeling and Correlation of

Although the choice of variables is known to be a major factor in the ... The illustrative examples of both effective and ineffective choices of varia...
0 downloads 0 Views 190KB Size
Ind. Eng. Chem. Res. 2002, 41, 3907-3929

3907

Critical Choice of Variables in the Modeling and Correlation of Mechanisms for Transport and Reaction Stuart W. Churchill† Department of Chemical Engineering, University of Pennsylvania, 311A Towne Building, 220 South 33rd Street, Philadelphia, Pennsylvania 19104

Although the choice of variables is known to be a major factor in the modeling of mechanisms of transport and reaction, as well as in the derivation of analytical solutions, the efficient execution of numerical solutions, and the development of correlating equations thereof, the analysis presented here may be the broadest and most detailed one yet attempted in this context. Although some general principles and methodologies are identified and described herein, the most effective techniques and choices are found to be somewhat specific to the particular technical topic. The illustrative examples of both effective and ineffective choices of variables are taken from fluid flow, chemical conversions, thermal conduction, thermal radiation, and thermal convection because of their general familiarity, but the results are presumed to be applicable to all aspects of chemical engineering and allied topics, both new and old. The techniques identified as most effective are dimensional analysis, the elimination of unimportant variables, the derivation of asymptotes, the introduction of virtual variables, the recognition of analogies, and the critical testing of results with experimental data or computed values. The most underutilized technique is that of speculation: asking the question “What if?”. This particular technique is applicable and often effective in conjunction with each of the others. Introduction The choice of variables has long been recognized to be of primary importance in the formulation of a model for a mechanism of transport or reaction, as well as in the formulation of algebraic correlating equations. The common objective of these two endeavors is to identify those variables and groups of variables that describe the mechanism most efficiently and to identify those variables that may be eliminated, replaced, or treated as minor parameters without the loss of significant functionality or accuracy. The primary difference in the selection of variables for a model and for a correlation is in the level of abstraction. Models generally consist of one or more differential or algebraic equations involving dependent variables such as the local and/or transient temperature, velocity, and concentration and independent variables such as the spatial coordinates and time. On the other hand, the variables in a correlation generally consist of dimensionless groupings such as the friction factor and the Reynolds number, which include bounding quantities such as the shear stress at the wall and space-averaged quantities such as the mixed-mean velocity. Despite these distinctions, the optimal variables for modeling and correlation, as well as the processes for their identification, are closely related, thereby justifying their joint consideration herein. The objective of this investigation has been to identify general principles and to formulate general guidelines for the selection of optimal variables for the modeling and correlation of mechanisms of transport and reaction. With one exception to be noted later in this introduction, that task does not appear to have ever been undertaken in a generalized sense, although it is usually an inherent step in a specialized sense in every theoretical and †

Fax: 215-573-2093. E-mail: [email protected].

experimental undertaking. Early in the course of this investigation, its organization in terms of technical topics such as thermal radiation, chemical conversions, and fluid flow was recognized to be more effective than its organization in terms of methodologies such as dimensional analysis, changes of variable, the formulation of analogies, and the speculative elimination of variables. This is because many of these techniques are uniquely effective with a particular technical topic. The very term illustrative example implies selectivity and omissions; although a number of topics in chemical engineering are examined, the examples herein are weighted in favor of fluid mechanics, heat transfer, and chemical reactions because of their greater familiarity to the author, and many utilize his own work for the same reason. The methodology of choosing optimal variables for single mechanisms or single combinations thereof is not to be confused with the quite different and more complex methodology of synthesizing an optimal situation for an extended set of coupled processes that is required for purposes of control (see, for example, work by Foss1 or Morari et al.2) or for the optimal operation of a chemical plant (see, for example, work by Edgar et al.3). Starting Point If a correlating equation is to be developed for a set of experimental or computed values for some limited process or mechanism, the starting point is usually the preparation of a list of variables. Although a differential model may provide some guidance for the variables to be listed, that guidance is usually insufficient because of the aforementioned difference in the level of abstraction. Dimensional analysis is then applied. The method of Rayleigh4 is infallible insofar as the list is complete and not redundant. Redundancy arises from the inclusion of variables that are not independent of one

10.1021/ie010776o CCC: $22.00 © 2002 American Chemical Society Published on Web 06/26/2002

3908

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

another, for example, the inclusion of both the axial pressure gradient and the mixed-mean velocity in a list intended to describe the velocity distribution in flow through a pipe. Another source of redundancy is the separate listing of variables that function only as a product, for example, w1 and c1 and w2 and c2 rather than w1, c1 and w2, c2 in a listing for the extent of heating in a double-pipe exchanger, which results erroneously in separate dimensionless groupings of w1/ w2 and c1/c2 rather than just w1c1/w2c2. Redundancy is ordinarily not a serious problem because it will usually be revealed when the correlation is tested. On the other hand, the omission of a significant variable will almost certainly have serious consequences. For example, the failure of a correlating equation to represent the values closely may then be attributed falsely to experimental error. Most of the countless subsequent “improvements” of the methodology of Rayleigh are either unnecessary or erroneous. The dimensional analysis of an algebraic, differential, or integral model, together with initial and/or boundary conditions, by the method of Hellums and Churchill5 is also infallible insofar as the model is complete and selfconsistent. However, as contrasted with the method of Rayleigh for a list of variables, the results are critically dependent on the particular choice of variables in which the model is formulated. For example, expression of the model in terms of the stream function rather than the velocity components alone may lead to a fundamentally different result. The speculative deletion of a possibly insignificant variable from a list or of a term from a model often leads to a greatly simplified and/or generalized representation. This technique provides a great opportunity for ingenuity. However, it is fraught with risk and must always be tested. For example, if the time dependence is omitted and the process turns out to be oscillatory, any solution obtained by virtue of that postulate is meaningless. Conceptual simplifications or introductions, such as a global reaction mechanism, an effective viscosity, a fictive film thickness, or the existence of a pseudo stationary state, generally lead to results of uncertain validity. One earlier analysis may be mentioned at this point. Klinkenberg and Mooy6 applied dimensional analysis to the partial differential equations of conservation for momentum, energy, and chemical components in general as well as in various reduced forms. They thereby identified and characterized physically a large number of dimensionless groups. Their analysis is instructive and generally sound, but not having available the aforementioned method of Hellums and Churchill or its equivalent, they were unable to identify by dimensional analysis alone the reduced groupings that result from similar transformations, except in a few cases by virtue of reference to the classical solutions that made use of such a transformation or implicitly defined it. Also, they did not consider the time-averaged form of the equations of conservation and were therefore forced to fall back on empirical correlating equations as a source of dimensionless groupings for turbulent flow and convection. They presented, without any theoretical justification, some relationships in the form of products of powers of dimensionless groups but did note that these relationships were inherently approximative except as asymptotes.

Attention is now turned to particular technical topics as a source of illustrative choices of variables. Thermal Conduction Thermal conduction is chosen as the first technical process to be examined with respect to the choice of variables because of its relative simplicity and its welldeveloped mathematical structure. Steady thermal conduction in a solid or stagnant fluid is often presumed to be represented together with appropriate boundary conditions by

∇ ‚ k∇T ) 0

(1)

thus introducing the temperature and the thermal conductivity as variables. Van Dusen7 suggested the reduction of eq 1 by replacing k∇T with k′∇T ′, where k′ is invariant, to obtain

∇ 2T ′ ) 0

(2)

This transformation reduces all problems in pure steady conduction to ones with an invariant thermal conductivity. If the value of T ′ is arbitrarily taken to be zero at some reference temperature T0, the solution in T ′ can be converted to one in T by means of the integral

T′)

∫TT(k′k ) dT 0

(3)

The application of this transformation to transient conduction, as represented by

∂T ∂t

(4)

(Fck )∇ T ′ ) ∂T∂t ′

(5)

∇‚k∇T ) Fc results in 2

Because k/Fc generally varies less with temperature than does k, this is a somewhat improved representation. For invariant k, F, and c, the application of the operator ∆ to eq 4, followed by reversal of the order of differentiation on the right-hand side, gives

(Fck )∇ (∇T) ) ∂t∂ (∇T) 2

(6)

which may be expressed as

(Fck )∇ j ) ∂t∂j 2

(7)

where j ) -k∇T is the heat flux density vector. This latter change of variable has two primary consequences. First, eq 7 is more convenient than eq 5 if the heat flux density rather than the temperature is specified initially and at the boundaries. Second, all of the solutions in the literature for eq 5 are applicable for eq 7 insofar as the boundary and initial conditions are congruent. The method of dimensional analysis of Hellums and Churchill5 reduces eq 5 to an ordinary differential equation that may be integrated formally for a semiinfinite region, initially at T0 and with a temperature Tw imposed and maintained at the surface at x ) 0 for t > 0, to obtain the well-known solution

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3909

{

T - T0 x ) erfc Tw - T 0 2(Rt)1/2

}

(8)

which implies that

j0 ) k(Tw - T0)/(πRt)1/2

(9)

where R ) k/Fc and j0 is the heat flux density at the surface. This solution has long been recognized to be faulty for very short times in that a finite temperature rise is implied for all distances and times. The rate of transport of energy is thereby implied to have a velocity exceeding that of light. The inherent postulate of an invariant density and therefore of no mass movement is the source of this anomaly. If eq 5 is solved simultaneously with the equation of conservation for momentum, the continuity equation, and an equation of state that provides for a change of density with temperature, a compressive, slightly supersonic pressure wave is predicted with an accompanying rise in temperature. Rayleigh8 in 1899 derived an approximate solution for the compressive wave motion in a gas heated instantaneously at a plane but did not associate his results with the anomalous prediction of eq 5. This behavior, which has now been completely resolved both theoretically and experimentally by Brown and Churchill9 and which occurs with liquids and solids as well as with gases, needs to be accounted for in a practical sense only for the very rapid imposition of a large temperature difference, such as with lightning, for which the pressure wave is evident as thunder. The final development of a successful numerical solution for this process revealed the dramatic effect of the choice of an optimal relationship between the grid size and the time step. The latter relationship may be considered to be a variable in the context of the numerical algorithm. All prior numerical solutions satisfied separate relationships between the grid size and time step, as imposed by the finitedifference models for diffusion and wave motion, respectively, but did not satisfy the combined relationship. This combination required for stability an increase of 104 in the number of grid spaces and also in the number of time steps. The combined increase of 108 in the required number of numerical steps for stability was overlooked by earlier investigators because of the choice of an insensitive criterion for convergence, namely, the amplitude of the pressure wave. A more critical criterion is the shape of the pressure wave. The unique lesson from this example is the possible critical dependence on variables such as ∆x and ∆t within the numerical algorithm. Important if not unique lessons are (1) the critical effect of the retention in the model of the density as a distinct variable even though it changes only slightly and (2) the testing of the numerical results with experimental measurements. This behavior provides a still further example of the erroneous choice of a model. Solutions based on the so-called hyperbolic equation for conduction

∂2T ∂T R ∂2T + R ) ∂t uw2 ∂2t ∂x2

(10)

with uw usually chosen arbitrarily to be the velocity of sound, appear throughout the literature. Equation 10 avoids the prediction of an infinite rate of transport of energy but is quite fallacious in its predictions other-

wise, as well as in its structure. The lesson here is that elegant models and solutions may have nothing to do with the real world. Transient radial conduction in a spherical region for the idealized conditions corresponding to eq 5 may be represented by

R ∂ 2 ∂T ∂T r ) 2 ∂r ∂r ∂t r

(

)

(11)

Substituting U/r for T reduces eq 11 to

R

∂2U ∂U ) ∂t ∂r2

(12)

Solutions for a planar region are thus applicable for a spherical region insofar as the boundary conditions can be matched. The analogous expression for a cylindrical region is

∂T ∂T R ∂ r ) r ∂r ∂r ∂t

( )

(13)

Substituting ez for r reduces eq 13 to

R ∂2T ∂T ) ∂t e2z ∂2z

(14)

Equation 14 is possibly advantageous over eq 13 for numerical integrations, but its primary advantage is for steady-state conduction for which it becomes congruent with the equivalent expression for a planar region. Churchill10 observed that the solution for the temperature distribution in the region outside a sphere subjected to a step in wall temperature can be expressed as

() {

T - T0 a r-a ) erfc Tw - T0 r 2(Rt)1/2

}

(15)

}

(16)

and outside a cylinder as

()

T - T0 a ) T - T0 r

1/2

{

erfc

r-a 2(Rt)1/2

Hence, a common exact solution for all three geometries is

() {

T - T0 a ) Tw - T 0 r

n

erfc

x 2(Rt)1/2

}

(17)

where n ) 0, 1/2, and 1 for a plane, cylinder, and sphere, respectively. Equation 9 is an approximation for the heat flux density for short times for all three geometries. Equations 12 and 17 illustrate generalities that can be achieved for different geometries by the proper choice of variables. Churchill and Chang11 noted that thermal conduction from bodies of finite surface area to infinite surrounding may be represented approximately by a numerical value of Nu ) hL/k, where L ) (4A/π)1/2 for planar surfaces and L ) (A/π)1/2 for solids. They also derived exact values of L for simple shapes and developed upper and lower bounds for more complex ones. The following solution for the heat flux density to the isothermally heated surface of a finite region of thickness δ with the second surface maintained at the initial

3910

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

temperature may be derived by separation of variables:

j)

( )

k(Tw - T0) δ2 δ

πRt



[1 + 2

∑ (-1)ne-n δ /Rt] n)1 2 2

(18) ce ) cs[(Τf - Τw)/(Τ0 - Τw) erf{λ}]2

whereas the following alternative solution may be derived by means of the Laplace transform: ∞

j ) k(Tw - T0)[1 + 2

e-n π Rt/δ ] ∑ n)1 2 2

2

(19)

Equation 18 converges rapidly for long times but not for short ones, and eq 19 has the converse characteristics. Thus, δ2/Rt and Rt/δ2 may be considered as alternative and complementary combinations of variables, and thereby as alternative effective dimensionless variables that are associated with the two methods of solution. This methodology is applicable with the same advantages for many problems in transient thermal conduction. A series solution of the following form can often be derived for transient thermal conduction in threedimensional bounded regions: ∞

T{x,y,z,t} )

fne-R ∑ n)1

2t

+ T{x,y,z,∞}

n

(20)

Schoenherr and Churchill12 suggested the following efficient method for solving the differential model numerically rather than summing this series. For sufficiently long times, the series may be approximated by one term, that is, by

T{x,y,z,t} = f{x,y,z} e-R1 t + T∞{x,y,z,∞} 2

namely, eq 9, could be converted to that of Neumann for freezing simply by introducing the following expression for the effective heat capacity of the frozen region:

(21)

where here λ is the eigenvalue of the exact solution of Neumann for freezing in the semi-infinite region. They confirmed the usefulness of this virtual variable by comparing the resulting predictions for the region outside a cylinder and in a two-dimensional quarter space with numerical solutions of the unreduced differential model. Churchill and Gupta15 concluded, however, that the solution of Neumann must be considered highly idealized when applied to a finely pored material because the latent heat of freezing is given up over a range of temperature owing to the capillary pressure and, furthermore, because liquid water is transported toward the freezing front owing to its removal by freezing as well as to the variation of its surface tension with temperature. Thus, physical variables not included in the idealized differential or algebraic model may have a significant influence. The validity of the concept of an effective thermal conductivity for porous media is questionable for the transient behavior, as indicated in the previous paragraph, but not for the steady state. Churchill16 reviewed the literature for the latter regime and devised several generalized correlating equations for this virtual variable, of which only a concept originating with Maxwell17 will be examined here. Using the principle of imbedding, Maxwell derived an exact solution for electrical conduction through a static dispersion of uniformly sized spheres sufficiently dilute so that the disturbance of the electric field by any one sphere does not disturb that of its neighbors significantly. That solution may be expressed in thermal terms as

Taking the first and second derivatives of eq 21 with respect to t and eliminating f{x,y,z} and R between these three expressions results in

T{x,y,z,∞} = T -

(∂T/∂t)2 ∂2T/∂t2

(22)

A finite-difference algorithm is used to solve the differential model, namely, eq 5 (possibly also including a source term, etc.), for T, ∂T/∂t, and ∂2T/∂t2. This calculation may be halted as soon as T{x,y,z,∞} ceases to change significantly. Convergence is very rapid and requires far less computational work than either the summation of the series in eq 20 or the execution of a conventional numerical solution. This procedure is also effective for steady-state problems, for which a false transient term with an arbitrarily large thermal diffusivity is simply added to the differential model. Despite the somewhat lengthy description here, this procedure is straightforward and can be conceived in the context of this paper as simply the substitution of the new variable (∂T/∂t)2/(∂2T/∂t2) for the infinite series of eq 20. Franz Neumann is credited by Carslaw and Jaeger13 with the derivation of an exact analytical solution for freezing in a stagnant semi-infinite region, such as a wet porous media with an initial temperature above the freezing point, following the imposition and maintenance of a temperature below the freezing point on the surface. Gupta and Churchill14 found that the solution for the heat flux density at the surface without freezing,

(23)

ke 1 + 2φ ) kc 1-φ

(24)

φ ) (kd - kc)/(kd - 2kc)

(25)

where

Here, ke is the effective conductivity of the dispersion, kc the thermal conductivity of the continuous material, kd that of the dispersed material, and  the volumetric fraction of the latter. Exact numerical solutions for orderly arrays and experimental data for randomly dispersed, uniformly sized spheres, when plotted as ke/ kc versus φ, deviate significantly from eq 24 only for  approaching max, corresponding to a packed bed, and then only upwardly on the mean, indicating that eq 24 is a lower bound. Furthermore, the deviations are essentially a function of φ only, with a minimal additional parametric dependence on kd/kc and . This same set of coordinates provides a good representative for nonspherical particles, nonuniformly sized spheres, and fluidized beds, again with eq 24 as a lower bound. A similar theoretical solution for an array of cylinders with axes parallel to each other and normal to the direction of conduction results in a good representation for nonuniformly sized cylinders with their axes in the plane normal to the direction of conduction but randomly oriented within the plane, such as with fibrous insulations. These theoretically based compound vari-

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3911

ables for high idealized conditions thus prove to be remarkably effective virtual variables for nonidealized conditions. It goes without saying that neither the quantity φ nor eq 24 would have ever been discovered in the absence of the solution of Maxwell. It should not be inferred that eq 24 and its analogue for cylinders are of sufficient accuracy for practical purposes for all of the cases mentioned (see, for example, work by Churchill16 and Molerus18). The unique merit of eq 24, together with eq 25, is in the insight provided by such a generalized representation and in the criterion provided for the quantitative determination of the effects of the various inherent idealizations and the associated parameters. A very useful if simplistic concept for the representation of resistances in series is that of Ohm’s law

δo ) k oA o

δi

∑kiAi

(26)

Here δi, ki, and Ai are the individual thicknesses, conductivities, and areas and δo, ko, and Ao the overall values. The use of ko as a virtual variable eliminates the interfacial temperatures from the model. Equation 26 is, of course, applicable for convection and/or radiation in series with conduction as well. Ohm’s law for resistances in parallel is also useful for conduction, convection, and radiation even though the elimination of interfacial temperatures is not involved. The objective of presenting these examples of effective charges of variable for conduction, some of which are well-known, is to suggest the possible applicability of the same or similar concepts for other more complicated processes. Thermal Radiation The rigorous description of radiative interchange between surfaces and between gases and surfaces is much more complicated than thermal conduction or, for that matter, than any of the other processes to be examined herein. This complexity has inspired many virtual variables and approximations out of necessity. A book by Hottel and Sarofim19 is in some ways a compilation of such expedients, only a few of which will be mentioned here because they are not generally applicable or adaptable for other processes. The simplest and perhaps most useful virtual variable for thermal radiative interchange is the radiative heattransfer coefficient

hr )

σ(T14 - T24) ) σ(T22 + T12)(T1 + T2) ) T1 - T2 4σTeffective3 (27)

The first term on the right-hand side of eq 27 is itself a simplified expression in that different values of the emissivity apply to the two surfaces or to a surface and a volume of gas (see, for example, the book by Hottel and Sarofim, pp 297-301). The principal value of this representation is the linearization of the dependence on T1 and T2, thereby allowing radiation to be combined with conduction or convection by means of Ohm’s laws. Insofar as the fractional variation in Teffective3 in absolute units is small, some arbitrary average may be used. Otherwise, Teffective3 can be calculated rigorously by iteration. The equivalent representation in terms of thermal conduction is

kr ) 4σTeffective3δ

(28)

where δ is some measure of the distance between the sources of radiation. The scattering, absorption, and reemission of radiation by a single particle, and hence radiative transfer through a dispersion of particles, such as in a cloud, depends on the diameter and the complex index of refraction of the particle as well as on the wavelength and degree of polarization of the impinging radiation. Mie20 solved the Maxwell equations exactly to obtain a series solution for this process for a single sphere in terms of Riccati-Hankel and Riccati-Bessel functions of half-integral order. Following a suggestion of Hartel,21 Chu and Churchill22 reexpressed this solution exactly for randomly polarized radiation by a series of Legendre polynomials, which makes computing of the angular distribution of the scattered radiation much simpler. In this instance a detailed process of reexpression rather than a simple substitution results in a reduction in the required number of tabulated functions and in the number of terms to be summed. Scattering by a dispersion of particles is very complicated, primarily because of the highly anisotropic distribution of scattering by single particles. This complexity has led to the introduction of many virtual variables and approximations in order to obtain solutions. For example, Schuster23 proposed modeling of the process in astronomical applications by dividing both singly and multiply scattered radiation into forward and backward components. The components of the two-flux method are uniquely defined, although the model for multiple scattering is approximate. Chu and Churchill24 improved upon this two-flux model by dividing the singly and multiply scattered radiation into six components, including four sidewise ones. The six discrete fluxes for single scattering as well as those within the dispersion are not uniquely defined and thereby are virtual variables. The approximation of scattering within a dense dispersion as a diffusional process is an obvious alternative to the use of discrete components, but difficulty is encountered with strongly peaked distributions and with the formulation of boundary conditions. Single scattering is never isotropic, but the angular distribution of the multiply scattered radiation within a dispersion approaches isotropy as the concentration of the particles increases. Chu et al.25 devised an exact higherorder diffusional model for the limiting case of isotropic scattering within a dispersion that incorporates the boundary conditions inherently. This methodology may be interpreted in the context of this investigation as the determination of an exact virtual diffusivity. Larkin and Churchill26 applied the two-flux method for heat transfer through fibrous and foamed thermal insulations, for which thermal radiation is a significant factor. They determined the effective scattering and absorption coefficients by direct transmission of blackbody radiation and reexpressed their results for convenience in terms of the effective thermal conductivity defined by eq 28. Chin and Churchill27 carried out similar measurements for packed beds of various granular materials and Cimini and Chen28 for fluidized beds. One last virtual variable for radiation will be mentioned. The monochromatic transmission τλ of radiation through a film of liquid of thickness δ may be represented by the Bouguer-Beer law

3912

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

τλ ) e-σλδ

(29)

where σλ is the monochromatic absorptivity. However, σλ varies almost discontinuously within an absorption band, whereas all current instrumentation has a finite bandwidth, putting experimental measurements of σλ into doubt. Elsasser,29 by idealizing an absorption band as an infinite array of equidistant lines, each of the Lorentz shape and the same half-width, defined an effective absorptivity in terms of the smoothed coefficient by means of the following expression:

erfc{γλδ} )

λ+(∆λ/2) exp{-σλδ} dλ ∫λ-(∆λ/2)

1 ∆λ

(30)

Chin and Churchill30 showed that any function of σλδ can be transformed into a function of γλδ by means of

1 ∆λ

F{γ x/z}

λ+(∆λ/2) λ 1 1 F{σλx} dx ) ∫0 1/2 dz ∫λ-(∆λ/2) π z (1 - z)1/2

(31)

which allows a solution in terms of σλ to be expressed in terms of γλ. Friedman and Churchill31 utilized this methodology to calculate the absorption of radiation by droplets of fuel and Chin and Churchill32 for the transmission of radiation through an atmospheric haze. The scattering of neutrons in nuclear reactors is closely analogous to the scattering of thermal radiation but with two significant differences; the scattering is much less anisotropic but the variable velocity of the neutrons must be taken into account. Even so, these virtual variables for light scattering may have applicability. Although the virtual variables of this section are not directly applicable for most of the other processes considered herein, they do suggest the potential of approximate representations for other very complex behavior. Chemical Conversions Chemical kinetics, as contrasted somewhat with thermal radiation and almost wholly with thermal conduction, is generally approached as a purely empirical subject. Consequently, the primary intent here is to identify and evaluate the virtual and approximate variables that are in present usage rather than to propose new ones. Chemical processing, which does have a theoretical structure, will only be examined briefly because the choice of variables is more constrained and less critical. Although a thermodynamicist might expect chemical reaction rates to be expressed in terms of the chemical potential or the fugacity of the individual reactants and products, in most textbooks on reaction kinetics and reaction engineering these rates are expressed in terms of molecular concentrations, presumably because of convenience and tradition. Convenience is perhaps a good excuse for this usage if the possible error is noted, but tradition is not. The classical simplification of chemical kinetics is the representation of a system of variable-order reactions in series and parallel by a global, first-order mechanism. In its simplest manifestation, a second-order reaction mechanism, kCACB, is represented by k′CA ) kCBCA. If component B is in large excess, a mean value may be chosen for CB with little uncertainty. For multiple

mechanisms, the formulation of an equivalent approximation is more difficult and the result more uncertain. Experimental data for the rate of reaction in gases at high temperature often exhibit fractional orders of reaction as well as a nonlinear reciprocal dependence on the reactants and products. Christiansen33 and others explained this somewhat surprising behavior by means of extended mechanisms involving free radicals with pseudo-steady-state concentrations, that is, with equal rates of formation and consumption. An observed dependence of the rate of reaction on pressure may be rationalized on the same basis, and analogous behavior usually is implied for reactions catalyzed by a solid surface. A global reaction mechanism in terms of fractional orders for the fuel and oxygen has often been postulated for combustion. Such a mechanism, with an appropriate virtual rate constant based on experimental data, may provide an adequate prediction of the rate of combustion and of the temperature profiles in the preflame, flame, and postflame zones because the bulk of the fuel is rapidly oxidized to H2O and CO, and the CO is then somewhat more slowly oxidized to a near-equilibrium concentration of CO2. However, this mechanism is fundamentally wrong in that no fuel burns directly but first must be pyrolyzed to molecular fragments that are then oxidized step by step through a long chain of reactions to CO, unburnt fuel, and partially burnt fuel. At the same time atoms of N2 are decomposed and oxidized to NO and NO2 through additional free-radical reactions. Zel’dovich34 proposed the supplementation of such global models with a psuedo-steady-state freeradical model for the postflame zone as suggested by the modeling of Christiansen. Although this model predicts the concentrations of NO, NO2, and CO in the burnt gases, which the global models are unable to do, the predictions are highly inaccurate. Pfefferle35 subsequently utilized a reaction mechanism involving 67 reactions and 31 components for the pyrolysis and combustion of ethane and thereby demonstrated that pseudo-steady-state concentrations of the important free radicals such as OH are not attained in the flame and postflame zone. This work explained the failure of the predictions of NO, NO2, and CO by the Zel’dovich model and produced far more accurate predictions. In the context of the current analysis, it may be noted that the method for local sensitivity analysis of Rabitz and Dougherty,36 together with a global trial-and-error method, was utilized by Pfefferle to select the 67 mechanisms from 115 initially proposed ones. Although many of the 67 rate constants and 67 equilibrium constants in that model are highly uncertain, the uncertainty of the overall model is fortuitously reduced by the multiple paths of reaction. All in all, the optimal selection of variables, including free radicals, together with the associated mechanisms, is a necessity for the calculation of accurate residual concentrations of the pollutants resulting from combustion. When the model of Pfefferle was applied to thermally stabilized combustion inside a ceramic tube, for which backmixing of the free radicals is virtually eliminated, an unexpected difficulty known mathematically as “stiffness” was encountered in the numerical integration of the differential component balances and the integrodifferential energy balance. This difficulty was resolved by White,37 who proposed using only the first few

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3913

reactions for thermal pyrolysis of the fuel and thereby deriving an analytical solution valid for the first few nanoseconds or nanometers. The full kinetic model and numerical integration are applied thereafter. Chemists ordinarily model chemical reactions at constant volume, as does occur with gases in enclosures of fixed size and with liquids. Hence, they commonly express the rate of change due to homogeneous reactions in Lagrangian forms such as

r)-

dCA dt

(32)

The latter expression is not applicable for a nonequimolar or nonisothermal reaction at constant pressure. The corresponding more general expression is

r)-

d(CA/F) 1 d(CAV) 1 NA )) -F V dt V dt dt

(33)

Engineers usually deal with homogeneous reactions in continuous flow through piping rather than in batch reactors. The corresponding rate of change in composition may be expressed as

r)-

dnA d(vCA) d(uCA) ))dV dV dL

(34)

where here nA is the molar rate of flow of component A, v the total volumetric rate of flow, and u the linear velocity. Chemists sometimes apply eq 32 for continuous reactors in terms of space times, that is, the time required for the fluid at its entering density to pass through the reactor. This expedient is twice wrong. First, flow reactors usually operate at nearly constant pressure throughout, and the density may therefore vary by virtue of changes in temperature due to the reactions or to external heating or cooling, as well as due to changes in the number of moles. Second, the correct time to be substituted is the residence time

tR )

∫0Ldtu

(35)

which cannot be evaluated in advance. The same objections apply to the space velocity, which is the reciprocal of the space time. Considerable confusion exists in the literature of chemical kinetics and chemical reaction engineering concerning the rate of reaction, which is due to a molecular or atomic mechanism and is represented by expressions such as r ) kCACB, and the rate of change in composition, which is due to one or more such mechanisms and is represented by expressions such as eqs 32-34 (see, for example, work by Kabel38). For some reason, such confusion between the definition of a rate mechanism and an expression for the rate of change due to that mechanism does not appear to occur commonly in the treatment of physical processes such as heat transfer. If diffusion is negligible, eq 34 is applicable for a filament of fluid, whose radial location varies during passage through the reactor owing to changes in density and/or viscosity. If this variation is also neglected, the overall conversion at any length may be determined by integrating the indicated conversion over the cross section. Such an integration is ordinarily avoided by postulating plug flow and thereby substituting the

mixed-mean velocity for u inside the integrand of eq 35. An attempt to justify this idealization by a vague reference to turbulent flow is always in error in one or more respects. First, the velocity distribution in turbulent flow differs considerably from uniformity. Second, as noted by Churchill and Pfefferle,39 turbulent flow is rarely achieved in tubular reactors. Third, the radial velocity distribution for laminar flow results in significant variations in conversion between the centerline and the wall, which leads to radial diffusion. With modern computing hardware and software, these idealizations are unnecessary even in the rare instance of turbulent flow (see, for example, work by Collins and Churchill40). Many of the above considerations relevant to idealized variables in models for homogeneous reactions have an analogue in catalytic conversions. However, the one example of the latter to be considered here was chosen because of its possible applicability to nonreactive systems rather than because of its analogy to homogeneous conversions. Thiele41 derived the following solution for the fractional reduction of the first-order irreversible rate of reaction in a porous slab with a catalytic inner surface, taking into account the resistance to diffusion:

η ) tanh{τ}/τ

(36)

Here τ ) δ(k′S′/D′v)1/2 is called the Theile modulus, δ is the half-thickness of the slab, k′ is a global first-order rate constant for the catalyzed reaction, S′ is the specific surface, that is, the area of the catalytic surface within the pore space per volume of the slab, and D′v is the effective diffusivity of the reactant in the porous slab. Here, k′ and D′v are virtual variables. Corresponding exact solutions for η have been derived for spheres and infinitely long cylinders. Aris42 showed that the approximation of the solutions for these three geometries for large values of τ can be generalized exactly. Churchill43 and others have extended this generalization for an external resistance and other reaction mechanisms, all by means of the choice of virtual variables. Curiously, eq 36 was derived independently and almost simultaneously by Damko¨hler44 and Zel’dovich,45 and even earlier by Hatta,46 without the use of virtual variables, in terms of its analogue for diffusion of a gas into a liquid film followed by chemical reaction. In retrospect, many of the virtual variables that have evolved for chemical conversions are no longer necessary or appropriate owing to the feasibility of using more exact representations in conjunction with numerical computations. Global reaction mechanisms, the postulate of a psuedo steady state for free radicals, and the use of such concepts as the quenching distance and the induction time is apt to persist for some time in the field of combustion because of the insecurity and unfamiliarity of most mechanical and aeronautical engineers with chemistry. At the same time, the use of such unnecessary simplifying concepts as plug flow, space time, and space velocity is apt to persist for flow reactors because of the insecurity and unfamiliarity of most chemists with fluid mechanics. A preference for elegance over accuracy appears to be pervasive in the field of chemical reactor design. As an example, the present author once attended a presentation entitled “The Characteristics of Large Sets of Second-Order, Reversible, Isothermal Reactions” (or the equivalent), which elicited from the aforementioned J. A. Christiansen the comment “the behavior you describe may be interesting from a math-

3914

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

ematical point of view, but it has nothing to do with chemistry”. Fluid Flow The flow of fluids not only is of direct interest for a wide range of conditions but also is a necessary complement for the analysis of convection, regeneration, mixing, transport of solids, separations, and chemical processing. This broad subject has inspired the formulation of so many idealizations, approximations, generalizations, and virtual variables that only a small sampling can be examined here. These examples were chosen on the basis of probable unfamiliarity rather than for their importance. Laminar Flow in Channels. Fully developed laminar flow in channels, for which exact solutions are known in a number of instances, may be utilized to illustrate the possibilities for and the merits of different dimensionless groupings. For example, the velocity distribution in fully developed laminar flow in a round tube may be expressed exactly as

u)

aτw r2 12µ a

[ ( )]

(37)

The corresponding mixed-mean velocity is

um ≡

∫01ud(ar )

2

)

aτw 4µ

(38)

Equations 37 and 38 indicate that fully developed laminar flow is independent of the density of the fluid, a condition inherent in their derivation. Of course, expressions involving the mass rate of flow w ) Fumπa2 are inherently dependent on the density. From eq 37 it follows that the velocity distribution in terms of µu/aτw is a function only of r/a and indeed that µu/aτw[1 - (r/ a)2] is invariant. Replacing τw in eq 37 by um from eq 35 indicates that u/um is also a function only of r/a and that u/um[1 - (r/a)2] is invariant. The choice between these two dimensionless expressions for the velocity distribution depends on whether τw or um is considered to be the independent variable. Equation 38 is often expressed as

f ) 16/Re

(39)

where f ≡ 2τw/Fum2 and Re ≡ 2aumF/µ, but this canonical form is misleading because these two dimensionless groups both include the density, which is not a variable. A better generalized form is simply

Po ) 16

(40)

where Po ) 4aτw/µum is the Poiseuille number. Equation 37 is not directly applicable for flow between parallel plates of unlimited breadth. However, if r is replaced by y, the distance from the wall of the tube, eq 37 can be reexpressed as

u)

aτw y 1 y 2 µ a 2a

[

( )]

(41)

If a is simply replaced by b, the half-spacing of the plates, eq 41 becomes directly applicable for a parallelplate channel of unlimited breadth. The dimensionless velocity distributions in terms of µu/aτw and µu/bτw are

then the same function of y/a and y/b, respectively, although those in terms of u/um are not. The Poiseuille number for a parallel-plate channel of unlimited breadth in terms of b is 12. For other channels Po has a fixed value or is a function of secondary variables, such as the radius ratio of a circular concentric annulus or the aspect ratio of a rectangular channel. As an example, consider flow through a helical coil of radius r0, tube radius a, and pitch b. Dean47 derived an analytical solution for the limiting case of r0/a f ∞ and b f 0 that may be expressed in the following form:

( )

Pos Dn2 ) 1 - 0.03058 Po 288

2

+ 0.00728

( ) Dn2 288

4

(42)

where Dn ) (2aumF/µ)(a/r0)1/2 is the Dean number and the subscript s indicates a straight tube. Equation 42 may also be written as

um ) 1 - 0.03058K2 + 0.01195K4 ums

(43)

where K ) (F2a7/2µ4r0)[-(1/48) (dP/dx)]2 and x is the axial distance. Because the applicability of this solution is limited to Dn < 20.45, for which Po/Pos < 1.02, its practical value is limited to the identification of Dn or K as the primary characteristic variable for all conditions, which is, however, a very useful discovery that does not follow directly from dimensional analysis of the general model. Truesdell and Adler48 generalized approximately the Dean number for flow through coils of finite pitch, for which there is no orthogonal set of coordinates, by replacing the radius of the coil, r0, by the radius of curvature, rc ) r0[1 + (b/2πr0)2] (see work by Manlapaz and Churchill49 for further details). Developing laminar flow in the inlet region of a round tube or parallel-plate channel provides an example of the error that may occur due to the wrong choice of a virtual variable in modeling. Dealy50 deduced from pure conceptual reasoning that in the regime of development the maximum in the velocity does not occur at the centerline or central plane, thereby negating the several approximate analytical solutions in the literature, all of which postulate a shrinking zone of uniform velocity. His reasoning was subsequently confirmed by essentially exact numerical solutions. Flow through Porous Media. An equation devised by Ergun,

-∆P3Dp2

Dpu0 ) 150 + 1.75 µ(1 - ) L(1 - ) µu0 2

(44)

where u0 is the superficial velocity,  the void fraction, and Dp the particle diameter, was found (see work by Ergun51 or Churchill52) to be very successful in representing experimental data for a wide range of conditions. For porous media not consisting of uniform spherical particles, Dp was represented by 6(1 - )/av, where av is the specific surface. Although it was not derived that way, the reduced form resulting from dropping of the rightmost term may be rationalized on the basis of eq 38 by purely geometric considerations. That is, it can be interpreted as a Poiseuille number. On the other hand, DeNevers53 rationalized the form of the reduced expression resulting from dropping of the laminar (constant) term on the basis of flow through a rough pipe. None of the early rationalizations for the

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3915

additivity of these two terms is very convincing, but Niven54 has recently shown that a model of the porous volume consisting of large and small channels in series leads to an additive dependence. This example wellillustrates the power of conceptual analysis. The two dimensionless variables of eq 44, together with the definition of Dp, incorporate all of the dimensional variables that have been quantified, and that expression itself predicts the complex dependence of the pressure drop on the superficial velocity, on the viscosity and density of the fluid, and on the void fraction and structure of the solid (the latter as represented by av). Just as with eqs 24 and 25 when applied to thermal conduction, eq 44 and the associated definition of av cannot be expected to provide precise predictions for all of the diverse structures encompassed by the term porous media. Indeed, in this instance there is some quantitative evidence concerning the magnitude of the uncertainty associated with the oversimplifications. MacDonald et al.55 tested the predictions at eq 44 with later experimental data for many diverse media, including uniformly sized fibers, unconsolidated granular particles, and even consolidated media. They concluded that the functional dependence was satisfactory but that a mean value of 180 rather than 150 for the leading coefficient and values ranging from 1.8 to 4.0, depending on the roughness, rather than 1.75 for the other coefficient, gave a better representation of the mean. Molerus and Schweinzer56 compared with experimental data for unconsolidated particles but not with eq 44 the prediction of an alternative expression of Molerus57 that differs functionally by virtue of an additive term in u0-1/2 to represent the contribution of the drag in the thin laminar boundary layer on the particles and by virtue of an additive term in u0-n to represent the contribution of turbulence. They also tested a second expression for very irregular particles in which the term in u0-n was dropped but a shape factor was included in all of the other terms. The predictions of these two expressions are presumably more accurate numerically and functionally but at the expense of six common and one differing arbitrary coefficients as compared to the two of eq 44. This marginal improvement may be possibly justified in an intrinsic sense by virtue of identifying the contributions of the laminar boundary layer and turbulent regimes, but probably not in practice. Flow through Orifices. The limiting expressions for the pressure drop through a sharp-edged orifice from some distance upstream to the vena contracta (the minimum in the diameter of the jetting flow), which is commonly used to measure the rate of flow, bear a surprising similarity to those for porous media in the laminar and inertial regimes. The following expression was derived by Roscoe58 for creeping flow through a circular orifice in the asymptotic limit of a0/a1 f 0:

u0 ) πa(-∆P)/3µ

(45)

Equation 45 may also be expressed as

C0 ≡ (Re0/12π)1/2

(46)

where C0 ≡ Fu02[1 - (a0/a1)]/2(-∆P) = Fu02/2(-∆P). Equation 45, and hence eq 46, has been found to represent experimental data well for Re0 < 5 and small values of a0/a1. The theoretical solution (see work by Flu¨gge59)

φ ) 0.5793

(47)

where φ ) uv/u0, the contraction ratio for the vena contracta, and C02 ) [(a1/a0)2 - 1]/[(1/φ2)(a1/a0)2 - 1] = µ2 represent the data equally well for Re0 > 400. On the other hand, most current textbooks and handbooks cite the erroneous value of π/(π + 2) ) 0.6110, which is the value of φ for a planar orifice for the equivalent limiting conditions of a0/a1 f 0, based on the solution for a free streamline. The numerical differences are small but are misleading from a theoretical point of view. The coefficients φ and C0 for a circular orifice go through a maximum between the domains of eqs 46 and 47, and this intermediate behavior does not yet seem to have been represented by a theoretically based expression. Equation 46 appears to refute the asserted representation of all laminar flows by Po. The explanation is that C0 is not a drag coefficient but simply a coefficient that yields a constant value for the orifice for purely inertial flow (Re0 f ∞). External Flows. Textbooks on fluid mechanics are generally focused on unconfined flow because of the more extensive theoretical structure. Only five examples will be mentioned here, and even then very briefly. Stokes60 derived an exact solution for creeping flow over a sphere, for which the total drag force was found to be given by

Ct ) 12/ReD

(48)

where Ct ≡ FD/πa2Fu∞2. This solution, whose applicability is limited to Re < 0.1, is unique in that corresponding ones for cylinders, flat plates, etc., do not exist. Blasius61 discovered a similarity transformation for unconfined flow along a flat plate in the thin-laminar-boundarylayer regime with the free-stream velocity postulated to exist outside the boundary layer and thereby derived a series solution for the velocity field and the following expression for the local drag coefficient due to friction:

Cf ) τw/Fu∞2 ) 0.332Rex1/2

(49)

Imai62 utilized the displacement of the free-stream velocity by the boundary layer to extend the Blasius solution for the regime near the leading edge, obtaining thereby for the integrated-mean coefficient

C hf )

1.168 0.664 + Rex Re 1/2

(50)

x

Van Dyke63 notes that “it is remarkable that 50 years were required to discover that one term of the boundary layer solution provides two terms for the friction drag”. However, it should be noted that the coefficient of 0.664 in eq 50 implies integration of Cf as given by eq 49, from the leading edge (x ) 0), where it does not apply, to x ) L. Newton64 prior to 1710 derived an expression of the resistance of a falling solid sphere that is equivalent to

C h t ) 1/4

(51)

He correctly deduced the functional dependence of the drag on the density, the velocity, and the diameter, but his quantitative reasoning is quite faulty and the choice of the essentially correct value of 1/4 for the coefficient undoubtedly benefited from his experimental observations.

3916

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

Goertler65 derived an integral change of variables for curved surfaces, including flat plates, cylinders, and spheres, that generalizes the differential model and was expected to result in a more rapidly converging series solution than that of Blasius.61 This last expectation has not been completely fulfilled, but the adaption of the transformation for free convection is, as will be shown, very effective. Rising Bubbles. The upward mention of single bubbles in a liquid provides a unique spectrum of dimensionless dependent variables, independent variables, and parameters, as well as of regimes of flow. The most explicit relationship is perhaps uT* ≡ φ{Dv*,ζ,ξ,M}, where

( ) ( )

1/3 FL2 1/3 2Re ) ) µLg∆F 3Cf dimensionless terminal (steady-state) velocity

uT* ≡ uT

( ) ( )

Dv* ≡ Dv

FLg∆F

1/3

3Re2Ct ) 2

1/3

) µL2 dimensionless volume-equivalent diameter ζ ≡ µd/µc ) viscosity ratio ξ ≡ Fd/Fc ) density ratio

and

M≡

µ4g∆F ) Morton number FL2σ3

or as Ct ≡ φ{Re,ζ,ξ,M}, where

Ct ≡

2Dvg∆F 2

)

3FLuT

2Dv* 3(uT*)2

) total drag coefficient

and

Re ≡

DvuTFL ) uT*Dv* ) Reynolds number µL

A number of alternative dimensionless groups have also been used to describe or correlate this behavior, including

(

)

DvuT2FL 2M(Re)4 ) We ≡ σ 3Ct

Eo¨ ≡

) M1/3Dv*(uT*)2 ) Weber number

Dv*g∆F 3 1/3 9 ) (We)Ct ) M(Re)4Ct2 ) M1/3 σ 2 4 (Dv*)2 ) Eo¨tvo¨s number

(

Fr ≡

Ar ≡ and

1/3

)

FLuT2 2 ) ) Froude number Dvg∆F 3Ct

Dv3FLg∆Ρ µι2

) (Dv*)3) Archimedes number

E ≡ Dmax/Dmin ) eccentricity ratio The choice between all of these alternatives is usually based on the intrinsic form of a theoretical solution or the simplicity of a correlating equation. The density ratio ξ is not a significant variable under ordinary conditions, and it may be dropped from the above functional relationships and ∆F replaced by FL. The factor g then serves as a marker of the role of FL in buoyancy as opposed to inertia. The viscosity ratio ζ is also not a significant factor even though it appears in some of the theoretical solutions mentioned below. As might be expected from this multiplicity of variables and parameters, a large number of solutions have been derived for various limiting conditions. These include Stokes’ law for small bubbles coated with a surfactant, the Hadamard-Rybczynski law for small bubbles with a free surface, the Levich-Ackeret equation for spherical bubbles with a free surface in the thinlaminar-boundary-layer regime, the Harper-Moore equation, which extends the Levich-Ackeret equation to lower values of Dv*, the Moore equation for a slightly deformed bubble in the laminar regime, the Davies and Taylor equation for a hemispherical (mushroom-cap) bubble, the Mendelson equation for wave motion along the surface of a malleable bubble, and the very similar Lehrer equation. References and details concerning these solutions are given by Churchill,66 and generalized correlating equations for Ct and uT* are also given by Churchill.67 Curiously, the Levich-Ackeret equation, as well as the Stokes and Hadamard-Rybczynski equations, predicts a constant value for CtRe, which implies no inertial effects in the thin-boundary-layer regime. The Moore equation for slightly distorted bubbles includes E, the eccentricity ratio, as a parameter. The latter quantity is actually a dependent variable for which experimental data have been correlated as a function of We or Eo¨ , implying independence of E from the viscosity. The quantities Ct and uT* bear a one-toone correspondence. Hence, separate solutions for the various regimes and separate overall correlating equations are redundant but convenient. Small rising bubbles sometimes follow a spiral path rather than a vertical line, resulting in a deviation from the theoretical expressions for both Ct and uT*. This multiplicity of dimensionless variables for rising bubbles is somewhat unique and potentially confusing and thereby provides an important example in the context of this investigation of the choice of both dimensional and dimensionless variables. Turbulent Flow. The structure that has been deduced for turbulent flow differs radically from that for thermal conduction, thermal radiation, chemical kinetics, and even the flows considered above in that it is almost wholly conceptual in origin rather than theoretical or empirical. In the interest of simplicity, the treatment here is primarily for fully developed flow in a round tube, but it is readily adapted for channels of other shapes and even for unconfined flow along a surface. The greatest advance of all time in the description of turbulent flow is that of time averaging by Reynolds.68 That process, which he actually carried out in terms of space averaging, replaces the time-dependent pressure and components of the velocity by time-averaged ones plus some additional quantities such as u′v′. The particular results for a round tube, fully developed flow,

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3917

invariant physical properties, and a steady state may be integrated once to obtain the following expression for the momentum balance in the negative radial direction:

du - Fu′v′ τ)µ dy

(52)

Here, in the interest of simplicity and generality, the term (-dP/dx) (r/2) has been replaced by the local shear stress τ. Equation 52 is the usual starting point for the analysis of all turbulent shear flows. Direct numerical simulation, which is based on the primitive timedependent representation, is the principal exception. This latter methodology, which is perhaps the second greatest advance ever in turbulent flow, will not be examined here because it does not involve the introduction of new variables. Prandtl,69 before turning to eq 52, conjectured that u ) φ{y,a,τw,µ,F} and then applied dimensional analysis to obtain

u

() F τw

1/2

{

}

y(τwF)1/2 a(τwF)1/2 , µ µ



(53)

for which he introduced the notation

u+ ) φ{y+,a+}

(54)

which has remained the standard symbolism to this day. He conjectured that the velocity profile near the wall might be negligibly dependent on the radius of the tube a, reducing eq 54 to

u+ ) φ{y+}

(55)

which has proven to be such a good approximation for all geometries that it is known as “the universal law of the wall”. He next conjectured that the velocity gradient near the centerline would depend primarily on the fluctuations in the velocity and negligibly on the viscosity, thereby leading from eq 53 to

{ay}

uc+ - u+ ) φ

(56)

which is known as the “law of the center”. The quantity uc+ - u+ is known as the dimensionless velocity defect. Millikan70 subsequently conjectured that a range of the radius might exist for which both eqs 55 and 56 provide a good approximation and recognized that the only expression with that mathematical functionality is

u+ ) A + B ln{y+}

(57)

As might be expected, eq 57 fails seriously near the wall, predicting negative values of u+ for y+ < e-A/B, and fails moderately near the centerline as implied by the finite value of du+/dy+ at y+ ) a+. Von Ka´rma´n,71 nevertheless, integrated eq 57 over the entire cross section of the tube to obtain

um+ ≡

∫0u+ d(ar )

2

)A-

3B + B ln{a+} 2

(58)

Equation 58 serves as an approximate expression for the Fanning friction factor because um+ ≡ (2/f)1/2. Prandtl72 finally turned to eq 52, conjectured a negligible variation in τ and a negligible value of u′v′

very near the wall, and obtained by integration the asymptote

u+ ) y+

(59)

To obtain more specific expressions than eqs 54-59, Prandtl73 subsequently conjectured that u′v′ in eq 52 might be modeled in terms of l2 (du/dy)2, where l is a mixing length analogous to the mean free path of gas molecules. Then, after replacing τ by τw(1 - y/a), neglecting the viscous term, and postulating an invariant value of l, he took the square root and integrated from the centerline to obtain

uc+ - u + )

( )( )

2 l0 y 13 a a

3/2

(60)

Many pairs of values of A and B in eq 58, differing for those in eq 57, have been proposed in order to compensate for the errors resulting from the application of eq 57 for the entire cross section. Churchill and Chan74 ultimately noted that insofar as the velocity profile in the boundary layer can be represented by interpolation between eqs 57 and 59 and insofar as the effect of the wake on um+ can be represented by a fixed value ∆ these errors can be compensated for quantitatively by the expression

um+ ) A -

( )

3B D C +∆- + + + 2 a a

2

+ B ln{a+} (61)

where C and D are constants depending on the interpolative expression. Because design calculations are invariably for specified values of w ) πa2Fum or -∆P rather than for a+ or um+, a whole literature has been generated in which eqs 58 and 61 and their counterparts are approximated by expressions explicit in Re ) 2umaF/µ ) 4w/πµ. This expedient is completely unnecessary because eqs 58 and 61 may readily be solved iteratively if Re/2um+ is substituted for a+. Colebrook75 correlated his own experimental pressure drops for naturally roughened pipes in terms of an effective roughness ec, defined so as to give the same value for the friction factor for asymptotically large values of Re as does an empirical expression in generic form for a uniformly roughened pipe, namely, the following correlating equation of Nikuradse:76

(2f )

1/2

) C + B ln{a+/e+}

(62)

He then proposed the equivalent of the following expression:

() 2 f

1/2

)A-

{

a+ 3B + B ln 2 1 + Eec+

}

(63)

which, with E ) exp{(∆/B) - (3/2) - (A/B)}, reduces to eq 58 for c+ ) 0 and to eq 62 for c+ . 1/E. The choice by Colebrook of c as a variable should be recognized as truly inspired. The objective of Prandtl73 and Boussinesq77 in introducing the mixing length and the eddy viscosity, respectively, into eq 52 was to replace u′v′ with a virtual variable that had a more constrained behavior and was perhaps more predictable or easier to correlate. How-

3918

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

ever, Churchill and Chan78 found that u′v′, when expressed dimensionlessly as (u′v′)+ ≡ -Fu′v′/τw, is actually better behaved than either of these virtual variables in every respect. Churchill79 subsequently showed that (u′v′)++ ) -Fu′v′/τ, which is simply the local fraction of the total shear stress due to the turbulent fluctuaions, is an even better dimensionless variable for round tubes in terms of predictability, correlation, mathematical manipulability, and understanding. One example of the latter is provided by the relationships between these variables, namely,

(l+)2 )

(u′v′)

(

+

++

)

y 1 - + [1 - (u′v′)++]2 a

(64)

and

µt (u′v′)++ ) µ 1 - (u′v′)++

(65)

Surprisingly, the mixing length and eddy viscosity are both revealed by their relationships to (u′v′)++ to be independent of their heuristic diffusional origins. However, because 0 < (u′v′)++ < 1 for all conditions except y+ ) 0, the mixing length is seen from eq 64 to be unbounded at y+ ) a+, undermining the whole concept. The eddy viscosity is, however, seen from eq 65 to be well-behaved in a round tube for all conditions. These same conclusions apply to flow between parallel plates, but in a circular concentric annulus the eddy viscosity, the mixing length, and (u′v′)++ are all singular at some location within the fluid and negative in an adjacent finite region for all conditions, while (u′v′)+ remains well-behaved. Equation 52, when reexpressed in terms of u+, y+, a+, and (u′v′)++, becomes

(

1-

)

y+ du+ [1 - (u′v′)++] ) + + a dy

([ ( ) ] 0.7

y+ 10

-8/7

|{

} ) |)

1 0.436y+ 1 6.95y+ 8/7 -7/8 1 + (67) 0.436a+ a+ + exp -

(

has been developed by Churchill.80 The unbounded value of l+ at y+ ) a+ refutes eq 60. On the other hand, the recognition of the existence of a finite limiting value of (u′v′)++ at that location permits the derivation of the correct asymptotic expression

uc+ - u+ )

r2 a+ [1 - (u′v′)0++] 2 a

()

(

u+ ) y+ 1 -

(68)

as a replacement. The equivalent of eq 68 was derived

y+ 2a+

)

(69)

is directly applicable for a parallel-plate channel if a+ is simply replaced by b+. He conjectured that this analogy might hold for turbulent flow and found that, indeed, experimental values of u+ ) φ{y+,a+} and u+ ) φ{y+,b+} are completely congruent in the two fully developed turbulent regimes, that is, for a+ and b+ greater that 145. However, the analogy fails between a+ ) b+ ) 45, which corresponds to the onset of turbulence in a parallel-plate channel, and a+ ) b+ ) 145. The point of onset of turbulence in a round tube, namely, a+ = 65, falls in this intermediate, nonanalogous regime. It follows from eq 52 that, insofar as the analogy of MacLeod holds for u+, it also holds for (u′v′)++. Although this analogy has been confirmed within the scatter of the best available experimental data for u+ and (u′v′)++, it has a speculative origin and may not be exact. Rothfus et al.83 proposed a somewhat complementary analogy between the velocity distribution and the mixed-mean velocity in the turbulent regime of a round tube based on the relationship of those variables in the laminar regime. The velocity u+ is replaced by um+ and y+ by a+/4, or conversely. When this is applied to eq 57, the result is

um+ ) A + B ln{4} + B ln{a+}

(70)

which is simply a crude approximation for eq 58. However, when applied to eq 63, the result is

(66)

from which it may be seen that, for y+ , a+, 1 (u′v′)++ is equivalent to du+/dy+. It follows that du+/ dy+ in expressions or derivations may be replaced by (1 - y+/a+)[1 - (u′v′)++] or simply by 1 - (u′v′)++ as appropriate. This is often convenient because a generalized correlating equation for (u′v′)++, namely

(u′v′)++ )

earlier by Hinze81 based on the equivalent postulate of a finite limiting value of the eddy viscosity at the centerline. MacLeod82 noted that the exact solution for the velocity distribution in laminar flow in a round tube, namely, eq 41, when expressed in terms of the notation of Prandtl for turbulent flow, namely, as

u+ ) A -

{

a+ 3B + B ln{4} + B ln 2 1 + Eec+

}

(71)

which provides a speculative but unique expression for the effect of natural roughness on the velocity distribution. For parallel plates, this analogy requires the replacement of y+ by b+/3 and u+ by um+. Expressions for the friction factor in a round tube are often conjectured to apply to other cross sections by means of an equivalent diameter. The most commonly chosen one is the hydraulic diameter Dh ) 4Ax/Pw, where Ax is the cross section of the fluid stream and Pw the wetted perimeter. A commonly used alternative is the laminar equivalent diameter, Dl, which yields the same relationship between the friction factor and the Reynolds number as that for a round tube, namely, eq 39. Insofar as the effects of the wake and the boundary layer are ignored, the prediction of these virtual variables may be compared as follows with the exact results for a parallel-plate channel. The expression for um+ for a parallel-plate channel, as obtained by integrating eq 57 over the cross section, is

um+ ) A - B + B ln{b+}

(72)

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3919

from which it follows that for a+ ) b+

(um+)P - (um+)R )

B 2

(73)

The analogue of eq 61 for a parallel-plate channel is

um+ ) A - B -

C′ + ∆′ + B ln{b+} a+

(74)

resulting, for b+ ) a+, in

(um+)P - (um+)R )

(

) ( )

B C - C′ D + + + 2 a a

2

- ∆ + ∆′ (75)

The hydraulic diameter, which is 4b, results in

3 um+ ) A - B + B ln{2} + B ln{b+} 2

(76)

and, for b+ ) a+, in

(um+)P - (um+)R ) B ln{2}

(77)

while the laminar equivalent diameter, which is 8b/3, results in

um+ ) A -

{}

3B 4 + B ln{b+} + B ln 2 3

(78)

The concept of a heat-transfer coefficient is perhaps the greatest advance of all time in terms of prediction and correlation for thermal convection. Newton85 noted the proportionality of the rate of heat transfer between a fluid stream and a surface to their temperature difference and to the area for transfer and is thereby generally credited with this concept, which replaces four variables, namely, j, Tf, Ts, and A, with one, namely, h ) j/A(Tf - Ts). The independence of h from Tf and Ts, which may now be rationalized on the basis of the linear dependence of the differential energy balance on temperature, has been adapted with reasonable success for mass transfer but not for momentum transfer owing to the nonlinear dependence of that rate on the velocity of the fluid. This illustrates both the possible utility and limitations of an analogy. Selective examples of the choice of specific variables for convection are presented below for (1) forced external convection, (2) free and natural convection, (3) film boiling and film melting, (4) combined free and forced convection, (5) forced laminar convection, (6) forced turbulent convection, and (7) the behavior of heat exchangers. Forced External Convection. Boussinesq86 derived perhaps the first theoretical solution for forced convection by postulating potential flow (as an approximation) and showing that the differential energy balance for two-dimensional planar flow could then be expressed as follows:

and, for b+ ) a+, in

{}

4 (um+)P - (um)R ) B ln 3

(79)

Langmuir84 devised a model for convection from curved surfaces in terms of conduction across an effective film thickness δL that may be adapted and expanded in series to obtain the following expression for the friction factor for a round tube in terms of that for a parallelplate channel:

[(Re)f]R )

{

-2

}

2 ln 1 [(Re)f]P

u∞

Convection The last general topic considered herein, namely, convection, proves to be a very rich one in terms of identifying special and optimal variables.

)

(81)

Here ψ and φ are the stream and potential functions, respectively. This model, and its counterpart in polar coordinates, has been solved by Boussinesq and others (see work by Galante and Churchill87) for both isothermal and uniformly heated surfaces for a great number of geometries. For asymptotically large values of Re, these solutions all take the form

= [(Re)f]P - 1 (80)

Equation 80 predicts a smaller value of [(Re)f]R, and therefore of f for a round tube than for a parallel plate, whereas eqs 73 and 77-79 all predict larger ones. This is not a contradiction but rather the consequence of comparing a result for equal values of Re with expressions for equal values of b+ and a+. The predictions of eqs 77, 79, and 80 are relatively poor and eq 72 is slightly inferior to eq 75, which is recommended herein on structural grounds as the most accurate expression, and hence as a standard. However, the corrections are all relatively small, which may account for the continued usage of eqs 77 and 79 in practice. In summary, analyses of turbulent flow in round tubes identify u+, um+, uc+ - u+, and (u′v′)++ as optimal independent variables, y+, a+, and ec+ as corresponding optimal dependent variables, and Dh, Dl, l, δL, and µt as inferior virtual ones.

(

∂T k ∂2T ∂2T + ) ∂φ Fc ∂ψ2 ∂φ2

h

(

x u∞Fck

)

1/2

)A

(82)

where A is a constant that depends on the geometry and thermal boundary conditions, u∞ is the free-stream velocity or the local velocity over the surface, and h is either the local or the integrated-mean heat-transfer coefficient. Equation 82 may also be expressed as

( )

u∞Fcx hx )A k k

1/2

) A(Pex)1/2 ) A

( )()

u∞Fx 1/2 cµ 1/2 ) µ k 1/2 1/2 A(Rex) Pr (83)

Conduction in the direction of flow, and in the case of spheroidal bodies in the transverse direction as well, results in higher-order terms as Pex decreases. Unfortunately, this beautiful set of solutions is in significant error in all cases owing to the presence in the real world of a boundary layer, and in the case of bluff bodies due to separation as well. This example illustrates not only the great possible gain from a simplification and its exploitation but also the possible risk. Despite their functional and numerical inaccuracy, these solutions for potential flow have two lasting merits. First, they are upper bounds for h for the real behavior and, second, they are valid asymptotes for Pr ) 0 and are thereby

3920

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

useful for the construction of correlating equations for the dependence on Pr. An exact solution is not possible for the other extreme of creeping flow over a cylinder, but Langmuir,84 as mentioned in connection with eq 80, derived an approximate expression that encompasses that regime by representing convection in terms of conduction across an effective annular film of thickness δL with a logarithmic-mean area. The result is

2πkL∆T(D + 2δL - D) q) D + 2δL δL ln D

{

}

{

2

2 ln 1 + Nu∞

}

= Nu∞ + 1

(84)

(85)

Here Nu∞ ) D/δL is an expression for convection for a sufficiently large cylinder or a value of Re such that curvature is not a significant factor. It may be noted that the virtual variable δLvanishes from the final forms (eqs 80 and 85). Equation 85 is equally applicable for free convection, and its approximate generalized form, namely

Nu ) Nu∞ + Nu0

(86)

has often been proposed on purely empirical grounds for all immersed objects. A theoretical expression or a correlating equation for the thin-laminar-boundarylayer regime is usually used in eq 86 for Nu∞, and a theoretical expression or an experimentally determined value for Nu0. For a sphere, the application of the Langmuir methodology with the geometric-mean area gives directly without the necessity of series expansion

Nu ) Nu∞ + 2

(87)

which conforms to eq 86 with an exact limiting value of 2 for Nu0. Creeping flow does occur for a sphere, which has allowed the derivation of the following two asymptotic expressions for Pe f 0 and Pe f ∞, respectively,

Nu ) 2 +

Pe + ... 2

(88)

and 1/3

Nu ) 0.9145Pe

Nux ) (RexPr/π)1/2

for Pr f 0

(89)

Equation 88 is of practical value and conforms to eq 85. On the other hand, eq 89 has no regime of applicability for heat transfer because creeping flow is restricted to Re < 0.1 and the highest value of the Prandtl number for any liquid is less than 100. The resulting maximum allowable value of Pe ) 10 is insufficient to predict values of Nu above the minimum of 2. This is an example of the possible error that may result from the misapplication of a model beyond its implicit range of validity. The characteristic dimension in eqs 85 and 8789, but not necessarily in eq 86, is implied to be the diameter.

(90)

and

Nux ) 0.33872Rex1/2Pr1/3

Equation 84 may be simplified to obtain the solution of Langmuir and that result expanded in series to obtain a secondary simplification as follows:

Nu )

Thin-laminar-boundary-layer theory, which utilizes the velocity distribution for potential flow outside the boundary layer, results in the following asymptotic solutions for a flat plate:

for Pr f ∞

(91)

Thin-laminar-boundary-layer theory is applicable for convection from cylinders and spheres only for the forward region ahead of the point of separation. In addition, the aspect ratio of a cylinder, L/D, is a significant variable for almost all values owing to the disturbance of the flow by end effects. These are again implicit limitations that are not apparent from the model itself. Natural Convection. Solutions for free and natural convection have almost all utilized the clever change of variables generally credited to Boussinesq (see work by Hermann88), namely, -g - (1/F) (∂F/∂y) = -g + (gF∞/F) = -g[1 - (T/T∞)] = gβ(T - T∞), together with the complementary idealizations that the variation in F with velocity is negligible, that elsewhere in the differential model F ) F∞, and that β(Tw - T∞) is negligible with respect to unity. For an ideal gas, β ) 1/T∞; for liquids, an average experimentally determined value is used. For the thin-laminar-boundary-layer regime, a similarity transform was first identified by Pohlhausen (see work by Schmidt and Beckmann89), who used it to derive

(

hx µ2k k gF 2β∆T ∞

) {} 1/4

cµ k



Nux ) Grx1/4φ{Pr}

or

(92)

For Pr f 0, eq 89 reduces to

(

hx k2x 2 k gF c2β∆T ∞

) () 1/4

)A

cu k

1/2

or Nux ) A(GrxPr2)1/4 (93)

and for Pr f ∞ to

(

hx kµx k gF 2cβ∆T ∞

)

1/4

()

)B

cµ k

1/4

or Nux ) B(GrxPr)1/4 (94)

where Grx ≡ gF∞2β∆Tx3/µ2. For uniform heating, jw ) h(Tw - T∞) is ordinarily specified rather than Tw. Substitution of jw/h in place of Tw - T∞ in eq 92 and rearrangement give

Nux ) (Grx*)1/5φ{Pr}

(95)

where Grx* ) gF∞2βjx4/µ2k. Even though eq 95 is more explicit, the commonality of the form of eq 94 for both modes of heating is perhaps an overriding factor in favor of the latter. Free and natural convection is sometimes expressed in terms of the Rayleigh number Ra ) Gr × Pr, but Gr is preferable because it, rather than Ra, characterizes the transition from laminar to turbulent motion.

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3921

Nusselt90 speculated that h might become independent of x in the turbulent regime, which implies that

Nux ) Grx1/3φ{Pr}

(96)

Computed and experimental results indicate that eq 96 is not representative immediately following the onset of turbulence but is eventually approached as x increases. On the other hand, the subsequent speculation of Frank-Kamenetskii91 of the independence of turbulent free convection from all diffusional mechanisms, and thereby from µ and k, which requires that

Nux ) AGrx1/2Pr

(97)

is clearly refuted by the computed and experimental results. Saville and Churchill92 adapted for free convection the integral transformation of variables by Goertler65 and achieved a single generalized and rapidly converging solution for horizontal cylinders and vertical axisymmetric bodies of fairly arbitrary shape in the thinlaminar-boundary-layer regime. The solutions for cylinders are adequate for the overall heat-transfer coefficient because, as contrasted with forced convection, the angular regime of separation and the width of the wake for free convection are quite constrained. Churchill and Churchill93 accomplished an alternative, less elegant generalization for both heat and mass transfer by means of an effective value for the characteristic length. One of the most famous fluid mechanicists in the world, in the role of a visitor, predicted that the fluid motion due to natural convection inside a horizontal cylinder, heated on one vertical wall and cooled on the other, would consist of a thin boundary layer coupled to a slowly rotating central core and advised that the process be modeled analytically in those terms. Fortunately, this proposal did not have its intended effect of discouraging the already planned, first ever, numerical solution for natural convection in an enclosure by Martini and Churchill94 because the calculations and the complementary experimental work revealed a boundary layer of rapidly changing and significant thickness and an oscillating but nonrotating central core. Modeling in terms of the conventional variables of boundarylayer theory would certainly have failed, as it indeed did in the hands of others. Aziz and Hellums95 achieved the first three-dimensional solutions for laminar natural convection in an enclosure by means of representing the equations of conservation in terms of the vector potential. The stream function provides a good visual representation for the fluid motion in two dimensions, but the vector potential does not serve that role in three dimensions. Ozoe et al.96 instead presented an isometric representation of particle paths computed from the numerical solution for the velocity field. Such particle paths provide a unique insight into the circulation in inclined, baffled, and nonuniformly heated enclosures and can be interpreted as a variable in that context. Samuels and Churchill97 added a false transient term to the stream-function equation for two-dimensional natural convection in enclosures, rendering it parabolic rather than elliptic and thereby expediting the numerical process of solution. They obtained meaningful results by the use of this virtual, vanishing variable but seriously misinterpreted them in part because of the

choice of an insensitive variable for the ordinate of a graphical representative. In a logarithmic plot of Nu versus Gr, the computed values of Nu for each value of Pr less than unity fall along a different, nearly straight line and appear to extrapolate to Nu ) 1 at different values of Pr, implying, falsely, a dependence of the critical Rayleigh number on the onset of motion on Pr. This misimpression was later corrected by Ozoe et al.,97 who first plotted Nu - 1 versus Pr and then achieved an even better extrapolation by plotting (Nu - 1)/ (Nu{Pr)10-3} - 1) versus Pr. Booker99 postulated a pseudo-steady-state fluid motion for the melt in Czochralski crystallization and obtained a convergent numerical solution, which was, however, shown to be invalid by virtue of experimental observations of transient oscillations by Yi et al.100 Nakano et al.101 subsequently confirmed the prevalence of temporal oscillations in natural convection in lowPrandlt-number fluids at high Grashof numbers in enclosures by carrying out numerical calculations for a rectangular region heated from below. These results raise doubts concerning the validity of the many numerical solutions for natural convection for low-Prandtlnumber fluids in which a steady state has been postulated. As indicated by the above examples, transformations of variables, asymptotics, and simplifications have had a great impact on the development of analytical and numerical solutions for free and natural convection. On the other hand, some seemingly reasonable postulates and simplifications have proven to be misleading or erroneous. Film Condensation and Its Analogues Nusselt102 derived the following expression for the mean heat-transfer coefficient for laminar film condensation of a pure saturated vapor on a vertical isothermal plate by neglecting the sensible heat and inertia of the liquid film and the drag of the vapor and postulating fully developed (parabolic), nonrippling flow:

(

)

hmL 4 λgF2L3 ) k 3 4µk(Ts - Tw)

1/4

(98)

Here L is the vertical height of the surface, λ is the latent heat of condensation, and k, µ, and F are the conductivity, viscosity, and density of the liquid. The slightly higher coefficient observed experimentally is ordinarily attributed to rippling. Churchill103 has evaluated the other idealizations quantitatively. Emmons104 recognized that this process was analogous, by virtue of a gravitationally generated motion, to free convection, film boiling, and film melting and derived the following generalized solution corresponding to eq 98:

(

hmL FFL3 ) k KµkR∆T

)

1/4

(99)

Here F is the gravitational force per unit volume of the film, K is a numerical factor characterizing the viscous shear, and R is the increase in the mass rate of flow of the film per unit of heat transfer. He approximated these three quantities as shown in Table 1 on the basis of a conceptual analysis. In Table 1, M is the mass of the cake of solid below which the solid is melting. The factors K, F, and R result

3922

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

Table 1. Estimated Factors in the Analogue of Emmons process

K

F

R

condensation free convection boiling melting

3 4 8 6

gF gFβ∆T g(FL - Fv) 3Mg/L3

1/λ 1/c∆T 1/λ 1/λ

in exact expressions functionally in the sense of the Nusselt solution for condensation but with slightly less accurate numerical coefficients. For example, the net numerical coefficient of eq 99 for condensation based on the factors in Table 1 is 3-1/4 ) 0.760 as compared to 43/4/3 ) 0.942 for eq 98. This analogy allows experimental data for all four processes to be plotted together as hmL/k versus FFL3/KµkR∆T. When L in eq 99 is replaced with w, the corresponding mass rate of flow of the film per unit breadth results in the following alternative form:

( ) ( )

hm

Kµ2 4k2FF

1/3

)

µ 4w

1/4

(100)

which also may be used for a common graphical representation with the advantage that 4w/µ, the Reynolds number for a film based on the hydraulic diameter, might define the upper limit of validity of eqs 99 and 100, as restricted by the postulate of laminar flow. However, capillary waves as defined by M ≡ gµ4/FR3 or We ≡ 4δum2F/σ, where σ is the surface tension and δ the thickness, may occur at even lesser lengths or rates of flow. In the section on Natural Convection, it was noted that hx is presumed to become independent of x as the latter increases. The same independence would then eventually apply to hm. It is therefore tempting to speculate that the analogue of eq 99 for the turbulent regime is

(

hmL FFL3 )A k KµkR∆T

)

1/3

(101)

It is encouraging but hardly a proof that taking K, F, and R from Table 1 for laminar convection converts eq 101 to the correct functional form for turbulent free convection, namely, that of eq 96. On the other hand, the limited and very scattered experimental data for film condensation, film boiling, and film melting in the turbulent regime do not appear to confirm independence of hm from L. Reasonably accurate solutions, in the sense of eq 98 of Nusselt, have been derived for free convection, film boiling, and film melting. The approximations corresponding to eqs 100 and 101, together with the numerical coefficients of Table 1, represent necessary sacrifices to achieve commonality and thereby insight. In such a qualitative but not a quantitative sense, these values represent one of the greatest generalizations ever achieved for heat transfer. It should be possible to construct similar conceptual generalities for other completely different topics. Assisting Free and Forced Convection. This behavior has inspired the choice of many combined or equivalent dimensionless variables with three curious consequences. First, all of these approaches lead to a result in the form of the generalized correlating equation proposed by Churchill and Usagi,105 namely

Nun ) NuNn + NuFn

(102)

Here NuN represents free (natural) convection, NuF represents forced convection, and n is an arbitrary exponent. Second, every value of n determined by seemingly reasonable conjecture has proven to be wrong. Third, this particular application of the ChurchillUsagi equation is the only one for which a value of the combining exponent n other than (1 has ever been derived theoretically. A skeletal description of the various conjectured combinations is as follows. McAdams106 proposed utilizing the higher values of NuN and NuF as a conservative approximation for combined free and forced convection in general. This is equivalent to letting n f ∞. Lemlich and Hoke107 proposed the postulate of an equivalent free-stream velocity to represent free convection from a horizontal cylinder but did not specify a combining rule. Kirscher and Loos108 proposed the use of Re + (Gr/2)1/2 for various bodies, which is the equivalent of n ) 2. Bo¨rner109 proposed adding Re and (Gr/2)1/2 vectorially, which is equivalent to n ) 4. Jackson and Yen110 chose n ) 4 on the basis of the additivity of the work of free and combined convection. Sparrow and Gregg111 derived the asymptotic expression

(

Grx Nu ) NuR 1 + φ{Pr} Rex2

)

(103)

for the perturbation of free convection on forced convection for a vertical plate. A combining exponent of n ) 4 can be inferred from eq 103. Eschgy112 derived the corresponding asymptotic expression for the perturbation of free convection due to forced convection

(

)

Rex Nu ) Nu 1 + φ{Pr} Grx1/3

(104)

from which it may be inferred that n ) 2. Acrivos113 chose n ) 4, without any rationalization, to fit his precise computed values for the forward point of a sphere. Finally, Churchill111 plotted precise computed values for isothermal and uniformly heated vertical plates in the thin-laminar-boundary-layer regime in the critical form of arithmetic split coordinates as recommended by Churchill and Usagi105 and thereby found that n ) 3 was not only better than n ) 2 or n ) 4 but also resulted in an essentially exact representation. Churchill114 then showed that experimental data for all geometries and conditions were well represented in terms of Nu - Nu0 and n ) 3. Ruckenstein115 was inspired by this result to derive an approximate analytical solution for a vertical plate, and by adding the velocity fields rationalized, the value of n ) 3. Earlier equivalent solutions for other processes were ultimately recognized as relevant. For example, even earlier, Martinelli and Boelter116 derived an analytical solution for free and natural convection in an open-ended vertical tube in the laminar regime, by the somewhat suspect technique of integral-boundary-layer theory, that can be interpreted to have the form of eq 102 with n ) 3, and Kitaura et al.117 derived a solution for a rotating and vibrating sphere in a forced flow in terms of the third-power mean of the individual expressions for Nu by adding the velocity fields, as did Ruckenstein later but independently. It may be concluded that assisting

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3923

convection of all types can be represented by such a combined independent variable. The lessons to be drawn are that (1) seemingly reasonable rationalizations may be wrong, (2) sensitive methods must be used to test such rationalizations, and (3) applicable solutions in the literature of indirectly related topics may be overlooked. Internal Forced Convection. Heat transfer to or from a fluid stream passing through a round tube is the most important thermal process in the chemical and petroleum industry other than distillation. The heattransfer coefficient for this process is usually expressed in terms of the temperature difference between the wall of the tube and the mixed-mean temperature of the fluid stream defined as



[

]

uFc Tm ≡ 0 T dR2 = (uFc)m 1



( )

u T dR2 (105) 0 um 1

The term Fc/(Fc)m is usually dropped for convenience even though it always varies somewhat with temperature as does also u/um. On the other hand, the neglect of the variation of u/um due to drag of the wall, just as the postulate of plug flow in chemical reactors, is never justifiable in convection. Graetz118 in 1883 and 1885 derived a solution in the form of an infinite series for developing convection in fully developed laminar flow in a round tube following a step in wall temperature. His solution may be expressed functionally as

Nu0 ) φ{Gz}

2Nu[1 - R2]ψ )

2Nu[1 - R2] )

(107)

This solution is based on the linearization of the velocity profile, the neglect of the effect of curvature on thermal conduction, and the replacement of the mixed-mean temperature with the inlet temperature, which is postulated to prevail across most of the tube and thereby to be applicable as a boundary condition for y f ∞. For x f ∞, Nu f 3.657, which is the solution for fully developed convection, for which x, F, c, um, and w are no longer variables. Fully developed convection, which was mentioned in the previous paragraph, is a more subtle concept than fully developed flow. The latter is described mathematically be setting the derivatives of the velocity to zero, while the former is described mathematically by setting the derivative of h and/or (Tw - T)/(Tw - Tm) to zero and neglecting longitudinal conduction. For a uniform wall temperature, setting (∂/∂x) [(Tw - T)/(Tw - Tm)] ) 0 leads to (∂T/∂x)/(∂Tm/∂x) ) (Tw - T)/(Tw - Tm), and together with the overall energy balance, jw ) (aumF/ 2)/(∂Tm/∂x) allows reductions of the time-averaged energy balance in laminar flow to

]

(108)

1 d dφ R R dR dR

[

]

(109)

where here φ ) (Tw - T)/(Tw - Tm). The solution of this model leads to Nu ) 48/11. The choice of Nu, R, and Ψ or φ as variables is critical to these formulations and solutions. For fully developed turbulent convection, the above conditions of uniform wall temperature and uniform heating again result in a similarity transformation, but this time it is more convenient to use the conventional notation of Prandtl,69 resulting in

(1 + γ)R )

dT+ + (T′v′)+ dy+

(110)

where T+ ) k(τwF)1/2(Tw - T)/µjw, (T′v′)+ ) FcT′v′/jw, and γ ) (j/jw)(τw/τ) - 1. For a uniform wall temperature

j 1 ) jw R

πD2Fu

NuD ) 1.167Gz1/3

[

where R ) r/a, Nu ) 2ajw/k(Tw - Tm), and ψ ) (Tw T)/(Tw - Tc). The centerline temperature is utilized here because ψ is then known at R ) 0. The value of 3.657 in the previous paragraph is a result of the solution of eq 108. For uniform heating, setting (∂/∂x) [(Tw - T)/ (Tw - Tm)] ) 0 results in ∂T/∂x ) ∂Tw/∂x ) ∂Tm/∂x and leads to

(106)

where Gz ) wc/kx ) mc/4kx. This dimensionless group is now known as the Graetz number. Dimensionless analysis of a simple listing of these variables would result in two groups such as D/x and wc/kD or umFcD/k. The combination in eq 106 is a consequence of constraints arising from the differential energy balance. The Graetz solution involves eigenvalues and eigenfunctions that must be evaluated separately and numerically and converges very slowly. As a consequence, Le´veque119 derived the following asymptotic solution for the inlet, where the series of Graetz converges most slowly:

1 d dψ R R dR dR

+

( ) +

∫0R uu + TT + 2

m

dR2

(111)

w

while for uniform heating

1 j ) jw R

∫0R

2

( )

u+ dR2 um

(112)

The heat flux density ratio j/jw, which is discussed in detail by Churchill and Balzhiser,121 was apparently first identified as a parameter by Reichardt,120 who suggested the use of the quantity γ as a perturbation representing the deviation of the heat flux density from the shear stress ratio. Equation 110 is the optimal starting point for theoretical modeling, but other historical approaches and the associated variables will be examined before turning back to this expression. Nusselt122 was apparently the first to apply dimensional analysis to a listing of variables for turbulent forced convection, but he erroneously postulated that h was proportional to the product of arbitrary powers of each of the variables and thereby obtained

NuD )

( )( )

DumF hD )A k µ

n

cµ k

m

) A(Re)nPrm (113)

rather than the less explicit but nevertheless correct expression for this set of variables, namely

NuD )

{

}

DumF cµ hD ) φ(Re,Pr) )φ , k µ k

(114)

The harm done by this pioneering effort is immeasurable because his postulate and result have been widely adopted, not only in heat transfer but also throughout chemical engineering. It has been known for some time that power dependences and their products occur only

3924

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

for asymptotic behavior but that recognition has not permeated sufficiently to prevent such expressions as eq 113 from still appearing in our textbooks, handbooks, and even current articles. Nusselt compounded the error by inferring from his own experimental data for turbulent convection that n and m had the same value, leading to

(

)

DumFc NuD ) A k

m

m

) A(ReDPr) ) A(PeD)

m

}

D(τwF)1/2 cµ hD ) φ{a+,Pr} (116) ) NuD ) φ , k µ k The postulate of independence of h from D for very large values of the Prandtl number reduces eq 116 to

(

)

{}

NuD µ cu h ) )φ ) φ{Pr} (117) 1/2 k (τ F)1/2 k ReD(f/2) w The postulate of large values of Pr further implies that the resistance to heat transfer is confined to the region near the wall where the turbulent shear stress is proportional to y3, leading to

Nu ) A(Re)(f/2)1/2Pr1/3

(118)

Equations 117 and 118 are well-confirmed by experimental data for large values of Pr, whereas the elimination of D from eq 114 leads to a relationship with no experimentally observed range of validity. Inclusion of τw in place of F rather than in place of um in eq 114 results in

{

}

τwD cµ hD , ) φ{ReD(f/2),Pr} (119) ) NuD ) φ k µum k which reduces, for the postulate of independence from D and c, to

NuD hµum ) )A kτw ReD(f/2)

(120)

Equation 120 is known to be valid only for the singular condition of Prt ) Pr = 0.87. This exercise in speculative dimensional analysis reveals that many different dimensionless groupings result from the initial choice of dimensional variables as well as from the postulate of independence from particular ones. Only experimental data or very general numerical solutions reveal the validity or optimality of these choices. The analogue of the eddy viscosity of Boussinesq77 is the eddy conductivity, which converts eq 110 to

(1 + γ)R ) +

dT dy+

( )

kt dT+ 1+ ) + k dy kt cµ µt µt Pr dT+ 1+ ) + 1+ (121) cµt k µ µ Prt dy

[

( )( )]

(1 + γ)R[1 - (T′v′)++] )

[ () ]

On the other hand, Churchill79 introduced (T′v′)++ )

dT+ dy+

(122)

From eqs 65, 121, and 122, it follows that

Prt (u′v′)++[1 - (T′v′)++] ) Pr [1 - (u′v′)++](T′v′)++

(115)

Other variables could have been chosen to represent turbulent convection with equal justification. For example, the choice of τw rather than um as an independent variable results in

{

pcT′v′/j ) (T′v′)+(jw/j) in eq 110, which results in

(123)

Also, eq 122 may be reexpressed in terms of (u′v′)++ and Prt/Pr as

(

(1 + γ)R ) 1 +

)

++ Pr (u′v′) dT+ Prt 1 - (u′v′)++ dy+

(124)

Abbrecht and Churchill123 found experimentally that Prt is a function only of (u′v′)++ and Pr and that it is the same function for both round tubes and parallel plates and for all modes of heating at the wall(s) that lead to fully developed convection. This same generalized functionality has been implied, at least for the turbulent core, by the expressions for Prt derived by Yahkot et al.124 and others by renormalization group theory. The relative simplicity of eq 124, which only requires a correlating equation for (u′v′)++ and a theoretical expression or correlating equation for Prt in order to determine NuD by numerical integration for all Pr and Re, is a direct consequence of the judicious choice of variables. Analogies are often effective in eliminating a variable common to both processes. Reynolds,125 in deriving the most famous and long-lasting of analogies, speculated that momentum and energy might be transferred from a stream of fluid to the wall of the tube through which it is flowing by the same net mass flux of eddies and then eliminated the virtual variable between the two expressions for the rates of transport to obtain

j w um cτw(Tm - Tw)

)1

or

Nu ) (Re)(Pr)(f/2)

(125)

It may be noted that the density, thermal conductivity, and viscosity of the fluid as well as the diameter of the tube are absent as explicit variables. Prandtl126 interposed a viscous boundary layer of effective thickness δ, across which the transport of momentum and energy was postulated to occur only by molecular diffusion, resulting in a relationship that can expressed as

Nu )

(Re)(Pr)(f/2) 1 + δ+(Pr - 1)(f/2)1/2

(126)

Jakob127 wrote that “there exists a whole literature about (the equivalent of δ+)”. Equation 126 represents a great advance over eq 125 in two respects. First, it reveals that the applicability of eq 125 is limited to Pr ) 1 because of the neglect of a diffusional resistance near the wall. Second, eq 126 predicts a coupled dependence of Nu on Pr and Re (including that of f ). Strangely, this coupling has been ignored or overlooked to the present day in many correlating equations. Equation 126 itself is in error at large values of Pr because of the postulate of negligible transport by

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3925

turbulence within the viscous boundary layer and at small values of Pr because of the neglect of thermal conduction in the turbulent core. Many complex analogies have been derived to correct for these latter deficiencies of the Prandtl analogy. Two will be examined here because of their consequences with respect to the objective of this investigation. Reichardt120 derived an analogy for a uniform wall temperature by eliminating dy between the exact differential energy and momentum balances and making a few seemingly reasonable idealizations in order to be able to integrate in closed form. Reichardt presented his analogy in terms of a number of algebraic expressions, tabulations, and graphs. An equivalent overall formula, corrected by substituting the exact algebraic expression based on direct numerical simulation for his erroneous expression for the turbulent shear stress very near the wall, is, per Churchill,128

Here the subscript mR4 indicates the integrated-mean value over R4 and wmR4 the integrated-mean value weighted by 1 - (u′v′)++. The implied independence of Prt ) Pr from y+ in the derivations of eqs 130 and 131 is supported by renormalization group theory. Churchill et al.132 were unsuccessful in developing correlating equations for intermediate values of Pr by combining eqs 129-131 in the form of eq 102. However, they subsequently recognized that eq 126 could be interpreted as

Nu )

1 Prt 1 Prt 1 + 1Pr Nu1 Pr Nu∞

( )

(

)

(132)

where here Nu is the analogue of eq 130 for uniform wall temperature. When eq 132 was reexpressed as

Nu ) 1 (1 + γ)m 2 Tm uc Prt 13.62 2 + Re f Tc um Pr Re f

( )( )( )

( ) ( )( ) ( 1/2

Tm Prt Tc Pr

1/3

1-

)

Prt Pr

(127)

Equation 127 is noteworthy in that it just, as the Reynolds analogy, is free of any explicit empiricism. (The coefficient 13.62 was evaluated on the basis of direct numerical simulations for the turbulent shear stress very near the wall.) The applicability of eq 127, just as eq 126, is obviously limited to Pr g1. Colburn129 followed an entirely different path. He noted the functional and numerical similarity of empirical correlating equations for the friction factor and the Nusselt number, eliminated their common leading coefficient, and thereby obtained

Nu .023 f ) ) 0.8 1/3 2 (Re)Pr Re

(128)

The exponent of 1/3 for the Prandtl number was not chosen for theoretical reasons but merely as a convenient compromise between experimentally determined values of 0.3 and 0.4 for cooling and heating of the fluid, respectively. He gave the grouping Nu/Re(Pr)1/3 the name j factor, which may have contributed to the continued widespread usage of eq 128, although, as may be noted by comparison with eq 127, it is in error functionally in almost every respect. Surprisingly, in view of its experimental antecedents, eq 128 has been shown by Churchill and Zajic130 and others to be in serious error numerically as well as functionally for all values of Pr. Starting with eq 124, Heng et al.131 derived and evaluated numerically the following solutions for turbulent convection with uniform heating for the only three conditions for which the uncertainty in the value of Prt can be avoided:

Nu0 ≡ Nu{Pr)0} ) 8/(1 + γ)mR4 Nu1 ≡ Nu{Pr)Prt} ) Re(f/2)/(1 + γ)wmR42 and

Nu∞ ≡ Nu{Prf∞} ) 0.07343(Pr/Prt)1/3Re(f/2)1/2 (131)

(129) (130)

Nu - Nu1 ) Nu∞ - Nu1

1 Nu∞ Pr 1+ / -1 Nu1 Prt

(

)

(133)

they recognized that it could be interpreted as a degenerate case of the staggered extension of eq 102 for three regimes by Churchill and Usagi.133 It follows that

Nui ) Nu1(Nu∞ - Nu1)

(

)

Pr - 1 /Nu∞ Prt

(134)

is a virtual intermediate asymptote. Without this latter identification, it is doubtful if the existence of a sigmoidal transition from Nu to Nu∞ with an intermediate point of inflection would ever have been discovered. Furthermore, this identification leads to the speculative derivation of an analogue of eq 132 for Pr < Prt. The absence of any allusion to geometry or the mode of heating in eq 132 and its analogue for Pr < Prt then leads to the conjecture that they might be directly applicable for all such conditions, and this has proven to be the case. Although these two expressions based in the analogy of Reichardt120 provide a better representation for experimental and computed values than any prior ones, some minor numerical discrepancies were noted by Churchill and Zajic,130 particularly for Pr = 0.01 and 10. These discrepancies were found to be, as expected, a consequence of the idealizations made by Reichardt. They further found that simply modifying eq 132 by substituting 1 - (Prt/Pr)2/3 for 1 - (Prt/Pr), on the basis of an analogy of Churchill134 (see work by Churchill and Zajic130 for an updated derivation), leads to an essentially exact representation for all values of Re and Pr, all geometries, and all thermal boundary conditions. In terms of the objectives of this investigation, expressing Nu as a function of the new variables Pr/Prt, Nu0, Nu1, and Nu∞ has proven to be a great advance over expressing it in terms of Re, Pr, geometry, and the mode of heating. The novelty and generality of this approach suggest that analogous applications might be found for other aspects of convection and for chemical engineering processes in general. Heat Exchange. Strenger et al.135 investigated both experimentally and computationally the behavior of a

3926

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

double-spiral heat exchanger with rectangular passages and with electrical heating and reversal of the flow at the inner core. When they plotted their results in the form of wc(T2 - T1)/q = (T2 - T1)/(T3 - T2) versus UA/ wc, the number of thermal transfer units, they were surprised by the appearance of a maximum in the ordinate and therefore in T2, as a function of the abscissa and therefore of w. Here, T1 is the temperature of the fluid entering the exchanger at the periphery, T2 is its temperature entering the core, and T3 is its temperature leaving the core. U is the mean value of the overall coefficient of heat transfer between the two fluid streams, A the area for heat exchange, q the external thermal input at the core, and w the mass rate of flow. The source of the maximum in T2 is the crossover in the temperature of the entering and outer exiting streams in adjacent passages. A maximum in T2 had never before been identified for double-spiral heat exchangers because, in the conventional plot of  ) (T2 - T1)/(T3 - T1), called the effectiveness, the maximum is too flat to be readily discernible, while the other common dependent variable Z ) wc(T2 - T1)/ U(∆T)lm, called the log-mean correction factor, does not go through a maximum. These latter choices are conventional because T1 and T3 are usually specified rather than T1 and q for nonreversed flow. The obvious lesson is that all seemingly equivalent groups of variables should be tested because they may uniquely reveal important aspects of behavior. Summary and Conclusions Most improvements in modeling and correlation are associated with the improved choice of individual variables or of dimensionless groups of variables. Although this association is generally recognized implicitly or after the fact, insufficient attention is given in many investigations to the possibilities inherent in the active pursuit of alternative choices and to comparisons of the consequences. Techniques for dimensional analysis of both models and lists of variables are straightforward and nearly infallible. Nevertheless, their misapplication in the past, as illustrated by eq 113, is a continuing legacy of our profession. Furthermore, these techniques are often misapplied or incompletely exploited even today. As amply illustrated herein, the deletion of individual variables is the most common source of improvement in the choice of dimensionless groups and thereby in modeling and correlation. Such deletions often have a conceptual origin, such as with the deletion of the density for slow flow on the basis of negligible inertia, but may also be undertaken on a purely speculative basis. Dimensionless analysis should be reapplied or else the prior result of dimensional analysis reduced appropriately following a deletion. The selective deletion of terms, as contrasted with variables, such as in the derivation of the model for a thin-laminar-boundary layer, is almost always done on a conceptual basis. The reapplication of dimensionless analysis is again essential. The possibilities and consequences of using alternative variables of equal validity is often overlooked. Equations 116 and 119, and thereby their reduced forms, such as eqs 118 and 120, are examples of the productive application of this technique. The deletion of a variable or the postulate of its invariance may be justifiable for most conditions but

nevertheless lead to the oversight of an important aspect of behavior for other conditions. Examples cited herein include the effect of the variation of the surface tension of water in a porous media with temperature and the effect of the variation of the density in the extended model for thermal conduction. The deletion of the viscosity was shown to lead to some valid simplifications for flow and some invalid ones for convection. The introduction of new variables with a simple physical sense, such as h, β, (u′v′)++, and Prt, is relatively straightforward, but that of virtual variables such as the eddy viscosity, the mixing length, the boundary layer thickness of Prandtl, the film thickness of Langmuir, and the effective roughness of Colebrook are not and for that reason remain a testament to the imagination of their originators. The grouping of variables in a solution for a limiting or idealized case may have extended utility. The solution of Maxwell for conduction in dispersions is a prime example. The analogies of Emmons, Langmuir, Reynolds, Prandtl, Reichardt, MacLeod, and Rothfus et al. provide generalized groupings of variables at the expense of some idealizations. This methodology has yet unexplored potential in many fields. Last, a somewhat differing procedure for the choice of variable may be noted for each of the several technical topics chosen herein for illustration. Thermal conduction is a mature field with a well-defined and well-explored mathematical structure. Even so, some aspects, such as the effect of variable physical properties, remain unresolved in part and therefore susceptible to improvement by the innovative choice of variables. Thermal radiation is so complex intrinsically that, despite considerable work by the giants of physics, approximate modeling and thereby the choice of new variables, either real or virtual, remain open subjects. The continued use of approximate models and variables such as global kinetics, plug flow, and space time in chemical reaction engineering, even though the more exact representations are known and computationably tractable, can only be attributed to technological inertia and a preference for elegance or extreme simplicity at the expense of accuracy. The opportunities for improvement in this respect are obvious. Fluid flow provides perhaps the best illustrations of the merits of the clever choice of appropriate reduced models for different regimes, the use of essentially exact models for limiting conditions, and the use of numerical solutions for intermediate behavior. However, opportunities remain in many aspects such as turbulent flow for which exact numerical simulations have slackened after a promising beginning. Convection provides perhaps even greater opportunities than fluid flow because of its critical dependence on flow and the superposition of parametric variability. Although the examples presented herein are all for particular classical topics with which the author is most familiar, not only are the concepts presumed to be relevant for new emerging topics but also the opportunities therein are unlimited by virtue of the absence of a historical structure. Literature Cited (1) Foss, A. S. Critique of Chemical Process Control Theory. AIChE J. 1973, 19, 209. (2) Morari, M.; Arkun, Y.; Stephanopolous, G. Studies in The Synthesis of Control Structures for Chemical Processes: Part 1.

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3927 Formulation of the Problem. Process Decomposition and the Classification of the Control Tasks. Analysis of the Optimizing Control Structures. AIChE J. 1980, 26, 220. (3) Edgar, T. F.; Himmelblau, D. M.; Lasdon, L. S. Optimization of Chemical Processes, 2nd ed.; McGraw-Hill: New York, 2001. (4) Rayleigh, Lord (J. W. Strutt). The Principle of Similitude. Nature 1915, 95, Mar, 66; 1915, 95, Aug, 644. (5) Hellums, J. D.; Churchill, S. W. Simplification of the Mathematical Description of Boundary and Initial Value Problems. AIChE J. 1964, 10, 40. (6) Klinkenberg, A.; Mooy, H. H. Dimensionless Groups in Fluid Friction, Heat, and Material Transfer. Chem. Eng. Prog. 1948, 44, 17. (7) Van Dusen, M. S. Note on the Theory of Heat Conduction. Bur. Stand J. Res. 1930, 4, 753. (8) Rayleigh, Lord (J. W. Strutt). On Conduction of Heat in a Spherical Mass of Air by Confined Walls of Constant Temperature. Philos. Mag. 1899, 34, 314. (9) Brown, M. A.; Churchill, S. W. Finite-Difference Computation of the Wave Motion Generated in a Gas by a Rapid Increase in the Bounding Temperature. Comput. Chem. Eng. 1999, 73, 351. (10) Churchill, S. W. Heat Transfer Rates and Temperature Fields for Underground Storage Tanks. Soc. Pet. Eng. J. 1962, 2 (No. 1), Mar, 28. (11) Churchill, S. W.; Chang, S. M. Approximations and Bounds for Conduction from Plane and Solid Surfaces to Infinite Surroundings. Lett. Heat Mass Transfer 1981, 6, 79. (12) Schoenherr, R. U.; Churchill, S. W. The Use of Extrapolation for the Solution of Heat Transfer Problems by FiniteDifference Methods. J. Heat Transfer, Trans. ASME 1970, 92C, 564. (13) Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solids, 2nd ed.; Clarendon Press: Oxford, 1959; p 282. (14) Gupta, J. P.; Churchill, S. W. A Model for the Migration of Moisture During the Freezing of Wet Sand. Chem. Eng. Prog. 1973, 69 (No. 135), 192. (15) Churchill, S. W.; Gupta, J. P. Approximations for Conduction with Freezing or Melting. Int. J. Heat Mass Transfer 1997, 20, 1251. (16) Churchill, S. W. The Thermal Conductivity of Dispersions and Packed BedssAn Illustration of the Unexploited Potential of Limiting Solutions for Correlation. In Advances in Transport Processing; Mujumdar, A. S., Mashelkar, R. A., Eds.; Wiley Eastern: New Delhi, India, 1985; Vol. 4, p 394. (17) Maxwell, J. C. A Treatise on Electricity and Magnetism; Clarendon Press: Oxford, 1873; Vol. 1, p 365. (18) Molerus, O. Heat Transfer in Moving Beds with a Stagnant Interstitial Gas. Chem. Eng. Sci. 1997, 40, 4151. (19) Hottel, H. C.; Sarofim, A. F. Radiative Transfer; McGrawHill Co.: New York, 1960. (20) Mie, G. Beitra¨ge zur Optik tru¨ber Medien. Ann. Phys. 1908, 25, 377. (21) Hartel, W. Zur Theorie der Lichtstreuung durch tru¨be Schichten, besonders Tru¨bgla¨ser. Licht 1940, 40, 141. (22) Chu, C.-M.; Churchill, S. W. Numerical Solution of Problems in Multiple Scattering of Electromagnetic Radiation. J. Phys. Chem. 1955, 59, 855. (23) Schuster, A. Radiation through a Foggy Atmosphere (in German). Astrophys. J. 1905, 22, 1. (24) Chu, C.-M.; Churchill, S. W. Numerical Solution of Problems in Multiple Scattering of Electromagnetic Radiation. J. Phys. Chem. 1955, 59, 855. (25) Chu, C.-M.; Churchill, S. W.; Pang, S. C. A Variable-Order Diffusion-Type Approximation for Multiple Scattering. In Electromagnetic Scattering; Kerker, M., Ed.; MacMillan Co.: New York, 1963; p 905. (26) Larkin, B. K.; Churchill, S. W. Heat Transfer by Radiation Through Porous Insulations. AIChE J. 1959, 5, 467. (27) Chin, J. H.; Churchill, S. W. Radiant Heat Transfer in Packed Beds. AIChE J. 1963, 9, 25. (28) Cimini, R.; Chen, J. C. Experimental Measurements of Radiant Transmission through Packed and Fluidized Beds. Exp. Heat Transfer 1987, 1, 145. (29) Elsasser, W. M. Heat Transfer by Infrared Radiation in the Atmosphere, 2nd ed.; Harvard Meteor Studies Series 6; Harvard University: Cambridge, MA, 1942. (30) Chin, J. H.; Churchill, S. W. Incorporation of the Error Function Absorption Coefficient in the Transport Equation. Q. Appl. Math. 1960, 18, 93.

(31) Friedman, M. H.; Churchill, S. W. The Absorption of Thermal Radiation by Fuel Droplets. Chem Eng. Prog. Symp. Ser. 1965, 61 (No. 57), 1. (32) Chin, J. H.; Churchill, S. W. Radiant Heat Transfer Through the Atmosphere. Chem. Eng. Prog. Symp. Ser. 1960, 56 (No. 30), 117. (33) Christiansen, J. A. On the Reaction between Hydrogen and Bromine. Mat.-Fys. Medd.sK. Dan. Vidensk. Selsk. 1919, 1, 1. (34) Zel’dovich, Ya. B. The Oxidation of Nitrogen in Combustion and Explosions. Acta Physicochim. URSS 1946, 21, 577. (35) Pfefferle, L. D. Stability, Ignition, and Pollutant Formation Characteristics Mathematically Stabilized Plug-Flow Burner. Ph.D. Thesis, University of Pennsylvania, Philadelphia, PA, 1983. (36) Rabitz, H.; Dougherty, E. D. A Computational Algorithm for Green’s Function Method of Sensitivity Analysis in Chemical Kinetics. Int. J. Chem. Kinet. 1979, 11, 1237. (37) White, C. W., III. The Integration of Stiff Differential Equations in Chemical Reactor Modeling. Ph.D. Thesis, University of Pennsylvania, Philadelphia, PA, 1983. (38) Kabel, R. L. Rates. Chem. Eng. Commun. 1981, 9, 15. (39) Churchill, S. W.; Pfefferle, L. D. The Refractory Tube Burner as an Ideal Stationary Reactor. Int. Chem. Eng. Symp Ser. 1984, No. 87, 279. (40) Collins, L. R.; Churchill, S. W. The Effect of Laminarizing Flow in Postflame Reactions in a Thermally Stabilized Burner. Ind. Eng. Chem. Res. 1990, 29, 456. (41) Theile, E. W. Relation between Catalytic Activity and Size of Particles. Ind. Eng. Chem. 1939, 31, 916. (42) Aris, R. On Shape Factors for Irregular Particles. 1. The Steady-State Problem. Diffusion and Reaction. Chem Eng. Sci. 1956/57, 6, 262. (43) Churchill, S. W. A Generalized Expression for the Effectiveness Factor of Porous Catalyst Pellets. AIChE J. 1977, 23, 208. (44) Damko¨hler, G. Einfluss von Diffusion, Stro¨mung und Wa¨rmetransport auf die Ausbeute Chemisch-Technischen Reaktionen. Chem.-Ing.-Tech. 1937, 3, 430. (45) Zel’dovich, Ya. B. On the Theory of Reactions on Powders and Porous Substances. Acta Physicochim. URSS 1939, 10, 583. (46) Hatta, S. On The Absorption Velocity of Gases by Liquids. II. Theoretical Consideration of Gas Absorption due to Chemical Reaction. Tech. Rep.sTohoku Imp. Univ. 1932, 10, 119. (47) Dean, W. R. Note on the Motion of Fluid in a Curved Pipe. Philos. Mag., Ser. 7 1927, 4, 208; Streamline Motion of Fluid in a Curved Pipe. Philos. Mag., Ser. 7 1928, 5, 673. (48) Truesdell, L. C., Jr.; Adler, R. J. Numerical Treatment of Fully Developed Laminar Flow in Helically Coiled Tubes. AIChE J. 1970, 16, 1010. (49) Manlapaz, R. L.; Churchill, S. W. Fully Developed Laminar Flow in a Helically Coiled Tube of Finite Pitch. Chem Eng. Commun. 1980, 7, 54. (50) Dealy, J. Private communication. (51) Ergun, S. Fluid Flow through Packed Columns. Chem. Eng. Prog. 1952, 48, 89. (52) Churchill, S. W. The Interpretation and Use of Rate Data. The Rate-Process Concept, revised printing; Hemisphere Publishing Corp.: Washington, DC, 1979; Chapter 10, Figures 6-8. (53) DeNevers, N. Fluid Mechanics; Addison-Wesley: Reading, MA, 1970; p 384. (54) Niven, R. K. Insight into the Ergun and Wen and Yu Equations for Fluid Flow in Packed and Fluidised Beds. Chem. Eng. Sci. 2002, 57, 527. (55) MacDonald, T. F.; El-Sayer, M. S.; Mow, K.; Dullien, F. A. L. Flow through Porous MediasThe Ergun Equation Revisited. Ind. Eng. Chem. Fundam. 1979, 18, 199. (56) Molerus, O.; Schweinzer, J. Resistance of Particle Beds at Reynolds Numbers up to Re ≈ 104. Chem. Eng. Sci. 1989, 44, 1071. (57) Molerus, O. A Coherent Representation of Pressure Drop in Fixed Beds and of Bed Expansion for Particulate Fluidized Beds. Chem. Eng. Sci. 1980, 35, 1331. (58) Roscoe, R. The Flow of Viscous Fluids Round Plane Obstacles. Philos. Mag. 1951, 40, 338. (59) Flu¨gge, S., Ed. Stro˜mungs-mechnik. Handbuch der Physik, 9th ed.; Springer-Verlag: Berlin, 1960; Vol. 3; p 425. (60) Stokes, G. C. On The Effect of The Internal Friction in the Motion of Pendulums. Trans. Cambridge Philos. Soc. 1851, 9, 8. (61) Blasius, H. Grenzschichten in Flu¨ssigkeiten mit kleiner Reibung. Z. Math. Phys. 1908, 56, 1.

3928

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002

(62) Imai, I. Second Approximations to the Laminar Boundary Flow over a Flat Plate. J. Aerosol Sci. 1957, 24, 155. (63) Van Dyke, M. Perturbation Methods in Fluid Mechanics, annotated ed.; Parabolic Press: Stanford, CA, 1975; p 138. (64) Newton, I. Principia, The Motion of Bodies; S. Pepys: London, 1686; English translation of 2nd ed.; Vol I, 1713; A. Motte, 1729; revised translation by F. Cajori, University of California Press: Berkeley, CA, 1966. (65) Goertler, H. A New Series for the Calculation of Steady Laminar Boundary Layer Flow. J. Math. Mech. 1957, 6, 1. (66) Churchill, S. W. A Theoretical Structure and Correlating Equation for the Motion of Single Bubbles. Chem. Eng. Process. 1989, 26, 269. (67) Churchill, S. W. Viscous Flows. The Practical Use of Theory; Butterworth Publishers: Stoneham, MA, 1988. (68) Reynolds, O. On the Dynamical Theory of Incompressible Viscous Fluid and Determination of The Criterion. Philos. Trans. R. Soc London 1895, 186, 123. (69) Prandtl, L. U ¨ ber die ausgebildete Turbulenz. Verhdl. 2. Kongr. Techn. Mech., Zurich 1926, 62. (70) Millikan, C. B. A Critical Discussion of Turbulent Flows in Channels and Circular Tubes. Proceedings of the Fifth International Congress on Applied Mechanics, Cambridge, MA, 1939; John Wiley & Sons: New York, 1939; p 386. (71) von Ka´rma´n, Th. Mechanische A ¨ hnlichkeit und Turbulenz. Proc. Third International Congr. Appl. Mech. Stockholm 1930, Part 1, 94. (72) Prandtl, L. Bericht u¨ber Untersuchungen zur ausgebildeten Turbulenz. Z. Angew. Math. Mech. 1925, 5, 136. (73) Prandtl, L. Neuere Ergebnisse der Turbulenzforschung. Z. Ver. Dtsch. Ing. 1933, 77, 105. (74) Churchill, S. W.; Chan, C. Improved Correlating Equations for the Friction Factor for Fully Developed Turbulent Flow in Round Tubes and between Identical Parallel Plates, both Smooth and Rough. Ind. Eng. Chem. Res. 1994, 33, 2016. (75) Colebrook, C. F. Turbulent Flow in Pipes with Particular Reference to the Transition Region Between Smooth and Rough Pipe Laws. J. Inst. Civ. Eng. 1938-1939, 11, 133. (76) Nikuradse, J. Stro¨mungsgesetze in rauhen Rohren. Ver. Dtsch. Ing. Forschungsh. 1933, No. 361. (77) Boussinesq. J. Essai sur la the´orie des eaux courantes. Mem. presents divers savants. Acad. Sci. Inst. Fr. 1827, 23, 1. (78) Churchill, S. W.; Chan, C. Theoretically Based Correlating Equations for the Local Characteristics of Fully Turbulent Flow in Round Tubes and between Parallel Plates. Ind. Eng. Chem. Res. 1995, 34, 1332. (79) Churchill, S. W. New Simplified Models and Formulations for Turbulent Flow and Convection. AIChE J. 1997, 142, 11125. (80) Churchill, S. W. The Prediction of Turbulent Flow and Convection in a Round Tube. Adv. Heat Transfer 2001, 134, 255. (81) Hinze, J. C. Turbulence; McGraw-Hill Book Co.: New York, 1959. (82) MacLeod, A. L. Liquid Turbulence in a Gas-Liquid Absorption System. Ph.D. Thesis, Carnegie Institute of Technology, Pittsburgh, PA, 1951. (83) Rothfus, R. R.; Monrad, C. C.; Senecal, V. E. Velocity Distribution and Friction Factor in Smooth Concentric Annuli. Ind. Eng, Chem. 1950, 42, 2511. (84) Langmuir, I. Convection and Conduction of Heat in Gases. Phys. Rev. 1912, 13, 471. (85) Newton, L. Scala graduum caloris. Philos. Trans. R. Soc. London 1701, 22, 824. (86) Boussinesq, J. The´ orie analytique de la chaleur, mise en harmonie avec la theorie me´ chanique de la lumiere; GauthierVillars: Paris, 1903; Vol. I. Boussinesq, J. Calcul du poivoir refroidissant des courants fluids. J. Math. Pures Appl. 1905, 69, 285. (87) Galante, S. R.; Churchill, S. W. Applicability of Solutions for Convection in Potential Flow. Adv. Heat Transfer 1990, 20, 353. (88) Hermann, R. Wa¨rmeu¨bergang bei freier Stro¨mung am waagerechten Zylinder in zweiatomigen Gasen. Ver. Dtsch. Ing. Forschungsh. 1936, No. 379. (89) Schmidt, E.; Beckmann, W. Das Temperatur- and Geschwindigkeitsfeld von einer wa¨rmeabgebenden senkrechten Platte bei naturlicher Konvektion. Technol. Mech. Thermodyn. 1930, 1, 391. (90) Nussett, W. Das Grundgesetz des Wa¨rmeu¨berganges. Gesundh. Ing. 1925, 38, 477, 490.

(91) Frank-Kamenetskii, D. A. On the Limiting Form of the Free Convection Law for High Values of the Grashof Criterion (Translated Title). C. R. Acad. Sci. USSR 1937, 17, 9. (92) Saville, D. A.; Churchill, S. W. Laminar Free Convection in Boundary Layers near Horizontal Cylinders and Vertical Axisymmetic Bodies. J. Fluid Mech. 1967, 29, Part 2, 391 (93) Churchill, S. W.; Churchill, R. U. A Comprehensive Correlating Equation for Heat and Component Transfer by Free Convection. AIChE J. 1975, 21, 604. (94) Martini, W. R.; Churchill, S. W. Natural Convection Inside a Horizontal Cylinder. AIChE J. 1960, 6, 254. (95) Aziz, K. Hellums, J. D. Numerical Solution of the ThreeDimensional Equation of Motion for Laminar Convection. Phys. Fluids 1967, 10, 314. (96) Ozoe, H.; Yamamoto, K.; Churchill, S. W.; Sayama, H. Three-Dimensional, Numerical Analysis of Laminar Natural Convection in a Confined Fluid Heated from Below. J. Heat Transfer, Trans. ASME 1976, 98C, 202. (97) Samuels, M. R.; Churchill, S. W. Stability of a Fluid in a Region Heated from Below. AIChE J. 1967, 13, 77. (98) Ozoe, H.; Ukeba, H.; Churchill, S. W. Numerical Analysis of Natural Convection of Low Prandtl Number Fluids Heated From Below. Num. Heat Transfer, Part A 1994, 26, 363. (99) Booker, V. B. An Experimental and Three-Dimensional Numerical Investigation of the Czochralski Solidification of Single Silicon Crystals. Ph.D. Thesis, University of Pennsylvania, Philadelphia, PA, 1997. (100) Yi, K.-W.; Booker, V. B.; Eguchi, M.; Shyo, T.; Kakimoto, K. Structure of Temperature and Velocity Fields in the Si Melt of a Czochralski Crystal Growth System. J. Cryst. Growth 1995, 156, 383. (101) Nakano, A.; Ozoe, H.; Churchill, S. W. Numerical Computation of Natural Convection for a Low Prandtl Number Fluid in a Shallow Rectangular Region Heated from Below. Chem. Eng. J. 1998, 71, 175. (102) Nusselt, W. Die Oberfla¨chenkondensation des Wasserdampfes. Z. Ver. Dtsch. Ing. 1916, 60, 541, 569. (103) Churchill, S. W. Laminar Film Condensation. Int. J. Heat Transfer 1986, 29, 1219. (104) Emmons, H. W. Natural Convection Heat Transfer Correlation. Studies in Mathematics and Mechanics Presented to Richard on Mises by Friends, Colleagues, and Pupils; Academic Press: New York, 1954. (105) Churchill, S. W.; Usagi, R. U. General Expression for the Correlation of Rates of Transfer and Other Phenomena. AIChE J. 1972, 18, 1121. (106) McAdams, W. H., Ed. Heat Transmission, 2nd ed.; McGraw-Hill Book Co.: New York, 1942; p 217. (107) Lemlich, R.; Hoke, R. A Common Basis for the Correlation of Forced and Natural Convection to Horizontal Cylinders. AIChE J. 1956, 2, 249. (108) Kirscher, O.; Loos, G. Wa¨rme- und Stoffaustausch by erzewungener Stro¨mung an Ko¨rpern verscheidener Form, Teil 1. Chem. Ing. Tech. 1958, 2, 249. (109) Bo¨rner, H. U ¨ ber den Wa¨rme- and Stoffu¨bertragung um umspu¨lten Einzelkorpern bei u¨berlagerung von freierund erzewungener Stro¨mung. Ver. Dtsch. Ing. Forschungsh. 1965, No. 512. (110) Jackson, T. W.; Yen, H. W. Combining Forced and Free Convection Equations to Represent Combined Heat Transfer Coefficients for a Horizontal Cylinder. J. Heat Transfer, Trans. ASME 1971, 93C, 242. (111) Sparrow, E. M.; Gregg, J. L. Buoyancy Effects in Forced Convection Flow and Heat Transfer. J. Appl. Mech., Trans. ASME 1959, 81E, 133. (112) Eschgy, S. Forced Flow Effects on Free Convection Flow and Heat Transfer. J. Heat Transfer, Trans. ASME 1964, 86C, 290. (113) Acrivos, A. On the Combined Effedct of Forced and Free Convection Heat Transfer in Laminar Boundary Layer Flows. Chem. Eng. Sci. 1966, 21, 343. (114) Churchill, S. W. A Comprehensive Correlating Equation in Laminar, Assisting, Forced and Free Convection. AIChE J. 1977, 23, 10. (115) Ruckenstein, E. Equation Between Two Limiting Cases for the Heat Transfer Coefficient. AIChE J. 1978, 24, 940. (116) Martinelli, R. C.; Boelter, L. M. K. An Analytical Prediction of Superimposed Free and Forced Viscous Convection in a Vertical Pipe. Univ. California Publ. Eng. 1942, 5 (2), 23.

Ind. Eng. Chem. Res., Vol. 41, No. 16, 2002 3929 (117) Kitaura, Y.; Tanaka, H.; Ueda G.; Kojima, N. Mass Transfer from a Vibrating Single Sphere and Packed Beds (in Japanese). Kagaku Kogaku 1969, 33, 175. (118) Graetz, L. U ¨ ber die Wa¨rmeleitfa¨higkeiten von Flu¨ssigkeiten. Ann. Phys. 1883, 18, 79; 1885, 25, 337. (119) Le´veˆque, M. A. Les lois de la transmission de la chaleur par convection. Ann. Mines Ser. 12 1928, 13, 201. (120) Reichardt, H. Die Grundlagen des Turbulenten Wa¨rmeu¨berganges. Arch. Wa¨ rmetechnik. 1951, 617, 129. English translation: The Principles of Turbulent Heat Transfer; Report TM 1408; National Advisory Committee for Aeronautics: Washington, DC, Sept 1957. (121) Churchill, S. W.; Balzhiser, R. E. The Radial Heat Flux. Chem. Eng. Prog. Symp. Ser. 1959, 55, No. 29, 127. (122) Nusselt, W. Der Wa¨rmeu¨bergang in Rohrleitung. Mitt. Forsch. Arb. Ing.-wes. 1910, 89, 1750. (123) Abbrecht, P. H.; Churchill, S. W. The Thermal Entrance Region in Fully Developed Turbulent Flow. AIChE J. 1960, 6, 286. (124) Yahkot, V.; Orszag, S. A.; Yahkot, A. Heat Transfer in Turbulent Fluids I. Pipe Flow. Int. J. Heat Mass Transfer 1987, 30, 15. (125) Reynolds, O. On the Extent and Action of the Heating Surface of Stream Boilers. Proc. Lit. Soc. Manchester 1874, 14, 7. (126) Prandtl, L. Eine Beziehung zwischen Wa¨rmeaustausch und Stro¨mungswiderstand der Flu¨ssigkeiten. Phys. Z. 1920, 11, 1072. (127) Jakob, M. Heat Transfer; John Wiley & Sons: New York, 1949; Vol. I, p 503. (128) Churchill, S. W. Critique of the Classical Algebraic Analogies between Heat, Mass and Momentum Transfer. Ind. Eng. Chem. Res. 1997, 36, 3566.

(129) Colburn, A. P. A Method for Correlating Forced Convection Heat Transfer Data and a Comparison with Fluid Friction. Trans. AIChE 1933, 29, 174. (130) Churchill, S. W.; Zajic, S. C. The Prediction of Turbulent Convection with Minimal Empiricism. AIChE J. 2002, 48, 927. (131) Heng, L.; Chan, C.; Churchill, S. W. Essentially Exact Characteristics of Turbulent Convection in a Round Tube. Chem. Eng. J. 1998, 71, 163. (132) Churchill, S. W.; Shinoda, M.; Arai, N. A New Concept of Correlation for Turbulent Convection. Therm. Sci. Eng. 2000, 8 (No. 4), 49. (133) Churchill, S. W.; Usagi, R. A Standardized Procedure for the Production of Correlations in the Form of a Common Empirical Equation. Ind. Eng. Chem. Fundam. 1974, 13, 39. (134) Churchill, S. W. New Wine in New Bottles; Unexpected Findings in Heat Transfer. Part III. The Prediction of Turbulent Convection with Minimal Explicit Empiricism. Therm. Sci. Eng. 1997, 3, 13-30. (135) Strenger, M. R.; Churchill, S. W.; Retallick, W. B. Operational Characteristics of a Double-Spiral Heat Exchanger for the Catalytic Incineration of Air. Ind. Eng. Chem. Res. 1990, 29, 1977.

Received for review September 17, 2001 Revised manuscript received April 15, 2002 Accepted April 16, 2002 IE010776O