MODELING OF TEMPERATURE FIELDS IN A SOLID HEAT ACCUMULLATORS

LLC «Teplotehnika», Yavornytskyi D. Ave., 102, Dnipro, Ukraine, 49000, tel./fax +38 (0562) 33 33 06, e-mail director@teplotehnika.dp.ua, ORCID 0000-0002-9935-4778 Dep. «Heat Engineering», Dnipropetrovsk National University of Railway Transport named after Academician V. Lazaryan, Lazaryan St., 2, Dnipro, Ukraine, 49010, tel./fax +38 (056) 373 15 76, e-mail ivatire@mail.ru, ORCID 0000-0002-5948-9483 Dep. «Heat Engineering», Dnipropetrovsk National University of Railway Transport named after Academician V. Lazaryan, Lazaryan St., 2, Dnipro, Ukraine, 49010, tel. +38 (056) 373 15 87, e-mail gabrin62@mail.ru, ORCID 0000-0002-6115-7162


Introduction
Currently, one of the priority areas of energyefficiency is to save costs on heating in industrial and residential buildings by the stored thermal energy at night time and its return in the daytime.As a result, savings are achieved due to the difference in tariffs for cost of electricity in the daytime and at night one.Change to «Night» tariff allows paying for electricity on an average three times cheap-er in comparison with the normal mode of operation [1].One of the most common types of devices that allow accumulating and giving the resulting heat obtained in different ways, are the heat accumulators (HA), heat pipes and thermosiphons [2][3][4][5][6].Thermal storage devices may be used to implement such principal tasks as performing the distribution of a source and receiver of thermal energy in space and time, as well as smoothing the tem-perature field on the surface or in the volume of the object.Thermal storage devices are most widely used in the energy, engineering, transportation, chemical industry, agriculture.Consequently, research and development of methods for determining the operating modes and the weight and dimensional parameters of HA is an important task of energy conservation, actual in the contemporary conditions of energy deficit.

Purpose
To date a large number of works about HA were published.The functioning of the HA in the process of heat storage can be realized by two main mechanisms: the first is due to changes of the physical parameters in the thermal storage solid (TSS); the second -through the use of the binding energy between atoms and molecules of substances.
Capacitance-type batteries are the most common and simple.Heat capacity of substance, heating without its aggregative state change is used in them.Typical HA structural scheme is shown in Figure 1.It shows that HA always consists of insulated and thermal storage solid (TSS), heater, cooling systems, safety, regulation of heat supply and removal.
For the weight and dimension calculations one limits with mass determination [1].In determining the HA modes, one considers the heat transfer processes using classical approaches of thermal fields analysis, as well as techniques based on mathematical modeling of heat transfer [7].Mathematical models of HA functioning are focused on the description of the HA thermal field [8][9][10] and cannot be directly applied for calculations of temperature field distribution, for example, when con-vective heat transfers on the HA charge and discharge mode.In order to determine temperature stresses one can use [9].However, proposed before calculation methods [10][11][12][13][14] do not reflect the picture of heat transfer at active convective transfer occurring at HA charging / discharging.
The main objective of the work is to develop a method for calculating the temperature field TSS in the process of heat accumulation and removal at the design stage on the basis of mathematical modeling of the temperature field in condition of strong convective heat transfer.
Solid HA is a complex of multiple systems connected in a single structure constructively.Heating system is a mandatory element of the HA, in our case it is tubular heating elements (THEs).Heat generated by them is accumulated in the thermal storage solid of -HA charging is made.To use the stored heat, HA has a cooling system, in our case there are air channels.With the active circulation of the coolant -air, heat is removed from the TSS and supplied to the consumer.Heatdistribution system within the heating object space does not include in to HA complex.
The design concept of solid HA with convective heat transfer is shown in Figure 2. HA consists of a jar 1 which can be fixed on any rigid support, the front jar is closed with battery cap 2, on the jar thermal insulation 3, 4 is fixed, in which TSS 5 is placed.On the front surface of TSS finger baffles 6 are mounted for the cooling air flow direction, which is fed to the bottom of TA through incoming louvers 7, then, passing through the HA channels, enters to the mixer 8 and go through the outlet louvers 9 falls within the scope of the object of heat supply.

Methodology
Design scheme for the analysis of the temperature field in the HA can be shown in Fig. 3, where the following notation is introduced: Li -heat insulation layer; C -channel; La - layer of ТSS, Tcu -the temperature of the upper boundary of the channel; Tcl -the temperature of the lower boundary of the channel Tau -the temperature of the upper boundary of ТSS; Tal -the temperature of the lower boundary of ТSS.

Fig. 3. Diagram of the temperature field analysis
If we neglect the change of heat fluxes along the coordinate x , which is directed perpendicular to the plan of ТSS 5 (Fig. 2), the temperature field will depend on three independent variables, namely: spatial coordinates y and z , and time t .Using ratio 1 t z V = , where 1 V there is air movement velocity on the channel C, one can reduce constitutive equations to the form, where two independent variables y и z will take place.Then we can write such heat transfer equations for the system shown in Fig.
where ρ , p C , λ − thermal and physical characteristics of the material: density, heat capacity ratio and conductivity coefficient (subscripts 1 and 2 are used respectively for air and TSS); T − temperature.
Each of the two equations will have two boundary conditions on the coordinate y and on one initial condition on the coordinate z .
The presence of the thermal insulation on the upper boundary of the channel, and the lower boundary of the TSS let neglect with heat flow out the heat accumulator boundaries, in other words one can record Two other boundary conditions can be represented as where 21 q − heat flow coming into the channel from the heated ТSS; 1b T -coolant temperature at the bottom surface of the channel; 2t T − temperature at the top surface of ТSS; 12 α -heat transfer coefficients between the cooling coolant and the top surface of the heated ТSS.Initial conditions correspondingly for equations (1) and (2) will be where 1 ( ) f y , 2 ( ) f y − temperature functional de- pendences from the coordinate y .
In the first approximation temperature functional dependences can be taken as constants.Then instead of ( 7) and (8) there will be To solve equations ( 1) and ( 2) we use the Laplace integral transformation [13,14].
Using the theorem about the differentiation of the original, we obtain the operator analogs of equations ( 1) and ( 2) in such form where L T -temperature image T , including appropriate indexes; s -Laplace transformation variable; ( ) ( ) Thus, using the Laplace integral transformation, the transition from partial differential equations (1) and (2) (in originals) to the differential equations in ordinary derivatives (in images), that are solved much easier.
Operator equations for the boundary conditions (3) -( 6) will look like this Solutions of equations ( 11) and ( 12) have the form For determining the integration constants 11 C ,

12
C , 21 C and 22 C it is necessary to differentiate the last two equations on coordinate y and substitute boundary conditions ( 13) -( 16).
Substituting boundary conditions ( 13) and (15) into equation ( 17), as well as − ( 14) and ( 16) into equation ( 18), we obtain (after determining the integration constants 11 C , 12 C , 21 C and 22 C ) such equations in the images for determining the temperature fields ( ) Taking into account an expression ( 16) and (20), one can write down such ratio at 0 y = Then ( ) In order to go from the temperature image to the original, write the hyperbolic functions through exponential x e e − = + , ( ) x e e − = − ).
After appropriate changes expression (21) can be presented as follow ( ) where [ ] ( ) By analogy with the expression (22) one can convert the equation (20), namely ( ) where [ ] Using the general formula of transition from the image to the original [2] ( ) and multiplication theory (Borel theorem) one can obtain from the expression (22) such original for temperature distribution in the solid plug along the y-axis ( ) where Using the same technique as in the obtaining of expression (25), we find from (23) the original for temperature field distribution in the TSS ( ) where To determine the heat transfer coefficient 12 α in equations ( 25 where l ε -a correction factor that takes into account impact of the ratio of the cooling cavity length 0 L to its equivalent size e b on the heat transfer coefficient; Re -Reynolds criterion; Pr -Prandtl number; Pr WT -Prandtl number at a wall temperature of the cooling cavity.For transitional regime calculation is recommended to carry out by the graph, shown in Figure 4, at this the value NP is determined by expression ( ) The following relationship is the most acceptable for laminar regime 0,33 0,43 0,15 Re Pr where Gr -Grashof number.
To determine the criteria, one can use such expressions Thus, to find the temperature field distribution in a two-layer system accordingly to Fig. 2 it is necessary to solve the equations ( 25) and (26).However, this system generally comprises two unknown quantities, namely: 1b T and 2t T .At this, given quantities in boundary conditions (6)  taken as constant ones.In reality, they will depend on the coordinate j z .To take into account the last remark and achieve the required accuracy of calculations, it should 1b T and 2t T to find on short segments along z axle.
The initial values 1n T and 2n T are also should be constantly changed at each segment.Thus, the final values of the temperature field distribution on the previous segment will correspond to the initial values at the next segment along the axis z .For determining the unknown boundary temperature values from these equations, we obtain the following system of equations ( ) In the last two equations index j characterizes the values of the corresponding ones on each segment j z at zero value for the second coordinate For convenience, the solution of equations ( 32) -( 33) are presented in a matrix form 1 , 0,0 0,1 0 where CV T = .
To solve the above problem program block in the mathematical MathCAD package was developed.Re-solving results are shown in Fig. 4, 5.In this case the initial values are following: The indices correspond to the following designations: 1-channel, 2-TSS.Designations correspond to Standards.One should take into consideration that the channel length L is determined by the number of baffles in HA.
Length temperature behavior of TSS under specified conditions is shown in Fig. 5.The number of partitions along the channel (iz) is 30, the number of partitions in channel depth and thickness of TSS (iy) is 20.The coordinate system corresponds to shown one in Figure 3.As can be seen from Fig. 5 the temperature at TSS, depth of 15mm increases from normal oneat the beginning of the channel and at a length of 2 m is already 570 о С.In the mid-plane of TSS (iy=10), the temperature will be higher and 670 о С.
Air temperature behavior in the channel along length under given conditions is shown in Fig. 6.The air temperature in the channel will change only slightly, at the outlet from TSS and input to the mixer will be 770 о С. Temperature behavior de- pending on selection of calculation fixed point will be also insignificant and vary within 1-2 о С.As seen from the graph (Fig. 7), depth temperature behavior of TSS has exponential nature, in depth of TSS varies within 50 о С.
Air temperature behavior in the channel in depth at a fixed channel length under given conditions is shown in Fig. 8.
As seen from the graph (Fig. 8), air temperature behavior in the channel in depth has a logarithmic character, by channel depth varies slightly within1-2 о С.
After analyzing the above data, one can draw the following conclusion: temperature behavior of TSS in depth and length has exponential nature, it is more essential along the length than depth.Air temperature behavior in the length and depth of the channel varies insignificantly.

Originality and practical value
Technical analysis shows that the proposed method of estimating the temperature field distribution of solid heat accumulator in different modes is effective, technically feasible and allows determining the operation modes of the solid heat accumulator at the specified weight and dimensional characteristics in the design stage of solid heat accumulators.

Conclusions
The method of calculation for temperature fields of solid heat accumulators on charging / discharging modes was proposed.

2
) and (26) one can use expression in common case[15]    where Nu -Nusselt criterion; e b -equivalent size of the channel.In general case, it is divided into three modes: turbulent of turbulent regime one can use the following expression to determine the Nusselt cri-

Fig. 4 .
Fig. 4. The graph to determine the Nusselt criterion for the transitional regime Equivalent size can be determined from the formula 4 g e S b P ⋅ = , where g S -square of effective cross-section; Pfull (wetted) perimeter, regardless of what part of the perimeter is involved in heat transfer.For heating liquid fluid one can take ( ) 0,25 Pr Pr 1 WT

Fig. 6 .
Fig. 6.Temperature curve along the channel length, at a fixed depth Temperature curve of TSS in depth at a fixed length and given conditions is shown in Fig. 7.

Fig. 7 .
Fig. 7. Depth temperature curve of TSS at a fixed length