Theoretical modeling and experimental study of auxiliary concrete accumulator for solar heating systems

This article presents the issue of the use of the control-volume finite-elements method (CVFEM) to solve transient heat conduction problem in the ceramic or concrete structure of heat storage system. The system can be used as auxiliary storage in solar based domestic hot water (DHW) and heating installations. The storage system consists of modular symmetric components forming parallel air channels. The modular design and symmetry of cross section enables to build a simple numerical model using a coarse mesh of finite volumes. It allows solving the problem using the simple algorithm. Analyzed method of modeling provides a short computation time while maintaining high calculation accuracy.


Introduction
The use of the thermal capacity of solid materials for heat storage is an alternative to heat storage systems in water tanks. Despite the lower thermal capacity, the accumulation material can be incorporated into the structure of a building or can be a part of its construction. It creates many possibilities for innovative solutions for heating systems. The idea of heat storage using the thermal capacity of solid materials has been known since ancient times (heating systems based on hypocaust). In solar heating systems, accumulators made of solids are often used in ACPS (active collectionpassive storage) systems [1]. In ACPS systems, heat obtained from solar collectors (air, liquid) is transferred directly to heated objects characterized by high thermal capacity. In systems with active charging and discharging, the air flow is forced.

Fig.1. Possible accumulator locations in the building
Adequate selection of geometric parameters of the accumulator, thermal properties of the storage material and convective heat transfer conditions enable heat utilization spread over the day. Ordinary concrete or concrete with improved accumulation properties, chamotte, ceramics, stone and even sand may be used to build the accumulator. The accumulator can be made in the form of a closed space, wall or a horizontal slab. Examples of locations of the heat accumulators in buildings are presented on Fig. 1. Another group of accumulators consists of accumulators using hollow units with ducted airflow [7][8][9][10][11]. One of the pioneers in that field was James Kachadorian. In 1978 he patented a heating system which used solar energy from the direct gains system or from air collectors and stored it in a gravel-filled or hollowed sub-floor space with air distributing system by natural or forced circulation [8]. Short and Kutscher [1] investigated a hybrid building heating system with a concrete accumulator with a channel structure in the form of a concrete floor slab The air heated in the collector was directed to the airwater heat exchanger, which, as a result, heated domestic hot water. In the research carried out by YuXiang and Athienitis [12][13][14][15], the use of a ventilated ceiling slab (VCS Ventilated Concrete Slab) was described as a heat accumulating element obtained from a solar installation equipped with a BIPV (Building Integrated Photovoltaic) collector integrated with the roof.

The concept and the model 2.1 The concept
The concept of the analyzed heating system is based on the simultaneous storage of thermal energy in a water tank and concrete or ceramic accumulator with a channel structure. The heat used to preheat hot water and support the heating of a building is obtained from a traditional solar installation equipped with liquid collectors. Firstly, the water tank is heated, and the surplus of acquired energy is accumulated in the structure of a concrete or ceramic exchanger -an accumulator which is charged and discharged by the air stream (flow in internal channels). During charging, the air circulates in a closed circuit and is heated in a plate-fin heat exchanger. At night, the accumulator is discharged by the stream of air blown directly into the heated space.

Experimental facility
For the system described above, a model installation was proposed. It consists of an accumulator made of concrete units forming two parallel air channels, air circulation system, two-speed coil plate-fin heat exchanger, and water tank heated by a solar collector. During the experiments, alternative heating with a gas boiler with an electric stabilizing heater was also used. The main element of the installation is a heat accumulator made of 16 concrete blocks, used in a fireplace building technology arranged in two parallel rows. Accumulating material fills a rectangular prism measuring 2240mm x 560mm x 280mm. Inlets and outlets from air ducts were profiled with steel sheet fittings. The accumulating material weights 689.6 kg.

Fig. 3. Experimental facility real view
The heat exchange surface of internal ducts is 2.75m 2 , and the hydraulic diameter of a single channel is 178mm. The density of accumulating material is equal ρ=2820kg/m 3 , specific heat capacity c=0.94 kJ/kgK and thermal conductivity λ=1.7W/mK.

Mathematical model
The mathematical model of the heating system based on the assumption that the heat for the storage system was supplied in a stream of water taken from the domestic hot water tank heated by the solar installation. The supply water temperature can, therefore, change according to the scheme of obtaining heat from the solar collectors and disassembling the hot water in the installation. The developed and verified the mathematical model of the plate-fin heat exchanger [17,18] had been used. For modeling transient heat transfer in the accumulator structure, it has been assumed that for the proper description it is sufficient to analyze the temperature distribution in multiple modules limited by the symmetry axes of the cross-section [11]. The concrete blocks used for the verification of the experimental model have a cross-section as presented in Fig. 4 characterized by four symmetry axes. An element corresponding to 1/8 of the cross-section of the fitting can be adopted for modeling the temperature field. In the experimental installation, two parallel rows were made of 8 blocks.  It was assumed that the temperature gradient in the material structure in the z-axis direction within individual modules is negligible. The heat conduction equation for steady-state, two-dimensional temperature field in the part of an accumulator cross-section determined by the symmetry planes ( Fig. 4) has the following form The outer surface of the air duct is thermally insulated. The convective boundary condition at the inner surface of a duct is as follows: The energy conservation equation for the air flowing through the internal channel is given by The initial temperature of concrete units and air are uniform | =0 = | =0 = 0 (4) The air temperature at the accumulator inlet Ta,in is also known | =0 = , ( ) (5) The inlet temperature Ta,in is equal the air temperature downstream a plate fin-and-tube heat exchanger. The initial-boundary value problem (5) was solved using Control-Volume Finite-Element Method (CVFEM) for determining the temperature in concrete segments and finite difference method (FDM) for determining the air temperature. The finite difference grid for modeling air flow is defined as follows: where Δz is equal the length of a concrete unit (Figures 4,and,5). The coordinate zi refers to the center of the unit.
To solve the problem of transient heat transfer in the solid accumulator, the following assumptions were adopted:  the channel was divided into n = 8 segments, with a length corresponding to the length of a single accumulation module.  the temperature of the storage mass was computed at the cross-section located in the middle of the length of each module (indices i = 1, ...,n)  the staggered grid was used in the air domain  section j = 1 corresponds to the inlet plane to the duct (inlet air temperature)  the cross-sections j, and j+1 correspond to the inlet and outlet of air and coincide with contact planes of the concrete modules  the cross-section j = n+1 corresponds to the outlet plane from the duct. Equation The air temperature at the point zi was assumed to be equal the arithmetic mean of the temperatures Ta,j and Ta,j+1 (Fig. 5) , = , + , +1 2 The superscript k+1in Eq. (8) = 1, … ,9, = 0,1, … Next, the determination of temperature of concrete unit will be discussed. The control volume finite element method (CVFEM) [18,19] is used to solve the problem of two-dimensional transient heat transfer in the planes of the accumulator cross-section, This method allows modeling the temperature field in a repeatable section of the accumulator section using a triangle mesh based on a relatively small number of nodal points. In the CVFEM method, an energy balance is created for the control area surrounding each node. The control area around the node 1 is formed by one-third of all of the triangles with sides coming out of the joint node 1. The energy balance equation is written for the 1-a-O-c area, in which the point O is the center of gravity of the triangle 123, and the points a and c are respectively the centers of the sides 1-2 and 1-3 (Fig. 6). The surface area of a figure 1-a-O-c is one-third of a surface area of a triangle 123. In the case of a single triangular element 123, the energy balance equation for the control volume associated with the node 1 can be written as follows: )] s 1a + q ṡ s 1c (14) Energy balance equations written for all control areas form a system of ordinary differential equations for node temperatures.
In a numerical solution, it is convenient to set a system of balance equations in a matrix form. For perfect insulation of the outer accumulator surface, the equation system is When the thermal insulation is not perfect, then heat loss occurs on the outer accumulator surface, and equation system has the following form where: D column matrix of temperature-time derivatives at nodes 1 ... n, A matrix containing surface areas of elements around nodes, B matrix of coefficients in terms including conduction, C matrix of coefficients regarding convection, E matrix of coefficients defining the length of the element boundary section with the heat flux set on the boundary, T column matrix of temperatures at the nodes and the air temperature in the channel. For the adopted cross-section with 16 The matrices A, B, C, and E, can be generated quickly entering the coordinates of individual nodes.
In the case of the analyzed concrete accumulator, the symmetrical fragment of the cross-section was first divided into triangular elements and then 16 control areas were created. Because of the symmetry six fictitious nodes were added (Fig.7). The derivative of temperature in the normal direction (heat flux) to the symmetry planes 1-2-3-4 and 13-14-16 are equal zero. If the outer surface of a concrete unit was insulated perfectly then, zero heat flux was assumed. In this case, the entire external surface of the analyzed area  Calculations were also carried out for the case of imperfect insulation when there were heat losses from the outer accumulator surface to the environment (nodes 4-5- 10-11-15-16).
The system of ordinary differential equations of the first order for the node temperatures at solid and air domains was solved by the explicit finite difference method. The number of equations for the whole accumulator is 16×(n-1), where n=8 is the number of modules forming the accumulator channel. The use of the explicit method is a significant advantage in online control systems based on the mathematical model of the accumulator. The explicit method does not require iteration or solving a system of algebraic equations at each time step when controlling the air temperature at the accumulator exit.
The Gnielinski formula (20) was applied in the analysis of heat exchange in concrete accumulator channels.

Theoretical calculations
At first, the mathematical model developed in the paper was used to analyze the process of heating and cooling the accumulator under test conditions. The constant temperature of 60 o C (corresponding to the maximum real temperature obtained from the solar source) was maintained to charge the accumulator. The analysis was performed for three different flow rates 1.1 m/s, 2.2 m/s, and 4.4 m/s (Figs. 8 and 9).

Experimental verification of the model
Based on the developed mathematical model, Fortran software was written. The software allows to simulate the system operation using input data (calculated or measured) ie. air and water temperatures, air flow velocity upstream the air-cooled heat exchanger. In the tests concerning the full charging and discharging cycle of the accumulator, a comparison between measured heat flow rates on the water and air side of the heat exchanger and the accumulator were done. Also the energy taken from the air and transferred to the concrete accumulator during given period were compared. In addition, the results of temperature measurements at selected points of the solid storage mass were compared with the results of modeling.

Fig. 10. Diagram of measurement points
To determine the heat flux transferred from the insulation to the environment in the charging period, the heat flux was measured, and compared with computations. The overall heat transfer coefficient between the outer surface of the concrete duct and environment was estimated based on experimental tests. This allowed to take into account heat losses by insulation to the environment in the accumulator mathematical model. The time variations of air temperatures and accumulated energy obtained by the mathematical model show good compliance with the measurements. A comparison of the temperature of the air, temperature of the concrete accumulator at two points and energy accumulated in the concrete is depicted in Figures 11,12, and 13 for wa=2.0m/s (Rea =21600).

Conclusions
The proposed method of modeling a heating system with a duct concrete heat accumulator allows simulation of any charging and discharging of the accumulator for any time changes of heating water and air supply parameters. The simple cross-section of the accumulation element adopted for the tests enabled effective laboratory tests and their comparison with numerical simulations. In practice, the analyzed heating system can be used in conditions of a large daily amplitude of ambient air temperature. In the practical implementation of the accumulator, it is recommended to use more extended heat transfer surfaces in the accumulator channels.