- Email: [email protected]
Dynamic response of a packed bed thermal storage system—a model for solar air heating
,Solar Energy, Vol. 29, No. 3. pp. 201-206, 1982 Printed in Great Britain.
0038~92X]82/090201-06503.0010 Pergamon Press Ltd.
DYNAMIC RESPONSE OF A PACKED BED THERMAL STORAGE SYSTEM--A MODEL FOR SOLAR AIR HEATING A. E. SAEZ and B. J. McCoY Department of Chemical Engineering,University of California, Davis, CA 95616, U.S.A.
(Received 18 September 1981; revision accepted 4 January 1982) Abstraet--A mathematical model for simulating the dynamic temperature response of a packed column to an arbitrary time-dependent inlet air temperature is developed. The model includes axial thermal dispersion as well as intra-particle conduction, features that have usually been neglected but can be important in solar energy applications. Solutions, presented in terms of moments of the temperature response to an impulse of heat at the inlet, can be evaluated by simple numerical quadrature. Results of the model compare favorably with experimental data found in the literature. The model is used to optimize heat storage in a rock bin system subject to a realistic transient inlet temperature.
INTRODUCTION
THEORY
A packed-bed thermal storage system, e.g. a rock bin, consists of a column packed with pebbles or some other heat-absorbing solid material through which air flows when energy is transferred. Such systems operate in one of three modes:
Our model is identical to that used by Sagara et al. [I] to analyze pulse response data for heat transfer in a column packed with porous spheres. Wakao et a/.[2] have recommended the same model for analyzing data to determine heat transfer coefficients in packed beds. Consider a column of length L packed with solid spheres of diameter 2R. We make the following assumptions: • the flow in the column is axial and incompressible • viscous heating is negligible • radiant heat transfer between particles is negligible • physical properties of the fluid are constant • surface temperature of a particle is uniform • heat losses from the column are negligible • temperature and velocity profiles in a column crosssection are uniform. Under these conditions, the equation for the interparticle fluid T(t, z) is
(i) Charge mode Air is heated in the collection system and enters the bed, increasing the bed temperature. The charge mode for solar heating occurs usually in the morning and early afternoon hours, when the energy collection is large. (ii) Discharge mode Cool air is heated by passing through the bed, usually at night. (iii) Stand-by mode No air flows through the bed. This mode occurs when further storage of heat cannot occur because of the decreasing inlet air temperature, and no air heating from the bed is required. This is the transition period between the charge and discharge modes. The accurate mathematical description of this timedependent system is the objective of the present work. While accuracy is required, we also desire simplicity in the use of final equations. In the following sections we present a model for an adiabatic (well-insulated) column operating under transient conditions. The model, which includes axial thermal dispersion and intra-particle heat conduction, is compared to more limited models in the literature. We show that the temporal moments of the response to an impulse of heat at the inlet can be used to construct the response to any time-varying inlet condition. An expression for the total energy stored in the bed as a function of time is developed, and used to calculate when the charge mode should end for a realistic solar input. This optimization provides a basis for improved system control.
02T
(1)
where e is the interparticle void fraction, p and cp are the density and heat capacity of the fluid, Vo= ev is the fluid volumetric flow rate per unit column cross-sectional area, as = 3(1- e)/R is the surface area of spheres per unit column volume, and hs is the fluid-to-particle heat transfer coefficient. The effective axial thermal dispersion coefficient, kz may account for longitudinal particle-to-particle conduction, as well as fluid conduction and hydrodynamic dispersion of heat axially in the column. The equation for the temperature inside a sphere of radius R as a function of time and radial position, Ts(t, r), is
OT~ 1 ~ ~r2 8T~ p.cps - ~ = k. ~ Tr \ ~ - / 201
(2)
A. E. SAEZand B. J. McCoy
202
where p, and cps are the mass density and specific heat capacity and ks is the thermal conductivity coefficient of the solid spheres. At the column entrance we have the arbitrary input temperature,
T(t, O) = To(t).
(3)
The moments may be derived from the Laplace transform of ~(t), m. = ( - I)"(d"~/ds")[~=o which follows directly from the definition of O(s). Letting
When the column is very long compared to the diameter of the spheres, we may apply the usual boundary condition [ 1],
T(t, oo) = finite.
(4)
t~= T - T, and
~,s: L - T,
At the surface of spherical particles we require the heat flux to be continuous,
h s ( T - T~(r : R)) = k~(dT, lOr)[r=R,
(5)
(12)
we may apply Laplace transformation to eqns (1)-(6), solve for 0(s, z), and then use eqn (12) to calculate the moment expressions [1],
and at the center of a sphere, ~rr [T~(t, r)],=o : 0.
(6)
We assume that the bed is initially at some reference temperature, T, T(0, z) = T,.
(7)
TA0,~ = ~.
(8)
too(z) = too(0)
(13)
#'t(z) = #I(0) + (1 + 8o)zlv
(14)
#2(Z): # 2 ( O ) + [ t $ 1 + ~ ( I L epcpv
MOMENTS
Sagara et a/.[1] have used moment methods to interpret heat transfer processes in packed columns. The present objective is to use moments to predict the dynamic response of a packed bed to an arbitrary temperature stimulus at the column inlet. The temporal moments are integrals over time of a distribution function, ¢J, f~
m, =Jo t"6(t)dt, n =0, 1 , 2 . . .
(9)
(15)
where 1 - ~p~cps 60- - - - -
e
From the coupling boundary condition (5), an order of magnitude analysis shows that when the Biot number, Bi = h~2Rlks, is not small, intraparticle effects cannot be ignored. Thus, for particles of small thermal conductivity or large diameter, or for high heat transfer between fluid and particle, which might occur for large fluid velocities, Equation (2) is necessary for the accurate description of packed bed heat storage. If B i ~ 1, intraparticle heat transfer is negligible, and at any crosssection of the column T = Ts = constant is a good approximation.
+t$o)212Z/V
(16)
pcp
8~ = 1 - E Ps c..._p~
pcp 15
1
5
+
"
(17)
The zeroth moment, the area under the curve of T - T, plotted versus time, can be interpreted as the amount of heat entering the column and moving conservatively along the column. The first moment is the average temporal position of the peak, and the second central moment indicates the width of the peak. For the present purposes we may consider that inlet fluid temperature pulse is a delta function,
T(O, t ) - T~ = 0o8(t)
(18)
so that too(Z)= 0o, the strength of the impulse, and #'~(0) = O, #2(0) = O. Because energy is stored in the solid spheres, we need the single-particle average temperature, defined as
3/?
(Ts) = ~--~
T~ P dr.
(19)
The normalized nth moment is defined by
#"
n = 0, 1 , 2 . . .
(10)
(t-#D"O(t)dt, n = 2 , 3 . . .
(11)
= m,]mo,
and the central moments are
#" =-moo
The moments of this temperature, relative to T, may be found by applying the single-particle averaging operation to eqn (2), finding the Laplace transform, and using eqn (12). We find the following expressions for the moments of the temperature pulse in the solid particles:
mos(Z) = too(Z)
(20)
Dynamic response of a packed bed thermal storage system /*[,(z) =/*[(z) + ~
(~~~ h - +- ~ )1l
2 4 /*=,(z)= /*=(z)+p,2 c;,R [ 13 2 ,17-Tk7+ sk-
(21)
I
+ hZ /• (22)
The first and second particle moments show delays relative to the fluid moments, and reduce to fluid moments for very large values of ks and hs, i.e. for large rates of intraparticle heat transfer. The dependence of the second moment on particle diameter, 2R, shows that broadening of a temperature profile within the column is decreased with smaller particles. Thus, since smaller particles are favorable for heat storage but have the disadvantage of increasing the pressure drop across the column, a trade-off of costs is required to determine the optimum particle size. The temperature response for the delta function input can be expressed as a series expansion in orthogonal polynomials[3]. Hermite polynomials, which are simple to use in this context, provide sufficient accuracy for most purposes. We write the series expansion as T - Ti = mo exp (- x 2) ~
a.(z)H.(x)
(23)
203
rates[3]. The temperature profile during the discharge mode may be expressed as the response to a step down in temperature if the bed temperature is uniform. When axial thermal dispersion and intra-particle effects are very important the time-response to a delta function input is not well represented by a Gaussian distribution. In that case, the use of Laguerre polynomial expansions is recommended, since the first term of the series will be a Poisson distribution function which will represent better the solution. In terms of Laguerre polymonials, L,, the expression is T-
T, -
mo [ M ] ~-'
(A--'l)! \/*~1
(-At)~__o
exp ~
(M)
k.L. klt,/ (28)
where t2
A ---.-/* J #2
(29)
The coefficients k, can be found by applying the orthogonality relationship: ko(z)= l, k~(z)=0, k2(z)= 0. Applying the convolution theorem gives the response for any kind of input
n=O
- t')T(t,z)- Ti= [,o' [To(t-~A'-'~.
where
T~] t~]) /At'\a-' exp ~--~) /At'\ dt'. (30)
t-/*~
x = V(2/.2)'
(24)
The first term of the expansion yields a Gaussian distribution in time centered at /*~ with variance /*2. The first three coefficients, a,(z), of the expansion can be found from the orthogonality relation for Hermite polynomials to be ao(z) = l/X/(21r/*2), al(z) = 0, and a2(z) = O. In Laplace transform space the solution to eqns (1)-(6) can be expressed as the product of the inlet temperature and the response to the delta function input; therefore, it follows from the convolution theorem that
For a step input, the solution is T - T~= (A - 1)! r A,
(31)
where F is the incomplete gamma function, defined by:
F( &-~2 ) = fo(~it/~2)x~-' e Xdx.
(32)
COMPARISON WITH EXPERIMENT 1
~* t
T(t, z)- T~- X/(2-1r/*2)Jo [To(t- t') - T~]exp [- (t' -/*~)2/2/z2] dt'.
(25)
For the special case of the step input, i.e. To(t)- T,=
O
, t< 0
(26)
-T,, t > 0 '
Equation (25) yields t ' t-/*l
/.1 (27)
This representation is known to be accurate when /*~X/(2/.2)~2.5, i.e. for long columns and low flow
We have compared the prediction of our model with the experimental data of Vanden Brock and Clark[4] for an insulated column 8.00 cm in dia. and 24.2 cm in length packed with steel spheres 0.556cm in dia. The void fraction of the bed was 0.386. Because the Biot number is small for the air/steel system of these experiments (Bi=0.01), intra-partiele conduction should not affect results substantially. Thus, the main advantage of the present model over that of Clark et al.[5] in describing these experiments is that our model includes axial thermal dispersion. Outlet temperature responses were measured for step inputs of the air temperature; experimental data for one run is shown in Fig. 1. For the present model, physical properties of air and steel are evaluated at the mean temperature (1/2) (To+ T~). To test this assumption we made calculations for properties evaluated at extreme temperatures in a rock bin, 26 and 65°C. The maximum difference in the
A. E. SAEZand B. J. McCoy
2O4 7O
v
i
,
l~O
2~I0
60
50
40
30
~;0 t,
4;0
~~0
8;0
7&O
800
sec
Fig. I. Response to a temperature step input: Solid line is experimental data, dashed line is the present theory, (7) Vanden Broek and Clark[4] theory, [] Riaz theory[6]. calculated temperature profiles was 2°C, even though the properties were evaluated at two temperatures almost 40°C apart. Heat transfer coefficients are estimated with the correlation[2], Nu = 2 + 1.1 Pr 1/3Re °'6.
(33)
Thermal axial dispersion coefficients are estimated from [2]
kz = 0.7k + epcpvR.
(34)
The results of the numerical calculations made by Vanden Broek and Clark[4], shown in Fig. 1, are based on heat transfer coefficients regarded as adjustable parameters to fit the data. Heat transfer coefficients for the Riaz model[6] were based on eqn (33). The results of the present model, with Hermite polynomial expansions, are plotted in Fig. 1. We also plot the predictions of Riaz[6], and of Vanden Broek and Clark[4]. The model of Riaz is a single-phase conductivity model. The model of Vanden Broek and Clark ignores axial dispersion and intra-particle conduction. The accuracy of the different theories with respect to experimental data for four runs was represented by an average absolute deviation, defined by
r ~' [Tcalc-Texpldt' 8=-~1 Jo
(35)
where the absolute difference between calculated and experimental temperatures is averaged over At, the duration of the experiment. The present model gave values of 6 between 2 and 4°C. The model proposed by Riaz yielded larger errors, sometimes twice as large, and tended to underestimate the experimental curve. The Vanden Broek and Clark model, while tending to overestimate experimental data, showed deviations similar in magnitude to those of the present model. Even though the mean absolute deviation 8 for the present model is 3°C, errors up to 8°C were observed. The cause of the errors is likely to heat losses from the column. Vanden Broek and Clark[4] reported radial
temperature variations across the column cross-section, which can be attributed to heat losses. The temperature variations are not due to velocity non-uniformities, because velocity was measured experimentally[4] and found to be uniform over most of the column length. The criterion for uniform radial temperature when heat losses occur is that the column diameter or overall column heat transfer coefficient be small, or that the radial thermal dispersion coefficient be large. An interesting extension of the present work would be to include the heat loss effect by adding a loss term to eqn (1). The overall heat transfer coefficient could be treated as an empirical parameter, whose value would be adjusted to fit experimental temperature profile data. Adjusting the value h~ in the present model did not improve the fit to the experimental data. Another possible cause for the deviation between the present model and experimental data is an overestimation of the axial thermal dispersion coefficient; the smaller the value of this coefficient, the sharper the temperature curve. The predicted curves are broader than those measured experimentally, indicating possibly that the correlation may provide a value of kz that is too large. Hughes et a/.[7], have modified the Schumann model with empirical corrections for intraparticle heat conduction and axial dispersion. When an energy loss term is added, a finite difference method is required to solve the system of partial differential equations. The present model may also be modified for heat losses, and once temporal moments have been determined, a single numerical integration provides the temperature response for an inlet function. Further work on the topic of energy losses will be reported at a future time. APPLICATION TO A SOLAR E N E R G Y S T O R A G E SYSTEM
As an illustration of the applications of the model to a solar energy storage system, a typical rock-bin unit was simulated. Consider a column of dia. 1.22 m and length 6.1 m packed with granite rocks 6 cm in dia. The objective is to find the hour of the day to discontinue the charge mode in order to store the maximum amount of energy. The hourly ambient temperatures and energy collection rates are taken from experiments[8] at the New Mexico State University Solar House on a sunny day, 7 February 1977. A typical air flow rate, 0.30kg/sm 2, is chosen. The physical properties of air are evaluated at 25°C, and the properties of the rock are found in Clark[9]. The basic data for the simulation are given in Ref. (7). The heat transfer coefficient and thermal axial dispersion coefficient are calculated, by using the correlations (33) and (34), to be 28.5 W/m2K and 9.64 W/mK respectively. We will assume that the bed is initially at a uniform temperature of T~ = 5.6°C, and the entering air is at ambient temperature. We consider that all the heat collected is transferred to the flowing air. With these assumptions the following heat balance applies to the collector where To is the outlet temperature of the collector and therefore the inlet to the rock-bin:
Dynamic response of a packed bed thermal storage system
q = Gocp(To- T~)dc21r[4
(36)
from which we obtain
To - T~ = 4qlcrGodc2cp + Te - T~.
(37)
Note that T, and q, and therefore To- Ti, vary with time. In Fig. 2, we have plotted To- T~ as a function of time. We approximated To- T~ by straight line segments over every hour interval to introduce To- T~ as data for the numerical integration of eqn (25). Since axial dispersion and intraparticle conduction are not large in this case, comparison with the Laguerre expansion showed the Hermite polynomial expansion to be the better approximation for the temperature response for air and particles. The total energy, Q, stored at a certain time can be expressed in terms of energy available relative to a reference level, which is taken to be the minimum night temperature, 5.6°C. We have
Q = ~R2~pcp f o L ( T - T,-)dz + IrR2(l - ~)p, cp~fo L (T~ - T~) dz,
(38)
where, since pscp~ >>pcp, the latter term dominates. The available energy Q can thus be calculated as a function
60
5C
4C tJ o 3C i 2C
5
tO
t,
,~
of time by numerically integrating results of eqns (37) and (25); Q is plotted versus time in Fig. 3. From Fig. 3, the time at which the thermal energy availability is maximum is t = 9.5 hr, i.e. 4:30 p.m. This is the time when the charge mode should end, and the stand-by mode should begin until energy for heating is required. Continuing the charge mode beyond this time will actually decrease the total stored energy. It is useful to note that the maximum available energy, 380M J, exceeds the space heating needs of an average house. However, since the process of withdrawing energy from the bed is not 100 per cent efficient, the total 380 MJ is not all usable. CONCLUSION
The model presented in this paper, advocated for applications of heat transfer in packed beds, has also been successful in its mass transfer form in describing chromatographic and breakthrough curve processes[3]. The model is based upon reasonable assumptions and contains transport parameters, principally the longitudinal thermal dispersion coefficient, kz, and fluid-to-particle heat transfer coefficient, hs, which have been well studied. Even though the model includes axial thermal dispersion and intraparticle effects, the solutions are rather simple and sufficiently accurate for practical application, so that their incorporation in designing procedures for solar energy is worthwhile. A computer program is available [10] that gives the air temperature and particle average temperature response for a packed bed subject to an arbitrary temperature input. The program calculates the temperatures by using Hermite and Laguerre expansions and includes up to the third term in the expansions. Although this program is in Fortran, the computations involve only numerical integration and parameter substitution, which are readily performed on a programmable calculator. Two kinds of polynomial expansions were used. Hermite expansions give better results for the case of low axial dispersion and low intra-particle effects. Laguerre expansions are better in the opposite case. The model was used to simulate a solar energy storage system undergoing real conditions to illustrate its manageability in the simulation of practical processes.
hr
Fig. 2. Rock-bin air temperature input. NOMENCLATURE
a, Bi c~ dc GO h, k
surface area of particles per unit volume of bed h,2RIk,, Blot number specific heat column diameter mass flux per unit of bed cross-sectional area fluid-to-particle heat transfer coefficient thermal conductivity kz thermal axial dispersion coefficient
400
300
200
IC
0
/
m.
5
t~ t',
hr
Fig. 3. Thermal energy availability for the rock-bin unit. SE VoL29, No, 3---C
205
nth temporal moment
Nu h,2R/k, Nusselt number Pr cpp,/k,Prandtl number Q thermal energy availability q rate of collection of thermal energy R particle radius r radial coordinate inside a particle
206 s T T, To T~ Ts v z
A.E. SAEZand B. J. McCoY Laplace parameter air temperature environment temperature inlet temperature initial temperature solid particle temperature fluid velocity axial coordinate in the column
Greek symbols • bed void fraction # fluid viscosity #n nth central moment p.'~ nth normalized moment p density Subscript s solid Acknowledgement--This research was supported by the U.C. Appropriate Technology Program under UCAT Project No. 80126SE-3998.
REFERENCES
l. M. Sagara, P. Schneider and J. M. Smith, Chem. Engng J. 1. 47 (1970). 2. N. Wakao, S. Kaguei and T. Funazri, Chem. Engng Sci. 34, 325 (1979). 3. M.-S. Razavi, B. J. McCoy and R. G. Carbonell. Chem. Engng J. 16, 211 (1978). 4. C. Vanden Broek and J. A. Clark, Von Karman Institute for Fluid Dynamics. Lecture Series Vol. 1, (1979). 5. J. A. Clark and V. Arpaci, Von Karman Institute for Fluid Dynamics. Lecture Series Vol. 1, (1979). 6. M. Riaz, Solar Energy 21,123 (1978). 7. P. J. Hughes, S. A. Klein and D. J. Close, J. Heat Trans[er 98, 336 (1976). 8. T. R. Mancini, J. L. Peterson and P. R. Smith, Heat Trans[er in Solar Energy Systems, 45. A.S.M.E. Proceedings (1977). 9. J. A. Clark, Von Karman Institute for Fluid Dynamics. Lecture Series Vol. 1 (1979). 10. A. E. Saez, M. S. Thesis, Dynamic heat transfer in packed beds--applications to solar energy storage, Department of Chemical Engineering, University of California, Davis, CA 95616 (1981).
0038~92X]82/090201-06503.0010 Pergamon Press Ltd.
DYNAMIC RESPONSE OF A PACKED BED THERMAL STORAGE SYSTEM--A MODEL FOR SOLAR AIR HEATING A. E. SAEZ and B. J. McCoY Department of Chemical Engineering,University of California, Davis, CA 95616, U.S.A.
(Received 18 September 1981; revision accepted 4 January 1982) Abstraet--A mathematical model for simulating the dynamic temperature response of a packed column to an arbitrary time-dependent inlet air temperature is developed. The model includes axial thermal dispersion as well as intra-particle conduction, features that have usually been neglected but can be important in solar energy applications. Solutions, presented in terms of moments of the temperature response to an impulse of heat at the inlet, can be evaluated by simple numerical quadrature. Results of the model compare favorably with experimental data found in the literature. The model is used to optimize heat storage in a rock bin system subject to a realistic transient inlet temperature.
INTRODUCTION
THEORY
A packed-bed thermal storage system, e.g. a rock bin, consists of a column packed with pebbles or some other heat-absorbing solid material through which air flows when energy is transferred. Such systems operate in one of three modes:
Our model is identical to that used by Sagara et al. [I] to analyze pulse response data for heat transfer in a column packed with porous spheres. Wakao et a/.[2] have recommended the same model for analyzing data to determine heat transfer coefficients in packed beds. Consider a column of length L packed with solid spheres of diameter 2R. We make the following assumptions: • the flow in the column is axial and incompressible • viscous heating is negligible • radiant heat transfer between particles is negligible • physical properties of the fluid are constant • surface temperature of a particle is uniform • heat losses from the column are negligible • temperature and velocity profiles in a column crosssection are uniform. Under these conditions, the equation for the interparticle fluid T(t, z) is
(i) Charge mode Air is heated in the collection system and enters the bed, increasing the bed temperature. The charge mode for solar heating occurs usually in the morning and early afternoon hours, when the energy collection is large. (ii) Discharge mode Cool air is heated by passing through the bed, usually at night. (iii) Stand-by mode No air flows through the bed. This mode occurs when further storage of heat cannot occur because of the decreasing inlet air temperature, and no air heating from the bed is required. This is the transition period between the charge and discharge modes. The accurate mathematical description of this timedependent system is the objective of the present work. While accuracy is required, we also desire simplicity in the use of final equations. In the following sections we present a model for an adiabatic (well-insulated) column operating under transient conditions. The model, which includes axial thermal dispersion and intra-particle heat conduction, is compared to more limited models in the literature. We show that the temporal moments of the response to an impulse of heat at the inlet can be used to construct the response to any time-varying inlet condition. An expression for the total energy stored in the bed as a function of time is developed, and used to calculate when the charge mode should end for a realistic solar input. This optimization provides a basis for improved system control.
02T
(1)
where e is the interparticle void fraction, p and cp are the density and heat capacity of the fluid, Vo= ev is the fluid volumetric flow rate per unit column cross-sectional area, as = 3(1- e)/R is the surface area of spheres per unit column volume, and hs is the fluid-to-particle heat transfer coefficient. The effective axial thermal dispersion coefficient, kz may account for longitudinal particle-to-particle conduction, as well as fluid conduction and hydrodynamic dispersion of heat axially in the column. The equation for the temperature inside a sphere of radius R as a function of time and radial position, Ts(t, r), is
OT~ 1 ~ ~r2 8T~ p.cps - ~ = k. ~ Tr \ ~ - / 201
(2)
A. E. SAEZand B. J. McCoy
202
where p, and cps are the mass density and specific heat capacity and ks is the thermal conductivity coefficient of the solid spheres. At the column entrance we have the arbitrary input temperature,
T(t, O) = To(t).
(3)
The moments may be derived from the Laplace transform of ~(t), m. = ( - I)"(d"~/ds")[~=o which follows directly from the definition of O(s). Letting
When the column is very long compared to the diameter of the spheres, we may apply the usual boundary condition [ 1],
T(t, oo) = finite.
(4)
t~= T - T, and
~,s: L - T,
At the surface of spherical particles we require the heat flux to be continuous,
h s ( T - T~(r : R)) = k~(dT, lOr)[r=R,
(5)
(12)
we may apply Laplace transformation to eqns (1)-(6), solve for 0(s, z), and then use eqn (12) to calculate the moment expressions [1],
and at the center of a sphere, ~rr [T~(t, r)],=o : 0.
(6)
We assume that the bed is initially at some reference temperature, T, T(0, z) = T,.
(7)
TA0,~ = ~.
(8)
too(z) = too(0)
(13)
#'t(z) = #I(0) + (1 + 8o)zlv
(14)
#2(Z): # 2 ( O ) + [ t $ 1 + ~ ( I L epcpv
MOMENTS
Sagara et a/.[1] have used moment methods to interpret heat transfer processes in packed columns. The present objective is to use moments to predict the dynamic response of a packed bed to an arbitrary temperature stimulus at the column inlet. The temporal moments are integrals over time of a distribution function, ¢J, f~
m, =Jo t"6(t)dt, n =0, 1 , 2 . . .
(9)
(15)
where 1 - ~p~cps 60- - - - -
e
From the coupling boundary condition (5), an order of magnitude analysis shows that when the Biot number, Bi = h~2Rlks, is not small, intraparticle effects cannot be ignored. Thus, for particles of small thermal conductivity or large diameter, or for high heat transfer between fluid and particle, which might occur for large fluid velocities, Equation (2) is necessary for the accurate description of packed bed heat storage. If B i ~ 1, intraparticle heat transfer is negligible, and at any crosssection of the column T = Ts = constant is a good approximation.
+t$o)212Z/V
(16)
pcp
8~ = 1 - E Ps c..._p~
pcp 15
1
5
+
"
(17)
The zeroth moment, the area under the curve of T - T, plotted versus time, can be interpreted as the amount of heat entering the column and moving conservatively along the column. The first moment is the average temporal position of the peak, and the second central moment indicates the width of the peak. For the present purposes we may consider that inlet fluid temperature pulse is a delta function,
T(O, t ) - T~ = 0o8(t)
(18)
so that too(Z)= 0o, the strength of the impulse, and #'~(0) = O, #2(0) = O. Because energy is stored in the solid spheres, we need the single-particle average temperature, defined as
3/?
(Ts) = ~--~
T~ P dr.
(19)
The normalized nth moment is defined by
#"
n = 0, 1 , 2 . . .
(10)
(t-#D"O(t)dt, n = 2 , 3 . . .
(11)
= m,]mo,
and the central moments are
#" =-moo
The moments of this temperature, relative to T, may be found by applying the single-particle averaging operation to eqn (2), finding the Laplace transform, and using eqn (12). We find the following expressions for the moments of the temperature pulse in the solid particles:
mos(Z) = too(Z)
(20)
Dynamic response of a packed bed thermal storage system /*[,(z) =/*[(z) + ~
(~~~ h - +- ~ )1l
2 4 /*=,(z)= /*=(z)+p,2 c;,R [ 13 2 ,17-Tk7+ sk-
(21)
I
+ hZ /• (22)
The first and second particle moments show delays relative to the fluid moments, and reduce to fluid moments for very large values of ks and hs, i.e. for large rates of intraparticle heat transfer. The dependence of the second moment on particle diameter, 2R, shows that broadening of a temperature profile within the column is decreased with smaller particles. Thus, since smaller particles are favorable for heat storage but have the disadvantage of increasing the pressure drop across the column, a trade-off of costs is required to determine the optimum particle size. The temperature response for the delta function input can be expressed as a series expansion in orthogonal polynomials[3]. Hermite polynomials, which are simple to use in this context, provide sufficient accuracy for most purposes. We write the series expansion as T - Ti = mo exp (- x 2) ~
a.(z)H.(x)
(23)
203
rates[3]. The temperature profile during the discharge mode may be expressed as the response to a step down in temperature if the bed temperature is uniform. When axial thermal dispersion and intra-particle effects are very important the time-response to a delta function input is not well represented by a Gaussian distribution. In that case, the use of Laguerre polynomial expansions is recommended, since the first term of the series will be a Poisson distribution function which will represent better the solution. In terms of Laguerre polymonials, L,, the expression is T-
T, -
mo [ M ] ~-'
(A--'l)! \/*~1
(-At)~__o
exp ~
(M)
k.L. klt,/ (28)
where t2
A ---.-/* J #2
(29)
The coefficients k, can be found by applying the orthogonality relationship: ko(z)= l, k~(z)=0, k2(z)= 0. Applying the convolution theorem gives the response for any kind of input
n=O
- t')T(t,z)- Ti= [,o' [To(t-~A'-'~.
where
T~] t~]) /At'\a-' exp ~--~) /At'\ dt'. (30)
t-/*~
x = V(2/.2)'
(24)
The first term of the expansion yields a Gaussian distribution in time centered at /*~ with variance /*2. The first three coefficients, a,(z), of the expansion can be found from the orthogonality relation for Hermite polynomials to be ao(z) = l/X/(21r/*2), al(z) = 0, and a2(z) = O. In Laplace transform space the solution to eqns (1)-(6) can be expressed as the product of the inlet temperature and the response to the delta function input; therefore, it follows from the convolution theorem that
For a step input, the solution is T - T~= (A - 1)! r A,
(31)
where F is the incomplete gamma function, defined by:
F( &-~2 ) = fo(~it/~2)x~-' e Xdx.
(32)
COMPARISON WITH EXPERIMENT 1
~* t
T(t, z)- T~- X/(2-1r/*2)Jo [To(t- t') - T~]exp [- (t' -/*~)2/2/z2] dt'.
(25)
For the special case of the step input, i.e. To(t)- T,=
O
, t< 0
(26)
-T,, t > 0 '
Equation (25) yields t ' t-/*l
/.1 (27)
This representation is known to be accurate when /*~X/(2/.2)~2.5, i.e. for long columns and low flow
We have compared the prediction of our model with the experimental data of Vanden Brock and Clark[4] for an insulated column 8.00 cm in dia. and 24.2 cm in length packed with steel spheres 0.556cm in dia. The void fraction of the bed was 0.386. Because the Biot number is small for the air/steel system of these experiments (Bi=0.01), intra-partiele conduction should not affect results substantially. Thus, the main advantage of the present model over that of Clark et al.[5] in describing these experiments is that our model includes axial thermal dispersion. Outlet temperature responses were measured for step inputs of the air temperature; experimental data for one run is shown in Fig. 1. For the present model, physical properties of air and steel are evaluated at the mean temperature (1/2) (To+ T~). To test this assumption we made calculations for properties evaluated at extreme temperatures in a rock bin, 26 and 65°C. The maximum difference in the
A. E. SAEZand B. J. McCoy
2O4 7O
v
i
,
l~O
2~I0
60
50
40
30
~;0 t,
4;0
~~0
8;0
7&O
800
sec
Fig. I. Response to a temperature step input: Solid line is experimental data, dashed line is the present theory, (7) Vanden Broek and Clark[4] theory, [] Riaz theory[6]. calculated temperature profiles was 2°C, even though the properties were evaluated at two temperatures almost 40°C apart. Heat transfer coefficients are estimated with the correlation[2], Nu = 2 + 1.1 Pr 1/3Re °'6.
(33)
Thermal axial dispersion coefficients are estimated from [2]
kz = 0.7k + epcpvR.
(34)
The results of the numerical calculations made by Vanden Broek and Clark[4], shown in Fig. 1, are based on heat transfer coefficients regarded as adjustable parameters to fit the data. Heat transfer coefficients for the Riaz model[6] were based on eqn (33). The results of the present model, with Hermite polynomial expansions, are plotted in Fig. 1. We also plot the predictions of Riaz[6], and of Vanden Broek and Clark[4]. The model of Riaz is a single-phase conductivity model. The model of Vanden Broek and Clark ignores axial dispersion and intra-particle conduction. The accuracy of the different theories with respect to experimental data for four runs was represented by an average absolute deviation, defined by
r ~' [Tcalc-Texpldt' 8=-~1 Jo
(35)
where the absolute difference between calculated and experimental temperatures is averaged over At, the duration of the experiment. The present model gave values of 6 between 2 and 4°C. The model proposed by Riaz yielded larger errors, sometimes twice as large, and tended to underestimate the experimental curve. The Vanden Broek and Clark model, while tending to overestimate experimental data, showed deviations similar in magnitude to those of the present model. Even though the mean absolute deviation 8 for the present model is 3°C, errors up to 8°C were observed. The cause of the errors is likely to heat losses from the column. Vanden Broek and Clark[4] reported radial
temperature variations across the column cross-section, which can be attributed to heat losses. The temperature variations are not due to velocity non-uniformities, because velocity was measured experimentally[4] and found to be uniform over most of the column length. The criterion for uniform radial temperature when heat losses occur is that the column diameter or overall column heat transfer coefficient be small, or that the radial thermal dispersion coefficient be large. An interesting extension of the present work would be to include the heat loss effect by adding a loss term to eqn (1). The overall heat transfer coefficient could be treated as an empirical parameter, whose value would be adjusted to fit experimental temperature profile data. Adjusting the value h~ in the present model did not improve the fit to the experimental data. Another possible cause for the deviation between the present model and experimental data is an overestimation of the axial thermal dispersion coefficient; the smaller the value of this coefficient, the sharper the temperature curve. The predicted curves are broader than those measured experimentally, indicating possibly that the correlation may provide a value of kz that is too large. Hughes et a/.[7], have modified the Schumann model with empirical corrections for intraparticle heat conduction and axial dispersion. When an energy loss term is added, a finite difference method is required to solve the system of partial differential equations. The present model may also be modified for heat losses, and once temporal moments have been determined, a single numerical integration provides the temperature response for an inlet function. Further work on the topic of energy losses will be reported at a future time. APPLICATION TO A SOLAR E N E R G Y S T O R A G E SYSTEM
As an illustration of the applications of the model to a solar energy storage system, a typical rock-bin unit was simulated. Consider a column of dia. 1.22 m and length 6.1 m packed with granite rocks 6 cm in dia. The objective is to find the hour of the day to discontinue the charge mode in order to store the maximum amount of energy. The hourly ambient temperatures and energy collection rates are taken from experiments[8] at the New Mexico State University Solar House on a sunny day, 7 February 1977. A typical air flow rate, 0.30kg/sm 2, is chosen. The physical properties of air are evaluated at 25°C, and the properties of the rock are found in Clark[9]. The basic data for the simulation are given in Ref. (7). The heat transfer coefficient and thermal axial dispersion coefficient are calculated, by using the correlations (33) and (34), to be 28.5 W/m2K and 9.64 W/mK respectively. We will assume that the bed is initially at a uniform temperature of T~ = 5.6°C, and the entering air is at ambient temperature. We consider that all the heat collected is transferred to the flowing air. With these assumptions the following heat balance applies to the collector where To is the outlet temperature of the collector and therefore the inlet to the rock-bin:
Dynamic response of a packed bed thermal storage system
q = Gocp(To- T~)dc21r[4
(36)
from which we obtain
To - T~ = 4qlcrGodc2cp + Te - T~.
(37)
Note that T, and q, and therefore To- Ti, vary with time. In Fig. 2, we have plotted To- T~ as a function of time. We approximated To- T~ by straight line segments over every hour interval to introduce To- T~ as data for the numerical integration of eqn (25). Since axial dispersion and intraparticle conduction are not large in this case, comparison with the Laguerre expansion showed the Hermite polynomial expansion to be the better approximation for the temperature response for air and particles. The total energy, Q, stored at a certain time can be expressed in terms of energy available relative to a reference level, which is taken to be the minimum night temperature, 5.6°C. We have
Q = ~R2~pcp f o L ( T - T,-)dz + IrR2(l - ~)p, cp~fo L (T~ - T~) dz,
(38)
where, since pscp~ >>pcp, the latter term dominates. The available energy Q can thus be calculated as a function
60
5C
4C tJ o 3C i 2C
5
tO
t,
,~
of time by numerically integrating results of eqns (37) and (25); Q is plotted versus time in Fig. 3. From Fig. 3, the time at which the thermal energy availability is maximum is t = 9.5 hr, i.e. 4:30 p.m. This is the time when the charge mode should end, and the stand-by mode should begin until energy for heating is required. Continuing the charge mode beyond this time will actually decrease the total stored energy. It is useful to note that the maximum available energy, 380M J, exceeds the space heating needs of an average house. However, since the process of withdrawing energy from the bed is not 100 per cent efficient, the total 380 MJ is not all usable. CONCLUSION
The model presented in this paper, advocated for applications of heat transfer in packed beds, has also been successful in its mass transfer form in describing chromatographic and breakthrough curve processes[3]. The model is based upon reasonable assumptions and contains transport parameters, principally the longitudinal thermal dispersion coefficient, kz, and fluid-to-particle heat transfer coefficient, hs, which have been well studied. Even though the model includes axial thermal dispersion and intraparticle effects, the solutions are rather simple and sufficiently accurate for practical application, so that their incorporation in designing procedures for solar energy is worthwhile. A computer program is available [10] that gives the air temperature and particle average temperature response for a packed bed subject to an arbitrary temperature input. The program calculates the temperatures by using Hermite and Laguerre expansions and includes up to the third term in the expansions. Although this program is in Fortran, the computations involve only numerical integration and parameter substitution, which are readily performed on a programmable calculator. Two kinds of polynomial expansions were used. Hermite expansions give better results for the case of low axial dispersion and low intra-particle effects. Laguerre expansions are better in the opposite case. The model was used to simulate a solar energy storage system undergoing real conditions to illustrate its manageability in the simulation of practical processes.
hr
Fig. 2. Rock-bin air temperature input. NOMENCLATURE
a, Bi c~ dc GO h, k
surface area of particles per unit volume of bed h,2RIk,, Blot number specific heat column diameter mass flux per unit of bed cross-sectional area fluid-to-particle heat transfer coefficient thermal conductivity kz thermal axial dispersion coefficient
400
300
200
IC
0
/
m.
5
t~ t',
hr
Fig. 3. Thermal energy availability for the rock-bin unit. SE VoL29, No, 3---C
205
nth temporal moment
Nu h,2R/k, Nusselt number Pr cpp,/k,Prandtl number Q thermal energy availability q rate of collection of thermal energy R particle radius r radial coordinate inside a particle
206 s T T, To T~ Ts v z
A.E. SAEZand B. J. McCoY Laplace parameter air temperature environment temperature inlet temperature initial temperature solid particle temperature fluid velocity axial coordinate in the column
Greek symbols • bed void fraction # fluid viscosity #n nth central moment p.'~ nth normalized moment p density Subscript s solid Acknowledgement--This research was supported by the U.C. Appropriate Technology Program under UCAT Project No. 80126SE-3998.
REFERENCES
l. M. Sagara, P. Schneider and J. M. Smith, Chem. Engng J. 1. 47 (1970). 2. N. Wakao, S. Kaguei and T. Funazri, Chem. Engng Sci. 34, 325 (1979). 3. M.-S. Razavi, B. J. McCoy and R. G. Carbonell. Chem. Engng J. 16, 211 (1978). 4. C. Vanden Broek and J. A. Clark, Von Karman Institute for Fluid Dynamics. Lecture Series Vol. 1, (1979). 5. J. A. Clark and V. Arpaci, Von Karman Institute for Fluid Dynamics. Lecture Series Vol. 1, (1979). 6. M. Riaz, Solar Energy 21,123 (1978). 7. P. J. Hughes, S. A. Klein and D. J. Close, J. Heat Trans[er 98, 336 (1976). 8. T. R. Mancini, J. L. Peterson and P. R. Smith, Heat Trans[er in Solar Energy Systems, 45. A.S.M.E. Proceedings (1977). 9. J. A. Clark, Von Karman Institute for Fluid Dynamics. Lecture Series Vol. 1 (1979). 10. A. E. Saez, M. S. Thesis, Dynamic heat transfer in packed beds--applications to solar energy storage, Department of Chemical Engineering, University of California, Davis, CA 95616 (1981).