Numerical modelling method of heat exchange in heat accumulators with many phases and free phase boundaries movement

A correct implicit unconditionally stable numerical method for modelling heat transfer in heat accumulators with many phases and free phase boundaries movement of has been developed and presented. The stability of the implicit numerical scheme for approximating the energy equation is studied. As an example, a hardening process of zinc casting with a thickness of 10 cm was considered and compared with the results of calculations by another author.


Introduction
The thermal energy accumulation problem via phase transition is one of the most important at present. Its solution is complicated by many different conditions, regime parameters and states of matter. Until now, a large number of scientific papers have been devoted to the problems of heat transfer during phase transitions and increase over time. This can be estimated by the number of international publications. The most significant article is considered by Stefan [1], which proposed a solution to the problem of ice melting in the Arctic Ocean, it was later called the "Stefan problem". The total number of international publications citing this article is about 900, approximately 50 publications were published in 2016 and 2017, and 60 each in 2018 and 2019. Thus, interest in the class of heat transfer problems in phase transitions is steadily growing. It should be noted that the solidification of a round drop of liquid was first investigated by Lame back in 1831 [2]. The solution of this class of problems allows one to predict not only ice melting in the ocean, but also to study heat transfer in the soil [3], during metal casting and 3D printing [4], as well as when creating heat accumulators. Freezing of the soil leads to the death of the shallow root system of trees, which affects the growth activity of forests [5]. In general, the class of tasks is extremely wide and their solution can significantly affect energy efficiency. The use of heat accumulators with energy consumption can have the greatest effect. It is known that the schedule of consumption of heat and electric energy is uneven due to the periodicity of seasons and daily activities of a person. The accumulation of energy at the consumer during periods of time with a minimum load of the energy source and its expenditure during periods of peak loads can reduce the unevenness of the load on the source, reduce its installed power and the price of the generated energy. As a stimulus measure, measures have already been taken for multi-tariff metering of electricity in many countries. This allows the consumer to efficiently use heat accumulators for heating, storing energy at night. A fairly voluminous modern review of practical applications, accumulation agents, and methods for calculating heat transfer during phase transitions was made in a review article [6].
Experimental studies are not very difficult, however, to predict and summarize the results of studies, it is necessary to use computational studies, which requires the search for a unified technique that is well suited for all cases. To calculate and design systems with a phase transition, it is necessary to determine the temperature field in a substance, the boundaries of regions with phase transitions, and their change in time. One of the correct methods for numerically solving the phase transition problem is described in [7]. Heat transfer during phase transformations is an unsteady process in which the phase state of a substance changes and the phase boundaries move. The latter greatly complicates the task, since it is necessary to solve the equation of displacement of these boundaries. They have a boundary condition for the energy equation. Non-stationary processes in practice do not have a constant boundary condition, i.e. the temperature at the boundaries of the computational domain can be at different instants of time both above and below the temperature of the phase transition. This means that in the general case there can be unlimited number of phase zones in the considered region and the solution of such a problem by the method described in [7] becomes impossible, as well as the search for a solution that is established over time with periodically changing boundary conditions. In addition, there can also be many phase states of a substance, which further complicates the task. The problems of phase transitions can be classified according to the following criteria: the number of phases of a substance (two, three or more), the number of phase boundaries in the computational domain (one, two or more), the coordinate system (rectangular, cylindrical, spherical), the number of measurements (one, two or three), the presence of internal sources of heat, boundary conditions. Boundary conditions can be classified by gender and given law of non-stationary dependence. Since the result of solving many non-stationary problems can be represented as a set of harmonic solutions, the formulation of the problem established in the period of interest is of the greatest interest. As one of the limiting cases, the transition problem should also be considered. Earlier in [8], the correct method for solving the problem of phase transitions of the first kind, called the enthalpy [9], was considered and it was shown that the application of the classical difference scheme allows solving the problem only by the explicit conditionally stable method, which is often inconvenient and waste. In the present work, we propose a correct implicit method for numerical modelling of the heat exchange established during the heat exchange period with a multitude of phase transition boundaries in the computational domain and a multitude of phases.

Problem statement
The heat transfer in the heat accumulator in the absence of convection in the one-dimensional case for the cylindrical and rectangular coordinate systems is described by a linear parabolic differential energy equation where T is the temperature, K; t is time, s; x is the coordinate, m; ρ is the density, kg/m 3 ; cp is the specific heat, J/(kg·K); is the coefficient that allows to choose a coordinate system. An analysis of the types of the differential energy equation for solving the problem of phase transitions was carried out in [8]. Based on this analysis, the energy equation must be written as follows: where h is the enthalpy, J/kg. This type allows one to determine the change in the enthalpy of a substance over time and does not require determining the boundaries of phase transitions during the solution. The solution of equation (1) does not take into account the specific features of the phonon transfer of heat through the molecular structure of matter and corresponds to the theory of a continuous medium; therefore, it can be used to solve the problems of crystal growth with a strict orientation without defects, problems with a purely random orientation of the molecules, tasks when the number of defects is large and the distance between them is much smaller than any other scales. Since in such a formulation it is not necessary to establish the position of the phase boundaries, it becomes possible to solve problems 3 with many such boundaries. This formulation of the problem makes it possible to simulate enantiotropic transitions between allotropic states of matter. With some modifications and additions to the mathematical model, it is possible to solve the problems of monotropic transformations, such as, for example, the destruction of the crystal lattice of white tin and the transition to gray. To close the system of equations, one should describe the dependence of temperature on enthalpy. A qualitative example of such a dependence with two phase transitions is shown in Fig. 1.
where b is the coefficient of the linear equation, K; Is the enthalpy of transition of a substance from one to another state, n is the number of transitions between phases. The value of the coefficient b can be determined by the value of the partial dependence for h = 0. In this case, the phase is understood not only as the allotropic state of the substance, but also the state of the substance in the process of transformation, i.e., during melting and crystallization, evaporation and condensation, etc. Obviously, in the state of transformation, the heat capacity of the substance p c , coefficient b takes on the value of the saturation temperature, and the heat of the phase transition r takes on a certain value. Outside the transformation state, the specific heat cp and coefficient b take some values, and the phase transition heat r = 0. Dependence (2) is monotonic, continuous and consists of linear segments; therefore, it can be represented as a monotonic continuous piecewise linear formula For a greater generality of the obtained solutions and a decrease in the number of influencing parameters, it is necessary to bring it to a dimensionless form.
where   To solve the established periodic problems, the energy equation (1) should be reduced to a dimensionless form.
where   ω is circular frequency, s -1 . The ratio 0  is a certain normalization coefficient acting within the differential increment. When calculating the dimensional coordinate x, the size of all increments from the origin should be taken into account. Such a determination of the dimensionless coordinate X allows us to exclude the influence of changes in the density of various phases on the coordinate axis of the computational domain. The characteristic physical quantities x0, ρ0 are taken depending on the conditions of the problem being solved. To solve transient problems, energy equation (1) can also be reduced to a dimensionless form.
where  

Numerical solution method
Consider the system of equations (3) and (4). It can be numerically solved by the finite difference method. The second derivative is approximated with a second order of accuracy, which gives an unconditionally stable scheme. The stability of the solution to the system was studied in detail in [8]. It is shown that it is possible to develop only a conditionally stable scheme, while the most economical is an explicit scheme. However, in the process of obtaining a difference scheme, it is possible to develop an unconditionally stable economical implicit scheme. To do this, write the classic implicit scheme on a four-point grid pattern: where i, j is the numbers of the grid points along the X and tω axes, Δtω is the steps along tω and ΔX is the step along the X axis.
The stability of the circuit with respect to the time coordinate was confirmed in [8]. It is known that the approximation of the second derivative by the central difference scheme also leads to an unconditionally stable solution. It is seen that at the central point of the difference equations we have , ij H and , ij  . It is possible to write down a piecewise-defined dependence of temperature on enthalpy, however, to record the dependence of enthalpy on temperature, one should resort to a fraction of the phase, which is not convenient. It is possible to search for the solution of equation (6) by fixing the value of the right-hand side iteratively. In this case, a non-stable solution can be obtained when the change in the left side is greater than in the right. A study of the stability of such an iterative solution was carried out in [8] and showed that the scheme is conditionally stable. However, the substitution of formula (3) into equation (6) makes it possible to develop a new scheme.

Conclusion
A numerically modelling method of phase transitions processes in heat accumulators is developed and presented. The stability of the developed two-layer numerical schemes for approximating the energy equation with the monotonous continuous piecewise-linear formula of the dependence of temperature on enthalpy is substituted for it. The unconditional stability of the developed scheme is shown. Based on the proposed method, it is possible to study the appearance, disappearance, and displacement of the boundaries of the phase zones in the processes of heat and mass transfer. The method is applicable for calculating heat transfer with many boundaries of phase transitions in the calculation domain and with many phases of the substance. The developed method allows you to create software tools for designing heat accumulators with phase transitions. The scope of the method is wide enough and, if applied, can significantly reduce the load on the generation, transmission and energy consumption systems, and will increase the safety and stability of many energy processes.