**6. Experimental Verification of the Developed Numerical Model of Accumulator Discharge**

The packing of the heat accumulator under test consisted of cylindrical ceramic elements and pipes, as well as other steel elements. For modelling, the following thermophysical properties and dimensions of the accumulator construction elements were assumed: Prandtl number for air, Pr*<sup>g</sup>* = 0.7; air thermal conductivity, *kg* = 0.236 W/(m·K); kinematic viscosity of the air, *<sup>ν</sup><sup>g</sup>* = 16.35·10−<sup>6</sup> m2/s; air density at 50 ◦C, *<sup>ρ</sup><sup>g</sup>* = 1.09 kg/m3; air specific heat at constant pressure, *cp*,*<sup>g</sup>* = 1000.0 J/(kg·K); air velocity before the packing, *wg* = 2.0 m/s; stainless steel thermal conductivity, *km* = 14.7 W/(m·K); stainless steel density, *<sup>ρ</sup><sup>m</sup>* = 7800 kg/m3; steel specific heat, *cm* = 519 J/(kg·K); total mass of steel elements, *mm* = 252.8 kg; heat accumulator casing cross-sectional area, *Ar* = 2.17 m2; length of heat accumulator, *Lr* = 2 m; thermal conductivity of ceramic cylinder, *kw* = 6.4 W/(m·K); density of ceramic cylinder, *<sup>ρ</sup><sup>w</sup>* = 2700 kg/m3; specific heat of ceramic cylinder, *cw* = 519 J/(kg·K); total mass of ceramic packing, *mw* = 252.8 kg; and ceramic packing outer surface area, *Aw* = 22.57 m2.

Air temperature measurements were carried out using K-type sheathed thermocouple sensors with a grounded junction. All thermocouples were pre-calibrated using a Pt100 platinum resistance sensor. The thermocouples were selected so that when measuring air temperatures below 100 ◦C, their readings would be the same with an accuracy of ±0.1 K. Uncertainty of temperature measurement in the air temperature range of −20 ◦C to +100 ◦C was assumed as for a resistance temperature sensor Pt100; i.e., ±0.35 K

Air velocity was measured using a digital vane anemometer with a maximum resolution of 0.01 m/s. The uncertainty of the sensor was ±0.5% of the upper limit of the measuring range or ±1.5% of the measured value.

Due to the lack of dependence in calculating the heat transfer coefficient between the air and the elements of the heat accumulator with such a complex packing structure, a procedure for determining the mean heat transfer coefficient was developed.

The accumulator we built and tested was characterised by a complex flow and heat system inside the packing. Part of the air stream flowed inside the tubes in which the cylindrical filling elements were placed, and the other part flowed outside between the steel tubes. The air-side heat transfer coefficient *ha* was approximated using the following formula:

$$h\_{\mathcal{S}} = a + b\overline{w}\_{\mathcal{S}}^{\varepsilon} \tag{71}$$

where *a*, *b*, and *c* are constants to be determined.

The heat transfer coefficient *hg* in Formula (71) is expressed in W/(m2·K) and the mean air velocity through the accumulator *w<sup>c</sup> <sup>g</sup>* in m/s.

The mean velocity of air flow through the whole cross-section of the accumulator is calculated as: .

$$
\overline{w}\_{\mathcal{S}} = \frac{\dot{m}\_{\mathcal{S}}}{\overline{\rho}\_{\mathcal{S}} A\_{\mathcal{S}}} \tag{72}
$$

The free area *Ag* of the cross-section of the accumulator through which air flows is given by:

$$A\_{\mathcal{S}} = \varepsilon \pi R\_w^2 \tag{73}$$

where *ε* is the porosity of the entire accumulator packing and *Rw* is the internal diameter of the accumulator casing.

All constants in Equation (71) were selected so that the measured and calculated air temperature at the outlet from the accumulator agreed as well as possible. Constants *a*, *b*, and *c* were determined from the condition:

$$S = S(a, b, c) = \sum\_{i=1}^{n\_m} \left[ T\_{\mathbb{g},out}^{m \text{cas}}(t\_i) - T\_{\mathbb{g},out}^{calc}(t\_i) \right]^2 \to \min \tag{74}$$

The air temperature *Tcalc <sup>g</sup>*,*out* at the outlet from the accumulator for selected time points *ti*, *i* = 1, . . . , *nm* was calculated using the developed mathematical model.

The constants *a*, *b*, and *c* in correlation (71) were determined from condition (74). The constants *a*\*, *b*\*, and *c*\* at which the function *S* reached its minimum were determined using the Nelder–Mead method [38]. A library program for optimisation using the Nelder–Mead method can also be found in the IMSL library [39]. Calculations were performed with the use of MATLAB [40] using the FMINSEARCH program. In this paper, the number of measurement points was 1956. The air temperature at the battery outlet was measured with a step equal to 5 s. The total air temperature measurement time was 9775 s.

Calculations of the heat accumulator under test were carried out for the following input data:


The time changes in air velocity (*wg*,*inlet*) and temperature at the accumulator inlet (*Tg*,*inlet*) are shown in Figure 11.

**Figure 11.** Time changes of air velocity and temperature at the accumulator inlet.

The following values of parameters appearing in Equation (71) were estimated using FMINSEARCH [39]:

$$a = 1.2192, \; b = 8.42 \cdot 10^{-5}, \; c = 13.1104 \tag{75}$$

The time changes in the heat transfer coefficient calculated using Formula (71) and the coefficients determined using Equation (74) are shown in Figure 12.

**Figure 12.** Changes in the heat transfer coefficient *hg* on the surface of the heat accumulator packing elements as a function of time.

When analyzing the results shown in Figure 12, it can be seen that the average heat penetration coefficient on all air-flown fill elements was low. This was due to the low air flow velocity inside the accumulator, as the free cross-section through which the air flowed was large.

When analyzing the results presented in Figure 13, it can be seen that at the beginning of accumulator cooling, the differences between the measured and calculated air temperature at the accumulator outlet were more significant. The maximum value of the relative difference *eT* between the measured *T g meas* and the calculated *T g calc* of the accumulator outlet air temperature, which was calculated using the relationship *eT* = *T g meas* <sup>−</sup> *T g calc* / *T g meas* · 100%, did not exceed 25%. For times greater than about 700 s, the relative difference *eT* between the calculated and measured values of air temperature decreased to values in the order of a few percent. This was mainly due to the approximate heat transfer model in the ceramic cylinders and steel filling elements. The mathematical model of the accumulator assumed that these elements were modelled as objects with concentrated heat capacity. This meant that the temperature difference inside

the ceramic cylinders and the walls of the steel tubes was not taken into account. In reality, the temperature of the elements changed with time and inside the elements as well.

**Figure 13.** Measured air temperature at the accumulator inlet and comparison of the measured and calculated air temperature at the accumulator outlet as a function of time.

### **7. Conclusions**

Based on the calculations and measurements, the following conclusions were reached: The investigated heat accumulator was heated at night during the period of low electricity demand when the cost of 1 kWh is low. The packing of the accumulator was heated to a temperature of about 600 ◦C so that the mass of the accumulator could be smaller at a given heat demand.

The developed mathematical models of the accumulator—an analytical model, a numerical model based on the explicit finite-difference method, and the implicit Crank– Nicolson method—gave similar results.

The comparison between the calculated and measured air temperatures at the accumulator outlet was satisfactory except for the initial air heating period, when the differences between the numerical solution and the measured results were more significant.

The differences in the calculated and measured values of the air temperatures at the accumulator outlet in the initial phase of the cooling of the accumulator packing were partly due to the assumption that the temperature of the filling elements was constant in the entire volume and changed only in time (model with lumped heat capacity). However, at the beginning of the cooling of the filling elements, the temperature difference inside ceramic cylinders or steel tubes could be significant. The second reason for the differences between the calculated and measured temperatures at the air outlet from the accumulator was the complex structure of the accumulator filling, which, for example, made it difficult to determine the velocity distribution in the cross-section of the packing, as well as to find appropriate correlations for calculating the heat transfer coefficient. Due to the small length of the accumulator, the heat transfer between the packing elements and the flowing air was strongly influenced by the inlet section with a developing flow, where the heat transfer coefficient was much higher than for the developed air flow.

The developed numerical model of the accumulator can be used to calculate the appropriate sizing so that it can be used to heat a building during a period of high electricity prices. The model can also be used to calculate the changes in the air temperature at the outlet of the accumulator during the entire period of its discharge (heating of a building). The elaborated dynamic mathematical model of the accumulator can also be used in an automatic temperature control system for a heated room (building) based on the mathematical model. The designed heat accumulator can replace a conventional fossil-fuel-fired boiler in a boiler room or operate as an independent unit in a heated room. In both cases, the mathematical model of the accumulator can be used in the system for automatic air temperature control in the heated space.

**Author Contributions:** Methodology, formal analysis, writing—original draft, supervision, D.T.; conceptualisation, investigation, software, validation, T.S.; conceptualisation, methodology, software, validation, supervision, writing—review and editing, J.T. (Jan Taler); investigation, data curation, J.T. (Jarosław Tokarczyk). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
