Numerical modeling of heat transfer in the fixed-matrix regenerator working in the Electric Thermal Storage heating system

The study presents the concept of Electric Thermal Storage (ETS) central heating system. Thermal Energy Storage (TES) is carried out in the fixed-matrix regenerator. The energy conservation equations, determined for the discharge period of the regenerator operation, are implemented in MATLAB numerical procedures based on the Finite Volume Method (FVM). In the model pressure drops within the system are calculated, both for the airflow through the inner tubes, and between the tubes. The flow distribution calculations show that the assumption of even air flow distribution would not be justified. Subsequently, the values of heat transfer coefficients are determined for the four distinct heat transfer surfaces, for the variable axial coordinate z and during the time of the system operation. The use of two different criterion equations is considered, for determining the mean Nusselt number Num for fluid flow through the concentric annular duct, as well as for the local Nusselt number Nuz calculated for the fluid flow through a circular or noncircular ducts. The most appropriate approach is selected by comparing the calculation results with experimental data. Taking into account the relative error, RMSE, and MAPE values calculated, it may be concluded that the Taler correlation – for non-circular ducts – gives results closer to the experimental data obtained.


Introduction
The presented study is related to sensible heat storage carried out in regenerative heat exchangers, also called as regenerators. Regenerative heat exchangers are classified as fixed-matrix and rotary regenerators. The fixed-matrix or fixed-bed regenerator is a periodic-flow heat transfer device with high thermal capacity matrix through which the hot fluid stream and cold fluid stream pass alternately [1,2]. Specific parameters determine regenerators performance, i.e., geometry, thermal and mechanical properties of the solid matrix, working fluid properties, as well as a charge and discharge time-length periods [3].
As shown by Ismail et al. [10], two broad groups of models of heat transfer between the working fluid and the regenerator bed are commonly used. The first group includes models in which the regenerator bed and the fluid is considered as one phase. Therefore they are called single-phase models. The second group comprises of models where the bed and the fluid domains are treated as a separate phase, hence the name -two-phase models. In the second group, three basic models may be distinguished, the Schumann's model, the continuous solid phase model, and the concentric dispersion model, that assumes a thermal gradient in the solid particles [10].
In the paper numerical investigation of heat transfer between fixed matrix and working fluid during cooling period is presented. Thus, heat balance equations for the regenerator's matrix and working fluid domains are derived for the concentric dispersion model. The heat balance equations are formulated for the novel design of fixed-bed regenerator that works as a Thermal Energy Storage (TES) unit in Electric Thermal Storage (ETS) central heating system. The governing system of differential equations with formulated boundary and initial conditions is solved using a Finite Volume Method. In the study, an appropriate correlation for determining heat transfer coefficients is also chosen. The numerical model is verified by comparing the computation results with experimental data.
The heating system with the regenerator that works as TES is shown in Fig. 1. Two periods may be distinguished during the TES unit operation. In the first stage, when the off-peak period of power system load occurs, the regenerators matrix is heated up with lowpriced energy. Twenty four electrical heaters, 500W of power output each, are installed in the TES unit bed what gives 12 kW of power output in total. The charging period is extended during the time to achieve uniform matrix temperature throughout the entire bed. The bed temperature may reach even 700°C and is measured in three various cross-sections of the regenerator bed. The test facility is equipped with the plate fin-and-tube heat exchanger, as shown in Fig. 1, being investigated and tested in a broad range of previous papers [11][12][13][14][15][16][17][18]. Subsequently, during the discharge period, when electricity price is high (daytime, and on-peak periods), the stored heat is gradually released and used for space heating. Thus, according to Fig. 1, the regenerator matrix is cooled by the air that is circulating in a closed circuit. The centrifugal fan, located at the regenerator inlet section, controls the airflow. When leaving the TES unit, hot air is cooled in the plate fin-and-tube heat exchanger and flows back through the air duct to the regenerator inlet. Hot water, heated up in the heat exchanger, feeds the central heating system which consists of two radiators located in the adjacent room. Thus, in the discharge period, the heat exchanger performs the function of a low-temperature water boiler.
The regenerator geometry is shown in Fig. 2. The bed consists of eight steel tubes with the length of 2.0 m being filled with ceramic cylinders arranged one by one in seven rows in each tube. There are 64 cylinders in one row. The total mass of the ceramic elements is 200,7 kg. The ceramic cylinders of diameter dc = 30 mm and 30 mm height, are made of refractory concrete, which thermal properties are as follows: specific heat, cc = 790 J/(kg·K); thermal conductivity, kc = 6.4 W/(m·K); density, ρc = 2700 kg/m 3 . Other bed parameters may be found in [19].
The experimental setup scheme along with measurement points marked is shown in Fig. 1. The air mass flow rate ṁf is determined using thermal mass flow meter installed in the air duct, while the water mass flow rate ṁw is measured using turbine flow meter installed at the plate fin-and-tube heat exchanger outlet. The air temperatures are measured at the regenerator inlet Tf,inlet , and outlet section Tf,outlet. The water temperature measurements are carried out both at the heat exchanger inlet and outlet, determining Tw,inlet , and Tw,outlet respectively. Temperature measurements are performed using pre-calibrated T-type thermocouples (Cu-CuNi).

Model formulation
In the previous studies, an analytical solution [8,20,21] and lumped-parameter model [19] of the regenerator bed cooling is presented. The models assumed that the bed is a porous medium, where radial heat conduction is neglected. The two-dimensional model, shown in this study, takes into account the bed geometry. The energy conservation equations are formulated for the solid as ceramic cylinders, steel tubes of the regenerator core and for the outer shell, as well as for the fluid -air flowing in the air gaps across the bed and between the tubes.
When taking into account the simplifying assumptions, that are listed below, the energy conservation equations for the fluid and solid domains may be derived. The following assumptions are made when formulating the model equations: 1. Heat conduction in the air in the axial and radial direction is neglected. 2. Contact between the neighboring solids is neglected so that there is no heat conduction between adjacent ceramics rows as well as between ceramics and steel tubes. 3. Air temperature is a function of time t and axial coordinate z and is uniform at a given regenerator cross-section. 4. Physical properties of the ceramics and steel are considered as temperature independent. 5. There is no heat transfer between the system and the surroundings.
According to the abovementioned assumptions, the heat balance equations for solid and fluid domains are written in the following form: where c subscript denotes the ceramics parameter and f,int stands for the air flowing in the inner tubes.  for the airflow inside the inner tubes, contacting with the ceramic cylinders ,, ,, where the air velocity inside the inner tubes is denoted as wf,int.  for the steel tubes of the regenerator core where subscript t indicates the tubes parameters, and st stands for steel.  for the air flowing between the tubes ,, , Air velocity in the space between the tubes is given as wf,out.  for the regenerator outer shell The system of partial differential equations, along with boundary conditions, is solved using Finite Volume Method [22,23]. The primary goal of the calculations is to determine the heat transfer fluid temperature at the regenerator outlet, for the axial coordinate z = L, where L is the regenerator length. Input parameters of the model are the total air mass flow rate ṁf , air temperature at the inlet section Tf,inlet, and initial temperature T0 of the whole bed.

Modeling and analysis
When determining the airflow distribution within the bed, pressure drop calculations in the system are performed. Basic principles of fluid mechanics state that the fluid flows at a higher rate through the channel, wherein encounters lower pressure drop. In this case, the flow resistance is determined, according to [3], for the subsequent parts of the regenerator, inlet and outlet sections, interior of the steel tubes, and space between the tubes. Based on the estimated pressure drops, it is possible to determine the air flow rates in the inner tubes ṁf,int and in the space between the tubes ṁf,out with the constraint that its sum must be equal to the total mass flow rate ṁf where p is the number of steel tubes in the regenerator core, for the considered regenerator p = 8. Acquiring information on air distribution is necessary when calculating heat transfer coefficients, and therefore modeling convective heat transfer in the wall-fluid system for the considered bed geometry. An exemplary experimental dataset, collected during the test stand operation, is shown in Fig. 3 and 4. As it is shown in Fig. 4, air flow rate is constant during the system operation. Therefore, for ṁf = 0.07891 kg/s values of ṁf,int and ṁf,out are calculated, and equals to ṁf,int = 0.00175 kg/s and ṁf,out = 0.06491 kg/s respectively. Thus, less than 17.5% of the total air flow ṁf goes through the inner tubes. In this case, the crucial role is played by linear pressure losses. The side surfaces roughness of the ceramic cylinders, approx. 0.3 mm, results in significantly higher pressure losses due to the liquid layers friction on the surface of the cylinder. For stainless steel, the roughness value is approx. 0.045 mm.  It can be seen in Fig. 4 that at the initial time outlet air temperature Tf,outlet is more than 15°C lower than for the next points marked. In this case, a dynamic temperature measurement error occurs which is related to the thermal inertia of the thermocouples [24][25][26].
After determining the flow distribution within the regenerator bed, calculation of heat transfer coefficients is performed. In the study, two approaches for determining the heat transfer coefficients for the air flow through the regenerator bed are compared. The first methodology is presented by Gnielinski in [3]. Herein a chapter devoted to heat transfer in concentric annular ducts may be found. In the second approach, the use of modified Petukhov and Kirillov correlation for the fluid flow through a circular or non-circular duct is considered. The modification is proposed by Taler in [27,28].

Heat transfer in concentric annular ducts
For the case of airflow through the tubes filled with ceramics, it is assumed that each of the cylinder rows is contacting with a ring-shaped air layer of equivalent diameter equal to 35.7 mm. While for the steel tubes in the space between the tubes it is 122.2 mm. The air mass flow rate through the single annular duct is an equal part of the flow rates ṁf,int and ṁf,out. So that, for the inside of inner tubes, the mass flow rate per a single row of ceramic cylinders is equal to ⅐ of ṁf,int, and for the single tube it is ⅛ of ṁf,out.
In the case of thermally and hydrodynamically developing flow in a concentric annular duct, the mean Nusselt number in the laminar region may be calculated from Martin correlation: where Num,1, Num,2 and Num,3 are the mean Nusselt numbers, calculated for Re < 2300, for boundary conditions specified in [3]. Moreover, K is a correction factor taking into account the effect of changing the fluid physical properties with the temperature where λ is a friction factor an annular duct, Pr is a Prandtl number, and Fann is a correction factor [3]. Hydraulic diameter dh is a calculated as  h o i d d d (11) where do is the outer diameter, and di is the inner diameter of the annuli. In the above equations, the ranges of validity are 10 4 < Re < 10 6 ; 0.6 < Pr < 1000; 0 ≤ dh/L ≤ 1.
For the transition region between laminar and fully developed turbulent flow, for 2300 ≤ Re ≤ 10 4  (13) The drawback of formulae (12) is the adoption of the Reynolds number limit value of Re = 2300, and Re = 10 4 , which marks the boundaries of the transitional flow region. So that, the first derivative dNu/dRe in the transition point is no longer continuous. Consequently, when changing the flow region, from the laminar to the transitional, at Re = 2300, and from the transitional to the turbulent region, at Re = 10 4 , a rapid change in the course of the function Nu = f (Re) occurs. For this reason, a new correlation is proposed by Taler [27,28] that is valid for the transitional and the turbulent flow region.

Heat transfer in non-circular ducts
As the second case, the correlation for determining the local value of the Nusselt number Nuz for the fluid flow in the plain tubes is used. The relationship is proposed by Taler [27,28] (15) where Nuz,1, Nuz,2 and Nuz,3 are the local Nusselt numbers, calculated for Re = 2300, for boundary conditions specified in [27]. In the equations mentioned above, λ is a friction factor a plain tube, calculated from Taler equation [29], z is the axial coordinate, changing from 0 to L, and dh is the hydraulic diameter calculated as 4 times the flow area divided by the wetted perimeter of the conduit. For z = 0 is the point from which heating or cooling of the fluid begins. So that, it is a coordinate of the regenerator inlet section.

Determination of heat transfer coefficients
The heat transfer coefficients are determined for the four different contact surfaces:  hc for the lateral surface of ceramic cylinders;  ht,int for the inner surface of the tubes filled with the ceramic cylinders;  ht,out for the outer surface of the tubes;  hsh for the internal surface of the outer shell.
Values of the relevant heat transfer coefficients are determined in the model for each time-step and for each of the N control volumes of the solid body being in contact with the fluid. Figures 5 and 6 show the time courses of the heat transfer coefficient values determined for various axial coordinates z, where z = 0 is the inlet, and z = 2 m is the outlet cross-section of the regenerator bed.

Experimental validation
Experimental verification of the proposed heat transfer correlations is carried out with comparing the temperature Tf,outlet, obtained from measurements, with mass-average temperature  (16) where Tf,int , and Tf,out are the air temperatures determined, respectively, within the steel tubes, and in the space between the tubes. In Figs 7 and 8, computation results of the air outlet temperature along with its experimental values are presented. To better compare the numerical and experimental results, the relative error is calculated, for each of the time-step number n, and for the both considered cases TT MAPE N T (19) where Nt is the total number of time steps. For the results shown in Fig. 7, when using Eqs (7) -(13), RMSE = 5,77°C, and MAPE = 5,43%. When considering Taler correlation, Eqs (14) - (15), RMSE = 2,71°C, and MAPE = 2,66%, see Fig. 8. As it is shown, the RMSE and MAPE values are more than 2 times lower for the Taler correlation. Based on the RMSE and MAPE values, as well as maximum values of the relative error, presented in Fig. 9, it may be concluded, that a good agreement between the experiment and calculation results is achieved when using Taler approach. A significant relative error for an initial time, exceeding 16% in both cases, is caused by a dynamic temperature measurement error. In the remaining time period, εf is not exceeding ±5%.

Conclusions
The study presents the concept of Electric Thermal Storage (ETS) central heating system. Thermal Energy Storage (TES) in the system is carried out in the fixedmatrix regenerator. In the article construction of the regenerator is presented. The test stand built in the laboratory of the Institute of Thermal Power Engineering at the Cracow University of Technology is used to test the proposed ETS central heating system.
The energy conservation equations, formulated for the two-dimensional model, takes into account the regenerator bed geometry. After the model validation, the developed heat balance equations are implemented in MATLAB numerical procedures based on the Finite Volume Method. In the model pressure drops within the system are calculated, both for the airflow through the inner tubes, and between the tubes. The flow distribution calculations show that in the tubes filled with ceramics flow resistance is much higher than for the flow between the tubes. Therefore, approx. 17.5% of the total flow rate ṁf goes through the eight tubes filled with the ceramic cylinders. Thus, the assumption of even air flow distribution would not be justified.
Subsequently, the values of heat transfer coefficients hc, ht,int, ht,out, and hsh are determined for the four distinct heat transfer surfaces, for the variable axial coordinate z and during the time of the system operation. The use of two different criterion equations is considered: first one, proposed by Gnielinski, for determining the mean Nusselt number Num for fluid flow through the duct of the annular cross-section; the second one is recommended by Taler for the fluid flow through a circular or non-circular ducts. The Taler correlation is further modified for determining local Nusselt numbers Nuz. The most appropriate approach is selected by comparing the calculation results with experimental data. Taking into account the relative error, RMSE, and MAPE values calculated, it may be concluded that better results are obtained in this case for the Taler correlation.
The numerical model may be used in the design and performance assessment of the presented accumulative heating systems. The development of mathematical procedures that allows determining the temperature of the air leaving the TES unit during the time, with the selected parameters of the system, i.e., mass flow rate, number of tubes in the core, shape, and material from which bed particles are made, etc. The numerical procedures may support the design and selection of devices at the planning stage of installation capacity. Moreover, the model may also be useful during the implementation of a model-based control system that is planned to be developed in the future.