Thermal performance of packed bed thermal energy storage units using multiple granular phase change composites

Applied Energy 86 (2009) 2704–2720

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Thermal performance of packed bed thermal energy storage units using multiple granular phase change composites Mohamed Rady * Mechanical Engineering Department, High Institute of Technology, Benha University, Benha 13 512, Egypt

a r t i c l e

i n f o

Article history: Received 13 February 2009 Received in revised form 17 April 2009 Accepted 21 April 2009 Available online 19 May 2009 Keywords: Energy storage Phase change Packed bed Multiple PCM

a b s t r a c t The present article reports on the utilization of multiple granular phase change composites (GPCC) with different ranges of phase change temperatures in a packed bed thermal energy storage system. Small particle diameter of GPCC allows simple mixing of two or three ranges of GPCCs in a packed bed for enhancement of storage unit performance. Experiments have been carried out to characterize the phase changing characteristics of two GPCCs chosen for this purpose. Packed bed column experiments have been carried out to provide basic understanding of the heat transfer process in the composite bed consisting of a mixture of GPCCs at different values of mixing ratio. A mathematical model has been developed for the analysis of charging and discharging process dynamics. Once validated, the model has been used to perform a parametric study to investigate the overall bed performance at different values of mixing ratio and Reynolds number. An optimization of the value of mixing ratio has been obtained based on the overall charging and discharging times as well as the exergy efficiency. It has been demonstrated that, as compared to the use of single GPCC, careful choice of the mixing ratio of GPCCs in a composite bed can result in a significant enhancement of the overall storage unit performance. As compared to the use of multiple sequential layers of GPCCs, using units composed of a mixture of GPCCs with an optimized mixing ratio results in a remarkable improvement of the unit performance without limitations on the charging and discharging directions during practical applications. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The use of phase change materials (PCMs) for latent heat thermal energy storage (LHTES) has gained considerable importance in a wide variety of applications such as solar energy utilization, industrial waste heat recovery, and electrical power load shifting. Traditional LHTES systems usually use a single PCM. However, the use of multiple PCMs for LHTES systems has acquired great significance [1–4] because of their potential for superior thermal performance. Studies on the utilization of multiple PCM can be classified according to the geometry of the storage system. A great number of studies have been concerned with shell and tube heat storage units. In these units, the PCM fills the annular shell space around the tube, while the heat transfer fluid flows within the tube and exchanges heat with the PCM. Different PCMs have been arranged either in the form of segments in the shell side of the heat exchanger along the flow direction [5–7] or in the form of multiple slabs perpendicular to the flow direction [8–11]. For such configurations the low thermal conductivity of PCM is a limiting performance parameter. For packed bed configurations, the PCM is * Present address: TREFLE, Bordeaux University, F33607 Pessac, France. E-mail address: [email protected] 0306-2619/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.apenergy.2009.04.027

contained in spherical or cylindrical capsules and the heat transfer fluid is in direct contact with the capsule surface [1,2,12]. This offers the advantage of high heat transfer rate between the fluid and PCM. Among the studies reported that used multiple PCM segments arranged in the shell side along the flow direction, Gong and Mujumdar [5] proposed a solar receiver store consisting of multiple PCMs and found that the fluctuation of the outlet temperature of the heat transfer fluid can be greatly damped as compared with a single PCM. Haiting et al. [6] proposed a PCM receiver model composed of three different phase change temperature materials for the NASA 2 kW solar dynamic power system. Their results showed that it is possible to improve the receiver performance and to reduce both the fluctuation of working fluid temperature and the weight of the heat receiver by using multiple PCMs. A model for the shell and tube LHTES unit using multiple PCMs is developed and solved numerically by Fang and Chen [7] to investigate the effects of multiple PCMs’ melting temperatures and fractions on the melted fraction, stored thermal energy and fluid outlet temperature of the LHTES unit. Numerical results indicated that PCMs’ fractions and melting temperatures play an important role in the performance of the LHTES unit. As a result, appropriate choosing of multiple PCMs is very significant for the performance improvement of the LTES unit.

M. Rady / Applied Energy 86 (2009) 2704–2720

2705

Nomenclature cp d E fl hsa H Hls k L MGPCC Mair Nu Pe Pr Re S t T Tm TGPCC Tair Ud Z

specific heat at constant pressure [J/(kg K)] Particle diameter (m) exergy (J) liquid fraction fluid-to-particle heat transfer coefficient [W/(m2 K)] enthalpy (J/kg) PCM latent heat of fusion (J/kg) thermal conductivity [W/(m K)] bed height (m) total mass of GPCC in the bed (kg) total mass of air in the bed (kg) Nusselt number Peclet number Prandtl number of gas Reynolds number volume specific surface area (m2/m3) time (s) temperature (K) phase change temperature (K) GPCC temperature, charging end (K) air temperature, charging end (K) Darcy velocity (m/s) distance along the bed height (m)

Using multiple PCM slabs, Gong and Majumdar [8] developed a 1-D numerical model to analyze cyclic energy storage systems and found that the instantaneous heat flux on the heat transfer surface could be enhanced using a composite PCM slab arrangement when energy is charged and discharged from the same side of the slab. Gong and Mujumdar [9,10] then investigated the charge and discharge kinetics of LHTES in composite PCM slabs and concluded that the magnitude of enhancement depends upon the arrangement of PCMs and their thermophysical properties. Shaikh and Lafdi [11] studied combined convection–diffusion phase change heat transfer process in varied configurations of composite PCMs slabs. As compared to using a single PCM slab, the numerical results from their parametric study indicated that the total energy charge rate can be significantly enhanced by using composite PCMs. For packed bed configurations, Farid and Kanzawa [1,2] proposed a heat storage module consisting of a number of cylindrical capsules filled with phase change materials with air flowing across them. They observed a significant improvement in the rate of heat transfer during energy charge and discharge when phase change materials with different melting temperatures were used. Adebiyi et al. [12] reported that the efficiency of a packed bed LHTES using five PCMs exceeded that of a unit using a single PCM by as much as 13–26%. One of the main purposes of the previous research on this topic is to enhance the charging and/or discharging rate(s) of LHTES systems by decreasing the freezing and/or melting time of PCMs. Wang et al. [13–15] proposed a 1D heat conduction model of LHTES system (spherical PCM, an infinite cylindrical PCM, and an infinite flat plate PCM) using materials of different phase change temperature (PCT) distributions and obtained the optimum linear and optimum parabolic PCT distributions of composite PCMs. The results obtained showed that the application of PCMs with particular PCT distributions to LHTES systems can enhance charging and discharging rates significantly. As noted by Wang et al. [15], the main objective of associating multiple PCMs is to realize a homogeneous phase change process (HPCP) in which any PCM starts and ends its phase change simultaneously, and the phase change rate of each component is identical and remains constant during the

Greek symbols g exergy efficiency a thermal diffusivity (m2/s) e porosity; volume fraction l air dynamic viscosity (kg/m s) m air kinematic viscosity (m2/s) q density (kg/m3) Subscripts a air c charging d discharging dis dispersion in inlet conditions init initial conditions l liquid phase m1 onset of melting m2 endset of melting avgm average melting temperature (K) o reference s solid phase

phase change process. HPCP is an ideal phase change process and has no obvious interface in the PCM. In real applications, finite numbers of PCMs are used to approximate the HPCP. Their phase-change temperatures must be the volumetrically averaged values of the optimum PCT distribution. In recent years, granular phase changing composites (GPCC) obtained by micro or macro-encapsulating PCM in highly porous solid structures with protecting envelopes have been developed using different encapsulation techniques [16]. The encapsulation process ensures that the PCM, when in the liquid form, does not leak out of granulate. Small particle diameters of GPCC open a wide range of applications and integration in energy systems. Their application in packed bed latent heat energy storage is advantageous as due to high heat transfer rate associated with large surface area of direct contact between the heat transfer fluid and the GPCC. Utilization of multiple GPCC with different ranges of PCT in a packed bed storage system represents the main interest of the present work. A possible traditional method includes the utilization of sequential layers of multiple GPCCs. Another interesting method, proposed in the present study, can be obtained by simple mixing of two or three ranges of GPCCs to offer a homogeneous phase change process along the column length. The phase change temperature range should be compatible with the operating temperature range of charging and discharging processes. Up to the author’s knowledge, the idea of simple mixing of GPCCs in packed bed units has not been studied in the literature. In the present study, two GPCCs (Rubitherm GR27 and GR41) are employed in a packed bed storage system. GR27 and GR41 are considered as two examples representative of GPCC materials and the results obtained in the present study are useful in evaluating the application of such materials in LHTES systems. The interest in these materials is motivated from the potential of their recent applications in desiccant cooling systems [17] and in air conditioning systems of buildings [18]. The composition of the granule is 65% ceramics and 35% paraffin wax by weight with particle diameter between 1 and 3 mm. The phase change temperature range of the two GPCCs (21–29 °C for GR27 and 31–45 °C for GR41) covers an interesting range for practical applications in low temperature

2706

M. Rady / Applied Energy 86 (2009) 2704–2720

LHTES systems. Packed bed column experiments have been carried out using different mixing ratios of GR27 and GR41 and the results are compared with the cases of using simple GR27 and GR41 packed beds. A mathematical model for the packed bed column consisting of separate energy equations for the solid GPCC particles and the heat transfer fluid has been developed. The model accounts for phase change in the bed including multiple GPCCs. Once validated by comparison with the experimental results, the present mathematical model has been used to perform a detailed parametric study on the energy storage and discharge characteristics of the packed bed using multiple GPCCs. Two different configurations of multiple PCMs have been considered in the present study. In the first configuration, a mixture of GR27 and GR41 with different mixing ratios is considered. The optimum mixing ratio is determined by maximizing the energy charge and discharge rates during charging and discharging processes. An optimization based on exergy analysis of the packed bed is also performed. For the purpose of comparison, the optimum configuration of packed bed unit consisting of a mixture of GPCCs is compared with the configuration using successive segments of GPCCs arranged along the bed length. 2. Characterization of materials Screen analysis of GR27 and GR41 shows that they contain 0.1% of particles with diameter dp > 4 mm, 7.5% with dp = 3.0–4.0 mm, 89.5% with dp = 1.0–3.0 mm, 2.3% with dp = 0.5–1.0 mm, and 0.6% with dp < 0,5 mm. Characterization results for GR27 and GR41 are shown in Table 1. These include the start and end temperatures of phase change (Tm1, Tm2), the specific heat of solid and liquid phases (cps, cpl), and the latent heat of phase change (Hls) which is the difference in enthalpy at temperature Tm2 and Tm1. These characteristics have been obtained using DSC and T-history methods as reported in [19,20]. GR27 and GR41 have approximately equal values of latent heat of phase change.

Table 1 Characterization of GR27 and GR41. Sample

Tm1 (°C)

Tm2 (°C)

cps (J/kg K)

cpl (J/kg K)

Hls (J/kg)

GR27 GR41

29 45

21 31

2044 2784

1921 2745

64,850 65,895

Other properties: ks = 0.2 W/m k, q = 1360 kg/m3, see http://www.rubitherm.com.

The evolution of latent heat during charging and discharging of LHTES unit as function of temperature is an important issue for modeling, optimization, and design of the unit. This evolution is expressed in terms of the variation of liquid fraction fl of GPCC with temperature as follows.

fl ðTÞ ¼

HðTÞ  Hs ðT m1 Þ Hls

ð1Þ

H(T) is the enthalpy at any temperature (T) and Hs(T) is the solid phase enthalpy. The functional dependence of liquid fraction on temperature for GR27 and GR41 during the phase change process obtained using the T-history method is shown in Fig. 1. It should be noted that, the evolution of latent heat with temperature during phase change differs for the two materials. As compared to GR41, GR27 is characterized by the evolution of a considerable amount of latent heat over a narrower temperature range (22–27 °C). The evolution of latent heat for GR41 is more distributed over the phase change temperature range (31–45 °C). These characteristics shows that, if mixed together or arranged in successive segments, these two materials covers a range of phase change temperature between 21 and 45 °C that is suitable for a practical range of low temperature energy storage applications with charging and discharging temperatures in the order of 50 °C and 15 °C, respectively. 3. Packed bed experiments Packed bed column experiments have been carried out in the present study. The main purpose of the experiments is to provide basic understanding of the heat transfer process in the composite bed and supply experimental data for validation of the numerical model. A schematic diagram of the experimental set-up is illustrated in Fig. 2. The packed column is made of cylindrical tube of PVC tube of 45 mm inside diameter and 2.5 mm wall thickness. The instrumented test section is 200 mm height. The bed test section is insulated with glass wool of 5 cm thickness. The packing (storage medium) consists of a homogeneous mixture of GR27 and GR41 granulates. The mixing of GR27 and GR41 has been performed on a volumetric basis. That is to have a mixture of mixing ratio 1:1, two equal volumes of GR27 and GR41 are thoroughly mixed before introduction into the column. Mixtures with mixing ratios of GR27/GR41 = 1.0, 0.5, and 0.25 have been prepared and tested for different values of Reynolds number. Air flow is produced by a centrifugal blower with variable mass flow rate. After passing thorough the PCM packed bed, it is ex-

Fig. 1. Liquid fraction–temperature relation for GR27 and GR41.

M. Rady / Applied Energy 86 (2009) 2704–2720

2707

Fig. 2. Schematic of the experimental set-up for packed-bed.

hausted from the top and then the flow rate is measured using Brooks R 215 B variable area float meter with measurement range of 0–130 l/h at standard conditions of 20 °C and 760 mm Hg. The flow rate is corrected using temperature measurements just before the flowmeter inlet. The uncertainty in flow rate measurement is ±5%. The inlet air is heated by electrical heaters and mixed so that the temperature profile at the inlet bed is uniform. Copper constantan thermocouples of 0.1 mm diameter are used to measure the temperature at specific locations inside the test section, as shown in Fig. 2. Air temperatures are also measured at the inlet and at the outlet of the test section. Thermocouples are connected to an HP data acquisition system for continuous monitoring and recording of temperatures in date files at specified time intervals for further analysis. The accuracy of temperature measurements is ±0.2 °C. The bed temperature is uniform at the beginning of each experiment. For charging process, reported in the present study, the experiment starts by switching on the electric heaters with the required power using the autotransformer and allowing hot air flow at the desired rate to enter into the column. Experiments have been carried out at different levels of inlet air temperature. This type of control results in a time varying inlet temperature at the bed inlet. Least square fitting of the data of inlet bed temperature with time is used for further employment in numerical simulations of a given experiment. It should be mentioned that, this time variation of inlet bed temperature is close to real situations where the inlet bed temperature is not constant and varies with time. It represents a realistic dynamic behavior of the bed that is useful for analysis and validation of numerical models. Fig. 3 shows the time variation of temperature of the inlet air (Tin), at 4 cm from the bed inlet (location T1 and T8 in Fig. 2), at 8 cm ((location T3 and T9 in Fig. 2), and at the bed top section (location T7 in Fig. 2) for different values of mixing ratios of GR27 and GR41. Due to the difficulty of repeating experiments for the same conditions, selected experiments are presented in Fig. 3. The differences in the time variation of inlet bed temperature with time are due to different values of electric heating power. The Reynolds number Re is calculated based on the Darcy velocity at the inlet of the test section and using an average particle diameter of 2 mm. The differences in measured temperature values between (T1 and T8) and (T3 and T9) give an idea about the uniformity of temperature in the radial direction. It should be noted that the experiments are terminated when a sufficient

length of the column is charged and all phases of the process are experienced within this column part. The completely charged column length corresponds to the positions where the bed temperature approaches the inlet air temperature. The time variations of temperature at these locations represent all the process dynamics during charging. The temperatures at the bed top section and outlet air temperature increase at a lower rate as can be seen in Fig. 3. Temperature variations within the length of charged column are used in the following analysis and validation. During charging of columns containing only GR27, see Fig. 3, the bed temperature increases relatively rapidly in the initial stage of the process below the melting temperature range. In the range between 22 °C and 27 °C, the rate of increase in bed temperature decreases as due to the phase change of GR27. After the end of melting process, the rate of increase in bed temperature is high. Fig. 3 also shows the results for GR41 column which exhibits a lower flattening of the temperature profile during phase change. This is compatible with GR41 latent heat evolution obtained during material characterization. The observed phase change temperatures from the experiments agree with the results from the characterization experiments. It is also interesting to note that the measured temperatures T1 and T8 are approximately equal during the initial sensible heat change period. They begin to exhibit a remarkable difference during phase change that becomes obvious during PCM melting at a given section. Small radial temperature difference during sensible heating and cooling periods are attributed as due to the heterogeneous nature of the bed particles that are manifested as local variations in the bed porosity. Minor changes in the PCM content of GPCC is an additional factor that may affect the radial temperature difference during phase change. The volume fraction of PCM at a given radial position is actually not constant. Different phase change times are experienced at different radial positions. Locations where complete melting is finished first (because of lower value of PCM fraction) begin to show a higher rate of increase in temperature with time as due to sensible heating. On the other hand, locations where the melting process is not yet terminated experience a small change in temperature as due to the absorption of latent heat. The difference in radial temperature is also function of the shape of variation of latent heat with temperature. Latent heat evolution over a low temperature range would increase the differences in temperatures. Thus the differences in temperatures along the radial direction are more remarkable for the case of GR27 column. This can be explained by the

2708

M. Rady / Applied Energy 86 (2009) 2704–2720

content in the mixture. Also, the difference in radial temperature variation decreases with the decrease of GR27/GR41 ratio in the mixture. From the point of view of latent heat energy storage, the charging and discharging rates depend on the phase change temperature range, the manner of evolution of latent heat with temperature, and the temperature difference between the heat transfer fluid and the phase change temperature (PCT). GR27 has a significant evolution of latent heat in a relatively low and narrow temperature range (21–27 °C) that may be advantageous during charging of the packed bed. However, the performance of GR27 packed bed during discharging is limited by the low temperature difference between the heat transfer fluid and the PCT. On the other hand, GR41 characteristics favour an improved packed bed performance during discharging as due to its high and relatively wider range of PCT. The overall bed performance during charging and discharging can be enhanced by associating GR27 and GR41 in a single packed bed. A parametric study is then needed to reach an optimized configuration of the packed bed with well defined charging and discharging conditions. Due to the difficulty of repeating experiments at the same working conditions for different column compositions, detailed experimental comparisons and parametric studies are not possible. Therefore, the present experimental data have been used to validate a mathematical model developed for this purpose. Once validated, the mathematical model is used to carry out a parametric study for charging and discharging conditions as explained in the following sections.

4. Mathematical modeling and validation Detailed modeling of the heat transfer and fluid flow processes that take place in such arrangement of GPCC is quite complex. The rate of heat transfer to or from the solid in the packed bed is a function of the physical properties of the fluid and the solid, the flow rate of the fluid, and the physical characteristics of the packed bed. The governing equations for heat transfer in the packed bed are obtained by considering a REV containing GPCC particles and air with volume fractions es and ea, (es + ea) = 1. The solid phase is composed of two GPCCs of volume fractions es1 and es2 where es = es1 + es2. They consist of the fluid and solid (GPCC) energy equations and take into account the thermodynamic non-equilibrium between the two phases of different temperatures (Ta, Ts). The air flow passing through the porous bed is considered to be steady and uniform with superficial (Darcy) velocity, Ud. Heat conduction in the wall material, of small thickness and thermal conductivity, and radiant heat transfer, for the present low differences in temperature, are assumed to be negligible. The air and GPCC energy equations are written as:

@ðqa ea cpa T a Þ @ qa U d cpa T a @ @T a þ ¼ ea ka þ hsa Ssa ðT s  T a Þ @t @z @z @z @ðqs es cps T s Þ @ @T s ¼ es ks  hsa Ssa ðT s  T a Þ @t @z @z   @fl1 @fl2  qs1 Hls1 es1 þ qs1 Hls2 es2 @t @t Fig. 3. Experimental variation of bed temperatures with time at different mixing ratios of GR27 and GR41, Re = 7.29.

functional dependence of latent heat evolution on temperature for GR27 which shows a considerable amount of latent heat evolution in a narrow temperature range. For packed bed columns consisting of a mixture of GR27 and GR41, Fig. 3 shows that flattening of temperature during phase change decreases with the decrease of GR27

ð2Þ

ð3Þ

The inter-phase heat transfer is described by the heat transfer coefficient (hsa). The value of cps and qs corresponds to the volume averaged specific heat and density of the solid phase, qs = es1qs1 + es2qs2 and cps = es1cPs1 + es2cPs2. The evolution of latent heat during phase change of the GPCC is accounted for by using the source term (qs1Hls1es1ofl1/ot + qs2Hls2es2ofl2/ot) in Eq. (3). Where fl1 and fl2 are the liquid fractions of GPCC undergoing phase change and Hls1 and Hls2 are the latent heat of fusion of the PCM material. The value of fl is related to the temperature Ts and determines

M. Rady / Applied Energy 86 (2009) 2704–2720

the rate of evolution of latent heat during phase change. An iterative technique has been adopted for the calculation of liquid fraction for each PCM, fl1 and fl2.

fl1iþ1 ¼ fl1i þ n1

and fl2iþ1 ¼ fl2i þ n2

ð4Þ

where i refers to the iteration index and n is the correction. The corrections n1 and n2 are written as:

n1 ¼

cps ðT s  T PC1 Þ Hls1

and n2 ¼

cps ðT s  T PC2 Þ Hls2

ð5Þ

At convergence n1 = n2 = 0 and the temperature of phase changing nodes is equal to the phase change temperatures TPC1 and TPC2. The present characterization of GPCC shows that the phase change process takes place over a temperature range (between Tm1 and Tm2). The value of TPC used in Eq. (5) is therefore function of the liquid fraction TPC = f(fl). This functional relationship is a characteristic of the phase change material as shown in Fig. 1. The following constraints have been also applied during iterations.

fl1iþ1 ¼ max½0:0; fl1iþ1 ;

fl2iþ1 ¼ max½0:0; fl2iþ1 

ð6Þ

fl2iþ1 ¼ min½1:0; fl2iþ1 

ð7Þ

And

fl1iþ1 ¼ min½1:0; fl1iþ1 ;

Under the conditions of equal values of latent heat for the two phase change materials Hls = Hls1 = Hls2, the global value of liquid fraction fl is the summation of liquid fractions of PCMs in the REV.

fl ¼ es1 fl1 þ es2 fl2

ð8Þ

The fluid-to-particle heat transfer coefficient hsa is calculated using the correlation developed by Galloway and Sage [21]: 1=2

Nu ¼ 2:0 þ c1 Red Pr1=3 þ c2 Red Pr1=2

ð9Þ

where Red = Udd/t and Pr = t/a. For randomly packed bed of spheres c1 and c2 are 1.354 and 0.0326, respectively. Beasley and Clark [22] used higher values of these coefficients (c1 = 2.031 and c2 = 0.049), resulting in an increase of about 50% in the heat transfer coefficient, and obtained better agreement between the numerical and experimental results. It should be noted that in Eq. (3), ks represents the thermal conductivity of the solid GPCC. It is a physical property determined by measurements. The value of ks (0.2 w/m K) reported by the manufacturer has been used in the present study, see Table 1. However, the effective thermal conductivity of the fluid phase (ka) consists of the stagnant and the dispersion conductivity.

ka ¼ kas þ kdis

ð10Þ

The dispersive component of the longitudinal thermal conductivity in the porous media was evaluated as function of the Peclet number [23,24] as follows:

kdis ¼ 0:022ka kdis ¼ 2:7ka

Pe2d 1e

Ped

e

1=2

ðPed < 10Þ

ð11Þ

ðPed > 10Þ

To complete the mathematical description of the present problem, the initial and boundary conditions are specified as follows: at t = 0: Ta = Ts = Tinit, at inlet z = 0: Ta = Tin, oTs/oz = 0, at outlet z = L: oTs/oz = oTa/oz = 0. Where Tinit is the initial bed temperature and Tin is the inlet air temperature. Eqs. (2) and (3) are discretized by using the finite volume method employing an implicit time scheme. The upwind scheme is used for discretization of the convection terms and the central difference scheme is used for discretization of the diffusion terms. The

2709

resulting equations are solved using the Gauss–Seidel iteration method. A time step of 0.001 s and a grid mesh of 0.1 mm are used in the simulations to ensure numerical stability and accuracy. The simulation results are independent of the grid size and time step used. Validation of the numerical model has been carried out by comparison with experimental results of packed bed column. Figs. 4 and 5 show complete comparisons of experimental (T1, T8, T3, T9, T7, Tout) and numerical simulation (T1 Sim, T3 Sim, T7 Sim, Tout Sim) results during charging of GR27 and GR41 columns. In these calculations, an average bed porosity of 0.42 is used, and the phase changing characteristics resulting from the T-history method, see Fig. 1, are employed for the functional relationship of liquid fraction and temperature. The numerical model reproduces the variation of bed temperatures with time during different phases of the process including sensible heating and latent heat evolution that are dependent on the bed material. The agreement between experimental and numerical results is reasonable. This agreement is acceptable in view of the uncertainty in model parameters that may affect the accuracy of agreement. These are related to uncertainties in the phase change characteristics of GPCC, bed porosity, fluid–particle heat transfer coefficient, particle diameter, and axial dispersion. Separate studies have been carried out to analyse the sensitivity of agreement between experimental and numerical results to the above model parameters. The results show an important sensitivity of agreement to the phase change characteristics of GPCC and estimation of bed porosity. The sensitivity of agreement to the estimation of particle-to-fluid heat transfer coefficient and the axial dispersion are found to be less important. Applications of different functions of fl  T results in different evolution of latent heat during phase change that affects the variation of bed temperature with time. The uncertainty in the functional form of fl  T is thus a source of error that leads to differences between numerical simulations and experiments. Another important source of uncertainty is related to packed bed porosity that has been the subject of extensive research [25,26]. Literature values vary for the bed voidage of a randomly packed bed of spheres at large column to particle diameter ratio. This variance is ascribed to the packing mode. Due to the stochastic nature of the packing process, accurate and reliable predictions of bed porosity are elusive. In the present study, comparisons of experimental and numerical results obtained using different values of bed porosity show that the effect of bed porosity is significant. The local changes in porosity can lead to large variations in predicted velocity profile and therefore non-uniform head loss along the packed bed. The accurate prediction of local porosity is therefore important for predicting heat transfer in packed beds. This is especially true for packed beds having a particle size and shape distribution as the present application. It seems that some form of multidimensional statistical description of the packing will be necessary for performance modeling. The average bed to particle diameter employed in the present study 25 provides a central region essentially free from the wall effect on void fraction. However, near the walls, the void fraction increases and may result in flow channeling. Comparisons of experimental and numerical results for different values of mixing ratio are shown in Fig. 6. As explained in Section 3, temperature variations in the completely charged portion of the column are used in these comparisons. Namely, the temperatures at points 1 and 8, see Fig. 2. Again, the present model adequately predicts the time variation of bed temperature for different values of GR27/GR41 mixing ratios. Different values of mixing ratio affect the time evolution of latent heat and temperature. In general, the present model predicts all the physical features

2710

M. Rady / Applied Energy 86 (2009) 2704–2720

Fig. 4. Comparison of experimental and numerical results (GR27, Re = 7.29).

Fig. 5. Comparison of experimental and numerical results (GR41, Re = 7.29).

of the process and can be used for further parametric studies reported in the following section. 5. Parametric study Numerical simulations have been carried out to investigate charging and discharging kinetics of packed bed column. The col-

umn length and diameter are taken as L = 200 mm and D = 45 mm. Different mixing ratios of GR27/GR41 (0.2, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.66, 2.5, 5.0) have been considered for different values of Reynolds number (Re = 10, 20, 30). Constant inlet air temperature is adopted during charging (Tin ch) and discharging (Tin dis) processes. For proper comparisons, the choice of Tin value is made in order to keep equal temperature difference between inlet air

M. Rady / Applied Energy 86 (2009) 2704–2720

2711

Fig. 6. Comparisons between numerical and experimental results for different mixing ratios, Re = 7.29.

temperature and phase change temperature during charging and discharging. That is (Tin ch  Tm2) = (Tin dis  Tm1), where Tm2 is taken as the endset temperature of phase change of GR41 and Tm1 is taken to be equal to onset temperature of phase change for GR27. This results into charging and discharging temperatures of 50 °C and 16 °C, respectively. Performance of the packed bed unit

is analyzed in terms of energy charging and discharging rates. Charging and discharging rates are important performance parameters that are usually used in the analysis and performance evaluation of LHTES systems [10,11]. These are estimated in terms of the rate of energy added or extracted from air and the accumulated energy stored or discharged from the unit as follows:

2712

M. Rady / Applied Energy 86 (2009) 2704–2720

Fig. 7. Variation of energy charging rate and accumulated energy charged during charging of packed bed at different values of GR27/GR41 mixing ratio (0.2, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.66, 2.5, 5.0) and Re = 10.

_ air cpa ðT in  T out Þj and Q ¼ q_ ¼ jm

Z

t

_ qdt

ð12Þ

0

Fig. 7 shows the variation of energy charging rate and accumulated energy charged with time during charging process of packed bed unit at different values of GR27/GR41 mixing ratio. It can be observed that, at the beginning of the process the value of q_ is constant and high. The value of q_ begins to decrease when the thermal front reaches the unit exit. After the arrival of thermal wave front to the unit outlet, the unit dynamics during charging can be divided into three main stages. The first stage is characterized by sensible heat exchange between the heat transfer fluid and the GPCC during which the value of q_ decreases rapidly. This is followed by a stage during which the value of q_ decreases slowly as due to latent heat energy storage. At the end of phase change, the value of q_ decreases rapidly and the charging process is terminated. For GR27 only unit, the first stage is very short and a substantial amount of energy storage occurs during phase change with high _ This can and approximately constant values of charging rate (q). be explained by high heat transfer rates as due to large temperature difference between the inlet air temperature and the phase change temperature of GR27. For GR41 only unit, the phase change stage occurs relatively later and the charging rates during this stage are relatively lower and decreasing with time. This can be explained by the low temperature difference between the charging air temperature and the phase change temperature of GR41 and also by the latent heat evolution characteristics of GR41 which occurs over a larger temperature range. Fig. 7 also shows the variation of energy charged with time. The accumulated energy charged reaches its maximum in a shorter time for GR27 only unit as compared to GR41 only unit. The maximum values of Q is equal for GR41 and GR27 units as due to equal values of specific heat and latent heat considered in the present analysis. Using units composed of mixtures of GR27 and GR41, one can observe that increasing the ratio of GR27/GR41 in the mixture increases the charging rates and decrease the charging process time.

Fig. 8 shows the variation of energy discharging rate and accumulated energy discharged with time during discharging process of packed bed unit at different values of GR27/GR41 mixing ratio. Similar process stages can be identified as noted during charging process. However, in contrast with the charging process, high values of discharging rate and lower discharging process time are obtained using GR41 only unit as compared to GR27 only unit. The increase of GR27 content in the mixture decreases the discharging rate and increases the process time. It should be mentioned that, one of the main requirements of LHTES units is to enhance the charging and/or discharging rate(s) of LHTES systems by decreasing the freezing and/or melting time of PCMs. Optimization of the overall unit performance should consider both charging and discharging times. The variation of charging, discharging, and total charging and discharging times for units composed of mixture of GR27 and GR41 at different values of mixing ratio (GR27/GR41) are shown in Fig. 9. The process time is calculated as the time at which the values of outlet and inlet air temperatures are equal which corresponds also to zero values of charging and discharging rates. The processes times are normalized with respect to the corresponding processes times of GR27 only units. Normalized charging time decreases and discharging time increases with the increase of GR27 ratio in the mixture. As expected the normalized values of process times are approximately independent of the value of Reynolds number. This can be explained by the fact that, for the same value of Reynolds number, the overall heat transfer coefficient is constant for all bed configurations. Increasing the Reynolds number decreases the charging and discharging times for all bed configurations in the same manner as due to equal changes in the overall heat transfer coefficient. The dependence on the Reynolds number can only arise from the relative importance of heat advection and diffusion. However, the present results show that, for the present range of Reynolds number, heat diffusion is of minor importance as due low value of thermal conductivity of air and bed material. Therefore, for a given bed material, the normalized charging and discharging

M. Rady / Applied Energy 86 (2009) 2704–2720

2713

Fig. 8. Variation of energy discharging rate and accumulated energy discharged during discharging of packed bed at different values of GR27/GR41 mixing ratio (0.2, 0.4, 0.5, 0.6, 0.8, 1.0, 1.25, 1.66, 2.5, 5.0) and Re = 10.

Fig. 9. Effect of mixing ratio on normalized charging and discharging times.

times are approximately independent on the value of Reynolds number. The variation of total charging and discharging time with mixing ratio has a minimum for a value of GR27/GR41 between 0.8 and 1.0. This minimum corresponds to an optimized working point taking into account both charging and discharging processes. Using packed bed units composed of a mixture of GR27 and GR41, the overall unit performance, expressed in terms of the amount of reduction of total

charging and discharging process times, can be enhanced by about 15% for GR27/GR41 ratio between 0.8 to 1 as compared to GR27 only unit. The composite unit benefits from favourable characteristics of GR27 during charging and GR41 during discharging. Fig. 10 shows the time variation of outlet air temperature from the unit. It can be observed that, using GR27/GR41 = 1, the outlet air temperature approaches the value obtained using GR27 only bed during charging and GR41 only bed during discharging.

2714

M. Rady / Applied Energy 86 (2009) 2704–2720

Fig. 10. Time variation of outlet air temperature during charging and discharging, Re = 10.

More insight into the thermal phenomena involved can be obtained by analyzing the temperature and liquid fraction variation along the bed height and making reference to the governing energy equations, Eqs. (2) and (3). For high values of heat transfer coefficient hsa, the difference in temperatures between the solid and fluid phases are small, this is made evident by comparing the calculated solid and fluid temperatures obtained in the present study, and a single energy equation with one temperature (T = Ts = Ta) can describe the system which takes the form:

@ðqcp TÞ @ qa U d cpa T @ @T þ ¼ keff þS @t @z @z @z

ð13Þ

where qcp = qaeacpa + qsescps, keff = eaka + esks and S is a source term representing the latent heat evolution. Neglecting heat diffusion, Eq. (13) indicates that high heat transfer rates are expected at high gradients of temperature in the flow direction. Variations of temperature and liquid fraction along the bed height at a process time of 6000 s during charging and discharging are shown in Figs. 11 and 12. At this time significant differences in the charging and discharging rates are observed for different values of mixing ratio as can be seen in Figs. 7 and 8. Thus comparisons of temperature variations at this time give an idea about the mechanism of enhancement of heat transfer rates. Fig. 11 shows that during charging, the temperature decrease along the flow direction is higher for GR27 and GR27/ GR41 = 1 beds as compared to GR41 bed. Keeping in mind the liquid fraction variation along the bed height, one can observe that these temperature drops corresponds to the melting temperature range of the bed material. During discharging, higher increase in temperature along the flow direction are obtained using GR41 and GR27/ GR41 = 1 beds as compared to GR27. In general, we notice that high rates of heat transfer are obtained when significant temperature decrease occurs along the flow direction for charging and significant increase for discharging. That is proper mixing of GR27 and GR41 results in a favourable temperature variation along the bed length that ultimately enhances the overall charging and discharging rates.

To conclude the present parametric study and demonstrate the advantages of mixing technique of GPCCs proposed in the present article, one should compare the performance of packed bed composed of a mixture of GR27 and GR41 with GR27/ GR41 = 1 with an equivalent packed bed composed of two successive equal segments of GR27 and GR41. Two additional column configurations are considered for this purpose. In the first configuration, the first half of the column is composed of GR27 and the second half is composed of GR41. Charging and discharging are made from the GR27 side. This configuration is referred to as GR27–GR41 unit. In the second configuration, the arrangement of GR27 and GR41 is reversed. Charging and discharging are made by the GR41 side. This configuration is referred to as GR41–GR27 unit. Figs. 13 and 14 show the charging and discharging dynamics for different configurations of the packed bed at Re = 10. Favourable charging dynamics are obtained for units using GR27, GR41– GR27, and GR27/GR41 = 1. On the other hand, favourable discharging dynamics are obtained for units using GR41, GR27– GR41, and GR27/GR41 = 1. An overall comparison in terms of the normalized total charging and discharging times is shown in Fig. 15. Included also in this comparison the configuration of a unit consisting of GR41–GR27 during charging and GR27– GR41 during discharging. This unit is refereed to as CH GR41– GR27 + DIS GR27–GR41. As compared to GR27 only unit, the overall unit performance is enhanced by 15%, 7%, and 7% using configurations of GR27/GR41 = 1, GR27–GR41, and GR41–GR27, respectively. It should be also noted that using column configurations of GR41–GR27 during charging and GR27–GR41 during discharging results in an enhancement of about 11% as compared to a value of 15% obtained using a mixture of GR27/GR41 = 1. In general using units composed of a mixture of GR27/GR41 = 1 results in an optimized unit performance without limitations on the charging and discharging directions during actual implementations.

M. Rady / Applied Energy 86 (2009) 2704–2720

2715

Fig. 11. Variation of temperature and liquid fraction along flow direction during charging, Re = 10, time = 6000 s.

Fig. 12. Variation of temperature and liquid fraction along flow direction during discharging, Re = 10, time = 6000 s.

6. Exergy analysis Further insight into the analysis of results is obtained by performing an exergy analysis based on the second law of thermodynamics. An exergy analysis is a powerful tool for the evaluation of the performance of latent heat storage systems [27,28]. It takes into account the quality and usefulness of the energy transferred. Exergy analysis is performed from a knowledge of numerically determined time dependent temperatures of the heat transfer fluid

(air) and GPCC obtained using the mathematical model and assumptions explained in Section 4. For a complete description of exergy analysis and references for equations used here, we refer to [27–29]. 6.1. Charging process The exergy destroyed during the charge process, Edc, can be obtained by making an exergy balance on the packed bed storage

2716

M. Rady / Applied Energy 86 (2009) 2704–2720

Fig. 13. Charging dynamics for different configurations of packed bed, Re = 10.

Fig. 14. Discharging dynamics for different configurations of packed bed, Re = 10.

unit, Eq. (14). Exergy destroyed is a positive quantity for actual processes and becomes zero for a reversible process. It represents the lost work potential.

Edc ¼ Ein  Eout  Ec

ð14Þ

The balance in Eq. (14) means that the exergy destroyed (Edc) is the equal to difference between the net exergy supplied (Ein  Eout)

and the charged exergy (Ec). Ein is the inlet exergy and Eout is the outlet exergy. The charged exergy is the sum of exergy charged in the GPCC (EGPCC) and air (Eair), Ec = EGPCC + Eair. They are written as follows:

Ein ¼

Z 0

t

  T _ pa T in  T o  T o ln in dt mc To

ð15Þ

M. Rady / Applied Energy 86 (2009) 2704–2720

2717

.

.

.

.

.

.

Fig. 15. Normalized charging and discharging times for different packed bed configurations.

Eout ¼

Z

t

0

  T _ pa T out  T o  T o ln out dt mc To

ð16Þ

where To is a reference temperature taken as 25 °C. In Eq. (14), it is assumed that exergy gained due to heat leakage to the surroundings are negligible. The exergy charged in the GPCC is written as:

EGPCC

  T GPCC ¼ M GPCC cpl ei T GPCC  T m2;i  T o ln T m2;i i   X T m1;i M GPCC cps ei T m1;i  T init  T o ln þ T init i   X To M GPCC Hls ei 1  þ T mavg;i i X

ð17Þ

where i refers to the GPCC material, i = 1 and 2 refers to GR27 and GR41 in the present study, MGPCC is the total mass of GPCC in the bed, TGPCC is the final GPCC temperature at the end of charging, and Tmavg is the average melting temperature of GPCC. The exergy charged in the air expressed as:

  T air Eair ¼ Mair cpa T air  T init  T o ln T init

ð18Þ

where Tair is the final air temperature in the bed and Tinit is the initial bed temperature. It should be noted that the exergy charged in the GPCC, Eq. (17), consists of contributions during the solid state, liquid state, and latent heat evolution during melting. The exergy efficiency during

Fig. 16. Variation of exergy destroyed with mixing ratio, Re = 10.

2718

M. Rady / Applied Energy 86 (2009) 2704–2720

the charging process is defined as the ratio of charged exergy to the net exergy supplied.

gc ¼

Ec Ein  Eout

ð19Þ

6.2. Discharging process The exergy destroyed during the discharge process, Edd, is equal to the difference between the initial charged exergy (Ec) and the sum of net discharged exergy (Eout  Ein) and remains of charged exergy (Erc). It is written as:

Edd ¼ Ec  ½ðEout  Ein Þ þ Erc 

ð20Þ

The value of Ec is obtained from the analysis of the end of charging process. Eout and Ein are calculated using Eqs. (15) and (16) with the temperature data of discharging process. Erc is the remains of charged exergy at the end of discharging process. The expression for calculation of Erc is similar to this given by Eqs. (17) and (18) but using the PCM state and temperature at

the end of discharging process. It is non-zero when the final temperature of the bed at the end of discharging is higher than the initial bed temperature at the beginning of charging process. In the present study, complete discharging of the bed is performed (TGPCC = Tair = Tin) with the inlet temperature (Tin) equal to the initial temperature during charging (Tinit). Thus the value of Erc is equal to zero. The exergy efficiency during the discharge process is defined as the ratio of the net discharged exergy to the initial charged exergy:

gd ¼

Eout  Ein Ec

ð21Þ

The exergy efficiency of a complete cycle consisting of a complete charging and discharging processes can be obtained by the product of exergy efficiency during charging and discharging.

gcd ¼ gc gd

ð22Þ

It represents the ratio of net discharged exergy during the discharging process to the net charged exergy during the charging process.

Fig. 17. Variation of exergy efficiency with mixing ratio, Re = 10.

Fig. 18. Variation of exergy destroyed and exergy efficiency for different bed configurations.

M. Rady / Applied Energy 86 (2009) 2704–2720

6.3. Results of exergy analysis The variations of exergy destroyed during charging and discharging at different values of GR27/GR41 mixing ratio are shown in Fig. 16. The exergy destroyed initially decreases with the increase of GR27/GR41 mixing ratio until reaching a minimum value between 0.8 and 1.0 and then increases. This corresponds to an optimized mixing ratio of the heat storage unit. It is interesting to note that this optimum point corresponds to the value obtained from the analysis of total charging and discharging times shown in Fig. 9. The lowest amount of exergy destroyed corresponds to the highest exergy efficiency as shown in Fig. 17. Comparisons of exergy destroyed and exergy efficiency for different bed configurations are shown in Fig. 18. Again, the results obtained using exergy analysis confirms the results obtained based on the optimization of total charging and discharging times. That is remarkable decrease in exergy destroyed and improvement in exergy efficiency using mixing ratio of GR27/GR41 = 0.8–1.0. This supports the idea of introducing this concept as a simple method for performance enhancement of LTES systems. Further studies, experimental and numerical, are planned in the near future to completely elaborate this potential.

7. Conclusions Utilization of multiple GPCCs with different ranges of PCT in a packed bed storage system represents the main interest of the present work. Small particle diameter of GPCC allows simple mixing of two or three ranges of GPCCs in a packed bed for enhancement of storage unit performance. Experiments have been carried out to characterize the phase changing characteristics of two GPCCs (Rubitherm GR27 and GR41) chosen to demonstrate the feasibility of this technique for low temperature LHTES applications. Packed bed column experiments have been carried out using different mixing ratios of GR27 and GR41 and the results are compared with the cases of using simple GR27 and GR41 packed beds. A mathematical model for the packed bed column consisting of separate energy equations for the solid GPCC particles and the heat transfer fluid has been developed. The model accounts for phase change in the bed including multiple GPCCs. Validation of the numerical model has been carried out by comparison with experimental results of packed bed column at different values of mixing ratio. The model has been used to perform a parametric study to investigate the overall bed performance at different values of mixing ratio and Reynolds number for constant inlet charging and discharging temperatures. The unit performance is evaluated in terms of the total charging and discharging times. It has been demonstrated that, as compared to the use of single GPCC, careful choice of the mixing ratio of GPCCs in a composite bed can result in a significant enhancement of the overall storage unit performance. The optimum mixing ratio is independent of the value of Reynolds number. As compared to GR27 only unit, the overall unit performance is enhanced by about 15% using a mixing ratio of GR27/ GR41 between 0.5 and 1.0. The optimum configuration of packed bed unit consisting of a mixture of GR27 and GR41 is also compared with configurations using two successive equal segments of GR27 and GR41 arranged along the bed length. As compared to GR27 only unit, the overall unit performance is enhanced by 7% using equal sequential layers of GR27 and GR41 in the unit. It should be also noted that using column configurations of GR41– GR27 during charging and GR27–GR41 during discharging results in an enhancement of about 11% as compared to a value of 15% obtained using a mixture of GR27/GR41 = 1. An exergy analysis confirms the above results and shows remarkable decrease in exergy

2719

destroyed and improvement in exergy efficiency using mixing ratio of GR27/GR41 = 0.8–1.0. In general, using units composed of a mixture of GPCCs with an optimized mixing ratio results in a remarkable improvement of the unit performance without limitations on the charging and discharging directions during practical applications. The present analysis supports the present idea of using the mixing concept of multiple GPCCs as a simple method for performance enhancement of LTES systems.

References [1] Farid MM, Kanzawa A. Thermal performance of a heat storage module using PCM’s with different melting temperatures: mathematical modeling. Trans ASME: J Sol Energy Eng 1989;111:152–7. [2] Farid MM, Kanzawa A. Thermal performance of a heat storage module using PCM’s with different melting temperatures: experimental. Trans ASME: J Sol Energy Eng 1990;112:125–31. [3] Gong ZX, Mujumdar AS. Enhancement of energy charge–discharge rates in composite slabs of different phase change materials. Int J Heat Mass Trans 1996;39:725–33. [4] Gong ZX, Mujumdar AS. Cyclic heat transfer in a novel storage unit of multiple phase change materials. Appl Therm Eng 1996;16:807–15. [5] Gong ZX, Mujumdar AS. A new solar receiver thermal store for space-based activities using multiple composite phase change materials. ASME J Solar Eng 1995;117:215–20. [6] Haiting C, Xiugan Y, Xinbin H. Thermal performance analysis for a heat receiver using multiple phase change materials. Appl Therm Eng 2003;23:2353–61. [7] Fang M, Chen G. Effects of different multiple PCMs on the performance of a latent thermal energy storage system. Appl Therm Eng 2007;27:994– 1000. [8] Gong ZX, Mujumdar AS. Cyclical phase change heat conduction in thin composite slabs. In: Wrobel LC, Brebbia CA, editors. Computational modeling of free and moving boundary problems, vol. 2. Computational Mechanics Publications; 1991. [9] Gong ZX, Mujumdar AS. Cyclical heat conduction in melting and freezing of composite slabs of different phase change materials. In: Proc. 10th Int. Heat Transfer Conf., vol. 6; 1994. p. 343–8. [10] Gong ZX, Mujumdar AS. Enhancement of energy charge–discharge rates in composite slabs of different phase change materials. Int J Heat Mass Trans 1996;39:725–33. [11] Shaikh S, Lafdi K. Effect of multiple phase change materials (PCMs) slab configurations on thermal energy storage. Energy Convers Manage 2006;47:2103–17. [12] Adebiyi GA, Hodge BK, Steele WG, Jaladeh A, Nsofor EC. Computer simulation of a high temperature thermal storage system employing multiple family of phase change storage materials. J Energy Resour Technol 1996;118:102–11. [13] Wang J, Chen G, Jiang H. Theoretical study on a novel phase change process. Int J Energy Res 1999;23:287–94. [14] Wang J, Chen G, Zheng F. Study on phase change temperature distributions of composite PCMs in thermal energy storage systems. Int. J. Energy Res 1999;23:277–85. [15] Wang J, Ouyang Y, Chen G. Experimental study on charging processes of a cylindrical heat storage capsule employing multiple-phase-change materials. Int J Energy Res 2001;25:439–47. [16] Zhang D, Zhou J, Wu K, Li Z. Granular phase changing composites for thermal energy storage. Solar Energy 2005;78:471–80. [17] Rady MA, Huzayyin AS, Arquis E, Monneyron P, Le Bot C. Study of heat and mass transfer in a dehumidifying desiccant bed with macro-encapsulated phase change materials. Renew Energy 2009;34:718–26. [18] Nagano K, Takeda S, Mochida T, Shimakura K, Nakamura T. Study of a floor supply air conditioning system using granular phase change material to augment building mass thermal storage-heat response in small scale experiments. Energy Build 2006;38:436–46. [19] Rady M, Arquis E. A comparative study of phase changing characteristics of granular phase changing materials using DSC and T-history methods. In: Proceedings of the fourth international conference on thermal engineering: theory and applications, Abu Dhabi (UAE); January 12–14, 2009. [20] Rady M. Study of phase changing characteristics of granular composites using differential scanning calorimetry. Energy Convers Manage; 2009. doi:10.1016/ j.enconman.2009.01.30. [21] Galloway TR, Sage BH. A model of the mechanism of transport in packed distended and fluidised beds. Chim Eng Sci 1970;25:495. [22] Beasley DE, Clark JA. Transient response of a packed for thermal energy storage bed. Int J Heat Mass Trans 1984;27:1659–69. [23] Kuwahara F, Nakayama A, Koyama H. A numerical study of thermal dispersion in porous media. J Heat Trans 1996;118:756–61. [24] Amiri A, Vafai K. Transient analysis of incompressible flow through a packed bed. Int J Heat Mass Trans 1998;41:4257–79.

2720

M. Rady / Applied Energy 86 (2009) 2704–2720

[25] Theuerkauf J, Witt P, Schwesig D. Analysis of particle porosity distribution in fixed beds using the discrete element method. Powder Technol 2006;165:92–9. [26] de Klerk A. Voidage variation in packed beds at small column to particle diameter ratio. AIChE J 2003;49:2022–9. [27] Bejan A. Second-law analysis in heat transfer and thermal design. Adv Heat Trans 1982;5:1–58.

[28] Wantanabe T, Kanzawa A. Second law optimization of a latent heat storage system with PCMs having different melting points. Heat Recov Syst CHP 1995;15:641–53. [29] Moran MJ, Shapiro HN. Fundamentals of engineering thermodynamics. 5th ed. Hoboken (MJ, USA): John Wiley and Sons, Inc.; 2004.