Next Article in Journal
Methanol Marker for the Detection of Insulating Paper Degradation in Transformer Insulating Oil
Previous Article in Journal
Colorimetric Method for Detection of Hydrazine Decomposition in Chemical Decontamination Process
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development of a Model for Performance Analysis of a Honeycomb Thermal Energy Storage for Solar Power Microturbine Applications

1
Department of Mechanical Engineering & Aeronautics, School of Mathematics, Computer Science & Engineering, University of London, London EC1V 0HB, UK
2
Institute for Thermal Power Engineering, College of Energy Engineering, Zhejiang University, 38 Zheda Road, Hangzhou 310027, China
*
Author to whom correspondence should be addressed.
Energies 2019, 12(20), 3968; https://doi.org/10.3390/en12203968
Submission received: 31 July 2019 / Revised: 29 September 2019 / Accepted: 8 October 2019 / Published: 18 October 2019
(This article belongs to the Section D: Energy Storage and Application)

Abstract

:
Solar power microturbines are required to produce steady power despite the fluctuating solar radiation, with concerns on the dispatchability of such plants where thermal energy storage may offer a solution to address the issue. This paper presents a mathematical model for performance prediction of a honeycomb sensible-heat thermal energy storage designed for application of concentrated solar power microturbine. The focus in the model is to consider the laminar developing boundary layers at the entry of the flow channels, which could have a profound effect on the heat-transfer coefficient due to large velocity and temperature gradients, an effect which has not been considered in the modelling of such storage systems. Analysing the thermal and hydrodynamic boundary layer development, the Nusselt number and the friction factor were evaluated using a validated conjugate heat-transfer method. The simulations results were used to develop accurate regression functions for Nusselt number and friction factor. These formulations have been adopted within a one-dimensional model to evaluate the performance of the storage under different operating conditions. The model was in good agreement with conjugate heat transfer results with maximum relative error below 2%. Two case studies are presented to demonstrate the applicability of the proposed methodology.

1. Introduction

The fluctuating and uncertain behaviour of renewable energy sources has led to concerns regarding the reliability of electric grids which incorporates them. This challenge has created an increased interest in the development of energy storage systems. However, this “need” for energy storage to improve renewable integration is an economic question, and the benefits from an integrated energy storage system may not justify the investment [1,2]. In comparison to other intermittent forms of renewable energy such as wind or photovoltaics, concentrated solar power (CSP) technologies have the potential for integration with cost-effective thermal energy storage to increase the dispatchability of the plant [3]. This is particularly important for solar dish CSP plants with relatively small size (1–25 kWe) characterised by a low thermal inertia when compared with large-scale technologies such as parabolic through or solar tower [4].
Solar dish CSP power plants, especially when operated by micro gas turbines, are also characterised by a relatively high maximum temperature of the working fluid [5,6]. This characteristic imposes limitations on the selection of the storage system and media. Moreover, one of the most promising applications of these plants is the possibility to work as a standalone system in rural areas. The system then needs to be compact and easy to transport, making the compactness of the thermal energy storage (TES) unit a major requirement [7].
Ceramic honeycomb structures have been widely used as heat-transfer media within sensible energy storage systems [8]. Their high heat-transfer surface per unit volume, large heat capacity, and good thermal shock resistance can ensure compactness, high performances, and high temperatures, major requirements for solar dish micro gas turbines, especially when operated in rural areas [9].
Wang et al. [9] performed an experimental investigation of a packed-bed honeycomb ceramic for a high-temperature thermal storage system to be operated by concentrated solar power. The study analysed the heat-transfer characteristics of the system and demonstrated that the honeycomb ceramic material can be used as the thermal storage material for high-temperature thermal storage systems. Andreozzi [1] investigated numerically a high-temperature thermal storage in a honeycomb solid matrix. The author performed a parametric analysis considering different porosity values and concluded that an optimal porosity value must be evaluated to identify the right trade-off between heat stored and pressure losses.
Modelling the performance of a sensible heat thermal energy storage system involves the solution of a system of two differential equations, one for the solid media and one for the working fluid. Luo et al. [10] proposed a one-dimensional model for a honeycomb ceramic thermal storage and performed a parametric analysis in order to evaluate the influence of the main geometrical parameters. The model was solved adopting the Euler implicit scheme and was validated with experimental results. Results of the validation demonstrated an outlet temperature difference between simulations and experimental results of about 5–20%. Cascetta et al. [11] experimentally and numerically investigated thermal behaviour of a packed-bed thermal storage system filled with solid alumina beads. Their one-dimensional model calculated fluid and solid temperature separately based on two LTNE (Local Thermal Nonequilibrium) equations, which predicted well the axial temperature distribution and could give exact energy stored in charging or discharging processes.
Given the small size of the plant, a TES system integrated within a solar dish CSP plant is usually characterised by relatively low mass-flow, rate which inevitably leads to low Reynolds numbers and, consequently, laminar flow [10]. For laminar flows, in the first portion of the duct, the development of the thermal and hydrodynamic boundary layers has a crucial impact on the heat transfer and pressure drop. The heat-transfer coefficient as well as the friction coefficient may vary significantly from the fully developed flow [12].
The aim of this paper is to develop a methodology to improve the accuracy of a one-dimensional model of a honeycomb TES characterising the heat transfer and pressure loss in order to understand the effect of a laminar developing flow on the charging and discharging processes. Both thermal and hydrodynamic development of the flow boundary layer will be analysed using conjugate heat-transfer simulations, with the objective of evaluating the heat-transfer coefficient and the friction factor with dissimilar inlet velocity and duct diameter. The results will then be used to develop accurate regression models for both heat-transfer and pressure-drop coefficients. These functions will be then integrated in the one-dimensional model to improve its prediction accuracy.

2. Materials and Methods

The boundary layer development problem was addressed for the first time by Graetz [13]. Graetz studied the development of the boundary layer and, under the hypothesis of uniform wall temperature, he derived Equation (2) to evaluate the Nusselt number in the entry region of a circular duct. The equation was obtained using an analytical solution of the governing Equation (1).
2 T f r 2 + 1 r T f r = u α T f x 2 T f x 2
where r and x denote the radial and the axial directions, respectively; u is the fluid velocity; and α is the thermal diffusivity of the fluid. The inlet parabolic fully developed velocity profile as well as negligible conduction in the axial direction can be assumed. Under these assumptions, Equation (1) can be transformed in its nondimensional form and solved by the method of separation of variables, resulting in Equation (2).
N u ( x * ) = n = 1 G n e 2 λ n 2 x * 2 n = 1 G n λ n 2 e 2 λ n 2 x *
x * = x D R e P r
where x * is defined by Equation (3) and where the constants G n and eigenvalues λ n for this infinite series are given in Table 1.
A similar equation can be derived considering constant heat flux as a boundary condition. Moreover, it is possible to derive a solution for any arbitrary surface temperature distribution by simply discretising the surface temperature into a number of constant-temperature steps and summing, or superposing, the constant-surface temperature and thermal-entry-length solutions for each step [14].
An alternative to the Graetz series solution is the relationship proposed by Léveque [15], reported in Equation (4).
N u = { 1.077   x * 1 / 3 0.7              x *   0.01   3.66 + 6.9 ( 10 3 x * ) 0.488 e ( 57.2 x * )     x * >   0.01
Similar expressions can be derived considering uniform wall heat flux. The Léveque solution for the thermal boundary layer development problem with constant heat flux is reported in Equation (4).
N u = 1.302   x * 1 / 3 1                       x *   5 × 10 5 1.302   x * 1 / 3 0.5            5 × 10 5 <   x *   10 3   4.37 + 8.7 10 3 x * 0.488 e 41 x *              x * >   10 3
The most realistic approach to analysing the laminar flow development problem considers also the development of the velocity profile along the duct. The hydrodynamic development of the fluid can be considered in Equation (1), introducing non-negligible complications in its solution. Churchil and Ozoe [16] derived the relationship in Equation (6) to evaluate Nusselt number simultaneously for the thermal and hydrodynamic flow development with constant temperature as a boundary condition.
N u = ( 0.6366 ( 4 x * π ) 1 2 [ 1 + ( P r 0.0468 ) 2 3   ] 1 4 )
Churcill and Ozoe [16] also derived Equation (7) for the Nusselt number considering constant heat flux as a boundary condition.
N u 4.364 [ 1 + ( π / 4 29.6 x * ) 2 ] 1 / 6   = [ 1 + ( π / 4 19.04 x * [ 1 + ( P r 0.0207 ) 2 3   ] 1 2 [ 1 + ( π / 4 29.6 x * ) 2 ] 1 3 ) 3 2 ] 1 / 3
Kristiawan et al. [17] numerically studied heat-transfer enhancement of nanofluids and verified the accuracy and reliability of the laminar flow model at entry by using the Shah–London equation in Equation (8) with pure water as the heat-transfer fluid.
N u   = { 4.364 + 0.0722 R e P r D x          ( R e P r D x 33.3 ) 1.953 · ( R e P r D x ) 1 3                   ( R e P r D x 33.3 )
The forced convection within a sensible heat thermal energy storage is a conjugate heat-transfer problem and results in neither uniform wall temperature nor uniform heat flux at the boundaries. The problem in this case must be formulated for the entire solid–fluid system rather than using a priori boundary conditions, and Equation (9) must be solved simultaneously with Equation (1).
2 T s r 2 + 1 r T s r + 2 T s x 2 = 0
According to Shah and London [18], the local Nusselt number for the conjugate problem falls in between the constant wall temperature boundary condition problem and the constant heat flux boundary condition problem of the conventional convection problem. In this situation, the Nusselt number is a function not only of the dimensionless length x * but also of the ratio between the wall thickness δ and the duct diameter and of the parameter R W defined by Equation (10). The contribution of R W can be neglected whenever its value falls below 0.05 [18]. In an air-fed honeycomb TES, the wall thickness and the thermal conductivity of the fluid are usually in the order of magnitudes of 10 2 to 10 3 . Consequently, a value of the R W parameter lower than 0.05 can be expected. Thus, the contribution of R W can be neglected.
R W = λ f δ λ s D
In Equation (10), k f and k s are the thermal conductivity coefficients of the fluid and the solid, respectively. In order to achieve a relatively high efficiency, micro gas turbines (MGTs) are usually operated by a recuperated Brayton–Joule cycle. This cycle is characterised by a relatively low-pressure ratio. Low pressure losses along the thermal energy storage ducts are then a major requirement to improve cycle efficiency. The pressure drop in a developing flow is due to both wall shear forces and fluid acceleration. According to Kays [14], the pressure drop for the hydrodynamic entrance region is not easy to evaluate. This can be either calculated by introducing an apparent Fanning friction factor or by considering the pressure drop to consist of two components: the loss due to shear stresses based on the fully developed flow and additional pressure drop as a result of momentum change and of kinetic energy loss. The second contribution to the pressure drop, in a dimensionless form, is designated as the incremental pressure drop number K(x) (Equation (11)).
Δ p ρ u m 2 / 2 g =   f a p p ( x d ) = f f d ( x d ) + K ( x )
where K ( x ) , the incremental pressure drop, increases monotonically from 0 at x = 0 to a constant value in the fully developed region. This coefficient is a function of the Reynolds number and the duct geometry, which can be calculated by solving the complete set of Navier–Stokes equations. In Equation (11), Δ p is the pressure loss, ρ is the density, u m is the average fluid velocity, x is the axial distance, d is the pipe diameter, and f f d is the friction factor for the fully developed flow. In general, this coefficient can be approximated as 16 / R e [19].
The investigated geometry is shown in Figure 1. Air flows through the honeycomb structure composed of N T circular ducts heating the surrounding material. As shown by the figure, each duct is characterised by a diameter d and a length L. Assuming that each tube exchanges heat with the surrounding storage material in an equivalent control area, it is possible to define the equivalent diameter using Equation (12).
D e q = 4 a 2 π N T
where a is the single storage element side length and NT is the number of elements in the entire storage unit.
To analyse heat transfer between the fluid and solid materials, the one-dimensional model adopted by Luo [10] was used. The model is based on the following assumption:
  • Zero conduction of the material along y and z axis;
  • Uniform initial temperature distribution of solid material ( T S ).
  • Adiabatic external wall of the thermal storage unit.
  • The influence of fluid viscosity on its temperature ( T f ) is neglected.
  • Uniform flow and temperature along y and z axis. The effect of hydrodynamic and thermal development of the boundary layer is taken into account only from the heat-transfer point of view through the heat-transfer coefficient.
The model considers the flow as one dimensional and uses the following governing equations. For the solid, one has the following:
ρ S C S T S t λ S 2 T S x 2 = h A ( T s T f )
For the fluid, the governing equation is as follows:
ρ f C P , f ( T f t u T f x ) = h A ( T s T f )
In the previous equations, ρ S and ρ f are the density of the solid material and the fluid, respectively; C S is the solid medium specific heat; C P , f is the specific heat at constant pressure of the working fluid; h is the heat-transfer coefficient; and A is the heat-transfer surface. In the study, conducted by Luo, the heat-transfer coefficient was evaluated using Equation (15), an empirical relationship available in the literature for porous media [20].
N u = h d λ f = 3.66 + 0.668 R e P r d x 1 + 0.04 ( R e P r d x ) 2 3
In Equation (15) λ f is the fluid thermal conductivity. This equation demonstrates low variability of the heat-transfer coefficient with the axial distance from the inlet not taking into account the effects of the laminar boundary layer development accurately. The considered material is a ceramic with the properties reported in Table 2.
The MGT operates with an open Brayton–Joule cycle; consequently, the heat-transfer fluid is air. The properties of air were evaluated as a function of temperature. The variation of air properties with the pressure were neglected.

Regression Model

To solve the abovementioned 1-D model, knowledge of the heat-transfer and pressure-loss coefficients is required. These coefficients are strongly influenced by the laminar boundary layer development, the operating conditions (i.e., mass flow rate), and the TES geometry. For these reasons, a function to take into account all these aspects is required.
The heat-transfer and pressure-loss coefficients were evaluated using a regression model obtained from the results of conjugate heat-transfer analysis completed using a commercial software (ANSYS CFX TM). According to Bejan [13], in general, for the entry region, the Nusselt number can be modelled using the following relationship:
N u = k 1 ( R e P r x / d ) k 2
where k 1 and k 2 are constants which must be evaluated by the regression model in order to fit computational or experimental results. The value of k 2 is usually around 0.5 [13] and is assumed as such, and k 1 was chosen to best fit the computational results.
A similar approach was adopted for the friction coefficient. As already mentioned, the friction coefficient can be evaluated by separating the contribution of the shear stresses from the kinetic energy loss related to the hydrodynamic development of the fluid. The regression model adopted in this study, reported in Equation (17), is based on the model described by Kays [14] in Equation (11). The model was modified in this study by inserting two exponents related to the Reynolds number and the ratio x / d to better fit the results.
f = a 1 R e a 2   + a 3 ( x / d ) a 4
In Equation (17), a 1 ,   a 2 ,   a 3 and a 4 are constants which will be evaluated by the regression model to fit the computational results. The first term represents the pressure loss due to shear stresses. It is inversely proportional to the Reynolds number and is prevailing for x . The second term represents the contribution of the hydrodynamic development and is prevailing for x 0 .
The numerical solution is obtained over the solid and fluid domains. Mapped mesh with sweeping method was adopted and includes an inflation in the near wall to resolve the heat transfer more accurately, as shown in Figure 2. Grid independence study was conducted to find a suitable grid density for the calculations. The computational domain was discretized using five different grid arrangements, as reported in Table 3. The heat-transfer coefficient was selected for comparison between the different grids. Results from conjugate heat-transfer analysis with different grid arrangements demonstrate acceptable variability of 0.3% passing from 4.5 × 10 6 to 5.8 × 10 6 nodes. Therefore, the grid arrangement of 4.5 × 10 6 for the computational domain has satisfactory grid independence and is adequate for the resolution of the conjugate heat-transfer problem.
The heat-transfer coefficient and pressure drop have been evaluated using Equations (18) and (19).
h = q w ( T ¯ f T w )
f = τ w 1 2 ρ ¯ u ¯ 2
where q w is the wall heat flux, T w is the wall temperature, and τ w = μ   u / y is the shear stress. ρ ¯ , u ¯ , and T ¯ f have been calculated as the mass-flow averaged value.
To estimate the accuracy of the conjugate heat-transfer model in the evaluation of the heat-transfer coefficient and friction factor, a comparison with the results obtained using the commercial software TEXSTANTM was performed. The model was validated for the case of uniform wall temperature. The results were also compared with the Graetz series in Equation (2) and with the relationship developed by Churchill and Ozoe [16] in Equation (6). As already mentioned, the Graetz series assumes a fully developed velocity profile at the duct inlet. The purpose of this comparison is then to evaluate the influence of the hydrodynamic flow development.
A validation case with a circular duct with constant wall temperature of 310 K and diameter of 0.07 m was considered. Air enters the duct with a velocity of 0.225 m/s uniformly distributed along the radius and a temperature of 300 K.
Figure 3 shows the results of the validation. Results show that the maximum relative error of the present model for the heat-transfer coefficient is always below 4% and an absolute value of the error in the order of magnitude of 10 3 , significantly lower when compared with the Graetz series and the Churchill and Ozoe [16] relationship, which presents absolute values of the error of 1 and 0.1, respectively. For what concerns the friction coefficient, despite the maximum relative error in the first part of the duct being relatively high, the error decreases significantly along the duct with an average absolute value around 10 3 . For this reason, the accuracy of the model was considered adequate for the purpose of this study.
The discrepancy in the evaluation of the friction factor between the TEXSTAN code and the present model can be justified by the fact that TEXSTAN assumes a quadratic velocity profile, as shown in Figure 4, and does not consider the influence of the velocity overshoot. The effects of axial diffusion of momentum and radial pressure gradient cause the velocity profiles to have peculiar behaviour in the entrance region. Due to this phenomenon, the velocity profiles have a local minimum at the tube centreline and symmetrically located maxima on either side near the walls. Both experimental and analytical results support the existence of overshoots in the velocity profile [18]. Near the entrance where the fluid particles first meet the wall, viscous friction rapidly decelerates the flow to zero velocity at the wall. The high velocity gradient at the wall results in higher shear stresses and pressure gradients.
As demonstrated in Figure 4 the present model takes into account this effect which causes major differences in the development of the velocity profile and, consequently, in the evaluation of the friction factor. Figure 4 shows the nondimensional radius as a function of the nondimensional velocity. rs is the radius of the duct, and vm is the mean velocity.

3. Results and Discussion

The storage unit is designed for integration with a parabolic dish concentrated solar power plant operated by a micro gas turbine. During high solar irradiation periods, the concentrator is designed to allow for additional thermal power for storage over the required power demand. In these periods, it is possible to store this energy and to reclaim it during low sun irradiation time frames. A plant integrating a thermal storage unit can be conceptualised as shown by Figure 5 as one of the possible arrangements. The thermal storage unit is located between the turbine and the receiver. During the charging phase, the inlet temperature of the fluid is the receiver outlet temperature. In the proposed design, due to receiver material limitations, its maximum value is around 1073 K [5]. For this reason, during the charging phase, the inlet temperature of the fluid is set constant to 1073 K and the solid material on the other side is at ambient temperature. During discharging phase, the heat stored by the TES can be reused during low solar irradiance periods. During this, the inlet temperature of the heat-transfer fluid may vary significantly depending on the Direct Normal Irradiation (DNI) and the state of charge of the thermal storage unit. For this reason, the more demanding operating condition was considered. In particular, during discharging phase, the inlet temperature of the fluid is set constant to 420 K, corresponding to the compressor exit temperature at design conditions.

3.1. Conjugate Heat Transfer

This section presents results from conjugate heat-transfer simulations with the aim of studying the influence of boundary layer development on the velocity, temperature, and pressure. Moreover, different duct diameters and mass flow rates were considered and their influence on the flow development was analysed.
Charging and discharging the thermal storage unit are inherently transient problems. Thus, transient simulations considering different diameters of the pipe and inlet velocity of the fluid were carried out. The computational domain is comprised from the solid domain and the fluid domain. Boundary conditions at the fluid–solid interface are exchanged within the solver using a conservative flux interface model. The wall boundary condition for the momentum equation was set to “no slip”.
Figure 6, Figure 7 and Figure 8 show the results of the conjugate heat-transfer analyses. Several simulations were carried out considering dissimilar diameters and mass flow rates. Figure 6 shows the velocity profile at different radii and channel lengths normalised by the inlet velocity. The results highlight that increasing the inlet velocity augments the effect of the velocity overshoot shifting the maximum toward the wall. The same effect can be observed for increasing the duct diameter. Moreover, increasing the duct diameter delays the development of the velocity profile along the duct. The same effect can be observed when increasing the inlet velocity.
Figure 7 shows the temperature profiles at different radii and channel lengths normalised by the inlet temperature. It can be observed that the thermal boundary layer development involves a larger part of the duct when the inlet velocity is increased. The same effect can be observed with a higher duct diameter. Furthermore, increasing the duct diameter delays the thermal development of the flow along the duct; an effect which can be achieved also by increasing the inlet velocity.
Although not shown in the figures, the influence of the equivalent diameter was also analysed, demonstrating negligible influence on the thermal boundary layer development. The equivalent diameter also demonstrated small variability of the heat-transfer coefficient. Varying its value from 0.02 m to 0.05 m, a typical range for this application, demonstrates a maximum variability always below 1% and an average variability of 0.1%. Thus, this parameter will not be considered in the following analysis.
Figure 8 shows the pressure gradient profile at different radii and channel length normalised by the inlet pressure. It is possible to distinguish between two different contributions; the first contribution is due to fluid velocity profile development and is relevant at the inlet of the duct where the effect of the velocity overshoot can be observed. In this section of the duct, the pressure gradient is significantly higher and is characterised by a maximum caused by the velocity overshoot. Increasing the inlet velocity and the duct diameter augments the pressure gradient caused by these effects. From this point of view, the inlet velocity has a major influence by significantly increasing the absolute value of the pressure gradient and the absolute value of the peak caused by the velocity overshoot.
The second effect is due to wall shear which occurs at the wall and is prevalent for x . In this section of the duct, the pressure gradient is significantly lower and constant along the radius. When the velocity is increased, the absolute value of the pressure gradient grows but much less if compared with the first section of the duct. For a larger diameter, a higher value of the pressure gradient is observed which, in this case, is due to the fact that the velocity profile is not fully developed as shown in Figure 6.

3.2. Regression Results

Results from the conjugate heat-transfer analysis described in Section 3.1 were used to create a relationship able to describe the behaviour of the heat-transfer coefficient and friction factor along the duct under different mass flow rates and duct diameters. The relationships obtained from the regression analysis are reported in Equations (21) and (22). The method of the least squares was adopted to calculate the coefficients for the model f ( x ) so the sum of squares is the smallest possible and Equation (17) is satisfied. The analysis starts with the guess and moves towards the minimum of the sum of squares.
min i = 1 n ( y i f ( x i ) ) 2
where i is the number of known values y i obtained from the conjugate heat-transfer analysis. The equations obtained from the regressions were compared with results obtained from the conjugate heat-transfer analysis. Figure 9 shows a comparison between the results obtained from the conjugate heat-transfer analysis for Nusselt number and friction coefficient and the results obtained from the relationships developed in this study. The regression analysis demonstrates good accuracy. Nusselt number regression shows a maximum relative error around 9%, an average relative error below 3%, and a coefficient of determination R 2 always above 0.998 in a wide range of diameters and Reynolds numbers. The maximum relative error calculated for the friction factor is around 25%, while the average relative error is below 9% and the correlation coefficient is always above 0.995. Despite a higher relative error, its maximum absolute value is relatively small in the order of magnitude of 10 3 . The average value instead falls below the order of magnitude of 10 4 .
N u = 0.41 ( P r   R e x d ) 0.5 + 2.25
f = 22.3 R e 1.2 + 0.025 ( x d ) 0.64
The results show that, for both Nusselt number and friction factor, the maximum error is very close to the entrance of the duct. This error decreases significantly along the duct and approaches the order of magnitude of 10 4 and a relative value of 0.5% for each trend line in the fully developed region. This is especially true for the Nusselt number. The friction factor error, on the other hand, seems to be dependent on the Reynolds number variation and its accuracy decreases slightly when the Reynolds number is increased. Nevertheless, the accuracy of the models was considered adequate for the purpose of this study.

3.3. One-Dimensional Model Validation

The equations developed so far were integrated within the one-dimensional model described in Section 2. In order to demonstrate the accuracy of the model, the charging phase of a single channel characterised by an internal diameter of 0.01 m, equivalent diameter of 0.015 m, and length of 0.2 m was analysed. Results from the model were compared with the already validated conjugate heat-transfer model in five different sections of the duct. These results were compared for five different Nusselt relationships: the one used by Luo [10] for porous media, the ones developed in this study and the Graetz series [13], and the two relationship derived by Churchill and Ozoe [16] for constant temperature boundary condition (T) and constant heat flux (H).
The results of the validation, shown in Figure 10, demonstrate good accuracy of the present model and a maximum relative error below 2%. On the other hand, the Nusselt relationship used by Luo and the Graetz series demonstrate a higher relative error, which can reach 10% for the Graetz series and almost 20% for the Nusselt relationship adopted by Luo. As expected, this discrepancy decreases along the duct as the influence of the developing boundary layer reduces and the Nusselt number gets closer to the fully developed value.
The relationship developed by Churchill and Ozoe [16] for constant heat flux boundary condition (H) demonstrates a maximum relative error of about 7%, showing good agreement when compared with the Graetz series and the relationship adopted by Luo. Nevertheless, its accuracy is lower than the relationship presented in this study. The second equation developed by Churchill and Ozoe [16], on the other hand, shows high relative error, demonstrating that it is not suitable for this application.
These results are generally in line with what was stated in Section 1 and by Shah [18]. The solution of the laminar boundary layer for a conjugate problem falls in between the constant wall heat flux boundary condition and the constant wall temperature boundary condition.
As stated in Section 1, the considered Graetz series reflects the behaviour of thermally developing flow with a fully developed velocity profile. Consequently, it can be concluded that the hydrodynamic development has a major influence on the heat transfer within a honeycomb TES system, leading to an error which approaches 10% in the first portion of the duct. It can be concluded that the relation proposed by Luo is not suitable for this application. This relationship leads to an error which approaches 20% in the first portion of the duct, and consequently, a honeycomb TES should not be modelled as a porous media.
All models demonstrate good accuracy in the prediction of the outlet temperature. Nevertheless, the present model is still the more accurate with a maximum relative error of about 3% lower when compared with the 4% and 4.5% of the Graetz series and the Luo model, respectively, and the 3.5% characterising Equation (7).

3.4. Case Study

To illustrate the functionality of the model, two case studies are presented. In the first case, a one-hour storage system has been studied. In order to achieve the target charging/discharging time with the given material and heat-transfer fluid, the total heat capacity of the storage media must be higher than or equal to the energy transferred by the working fluid for the target time and the following relationship must be satisfied [21].
V S C P , f   m ˙   Δ t C S   ρ S
where V S is the storage media volume; m ˙ is the mass flow rate; C P , f and C S are the specific heats of the fluid and solid, respectively; ρ S is the storage media density; and Δ t is the target charging/discharging time. The overall volume of the system can be calculated using Equation (24).
V S = N T π ( D e q 2 d 2 ) L 4
Results from the regression analysis from Section 3.2 demonstrate minimal pressure losses with the maximum duct diameter and the maximum heat-transfer coefficient with the highest Reynolds number possible. For this reason, the values d = 0.02 m and Re = 1500 were selected as they are the extreme values considered in this study. Adopting these values and considering an equivalent diameter of 0.025 m, the sizing of the storage leads to an overall number of ducts of about 160 and an overall length of the storage of about 0.2 m.
As mentioned in Section 3, during the charging phase, the inlet temperature of the fluid is set to a constant value of the receiver outlet temperature (1073 K). The solid material starting temperature instead is set at 300 K. During the discharging phase, the inlet temperature of the heat-transfer fluid may vary significantly depending on the DNI, the state of charge of the thermal storage unit, and the control strategy adopted by the plant. The more demanding operating condition was considered and set constant to 420 K, corresponding to the compressor exit temperature at design conditions.
Figure 11 shows the charging and discharging behaviour of the storage. The solid media achieves the required characteristics, ensuring a discharging time longer than 1 h. On the other hand, the fluid temperature cannot guarantee a temperature high enough to sustain the system using the stored energy only but seems to be suitable for peak shaving and smooth solar fluctuation occurring along the day. Moreover, during the charging phase, the outlet temperature of air is high enough to ensure an appropriate turbine inlet temperature to produce electricity and to charge the storage at the same time. The maximum value of pressure loss is around 0.05% during the charging phase and 0.02% during the discharging phase.
The second case of the performance of a solar power microturbine plant of Figure 5 has been studied for duration of a typical day in Rome. The variation of the DNI and the output power, in the absence of thermal storage, of the plant is presented in Figure 12. The DNI presents a drop during the afternoon, most probably caused by a cloud passage with consequent output power drop, which causes the shutdown of the plant between 12:08 and 2:38. This case study aims to demonstrate the capability of the thermal storage to reduce the effects of solar fluctuations. The storage is operated to maintain constant output power.
During the discharging phase, the inlet temperature of the heat-transfer fluid may vary significantly depending on the DNI, the state of charge of the thermal storage unit, and the control strategy adopted by the plant. Figure 13 presents the performances of the plant when the thermal storage is integrated. The operational strategy applied during the charging phase is at constant turbine inlet temperature (TIT) and during the discharging phase is at constant power output. Specifically, in this situation the operations are characterised by two charging and two discharging phases. The daily variation of the solar irradiation is such that a second charging phase can be included during the day. The details of the cycle calculation and optimisation are given in References [22,23,24].
During the first discharging phase, the plant is operated to maintain the power output constant to the value 8 kWe. Nevertheless, the irradiation drop is such that the plant cannot produce constant power output. As shown by the chart in Figure 13, the storage can provide enough thermal energy to smooth the solar fluctuation and produces at least 7 kWe during the low irradiation period. The second discharging phase starts when the plant power output is 6 kWe. The output is then kept constant or at its maximum value when the solar irradiation and the storage are not able to sustain the heat demand. The overall energy generated by the plant during the day is 61,164 kWhe and showed an increase of 11.7% with respect to the energy generated by the plant operation without thermal storage (54,027 kWhe). The maximum efficiency achieved during the day is around 16%. Moreover, the plant efficiency is generally high during operations and above 15% for more than 7 h.

4. Conclusions

In this paper, a model for a honeycomb thermal energy storage for solar power applications was presented. The storage is intended for integration with a micro gas turbine power cycle, and the optimisation objectives were chosen in order to meet the engine performance requirements. The model is based on the solution of the one-dimensional Navier–Stokes equations in the fluid domain. Conjugate heat-transfer analyses using a commercial software were performed in order to develop an accurate regression model which takes into account the effect of the laminar flow development at the duct entry on the heat transfer and pressure loss along the storage pipes. The conjugate heat-transfer model was validated against data available in the literature. Different Nusselt number relationships were then compared to verify the accuracy of the model. The results demonstrated higher accuracy of the present model when compared with the relationship adopted by Luo [10] and the Graetz series, with a maximum relative error below 2%, significantly more accurate than the Nusselt relationship adopted by Luo and the Graetz series, which demonstrate a relative error of about 10% for the Graetz series and almost 20% for the Nusselt relationship adopted by Luo. As expected, this difference reduces along the duct as the influence of the developing boundary layer reduces and as the Nusselt number gets close to the fully developed flow value. Finally, two case studies have been presented indicating that, during the charging phase, the outlet temperature of the air is high enough to guarantee an appropriate turbine inlet temperature to produce electricity and charge the storage at the same time.

Author Contributions

Conceptualization, D.I. and Q.Z.; methodology, D.I. and Q.Z.; software, X.Z. and D.I.; validation, X.Z. and D.I.; formal analysis, D.I. and X.Z.; investigation, D.I., Q.Z., J.A.Z. and A.S.; resources, A.S. and G.X; data curation, D.I. and J.A.Z.; writing—original draft preparation, D.I.; writing—review and editing, J.A.Z.; visualization, D.I. and J.A.Z.; supervision, Q.Z., A.S. and G.X.; project administration, A.S., J.A.Z. and G.X.; funding acquisition, A.S., J.A.Z. and G.X.

Funding

This research was part of the SolGATS (Concentrated Solar Power Micro Gas Turbine with Thermal Energy Storage, www.city.ac.uk/solgats) project which is funded by Innovate UK and Ministry of Science and Technology of the People’s Republic of China, grant number 102879.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

CHTConjugate Heat Transfer
CSPConcentrated Solar Power
DNIDirect Normal Irradiation
MGTMicro Gas Turbine
TESThermal Energy Storage
TITTurbine Inlet Temperature
LTNELocal Thermal Nonequilibrium
Greek symbols
αFluid thermal diffusivity (m2/s)
λsSolid thermal conductivity (W/mK)
λnEigenvalues
ρsSolid density (kg/m3)
ρfAir density (kg/m3)
wShear stress
Nomenclature
AHeat-transfer area (m2)
aSingle storage element side length (m)
a0, a1, a2, a3, a4Constants
CConstant
CsSolid specific heat (J/kgK)
Cp,fAir specific heat at constant pressure (J/kgK)
dDuct diameter (m)
DeqEquivalent diameter (m)
fFriction factor
fappApparent friction coefficient
f f d Fully developed flow friction coefficient
GnEigenvalue constants
hHeat-transfer coefficient (W/m2K)
K ( x ) Fully developed incremental pressure drop
k1,k2Regression constants
LStorage length (m)
m ˙ Mass flow fate (kg/s)
N T Number of elements
NuNusselt Number
PePeclet Number
PrPrandtl Number
ReReynolds Number
TfAir temperature (K)
TsSolid temperature (K)
TIT Turbine inlet temperature (K)
uAir velocity (m/s)
VsStorage material volume (m3)
xAxial position (m)
rRadial position (m)
x*Nondimensional position
RWGeometry ratio of duct as defined in Equation (9)
kfThermal conductivity of the fluid
ksThermal conductivity of the solid
pPressure drop
DeqEquivalent diameter
TwWall temperature
rsThe radius of the duct
vmThe mean velocity
PinInlet pressure
tTime
TTSThermal storage temperature

References

  1. Andreozzi, A.; Buonomo, B.; Manca, O.; Tamburrino, S. Thermal energy storages analysis for high temperature in air solar systems. Appl. Therm. Eng. 2014, 71, 130–141. [Google Scholar] [CrossRef]
  2. Denholm, P.; Ela, E.; Kirby, B.; Milligan, M. The Role of Energy Storage with Renewable Electricity Generation; Technical Report NREL/TP-6A2-47187; National Renewable Energy Laboratory: Golden, CO, USA, 2010.
  3. Blanco, M.; Santigosa, L.R. Advances in Concentrating Solar Thermal Research and Technology; A Volume in Woodhead Publishing Series in Energy 2017; Woodhead Publishing: Duxford, UK, 2016; pp. 213–246. [Google Scholar]
  4. Schiel, W.; Keck, T. Parabolic dish concentrating solar power (CSP). In Concentrating Solar Power Technology: Principles, Developments and Applications; Lovegrove, K., Stein, W., Eds.; Elsevier: Amsterdam, The Netherlands, 2012; pp. 284–322. [Google Scholar]
  5. Iaria, D.; Al Zaili, J.; Sayma, A. Solar Dish Micro Gas Turbine Technology for Distributed Power Generation. In Sustainable Energy Technology and Policies; Springer: Berlin/Heidelberg, Germany, 2018. [Google Scholar]
  6. Ellingwood, K.; Safdarnejad, S.; Rashid, K.; Powell, K. Leveraging energy storage in a solar-tower and combined cycle hybrid power plant. Energies 2019, 12, 40. [Google Scholar] [CrossRef]
  7. Lovegrove, K.; Stein, W. Concentrating Solar Power Technology; Woodhead Publishing: Cambridge, UK, 2012. [Google Scholar]
  8. Srikanth, O.; Khivsara, S.D.; Aswathi, R.; Madhusoodana, C.D.; Nath Das, R.; Srinivasan, V.; Dutta, P. Numerical and experimental evaluation of ceramic honeycombs for thermal energy storage. Trans. Indian Ceram. Soc. 2017, 76, 102–107. [Google Scholar] [CrossRef]
  9. Wang, Y.; Bai, F.W.; Wang, Z.F.; Kiriki, H.; Han, M.X.; Kubo, S. Experimental research of the heat transfer characteristics using a packed-bed of honeycomb ceramic for high temperature thermal storage system. Energy Procedia 2015, 69, 1059–1067. [Google Scholar] [CrossRef]
  10. Luo, Z.; Wang, C.; Xiao, G.; Ni, M.; Cen, K. Simulation and experimental study on honeycomb-ceramic thermal energy storage for solar thermal systems. Appl. Therm. Eng. 2014, 73, 622–628. [Google Scholar] [CrossRef]
  11. Cascetta, M.; Serra, F.; Arena, S.; Casti, E.; Cau, G.; Puddu, P. Experimental and numerical research activity on a packed bed TES system. Energies 2016, 9, 758. [Google Scholar] [CrossRef]
  12. Shah, R.K.; London, A.L. Laminar Flow Forced Convection Heat Transfer and Flow Friction in Straight and Curved Ducts—A Summary of Analytical Solution; Technical Report #75; Stanford University: Stanford, CA, USA, 1971. [Google Scholar]
  13. Bejan, A. Convection Heat Transfer, 4th ed.; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2013. [Google Scholar]
  14. Kays, W.M.; Crawford, M.E. Convective Heat and Mass Transfer; McGraw-Hill Series in Mechanical Engineering; McGraw-Hill Education: New York, NY, USA, 1966. [Google Scholar]
  15. Leveque, A. Les lois de la transmission de chaleur par convection. Ann. Mines Mem. Ser. 1928, 13, 381–412. [Google Scholar]
  16. Churchill, S.W.; Ozoe, H. Correlations for laminar forced convection with uniform heating in flow over a plate and in developing and fully developed flow in a tube. J. Heat Transf. 1973, 95, 78–84. [Google Scholar] [CrossRef]
  17. Kristiawan, B.; Santoso, B.; Wijayanta, A.; Aziz, M.; Miyazaki, T. Heat transfer enhancement of TiO2/water nanofluid at laminar and turbulent flows: A numerical approach for evaluating the effect of nanoparticle loadings. Energies 2018, 11, 1584. [Google Scholar] [CrossRef]
  18. Shah, R.K.; London, A.L. Laminar Flow Forced Convection in Ducts: A Source Book for Compact Heat Exchanger Analytical Data; Academic Press: Cambridge, MA, USA, 2014. [Google Scholar]
  19. Brkić, D.; Praks, P. Unified friction formulation from laminar to fully rough turbulent flow. Appl. Sci. 2018, 8, 2036. [Google Scholar] [CrossRef]
  20. Dai, G. Heat Transfer, 2nd ed.; Higher Education Press: Beijing, China, 1998. [Google Scholar]
  21. Li, P.; Chan, C.L. Thermal Energy Storage Analyses and Designs; Elsevier Inc.: Amsterdam, The Netherlands, 2017. [Google Scholar]
  22. Iaria, D.; Khader, M.; Al Zaili, J.; Sayma, A. Multi-objective optimisation of a centrifugal compressor for a micro gas turbine operated by concentrated solar power. In Proceedings of the Global Power and Propulsion Society Forum, Zurich, Switzerland, 16–18 January 2017. [Google Scholar]
  23. Iaria, D.; Al Zaili, J.; Sayma, A. Reducing Levelised Cost of Energy and Environmental Impact of a Hybrid Microturbine-Based Concentrated Solar Power Plant. In Proceedings of the Global Power and Propulsion Society Forum, Shanghai, China, 30 October–1 November 2017. [Google Scholar]
  24. Iaria, D.; Nipkey, H.; Al Zaili, J.; Sayma, A.; Assadi, M. Development and Validation of a Thermo-Economic Model for Design Optimisation and Off-Design Performance Evaluation of a Pure Solar Microturbine. Energies 2018, 11, 3199. [Google Scholar] [CrossRef]
Figure 1. Schematic of the thermal energy storage system.
Figure 1. Schematic of the thermal energy storage system.
Energies 12 03968 g001
Figure 2. Detail of the fluid domain mesh.
Figure 2. Detail of the fluid domain mesh.
Energies 12 03968 g002
Figure 3. Results of the conjugate transfer model validation: heat-transfer coefficient (right) and friction coefficient (left).
Figure 3. Results of the conjugate transfer model validation: heat-transfer coefficient (right) and friction coefficient (left).
Energies 12 03968 g003
Figure 4. Comparison between TEXSTAN and the present model for velocity profiles at different.
Figure 4. Comparison between TEXSTAN and the present model for velocity profiles at different.
Energies 12 03968 g004
Figure 5. Concentrated solar power (CSP) dish and micro gas turbine (MGT) power plant scheme.
Figure 5. Concentrated solar power (CSP) dish and micro gas turbine (MGT) power plant scheme.
Energies 12 03968 g005
Figure 6. Results of the conjugate heat-transfer analysis: velocity profiles at different locations.
Figure 6. Results of the conjugate heat-transfer analysis: velocity profiles at different locations.
Energies 12 03968 g006
Figure 7. Results of the conjugate heat-transfer analysis: temperature profiles at different locations.
Figure 7. Results of the conjugate heat-transfer analysis: temperature profiles at different locations.
Energies 12 03968 g007
Figure 8. Results of the conjugate heat transfer analysis Pressure profiles at different locations.
Figure 8. Results of the conjugate heat transfer analysis Pressure profiles at different locations.
Energies 12 03968 g008
Figure 9. Comparison between conjugate heat-transfer analysis (dots) and regression equation (lines) under variable diameter and Reynolds number. Nusselt number (left); Friction factor (right).
Figure 9. Comparison between conjugate heat-transfer analysis (dots) and regression equation (lines) under variable diameter and Reynolds number. Nusselt number (left); Friction factor (right).
Energies 12 03968 g009
Figure 10. Results of one-dimensional model validation: Comparison between different Nusselt relationships and conjugate heat-transfer analysis. The temperature variation of the solid with time (a) at x = 0; (b) at x = 0.05 m; (c) at x = 0.1 m; (d) at x = 0.15 m, (e) at x = 0.2 m; (f) the temperature variation of air with time.
Figure 10. Results of one-dimensional model validation: Comparison between different Nusselt relationships and conjugate heat-transfer analysis. The temperature variation of the solid with time (a) at x = 0; (b) at x = 0.05 m; (c) at x = 0.1 m; (d) at x = 0.15 m, (e) at x = 0.2 m; (f) the temperature variation of air with time.
Energies 12 03968 g010
Figure 11. Charging and discharging of the designed thermal energy storage (TES).
Figure 11. Charging and discharging of the designed thermal energy storage (TES).
Energies 12 03968 g011
Figure 12. Daily variation of Direct Normal Irradiation (DNI) (Blue) and output power (Yellow) for a typical cloudy day in Rome, Italy.
Figure 12. Daily variation of Direct Normal Irradiation (DNI) (Blue) and output power (Yellow) for a typical cloudy day in Rome, Italy.
Energies 12 03968 g012
Figure 13. Performances of the plant during day 2: (a) Power output and (b) turbine inlet temperature (TIT) (dashed line) and storage temperature (continue line).
Figure 13. Performances of the plant during day 2: (a) Power output and (b) turbine inlet temperature (TIT) (dashed line) and storage temperature (continue line).
Energies 12 03968 g013
Table 1. Infinite-series-solution functions for the circular tube, constant surface temperature, and thermal-entry length [14].
Table 1. Infinite-series-solution functions for the circular tube, constant surface temperature, and thermal-entry length [14].
n λ n 2 G n
07.3130.749
144.610.544
2113.90.463
3215.20.415
4348.60.383
>4 4 n + 8 / 3 1.01276 λ n 1 / 3
Table 2. Solid material properties.
Table 2. Solid material properties.
ρ S C S λ S
5000 kg/m31000 kJ/kgK5 W/mK
Table 3. Mesh independence study results.
Table 3. Mesh independence study results.
Number of Nodes 9.7   ×   10 5 1.8   ×   10 6 3.0   ×   10 6 4.5   ×   10 6 5.8   ×   10 6
Variation in the SolutionBaseline case8.7%2.3%1.8%0.3%

Share and Cite

MDPI and ACS Style

Iaria, D.; Zhou, X.; Al Zaili, J.; Zhang, Q.; Xiao, G.; Sayma, A. Development of a Model for Performance Analysis of a Honeycomb Thermal Energy Storage for Solar Power Microturbine Applications. Energies 2019, 12, 3968. https://doi.org/10.3390/en12203968

AMA Style

Iaria D, Zhou X, Al Zaili J, Zhang Q, Xiao G, Sayma A. Development of a Model for Performance Analysis of a Honeycomb Thermal Energy Storage for Solar Power Microturbine Applications. Energies. 2019; 12(20):3968. https://doi.org/10.3390/en12203968

Chicago/Turabian Style

Iaria, Davide, Xin Zhou, Jafar Al Zaili, Qiang Zhang, Gang Xiao, and Abdulnaser Sayma. 2019. "Development of a Model for Performance Analysis of a Honeycomb Thermal Energy Storage for Solar Power Microturbine Applications" Energies 12, no. 20: 3968. https://doi.org/10.3390/en12203968

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop