.J
..
:FLUID M~CHANICSIN CHEWICA~. ~NGIN~ER;IRIO~
I
STUART W. CHURCHILL, PETER
H.
ABBRECHTl, and CHIAO-MIN CHU
The University of Michigan, Ann Arbor, Mich.
Regenerative Heat Transfer in Two- and Three-Dimensional Flow through Porous Media A solution of striking simplicity, obtained by assuming local thermal equilibrium between solid and fluid, may be applied for calculating the regenerative extraction of energy from the earth with an air well
A investigation was made of transient heat transfer from porTHEORETICAL
ous media to a fluid in a general twoor three-dimensional flow field. The fundamental problem has received considerable attention in terms of both heat and mass transfer, but only for onedimensional flow. I n 1946 Theile (76) reviewed work in both heat and mass transfer. Recent contributions have dealt primarily with techniques of numerical evaluation, and with complexities characteristic of ion exchange, chromatography, and chemical reactions but not of heat transfer. Accordingly, only developments directly pertinent to this investigation are noted here. Two- and three-dimensional flow fields in porous media have received attention in conjunction with oil well operation. Although only highly symmetrical geometries have been considered, the general techniques of potential theory are applicable (73). The general heat transfer problem is first formulated mathematically. Methods of determining the flow field are then examined. Solutions for one-dimensional flow are next considered, then their extension to general flow fields. Finally the results are illustrated for the regenerative extraction of energy from the earth by means of an air well. Mathematical Formulation
The rate of flow of a fluid through porous media at low velocities can be 1 Present address, Minnesota Mining and Manufacturing Co., Detroit, Mich.
correlated by the cxpression (2) :
and the material balance for the fluid phase can be written
Two important assumptions are made in formulating the equations for heat transfer : Temperature gradients within the individual particles composing the solid phase are neglected. The bulk and surface temperatures of the solid phase at any point are thus taken to be the same. This assumption is most valid for small particles of high conductivity and for low rates offlow. Thermal conduction from particle to particle or through the fluid phase is neglected. The validity of this assumption can be checked by direct computation with measured effective condiictivities for porous media through which fluids are flowing ( 7 7 ) . The other extreme of infinite conductivity is subsequently considered. With these assumptions the energy balances can be written in terms of a coefficient for convection. For the solid phase
+ - v ( V p c T ) f ha(t
- T)= ' 0 (4) be
Equations 1, 2, 3, and 4 define the general problem and can be solved by numerical methods without further simplification, if the properties and coefficients are known as functions of pressure, temperature, and velocity. However, useful results can be obtained with a great deal less effort by assuming mean values for all properties and coefficients. Equation 2 can then be written as
and Equations 3 and 4 reduce to
and
where M = ha/p,c,(l
- X)
N = ha/pcX +
u
=
+
v/x
Equations 6 and 7 define the local temperature histories of the solid and fluid in terms of the effective velocity, + V, and the parameters, M and N , made up of mean properties and coefficients. Equations 1 and 5 define the flow and VOL. 49, NO. 6
JUNE 1957
1007
Because the Laplace equation is linear, a solution for a point source and an equal point sink a t a distance, 221, as illustrated in Figure 1, can be found by adding algebraically the potentials for the source and sink alone, giving I I I
PLANE OF SYMMETRY
i
i
x
-21
Figure 1. Representation of point source and point sink in an infinite region
pressure fields. Equation 5 can be solved and the velocity evaluated from Equation 1 independently of Equations G and 7 . Methods of solving Equation 5 are accordingly examined before the solution of Equations 6 and 7 is considered.
Determination of Flow Field Equation 5 , the Laplace equation, probably has been the subject of more extensive analysis than any other equation. The applicability of analytical, analog, and numerical methods is examined. Analytical Methods. Analytical solutions have been published for only a few boundary conditions pertinent to twoand three-dimensional flow through porous media. One technique for obtaining and extending such solutions is illustrated. The solution for radial flow a t a volumetric rate, q, from a point source in a three-dimensional infinite medium is
3
0
2
2'
3
The equation for the streamlines for this case is
This solution can also be interpreted as flow from the infinite plane equidistant from the source and sink to the sink. The flow pattern-Le., the configuration of the equipotential lines and streamlines-is independent of q, K, and p, but the total loss in potential is proportional to qp/K. Although the above solutions are for point sinks and sources, it can be inferred from Equations 8 and 9 that the total potential drop would be inversely proportional to the diameter of small, finite sinks and sources. The method of images next can be used to construct a solution for flow from a point source to a symmetrically located point sink in a slab finite in one dimension, as illustrated in Figure 2. To counterbalance the flow over the boundary of the slab a t B an image source is located at point 2 and an image sink at point 3. T o restore symmetry over plane A an image sink and source are then located a t 2' and 3', respectively. This process can be continued indefinitely. The solution is composed of the algebraic sum of the solutions for the individual sinks and sources. Only a few terms of the infinite series contribute appreciably in most cases. Similar constructions are possible in other highly symmetrical systems, but become unwieldy if not impossible for
IMAGE SINK
trical circuitry, ideal fluid flow, etc. Numerous techniques have been developed for this purpose. Most attention has been given to two dimensional fields, but a fluid mapper (72) or an electrolytic model ( 9 ) can be constructed for almost any boundary conditions yielding three-dimensional flow with axial symmetry. Fluid mapping yields a photograph of the streamlines. Determination of the potential and the velocity from the photograph involves considerable work and some practice. An electrical analog offers the possibility of a direct indication of the velocity on a recorder. Numerical Integration, I n principle the Laplace equation can be solved for any boundary conditions to any desired degree of accuracy by numerical integration. Equation 5 is approximated by a finite difference equation at an arbitrary number of points in the region. This set of equations is then solved, usually by trial and error. The required number of points for a given accuracy depends upon the complexity of the boundary conditions and the conformity of the grid to the actual strearnlines and equipotential lines. The process is too tedious to carry out by hand, but automatic digital computers are applicable. BLOWER
r 3TO HEAT EXCHANGER
I M A G E SOURCE
0
IMAGE SINK
IMAGE SOJRCE
Figure 2. Method of images applied to point source and point sink in a bounded region
1008
complex boundaries. Attention is therefore turned to more general methods. Analog Devices. As the Laplace equation applies to any potential function, flow through porous media can be represented by its counterpart in elec-
INDUSTRIAL AND ENGINEERING CHEMISTRY
Figure 3.
Air well
FLUID M E C H A N I C S Heat Transfer Solutions
Ns
The temperature history of the fluid and the porous media are obtained by solution of Equations 6 and 7, by using a previously determined velocity. field, the inlet temperature history of the fluid, and the initial temperature distribution in the porous media. The problem could be solved by numerical integration, but such transient problems are much more difficult to solve than the Laplace equation. So much success is achieved by analytical methods that these involved numerical techniques are not discussed. One-Dimensional Solutions. For oneflow at 'Onstant Equation 7 reduces to bT -u -F N ( t - T ) = (11) b8 *
where s is the distance from the inlet. Nusselt ( 7 4 solved Equations 6 and 11 for an arbitrary initial temperature distribution and a constant inlet temperature. Theile (76) adapted Nusselt's results for mass transfer with a variable inlet concentration and a uniform initial concentration and also corrected for accumulation by the void space, which Nusselt neglected. However, Equations 6 and 11 can be
T ( s ,e) = t(s, e ) = m(s) for
6
e
(12)
and an arbitrary inlet temperature history T ( s , e) = n(0) for
s =
0
E ( ~e,) = M
x
e-
x -M(B
6 - .)
-
e
and F(s, e) = Ne-M(O
- 3)
x
[m(r)
J e - ~ r -
7)
x
{
- m(o>~
1m(s)
- m ( o ) ~x
11
[2
T ( s ,e) = A(s, e) + B(s, e ) +
t(s,
e)
= D(S,e)
- _N s
(14)
+ E(s, 0) + F($, e)
where
+ m(0)
+ m(o)
(15)
NMs
(8
(e -
T)]
X d r (20)
x
dy (e - $) - ?-)Ix -
As implied by the initial condition (721, the heat capacity of and hence heat transfer to the fluid initially in the pore space are neglected. The slight error arising from this simplification will decrease rapidly with time. The process of solution is outlined at the end of this article. The results are as f o h w s for s < UO:
C(S, e)
47 - 6 + 4; 6)
x lo [2
(s
dr
Nusselt's solution is actually equivalent to the special case n(0) = n(0) = m(0); hence A(s, e), B(s, e), D(s, e), and E(s, 0) are all zero. The thermal counterpart of Theile's results corresponds to the special case m(s) = m(0) = n(0); hence A($, e), C(s, e), D(s, e), and F(s, e) are all zero. For the special case of uniform initial temperature to = m(s) = m(O), and constant inlet temperature To = n(0) = n(0); B(s, e), C(s, e), E(s, e), and F(s, e) are zero; and as noted by Furnas (5)
or functions from which the integrals can be obtained have been evaluated and graphed by many investigators, including Schumann (75), Furnas ( 5 ) ,Brinkley and Brinkley (I), Hougen and Marshall (8), and Hiester and Vermeulen (7). Klinkenberg (70) developed empirical equations
VOL. 49,
NO. 6
JUNE 1957
1009
I
TD.
hr
The effect of thermal conduction thmugh the bed, which bas bem neglected throughout, is also to round OE the tempUatirre profile. For the exheme case of infinite conductivity or perfect mixing in the bcd, and thermal quilibrium between the solid and thc fluid, the tempvature of the solid and fld,is d o r m and qual to the exit fluid temperature.. This tempaacurqis defined by the energy balance
ze YI 32 34 36 3 DIST4NCE FROM 4XlS OF SVMYETRI-K
&,-tion 23 5 mdjy .so] ad and boundary conditionsapreared by Equations 12and 13 tolield T(0, I)
and ~(e,
-
- U d )for
I(*
n
.
- &]-fo;
w
.
.i-arm~-
Ondbll&Onal flow be adaPtcd This Very Sbllph sohItion lhould be a cHro-.or. three-dimflow as reasonable approximation for the con. . ditiona bounded by Van Dermter (77). i ' the WO*te sysmn is dlThe effect of a finite rate ofheat transfer pect to the s e a , between the solid and the fluid is to,' m n d ORthe stepfunction and to cause : . , * bT (a the temperature wave in the solid to V V T - IVI lag behind the trmperature wave in the, f70oI, no& that thisx . whclr I is the distance m & d a l o q a fluid. &km time lag is approximately e q d to s e e . S a d , let '
I
.
> U&
(24)
. 3
< U.O
. (25)
W6ae
'
''
_.
dr
whae
~
,
0
6
0:
8000
.4,
a;
4000.
z
0.
0.
a
I
,.
Figure 5. Flow funaims for air well 15 feet below surface
.
'
,
1010
w-1
(
bT
bT
sional solutions for zao cnndwtivity are
o>
lM00
-
-
(31 Irr + N(L - T ) which has the Same form as Equation 11. Thus it is apparent that the onedimen-
-Uo-
DO 16000
Ve is the absolute velocity at J
Then
I I
(29)
AND
B#Mpwo
&Y
applicable to two- or three-dimensional flow in t a m s of variablca S and UOinstcad of s and U. Equation 18 is similarly applicable for two- or threedimensional flow only if conauction or mixing is limited to the dimtion of flow. The fluid teaching a sink along different s t r e a d ~ ~ ewill a not all have the same tern-. Thc mixed temperature en'tuing a sink must therefore be found by integration over the sink. The volumetric rate of flow between any
two streamlined..
whae xp is the distandc from the axis of
ptrahue at the sink along the srreamline *gh XU. : . ,
- For
stepfuncti& . solution, the
'
temperature of the individual .-n lines at .$&will either be toor To and
T m ~Lo=+f(d (Te- ~ f o )
"(33)
.
'
yhem f(xo) is the fraction of the total 9ow entering the sink inside the stream. lioethroughxo,andx~ismtheU'h%d , 'atreamhe whoae total length, $9, is
.
,
E
'
.
.
3 6 b
'
exactly equal t? Ud. working ploe off, q/Va and S(9 os. xa can readily he pnpar~dfor any given flow pattern. In turn a general plot of f(xcJ DS. ~
of T, w. B for a gim flow rate, Pair of amperaturn, and set of physical p m p ' t k can he Prepared. . r n O 9 expaegions have been pmrnted and diuxlascd in turns of q and 0, but it is a p p p t t h a t . & t e n p a h u e at the sink in a funetioir only of the total vol-,
.
MlLLW CUBIC FEET C.
'
.
. Ftgure.6. Mixed.tempr,atureatcairwell 15 feet below surface, ,
. ,
,
tae
some.di.mce MOW ground. he warin air ontering the pipe is thv.used an an energy source i n , a heat.pump
S d ' imensionii solution
problem defined by EguatiOna 6,.H,12,and13canhesimplifiedbythe'
point sink, the flow field given by Q u a tion8 9 and 10 is applicable: The following pmpcrtiw and condition8 are
asurndl
=,
=,is&
to
= 520 .
,
~.
I
9-.9-.
"
mbsti?tion
'
'
J = ,
. \
yi
137)
-
, . . . . . , aT. U x+ N ( f - T) "
.. . ' r a t h a than an ind of both q and 83 in sofaran the wliea and tned"p"y of 9. The mixed temperature at.& *nk for t h e m of infinite conductivitfor mixing & & dirstian of flow only in given hy the integral
Repramtstve meambn ' a and qui: potentiaUne8 ire &etched in Figure 4. . velocit*s for ,a flow of 36,000 eubic fcet per hour are shown at &e .intemections. A working plot of f 'V& and S(0 os. xb is shown in' Figure temperature
i.n;,
,
(36)
U
--M(f.~
,
'
NO,
T)=
-
0
at
-
(38) (39)
ar
(40)
m(rY
and
Tb, 0 ) = n(r)
'
(41)
Taking the L a p h transform with respect t o y then yields
-
& d is plomd total air & in Figun 6.' The tanperaturr that would u ar x N(Z .- TI 0 '* (42) he obtained for iniinite conductivity or TA = TI (Te to) mixing in the direction of flowin .included p 2 -k M(2 T)-=m b ) (43) for Cornpariaon. ' QWd For 'no conduction, the air entering and the sink r e m a i n s MOBM at the initial s ["' X I , w o ) ii@) . (4 t e m F a t u r e of the &rth until the temp a t u r e wave first tache the SidL along The solution of the &ormed problem (35) ' the axis of sybautiy. The mjxed tem- ,is " . . , M dmmce asymptoticauy HM)i Extraction of Enby an Air Well . to.&a-en- then air temperatun. Far in-' P+M U+ TW J)'= e, finite conductivity or mixing along the AS an exam+, 'the extiacted areamhrd , p l y , the temperature c a t e r a from the ground by an air Well sug(N - ")( gratid by D q (4)is examined. Such a . hg the sink is always 16s than that for P+M B',-,) in*rltetehed in P i 3. 'Cold no conducyion, buf tae W q c e heNJ m ( r ) d r -air'entcmthe's;rface o f ~ + e twear.the,& decwith the total P+M (45) ground a+ fl&;th-ph the,1yarmcr flow. .The difference in even l e s s ' f a pornus ea& to the open end of:a pipe .finite conductidties. and
-
-
-
+
...
-
>-
-
+
pa
Arg
'
-
I
.VOL
'\
4, NO. 6 .*
JUW 1987
1011
equivalent procedure, but the results are considerably more detailed.
The inverse transform of Equations 44 and 45 can be expressed in terms of the following functions, defined for s 2 0. F d y , s) = 1
+
21 (2
e)
and
F3Cy, s ) = e
- MY ZO(2
47)
For s < 0, all these functions are zero. Using formulas for the inverse Laplace transform given by Carslaw and Jaeger ( 3 ) ,it can be shown that
From these relationships, the inverse transforms of 45 and 46 are, respectively,
and
tb, s)
= m ( s ) e - MY
+M e
- _N S
Conclusions Transient heat transfer in one-dimensional flow through porous media can be computed from the analytical solution developed herein for any initial temperature distribution and any inlet temperature history. The much simpler solution obtained by assuming local thermal equilibrium is a reasonable approximation under many circumstances. These solutions can be utilized for any two- or three-dimensional flow field by a simple change of variables. An illustrative calculation suggests that thermal conduction through the porous media may often be negligible. Nomenclature = specific surface of solid, sq. ft./cu. ft. = heat capacity of fluid, B.t.u./ (lb.) (” F.) = heat capacity of solid, B.t.u./ (lb.) (” F.) = particle diameter, feet = fraction of flow inside streamline through x o = acceleration due to gravity, ft.,/ (hr.)2 \ , = conversion factor, (lb.) (ft.)/ (lb.F) (hr.)2 = convection coefficient, B.t.u./ (hr.) (sq. ft.) (” F.) = modified Bessel function of first kind and zero order = modified Bessel function of first kind and first order = permeability of solid, sq. feet = total distance through bed, feet = inverse Laplace transform of = ha/p,c,(l - X), (hr:) = initial temperature distribution in solid, F. = ha/pcX, (hr.)-I = inlet temperature history of fluid, O F . = Laplace transform of n ( y ) = pressure, lb./(sq. ft.) = variable of transformation = total flow rate, CLI. feet per hour = total flow, cu. feet = radial distance from point source, feet = radial distance from point sink, feet = distance along streamline, feet I
,
I
Equations 53 and 54 can readily be rearranged in the form of Equations 14 through 21. It is possible to obtain a solution, taking into account the heat capacity of the fluid initially in the void space by an
1 0 12
.
I
= temperature of solid, O F. = constant initial temperature of solid, ” F. j @ , s) = Laplace transform of t(y, s) T = temperature of fluid, ” F. ~ ( ps), = Laplace transform of T(y, s) TO = constant inlet temperature of fluid, ” F. TI = temperature at sink along any streamline, F. T,,, = mixed temperature at sink,
INDUSTRIAL AND ENGINEERING CHEMISTRY
O F .
-+
U
UO
uw
--c
= = =
-.,
V
=
vo
=
X
=
XO
=
XOO
=
X
=
Y
=
2
=
21
=
Z
=
P
=
P S
=
e 4
= =
w
=
li/
=
VIX, feet per hour Vo/X,feet per hour U/(1 [ P d l - X)/CPXI) = velocity of temperature wave, feet per hour
+
superficial velocity, feet per hour absolute velocity at s = 0, feet per hour distance from axis of symmetry, feet distance from axis of s p metry in source plane, feet distance from axis of symmetry to critical streamline in source plane, feet porosi ofsolid e distance from plane of symmetry, feet half distance between point sources, feet vertical distance, feet density of fluid, pounds per cu. foot density of solid, pounds per cu. foot time, hours Pg. g p Z , lb./(ft.) (hr.)2 viscosity, lb.,’(hr.) (ft.) stream function, \lb./(ft. (hr.)2
+
Literature Cited (1 ) Brinkley, S. R., Jr., Brinkley, R. F., J . Appl. Phys. 18, 582 (1948); Math. Tables, Aids Comp. 2, 221 (1947). (~. 2 ) Carman. P. C., “Flow of Gases through Porous Media,” Academic Press, New York, 1956. ( 3 ) Carslaw, H. S., Jaeger, J. C., “Conduction of Heat in Solids,” Oxford University Press, London, 1947. (4) Dana, H. J., “Typical Design Problems an the Residential Heat Pump,” State College of Washington, Eng. Expt. Sta., Pullman, 73 139491. ( 5 ) Fumai, C. C., Trans. A m . Inst. Chem. Engrs. 24, 142 (1930). ( 6 ) Gamson, B. W., Thodos, G., Hougen, 0.A., Ibid., 39, 1 (1943). (.7 .) Hiester, N. K., Vermeulen, T., Chem: Eng. Progr. 48, 505 (1952). (8) Hougen, 0. A., Marshall, W. R., Jr., Ibid., 43, 197 (1947). (9) Jakob, M., “Heat Transfer,” vol. I, p p . 401ff., Wiley, New York, 1949; Klinkenberp, A., IND. ENG. CHEM. 40, 1992 (1948). Mo_lino,_D. F., Hougen, J. O., Chem. mzg. rrogr. 48, 147 (1952). Moore, A. D., J . Appl. Phys. 20, 790 (1949). Muskat, Morris, “Physical Principles of Oil Production,” McGrawHill, New York, 1959. Nusselt, W., Z. Ver. deut. Ing. 71, 85 (1927). Schumann, T. E. W., J . Franklin Inst. 208, 405 (1929). Theile. E. W.. IND. ENG. CIIEM. ’ 38. G46 (194G). (17) Van’Deemt‘er, J: J., Zbid., 45, 1227 (1953); 46, 2300 (1954). RECEIVEDfor review January 2, 1957 ACCEPTEDApril 24, 1957 Division of Industrial and Engineering Chemistry, ACS, Symposium on Fluid Mechanics in Chemical Engineering, Lafayette, Ind., December 1956.
