An effective method for calculating the elements of thermal power plants, which are reduced to solving systems of partial differential equations

. Calculations of dynamic processes in the elements of thermal power plants (TPP) (heat exchangers, combustion chambers, turbomachines, etc.) are necessary to justify permissible and optimal operating modes, the choice of design characteristics elements, assessing their reliability, etc. Such tasks are reduced to solving partial differential equations. At present time for such calculations are mainly used finite-difference method and finite element method. These methods are cumbersome and complex. The article proposes a method, the main idea of which is to reduce the solution of equations to solving linear programming problems (LP) is demonstrated by the example heat exchanger of periodic action. The mathematical description includes the following energy balance equations for gas and ceramics, respectively, on the plane, where - indicates the length of the heat exchanger, and - the operating time. Also provides a more complex model, taking into account the spread of heat inside the balls of the ceramic backfill.


Introduction
Calculations of stationary and non-stationary modes of operation of a number of elements of thermal power plants (heat exchangers of various types, furnaces, combustion chambers, turbine grids, etc.) are reduced to solving systems of partial differential equations (PDE). The main methods for solving such systems are finite difference methods (FDM), finite volume methods (FVM) and finite element methods (FEM). When using FDM, a grid is constructed on the computational domain and for each of its nodes, based on the initial differential equations, a subsystem of algebraic equations is formed [1][2][3][4]. In these equations, the partial derivatives are replaced by the corresponding finite differences. Subsystems of algebraic equations of individual grid nodes are combined into a single system of algebraic equations, to which boundary conditions are added. It should be noted that, in this case, the accuracy of solving the system of PDE strongly depends on the values of the grid steps in spatial coordinates and in time. The desire to increase the accuracy of the solution leads to a reduction in the size of steps and, accordingly, to an increase in the number of grid nodes and the dimension of the system of algebraic equations. In many cases, this dimension becomes so large that the system cannot be solved as a whole without using one or another decomposition method. Reducing the dimension of the system can be achieved by using a grid with variable steps, but this greatly complicates the algorithm for solving the problem, which is especially noticeable for a computational space with complex geometry.
FVM is applicable to problems in which differential equations reflect the laws of conservation of mass (total or individual chemical elements), energy and momentum [5,6]. Most of the problems of heat and mass transfer belong to such tasks. Therefore, this method is most widely used in computational fluid dynamics. In accordance with FVM, the computational area is divided into control volumes for which an irregular geometric shape is permissible. For each volume, balance equations are formed that take into account the exchange of a given volume with adjacent volumes of mass, energy and momentum. These equations are algebraic, in which the derivatives are replaced by finite differences determined from the values of the corresponding parameters at the geometric centers of adjacent control volumes. In addition, the equations include the areas of the boundary surfaces between adjacent control volumes. Moreover, the balances of mass, energy and momentum are observed for the control volumes, regardless of the location of the surfaces dividing adjacent volumes. FDM allow more accurate and simpler than FVM to represent a complex computational domain. The disadvantages of both FDM and FVM include the impossibility of calculating the sought variables at points that are not grid nodes or centers of control volumes. FEM were originally intended for static calculations of building structures [7][8][9]. They are based on dividing the computational domain into a sufficiently large number of finite elements of simple shape, as a rule, polyhedra. Nodes are highlighted on each element. First of all, these are the vertices of polyhedra; however, it is possible to choose other points as nodes. For each element, for all functions sought from the system of differential equations, linear combinations of predetermined basis functions are sought, connecting the spatial coordinates and time with the corresponding sought variable. The set of such combinations for all elements must meet the following conditions: the minimum of the sum of the squared residuals for all nodes of all finite elements is achieved (the residuals are obtained by substituting the required derivatives of the corresponding linear combinations of basis functions in the differential equations); equality of the sought variables at the vertices of adjacent elements when determining them from linear combinations of the basis functions of these elements; equality of the calculated boundary conditions when determining them on the basis of the corresponding linear combinations of basis functions to the given boundary conditions. It should be noted that with a coordinated selection of the number of finite elements, the number of nodes in the elements and the number of basis functions, it is possible to achieve that the residuals at the nodes of elements, subject to the specified conditions, will be equal to zero. These conditions generate a system of algebraic equations, the solution of which gives linear combinations of basis functions that allow one to determine the desired variables at any point in the computational domain, which is an undoubted advantage of the FEM. It should be noted that if the initial system of PDE is linear, then the systems of algebraic equations to which the approximate solution of the system of PDE is reduced will be linear. When solving non-stationary problems using FDM, FVM, and FEM, the resulting systems of algebraic equations become extremely large and decomposition methods are used to solve them, which usually consist in dividing the solution in spatial coordinates and in time. A subsystem of equations related to one moment in time is distinguished. After solving it, the partial derivatives of the required quantities with respect to time are found. Using these derivatives, the values of the corresponding quantities are determined at the next time instant (at the next time layer). In this case, various explicit and implicit difference schemes are used [10]. In all the considered methods, the condition for a small deviation of the approximate solution of system of PDE from its exact solution is the smallness of the values of the characteristic geometric dimensions (grid steps, maximum dimensions of control volumes and finite elements). The most reasonable numerical criterion for such a deviation (the quality of an approximate solution) is the value of the residual maximum in absolute value at all considered (control) points of the computational domain; however, none of the considered methods uses this criterion. Taking into account the indicated FDM, FVM and FEM defects, a more effective method for solving the system of PDE is proposed. It is based on the search for such values of the coefficients of linear expansions of the basis functions, which represent the dependences of the functions sought from the system of PDE on the spatial coordinates and time, at which the residual maximum in absolute value reaches the minimum value, determined among all residuals at the given control points of the computational domain. The transition from minimizing the sum of the squared residuals to minimizing the maximal residual in absolute value significantly improves the quality of the approximate solution and allows one to go from small finite elements to sufficiently large blocks, within each of which their own linear expansions of basis functions are sought. The method is based on the assignment within the computational domain of control points, at each of which residuals are determined. All control points of the computational domain are divided into three groups. The first group is the internal checkpoints of the blocks. At these points, only the residuals of the original differential equations are calculated, which are obtained after substituting the desired functions and their partial derivatives in them, determined from linear expansions of the basis functions. The second group is the points lying on the boundaries of the blocks. At these points, the residuals of the differential equations are calculated for each adjacent block using its linear expansions. In addition, discrepancies between the sought-for values, as well as those partial derivatives that are included in the differential equations, calculated using linear expansions of adjacent blocks, are determined. The third group includes control points that lie on the boundaries of the computational domain. At these points, the composition of the residuals is supplemented by residuals that determine the accuracy of the approximation of the obtained solution to the initial and boundary conditions. In particular, discrepancies between the specified values of quantities at the boundaries of the computational domain and the values of these quantities calculated from linear combinations of basis functions are determined. If the computational domain is divided into subdomains, each of which is described by its own system of differential equations, then the indicated subdomains are divided into blocks. At the points lying on the boundaries of adjacent subdomains, discrepancies between the values of those sought-for functions that are included in the systems of differential equations of both subdomains are determined. It should be noted that when minimizing the residual maximum in absolute value, it is necessary to compare residuals having different dimensions and different physical meanings. Therefore, it is advisable to make such a comparison between the relative residuals obtained by dividing the absolute residuals by their maximum allowable values. If the initial system of PDE is linear, then the proposed method, which can be called the method of control points, is reduced to solving a linear programming problem [11][12].

Mathematical problem statement
There are N independent parameters N x x ,..., 1 . In the space of these parameters, the computational domain Q is specified. This domain is divided into L subdomains N Q Q ,..., 1 . Points lying on the boundary between two It should be noted that the k -th differential equation of the l -th subsystem of differential equations kl D in the general case includes not all l K of the sought functions, not all N K l  first derivatives of the sought-for functions with respect to independent parameters, not all Similar sets of numbers and indexes in these sets can be introduced for derivatives of higher orders. Taking into account these definitions, the system of differential equations can be written as: of residuals corresponding to differential equations are determined in it, residuals corresponding to the boundary conditions of the indicated types are determined at some points lying on the boundaries between blocks, between subdomains and on the outer boundaries of subdomains. As already indicated, the residuals depend on the coordinates of the corresponding control point and the expansion coefficients of the basis functions. As shown earlier, the boundary conditions are represented as equalities. Each such equality is associated with the residual obtained after substituting into it linear expansions of the basis functions and subtractions from the left side of the equality of the right side. Such residuals are determined for all control points lying on the boundaries between blocks, between subareas and on the outer boundaries of subareas. If the systems of differential equations are linear, then the residuals will be linear functions of the expansion coefficients. In general case:  The total number of residuals in the problem is determined by the number of control points, their distribution into internal and boundary, subsystems of differential equations and boundary conditions. The number of parameters to be optimized is also determined by the number of blocks and degrees of polynomials, with which the sought functions are described. It should be noted that using the polynomials found as a result of optimization, which approximate the sought functions, one can easily calculate in a large number of check points that do not coincide with the control points. If such a calculation shows unsatisfactory residual values, then additional control points can be selected from the check points with the maximum residual values, which are included in the original problem.
If the optimal value of the objective function is greater than 1 , i.e. The required accuracy of the solution of the problem is not ensured, then the number of blocks into which the subdomains of the problem or the degrees of polynomials are divided into which the required functions are approximated should be increased.

Heat exchanger mathematical model
As an example, consider the calculation of a periodic ceramic heat exchanger [13]. Such a heat exchanger is a cylindrical volume filled with a ball backfill. Consisting of ceramic balls of the same small radius R . Heating gas (products of combustion of ceramic fuel) is periodically supplied to the heat exchanger, which heats up the ball charge and heated air, which cools the ball charge. In this case, the movement of gas and air flows occurs along the axis of the cylinder in opposite directions. The heat exchanger is intermittent. First, the stage of heating the ball charge follows, when gas is supplied to the heat exchanger. This is followed by the stage of cooling the ball filling, when air is supplied to the heat exchanger. In the calculation, a simplifying assumption is made, which consists in the fact that the heat transfer between all balls located in one section perpendicular to the cylinder axis proceeds in the same way. Let introduce the following notation: , the temperature of hot gases, air, ceramics when it is heated, ceramics when it cools, respectively. At the stage of ceramic heating, the differential equation for the heat balance of the gas has the following form: Heat transfer between gas and ceramic is described by the following differential equation The differential equation for heat propagation in ceramics takes into account the heat propagation inside the ceramic ball. Radial heat propagation in a uniform ball of radius R is considered. It is assumed that at any time moment t the temperature at points located at the same distance r from the center of the ball will be the same. This means that the temperature in each ball depends only on r and on t . If spherical coordinates are introduced for each ball, then the process of heat propagation is described by the differential Fourier equation [14].
At the stage of ceramic cooling, the differential equation for the heat balance of the gas has the following form Heat transfer between gas and ceramics during ceramic cooling is described by the following differential equation.
And for the propagation of heat inside the ceramic backfill, the differential Fourier equation will take the following form: The model takes into account the boundary conditions (reversibility conditions), requiring that the temperature of the ceramics at all control points at the beginning of the heating stage be equal to the temperature of the ceramics at the corresponding control points of the end of the cooling stage and the temperature of the ceramics at the control points at the end of the heating stage equal to the temperature of the ceramics in the control points at the beginning of the cooling stage. Since the coefficients are functions of the temperatures of gas, air and ceramics calculated as a result of solving this problem, the calculation of the heat exchanger is carried out iteratively. At the first iteration, approximate values are set, after which problem (10) -(14) is solved. As a result of the solution, we obtain the temperatures. The coefficients are refined from the obtained temperatures, after which problem (10) -(14) is solved again. The refinement process goes on until the temperature difference obtained at adjacent iterations is less than a certain specified value.

Calculation results
The calculation of the heat exchanger is carried out with the following initial data. The radius of the cylindrical heat exchanger 2 m. The height of the heat exchanger will be 4 m. The duration of the heating stage and the duration of the cooling stage 300 sec. The radius of the backfill ball: 01 . 0 m. Throughout the entire heating stage, the gas temperature at the inlet to the heat exchanger is 1800 K. Throughout the entire cooling stage, the air temperature is 778 K. Gas and air consumption is 74 . 8 kg/sec and 427 . 11 kg/sec, respectively. During the calculations, 2 subdomains were identified subregion of the heating stage and the subregion of the cooling stage. Each sub-area was divided into 36 blocks using time division into 6 equal intervals and the length of the heat exchanger into 6 equal intervals.
The functions were approximated by 3rd degree polynomials. For gas temperatures, a polynomial is in two variables (time and length), and the temperature of ceramics is in 3 variables (Length, time, and radius). There were 10 control points along the radius of the ball. The total number of control points at the heating stage and at the cooling stage is 45760 for each of the subregions. Total residuals: 91520. The total number of constraints in the problem is: 183040. The number of parameters to be optimized is: 2161. Each of the polynomials included in 36 blocks of the subdomain is described by 10 (for gas or air) and 20 (for ceramics) parameters. The additional parameter Z is added to the total number of parameters. Taking into account the need to clarify the coefficients included in the differential equations and depending on the temperatures of gas, air and ceramics, 3 iterations of solving the linear programming problem were carried out. Optimal task value: 0.793. Calculation of residuals at check points (which were 3 times more than control points) showed that the maximum residual is: 0.832. Fig. 1 -Fig. 4 show the calculated values of the temperatures of gas, air and ceramics for some points of the heat exchanger.   It should be noted that the problem being solved is very difficult for the FDM, FVM and FEM methods. It is solved as follows. The temperature of the ceramic at the initial heating stage is set. Then, the processes of heating and cooling the ceramics are simulated many times until the temperature at the beginning of the heating stage stabilizes.
For ceramic values of thermal diffusivity coefficients, the considered system of differential equations is "rigid". This requires the specification of small characteristic sizes in spatial coordinates and in time, which leads to a significant consumption of computing resources. The use of this approach significantly reduces the resources required to solve the problem.

Conclusions
The analysis of existing methods for solving systems of partial differential equations, including their advantages and disadvantages. Based on the analysis of disadvantages of the existing methods, an effective numerical method for solving these problems (the method of control points) is proposed. Method is based on reducing the solution of systems of PDE to solving linear programming problems. The method is illustrated by calculation periodic heat exchanger with spherical filling.