Energy Storage in Concrete Bed

Abstract

The energy storage ability and temperature arrangement of a concrete bed which was charged and discharged at the same time was examined mathematically in this research. This was carried out by modeling a single globe-shaped concrete which was utilized to simulate a series of points along the concrete bed axis. Charging and discharging mode of the system were compared for 0.0094, 0.013, and 0.019 m3/s air flow rates. Higher change in temperature response was detected between the charging and fluid to solid heat transfer process at the inception of the concrete bed and the heat gain by the cool air flowing inside the copper tube was fairly high. The analysis of energy storage efficiency was also carried out and it was noticed that the globe-shaped concrete of 0.11 m diameter has the highest storage efficiency of 60.5% at 0.013 m3/s airflow rate.

1. Introduction

Heat can be transferred in a concrete bed through the following means: (i) heat transfer by convection from the bed wall to the fluid flowing inside the bed; (ii) heat transfer by convection from the globe-shaped concretes to fluid flowing in the bed, and this is known as fluid to particle mode; (iii) heat transfer by conduction from the bed walls to the globe-shaped concretes; (iv) heat transfer by conduction from one globe-shaped concrete to another, and this is known as particle to particle mode; (v) heat transfer by radiation; and (vi) heat transfer through fluid mixture [1]. The modes are as shown in Figure 1.

Figure 1.

Schematic showing the modes of heat transfers in concrete bed.

The particle to particle conduction mode can be further analyzed in axial and radial directions. Heat transfer through radiation mode will actually be important at higher temperatures.

Practically, it is found that two or three of the modes discussed earlier on can occur simultaneously. For example, the conduction between the particles may be affected by the convection between the particles and the fluid. This interaction among the different modes is one of the main reasons for the difficulty in correlating the total heat transfer and analyzing the experimental data in this field [2].

This research carried out numerically the temperature distribution in a concrete bed and also looked into the ability of a concrete bed to store energy simultaneously.

2. Review of literature

Anzelius [3] is the first author to publish a paper in heat transfer through packed beds but Schumann [4] is usually the first reference cited in most literature [1]. Both authors made some assumptions in order to find solution to equations that guide heat transfer for an incompressible fluids passing uniformly through a bed of solid particles with perfect conductivity. The following are the heat transfer equations derived for the system:

$$\frac{\left(T_s–T_{s,0}\right)}{\left(T_s–T_{s,0}\right)}=1–e^{-Y–Z}\mathrm{\Sigma }{Y}^n{M}^n\left(\mathit{yz}\right)=e^{-Y-Z}\mathrm{\Sigma }{Z}^n{M}_n$$ E1

$$\frac{\left(T_f–T_s\right)}{\left(T_{f,0}–T_{s,0}\right)}=1–e^{-Y–Z}\mathrm{\Sigma }\left[Y^n{M}^n\left(\mathit{yz}\right)\right]=e^{-Y-Z}\mathrm{\Sigma }\left[Z^n{M}_n\left(\mathit{yz}\right)\right]$$ E2

where T_s is the solid temperature, T_f is the temperature of fluid and Y and Z are dimensionless quantities. Eqs. (1) and (2) was solved simultaneously and the result were analyzed graphically to form Schumann curves. Knowing the outlet air and bed temperature, these curves can be used to determine the volumetric heat transfer coefficients and also the heat transfer coefficients of a packed bed undergoing heat exchange with a fluid provided the following conditions which were the simplifying assumptions made by Schumann were satisfied:

1. Due to the infinitesimal nature of the solid particles, resistance to heat transfer was so small.

2. Resistance to conduction heat transfer in the fluid was so small.

3. At any section of the bed, the heat transfer rate from fluid to solid or from solid to fluid was directly proportional to their average difference in temperature between the solid and fluid within the bed.

4. The transport properties of solid and fluid were not dependent on temperature, for example the density.

Furnas [5] utilized and expanded the Schumann curves to cover more range of temperatures and suggested an empirical relation for determination of heat transfer coefficient as shown in Eq. (3):

$$h_v=\frac{{\mathit{BG}}^{0.7}{T}^{0.3}{10}^{1.68\epsilon –3.56\epsilon 2}}{{d_p}^{0.9}}$$ E3

where, h_v is the volumetric heat transfer coefficient. B is a constant dependent on the bed material, G is the mass velocity of the fluid, T is the average air temperature, d_p is the particle diameter and ε is the porosity.

Saunders and Ford [6] utilized dimensional analysis to derive correlations to calculate heat transfer coefficient. The research was for spherical shaped alone and the application is not suitable for other solid particle geometries.

Another correlation for the determination of heat transfer coefficient between gases and randomly packed solid spheres was postulated by Kays and London [7]. Using the Colburn j-factor, the relationship was written as:

$$j_h=\frac{0.23}{{{\mathit{Re}}_p}^{0.3}}$$ E4

$$\text{where},j_h=\mathit{St}.{\mathit{Pr}}^{2/3}$$

$$200<{\mathit{Re}}_p<50,000\mathit{\text{and}}0.37<\epsilon <0.39$$

Löf and Hawley [8] studied heat transfer between air and packed bed of granitic gravel. Unsteady state heat transfer coefficients were correlated with the air mass velocity and particle diameter to obtain the equation:

$$h=0.652{\left(G/d_p\right)}^{0.7}$$ E5

This was calculated for 8 mm < d_p < 33 mm; 50 < Rep 500 and temperature range of 311–394 K. The author reached a conclusion that the temperature of the entering air had no appreciable effect on the hat transfer coefficient.

Leva et al. [9] studied and analyzed heat transfer coefficient between smooth spheres of low thermal conductivities and fluids (air and carbon dioxide) in packed beds and tubes of 50.8 and 6.4 mm diameters, respectively. The ratio of particles to tube diameters was varied from 0.08 to 0.27; gas flow rate was of Reynolds number range 250 to 3000. Correlation of film coefficient was determined as:

$$h=3.50\left(\frac{k}{D_t}\right)e^{-4.6\frac{\mathit{Dp}}{\mathit{Dt}}\left(\frac{\mathit{DpG}}{\mu }\right)0.7}$$ E6a

which is approximately:

$$h=0.40\left(k/D_t\right)\left(\mathit{DpG}/\mu \right)0.7$$ E6b

$$\text{Or},N=0.4{\mathit{Re}}^{0.7}$$ E6c

Maximum film coefficient was predicted and verified at a value of Dp/Dt equal 0.153.

Riaz [10] and Jefferson [11] studied the dynamic behavior of beds undergoing heat exchange with air using single and two phased modes. By incorporating factors of axial bed conduction and intra-particle resistance, which Schumann ignored, the heat transfer coefficients were evaluated and found to be 1 + Bi/5 times smaller than those predicted using Schumann curves.

Ball [12], Norton [13], Meek [14], Bradshaw and Meyers [15], Harker and Martyn [16] and also, Bouguettaia and Harker [17] have all researched on various packed beds using air and other gases as fluids and have developed correlations involving the heat transfer coefficient.

3. Methodology

3.1. Heat transfer model for a globe-shaped concrete bed

The modeling of heat transfer in a concrete bed was carried out mathematically. It was done through a single globe-shaped concrete which was simulated mathematically to represent series of points along the concrete bed axis.

A one dimensional finite difference formulation was used in modeling the single globe-shaped concrete material, where heat conduction to neighboring globe-shaped concrete was ignored.

Using this assumption reduced the globe-shaped concrete model to that of an isolated sphere in cross flow, where the total surface area of the sphere was exposed to convection. Also, the thermal properties of the materials within the bed accounted for temperature dependence.

3.2. Finite difference formulation of a single spherical shaped concrete material

Since conduction to other globe-shaped concrete has been neglected, the geometry allows the concrete to be reduced to one dimension along its radius.

A finite difference method was utilized to model this mathematically, [18]. For this approach, the globe-shaped concrete can be characterized by three different nodal equations:

1. a general, interior node

2. the center node

3. the surface node

All exposed to convection as shown in Figure 2.

Figure 2.

Schematic showing the cross-section of the globe-shaped concrete materials imbedded with copper tube.

For the general and interior node within the globe-shaped concrete model, the conduction equation for T_(r,t) is:

$${\rho }_c{C}_c\frac{\partial T}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(K_c{r}^2\frac{\partial T}{\partial r}\right)+{\dot{q}}_{\left(r,t\right)}$$ E7

where C_c = specific heat of concrete.

$K_c=$ thermal conductivity of concrete.

$\dot{q}=$ heat generation.

And this equation was represented in finite difference form.

The specific heat, thermal conductivity, and the heat generation, are temperature dependent and varied with the temperature along the radial direction.

Because the thermal properties are functions of temperature, and consequently functions of the globe-shaped concrete radius, the finite difference equations are derived by the volume integration over a finite difference node.

Multiplying Eq. (7) by r² and integrating both sides of the equation from r_n – Δr/2 to r_n + Δr/2 resulted to:

$${\rho }_c{C}_c\frac{\partial T}{\partial t}\underset{r_{n-\frac{\mathrm{\Delta }r}{2}}}{\overset{r_{n+\frac{\mathrm{\Delta }r}{2}}}{\int }}{r}^2\mathit{dr}=\underset{r_{n-\frac{\mathrm{\Delta }r}{2}}}{\overset{r_{n+\frac{\mathrm{\Delta }r}{2}}}{\int }}\frac{\partial }{\partial r}\left(K_c{r}^2\frac{\partial T}{\partial r}\right)\mathit{dr}+\underset{r_{n-\frac{\mathrm{\Delta }r}{2}}}{\overset{r_{n+\frac{\mathrm{\Delta }r}{2}}}{\int }}{r}^2{\dot{q}}_{\left(r,t\right)}\mathit{dr}$$ E8

The specific heat was assumed constant with respect to r, and therefore brought outside the integral.

By evaluating the integrals in Eq. (8) and representing the derivatives in finite difference form using the fully implicit method gives:

$${\rho }_c{C}_c\left(\frac{T_n^{Z+1}-T_n^Z}{\mathrm{\Delta }t}\right)\left(\frac{r_{n^{+}}^3-r_{n^{-}}^3}{3}\right)=\underset{r_{n-\frac{\mathrm{\Delta }r}{2}}}{\overset{r_{n+\frac{\mathrm{\Delta }r}{2}}}{\int }}2K_c\mathit{rdT}+\dot{q}\underset{r_{n-\frac{\mathrm{\Delta }r}{2}}}{\overset{r_{n+\frac{\mathrm{\Delta }r}{2}}}{\int }}{r}^2\mathit{dr}$$ E9

$${\rho }_c{C}_c\left(\frac{T_n^{Z+1}-T_n^Z}{\mathrm{\Delta }t}\right)\left(\frac{r_{n^{+}}^3-r_{n^{-}}^3}{3}\right)={\left[\frac{2K_c{r}^2}{2}\frac{\mathit{dT}}{\mathrm{\Delta }r}\right]}_{r_{n-\frac{\mathrm{\Delta }r}{2}}}^{r_{n+\frac{\mathrm{\Delta }r}{2}}}+{\dot{q}}_n{\left[\frac{r^3}{3}\right]}_{r_{n-\frac{\mathrm{\Delta }r}{2}}}^{r_{n+\frac{\mathrm{\Delta }r}{2}}}$$ E10

$$\begin{array}{l}{\rho }_c{C}_c\left(\frac{T_n^{Z+1}-T_n^Z}{\mathrm{\Delta }t}\right)\left(\frac{r_{n^{+}}^3-r_{n^{-}}^3}{3}\right)=K_{c\left(n^{+}\right)}{r}_{n^{+}}^2\left(\frac{T_{n+1}^{Z+1}-T_n^{Z+1}}{\mathrm{\Delta }r}\right)-\\ K_{c\left(n^{-}\right)}{r}_{n^{-}}^2\left(\frac{T_n^{Z+1}-T_{n-1}^{Z+1}}{\mathrm{\Delta }r}\right)+q_n{\left[\frac{r_{n^{+}}^3-r_{n^{-}}^3}{3}\right]}_{r_{n-\frac{\mathrm{\Delta }r}{2}}}^{r_{n+\frac{\mathrm{\Delta }r}{2}}}\end{array}$$ E11

where $T_n^Z$ and $T_n^{Z+1}$ indicate temperatures for an arbitrary node at times $t^Z$ and $t^{Z+1}$ respectively.

Also,

$${\mathrm{W}}_n=\mathrm{W}\mathrm{at}{T}_n^{Z+1}$$ E12

and,

$$K_{n^{+}}=\frac{K\left[T_{n+1}^{Z+1}+K\right]}{2},\mathrm{at}{T}_n^{Z+1}$$ E13

also,

$$r_{n^{+}}=r_n+\frac{\mathrm{\Delta }r}{2}$$ E14

$$r_{n^{-}}=r_n-\frac{\mathrm{\Delta }r}{2}$$ E15

Eq. (11) can be rearranged and solved for $T_n^Z$ plus a known quantity which resulted to:

$$\begin{array}{l}{\rho }_c{C}_c\left(\frac{T_n^{Z+1}-T_n^Z}{\mathrm{\Delta }t}\right)=\frac{3K_{c\left(n^{+}\right)}{r}_{n^{+}}^2}{r_{n^{+}}^3-r_{n^{-}}^3}\left(\frac{T_{n+1}^{Z+1}-T_n^{Z+1}}{\mathrm{\Delta }r}\right)-\\ \frac{3K_{c\left(n^{-}\right)}{r}_{n^{-}}^2}{r_{n^{+}}^3-r_{n^{-}}^3}\left(\frac{T_n^{Z+1}-T_{n-1}^{Z+1}}{\mathrm{\Delta }r}\right)+{\dot{q}}_n\left[\frac{3r_{n^{+}}^3-r_{n^{-}}^3}{3r_{n^{+}}^3-r_{n^{-}}^3}\right]\end{array}$$ E16

$$\begin{array}{l}{\rho }_c{C}_c\frac{T_n^{Z+1}}{\mathrm{\Delta }t}-{\rho }_c{C}_c\frac{T_n^Z}{\mathrm{\Delta }t}=\frac{3K_{c\left(n^{+}\right)}{r}_{n^{+}}^2}{r_{n^{+}}^3-r_{n^{-}}^3}\left(\frac{T_{n+1}^{Z+1}-T_n^{Z+1}}{\mathrm{\Delta }r}\right)-\\ \frac{3K_{c\left(n^{-}\right)}{r}_{n^{-}}^2}{r_{n^{+}}^3-r_{n^{-}}^3}\left(\frac{T_n^{Z+1}-T_{n-1}^{Z+1}}{\mathrm{\Delta }r}\right)+{\dot{q}}_n\end{array}$$ E17

Multiply Eq. (17) by Δt and divide by ρ_cC_c resulted to:

$$\begin{array}{l}{T}_n^{Z+1}-T_n^Z=\frac{3K_{c\left(n^{+}\right)}{r}_{n^{+}}^2\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\left[T_{n+1}^{Z+1}-T_n^{Z+1}\right]-\\ \frac{3K_{c\left(n^{-}\right)}{r}_{n^{-}}^2\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\left[T_n^{Z+1}-T_{n-1}^{Z+1}\right]+\frac{3\mathrm{\Delta }t{\dot{q}}_n}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}}\end{array}$$ E18

$$\begin{array}{l}\therefore T_n^Z+\frac{3\mathrm{\Delta }t{\dot{q}}_n}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}}=T_n^{Z+1}-\left[\frac{3K_{c\left(n^{+}\right)}{r}_{n^{+}}^2\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right]\left[T_{n+1}^{Z+1}\right]+\\ \left[\frac{3K_{c\left(n^{+}\right)}{r}_{n^{+}}^2\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right]\left[T_n^{Z+1}\right]+\left[\frac{3K_{c\left(n^{-}\right)}{r}_{n^{-}}^2\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right]\left[T_n^{Z+1}\right]+\\ \left[\frac{3K_{c\left(n^{-}\right)}{r}_{n^{-}}^2\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right]\left[T_{n-1}^{Z+1}\right]\end{array}$$ E19

Collecting the like terms from Eq. (19) yielded:

$$\begin{array}{l}\therefore T_n^Z+\frac{3\mathrm{\Delta }t{\dot{q}}_n}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}}=\left[\left(\frac{3\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r}\right)\frac{K_{c\left(n^{-}\right)}{r}_{n^{-}}^2}{\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right]\left[T_{n-1}^{Z+1}\right]+\\ \left[1+\left\{\left(\frac{3\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r}\right)\frac{\left(K_{c\left(n^{-}\right)}{r}_{n^{-}}^2+K_{c\left(n^{+}\right)}{r}_{n^{+}}^2\right)}{\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right\}\right]\left[T_n^{Z+1}\right]-\\ \left[\left(\frac{3\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }r}\right)\frac{K_{c\left(n^{+}\right)}{r}_{n^{+}}^2}{\left(r_{n^{+}}^3-r_{n^{-}}^3\right)}\right]\left[T_{n+1}^{Z+1}\right]\end{array}$$ E20

This resulting equation is valid for any general, interior node within the globe-shaped concrete 0 < r_n < R.

At the center node, where r_n = 0 the temperature profile is axisymmetric, and $\frac{\partial T}{\partial r}=0,$ when r = 0 thus, the temperature on either side of the node is equal.

$$T_{n-1}^{Z+1}=T_{n+1}^{Z+1}$$

$$\mathit{\text{likewise}},K_{c\left(n^{-}\right)}=K_{c\left(n^{+}\right)}$$

$\therefore$ Eq. (20) simplified to:

$$\therefore T_n^Z+\frac{\mathrm{\Delta }t{\dot{q}}_n}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}}=T_n^{Z+1}-\left(\frac{6\mathrm{\Delta }t}{{\rho }_{c\left(m\right)}{C}_{c\left(m\right)}\mathrm{\Delta }{r}^2}{K}_{c\left(n^{+}\right)}\right)T_{n+1}^{Z+1}$$ E21

This occur at r_n = 0.

This simplified form of Eq. (20) was used to represent the center node.

The conduction through the surface of the globe-shaped concrete is equal to the convection at the surface.

$$\therefore -K{\frac{\partial T}{\partial r}}_{\mathit{at},r=R}=-U_C\left(T_{r=R}-T_{\infty }\right)$$ E22

However, this boundary condition cannot be directly represented in finite difference form, since such formulation requires a volume element and Eq. (22) applies at a point.

Instead a first law energy balance was utilized to obtain the nodal equation for the surface of the globe-shaped concrete. This energy balance can be written as:

$${\dot{E}}_{\mathit{\text{in}}}-{\dot{E}}_{\mathit{\text{out}}}+{\dot{E}}_{\mathit{gen}}={\dot{E}}_{\mathit{st}}$$ E23

where,

$${\dot{E}}_{\mathit{\text{in}}}=-\mathit{KA}\frac{\partial T}{\partial r}$$ E24

$${\dot{E}}_{\mathit{\text{out}}}=U_{C}A\left(T-T_g\right)$$ E25

$${\dot{E}}_{\mathit{gen}}=\dot{q}V$$ E26

$${\dot{E}}_{\mathit{st}}=\mathit{\rho CV}\frac{\partial T}{\partial t}$$ E27

Representing Eq. (23) in a finite difference form consistent with Eq. (20) and (21) resulted to:

$$-\mathit{KA}\frac{\partial T}{\partial r}-U_{C}A\left(T-T_g\right)+\dot{q}V=\mathit{\rho CV}\frac{\partial T}{\partial t}$$ E28

This can be written in finite difference form to give:

$$\begin{array}{l}-4\pi r_{n^{-}}^2{K}_{n^{-}}\frac{T_n^{Z+1}-T_{n-1}^{Z+1}}{\mathrm{\Delta }r}-4\pi r_n^2{U}_C\left(T_n^{Z+1}-T_{\infty }\right)+\frac{4}{3}\pi \left(r_n^3-r_{n^{-}}^3\right){\dot{q}}_n=\\ \frac{4}{3}\mathit{\pi \rho }{C}_n\left(r_n^3-r_{n^{-}}^3\right)\frac{T_n^{Z+1}-T_n^Z}{\mathrm{\Delta }t}\end{array}$$ E29

where, U_c = convection coefficient.

Solving for $T_n^Z$ plus known quantities involving $\dot{q}$ and U_c, in a similar manner to Eq. (20) and (24) resulted to:

$$\begin{array}{l}\frac{-4\pi r_{n^{-}}^2{K}_{n^{-}}{T}_n^{Z+1}}{\mathrm{\Delta }r}+\frac{-4\pi r_{n^{-}}^2{K}_{n^{-}}{T}_{n-1}^{Z+1}}{\mathrm{\Delta }r}-4\pi r_n^2{U}_C{T}_n^{Z+1}-4\pi r_n^2{U}_C{T}_{\infty }+\\ \frac{4}{3}\pi \left(r_n^3-r_{n^{-}}^3\right){\dot{q}}_n=\frac{4}{3}\mathit{\pi \rho }{C}_n\frac{\left(r_n^3-r_{n^{-}}^3\right)}{\mathrm{\Delta }t}{T}_n^{Z+1}-\frac{4}{3}\mathit{\pi \rho }{C}_n\frac{\left(r_n^3-r_{n^{-}}^3\right)}{\mathrm{\Delta }t}{T}_n^Z\end{array}$$ E30

Multiply Eq. (30) by Δt and divide by $\frac{4}{3}\mathit{\pi \rho }\left(r_n^3-r_{n^{-}}^3\right)$ resulted to:

$$\begin{array}{l}{T}_n^Z-T_n^{Z+1}=\frac{3\mathrm{\Delta }t}{\mathrm{\Delta }\mathit{r\rho }{C}_n}\left(\frac{r_{n^{-}}^2{K}_{n^{-}}}{r_n^3-r_{n^{-}}^3}\right)T_n^{Z+1}-\frac{3\mathrm{\Delta }t}{\mathrm{\Delta }\mathit{r\rho }{C}_n}\left(\frac{r_{n^{-}}^2}{r_n^3-r_{n^{-}}^3}\right)T_{n-1}^{Z+1}+\\ \frac{3\mathrm{\Delta }{\mathit{tU}}_C}{\rho C_n}\left(\frac{r_{n^{-}}^2}{r_n^3-r_{n^{-}}^3}\right)T_n^{Z+1}+\frac{3\mathrm{\Delta }{\mathit{tU}}_C}{\rho C_n}\left(\frac{r_n^2}{r_n^3-r_{n^{-}}^3}\right)T_{\infty }-\frac{{\dot{q}}_n\mathrm{\Delta }t}{\rho C_n}\end{array}$$ E31

$$\begin{array}{l}{T}_n^Z+\frac{\mathrm{\Delta }t}{\rho C_n}{\dot{q}}_n+\frac{3\mathrm{\Delta }{\mathit{tU}}_C}{\rho C_n}\left(\frac{r_{n^{-}}^2}{r_n^3-r_{n^{-}}^3}\right)=-\left(\frac{3\mathrm{\Delta }t}{\mathrm{\Delta }\mathit{r\rho }{C}_n}\frac{r_{n^{-}}^2{K}_{n^{-}}}{r_n^3-r_{n^{-}}^3}\right)T_{n-1}^{Z+1}+\\ \left[1+\frac{3\mathrm{\Delta }t}{\rho C_n\mathrm{\Delta }r}\left(\frac{r_{n^{-}}^2{K}_{n^{-}}}{r_n^3-r_{n^{-}}^3}\right)+\frac{3\mathrm{\Delta }{\mathit{tU}}_C}{\rho C_n}\left(\frac{r_n^2}{r_n^3-r_{n^{-}}^3}\right)\right]T_n^{Z+1}\end{array}$$ E32

$$\begin{array}{l}{T}_n^Z=-\left(\frac{3\mathrm{\Delta }t}{\mathrm{\Delta }\mathit{r\rho }{C}_n}\frac{r_{n^{-}}^2{K}_{n^{-}}}{r_n^3-r_{n^{-}}^3}\right)T_{n-1}^{Z+1}+\\ \left[1+\frac{3\mathrm{\Delta }t}{\rho C_n\mathrm{\Delta }r}\left(\frac{r_{n^{-}}^2{K}_{n^{-}}}{r_n^3-r_{n^{-}}^3}\right)+\frac{3\mathrm{\Delta }{\mathit{tU}}_C}{\rho C_n}\left(\frac{r_n^2}{r_n^3-r_{n^{-}}^3}\right)\right]T_n^{Z+1}-\frac{\mathrm{\Delta }t}{\rho C_n}{\dot{q}}_n-\frac{3\mathrm{\Delta }{\mathit{tU}}_C}{\rho C_n}\left(\frac{r_{n^{-}}^2}{r_n^3-r_{n^{-}}^3}\right)\end{array}$$ E33

Eqs. (20), (21) (32) and (33) constitute a system of algebraic equations for heat transfer modeling in globe-shaped concrete.

4. Result and discussion

The values of Eq. (33) are obtained from the values in Table 1. Since the thermal properties are constant, average temperatures could therefore be used to determine thermal properties of bed materials.

Parameters	Values
Airflow rate	0.01316 m^3/s (28 cfm)
Air—density	1.07154 Kg/m^3
Air—specific heat capacity	1008 J/Kg K
Concrete—density	2400 Kg/m^3
Concrete—specific heat capacity	1130 J/Kg K
Copper tube—density	8900 Kg/m^3
Copper tube—specific heat capacity	384 J/Kg K
Area of globe-shaped concrete	0.013 m^2
Area of copper tube + header	0.664 m^2
Volumetric heat transfer coefficient	106.5 W/m^3 K

Table 1.

Modeling parameters.

The following data were obtained from the theoretical/mathematical modeling carried out on thermal performance of packed bed energy storage system as shown in Figure 3.

Figure 3.

Schematic of the storage tank systems.

The following are the definitions of the symbols:

Time = the interval time of measurements, in minutes.

T_s-in = the inlet air temperature to the packed bed storage tank in °C.

T_s-out = the outlet air temperature from the packed bed storage tank in °C.

T_t-in = the inlet air temperature to the copper tube in °C.

T_t-out = the outlet air temperature from the copper tube in °C.

T_A1, T_A2, T_A3, and T_A4 = the air stream temperatures (°C) through the bed at different heights of the storage tank 117.5, 235, 352.5, and 470 cm, respectively.

T_ci1, T_ci2, T_ci3, and T_ci4 = the core temperatures of the globe-shaped concrete (°C) through the bed at different heights of the storage tank 117.5, 235, 352.5, and 470 cm, respectively.

T_ti1, T_ti2, T_ti3, and T_ti4 = the temperatures of air flowing inside the copper tube (°C) through the bed at different heights of the storage tank 117.5, 235, 352.5, and 470 cm, respectively.

T_ct1, T_ct2, T_ct3, and T_ct4 = the temperatures of the contact made between globe-shaped concrete and imbedded copper tube (°C) through the bed at different heights of the storage tank 117.5, 235, 352.5, and 470 cm, respectively.

T_t1, T_t2, T_t3, and T_t4 = the surface temperatures of the copper tube (°C) through the bed at different heights of the storage tank 117.5, 235, 352.5, and 470 cm, respectively.

The results of the experimentation were shown in Figures 4, 5, 6 for globe-shaped concrete of size 0.11; 0.08 and 0.065 m diameter respectively while the discharging only temperature measurements were shown in Figures 7–9 respectively for air flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 4.

Average temperature measurement of charging packed bed storage system for globe-shaped concrete of size 0.11 m diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 5.

Average temperature measurement of charging packed bed storage system for globe-shaped concrete of size 0.08 m diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 6.

Average temperature measurement of charging packed bed storage system for globe-shaped concrete of size 0.065 m diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 7.

Average temperature measurement of discharging packed bed storage system for globe-shaped concrete of size 0.11 m diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 8.

Average temperature measurement of discharging packed bed storage system for globe-shaped concrete of size 0.08 m diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 9.

Average temperature measurement of discharging packed bed storage system for globe-shaped concrete of size 0.065 m diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 10 presents the comparison of the temperature variations with time at Ts-in, T_s-out, T_t-in, T_t-out, T_A1, T_A2, T_A3, T_A4, T_ci1, T_ci2, T_ci3, T_ci4, T_ti1, T_ti2, T_ti3, T_ti4, T_ct1, T_ct2, T_ct3, T_ct4, T_t1, T_t2, T_t3, and T_t4 during the simultaneous charging and discharging while Figure 11 presents for discharging only. The comparisons were presented for air flow rates of 0.0094, 0.013, and 0.019 m³/s.

Figure 10.

Comparison of average temperature measurement of charging packed bed storage system for globe-shaped concrete of size 0.065, 0.08, 0.11 m in diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

Figure 11.

Comparison of average temperature measurement of discharging packed bed storage system for globe-shaped concrete of size 0.065, 0.08, 0.11 m in diameter and flow rate of 0.0094, 0.013, and 0.019 m³/s.

These figures show that the difference of the temperature response between the charging and fluid to solid heat transfer process at the initial period (<30 min) of the packed bed was large (large inlet–outlet temperature difference means large heat supply), and the heat recovered by the cool air (approximately 27°C) flowing inside the copper tube was fairly high (larger inlet–outlet temperature difference compared with the later period indicates larger heat recovery).

Therefore, a relatively large part of the heat supplied by the simulated air heater was used to heat the air flowing inside the copper tube through conduction and convection and also stores the rest for continuous usage.

The following are the storage efficiency for globe-shaped concrete of size 0.11 m, 0.08 m and 0.065 m diameter at airflow rate of 0.0094, 0.013 and 0.019 m³/s (Figure 12):

Figure 12.

Storage efficiency of simultaneous charging and discharging packed bed storage system for globe-shaped concrete of diameter 0.065, 0.08, 0.11 m and air flow rate 0.0094, 0.013, and 0.019 m³/s.

For 0.11 m diameter globe-shaped concrete:

1. Storage efficiency at air flow rates of 0.0094 m³/s = 40.7%

2. Storage efficiency at air flow rates of 0.013 m³/s = 60.5%

3. Storage efficiency at air flow rates of 0.019 m³/s = 57.5%

For 0.08 m diameter globe-shaped concrete:

1. Storage efficiency at air flow rates of 0.0094 m³/s = 23.5%

2. Storage efficiency at air flow rates of 0.013 m³/s = 51.3%

3. Storage efficiency at air flow rates of 0.019 m³/s = 50.2%

For 0.065 m diameter globe-shaped concrete:

1. Storage efficiency at air flow rates of 0.0094 m³/s = 14.8%

2. Storage efficiency at air flow rates of 0.013 m³/s = 35.06%

3. Storage efficiency at air flow rates of 0.019 m³/s = 40.3%

5. Conclusion

The study led to the following findings and conclusions:

1. The mathematical model developed can accurately predict the temperature within the concrete bed for energy storage purpose.

2. The steady intermittent input temperature variation actually led to continuous discharge temperature at the copper tube outlet.

3. The mathematical model may be extended to specify the packed bed storage system dimensions.

4. Globe-shaped concrete of 0.11 m diameter has the highest storage efficiency of 60.5% at 0.013 m³/s airflow rate.