Introduction

Boundary value problems, which reflect technological or natural processes, have a high computational complexity. To solve them, various analytical or numerical approaches are used, neglecting the accuracy of the solution. Such boundary value problems can also include problems with heat conduction equations.

Solving this problem requires significant CPU time to perform calculations. This is especially true for modeling tasks, where the requirements for quality and detail of three-dimensional models that describe the processes of heat transfer, which, in turn, leads to increased computational costs.

The use of parallel computing is extremely important for obtaining an efficient and accurate forecast, as it provides the necessary acceleration of algorithm calculations when performing them on powerful multiprocessor systems. However, the development of effective parallel algorithms for multiprocessor systems is a difficult scientific and technical problem and requires solving interdisciplinary research using knowledge of mathematical, architectural and software models of parallel computing.

Literature review and problem statement

Equations of thermal conductivity are quite widely used. In work [1], the flow of viscous dissipation, chemical reactions, and convection heat on the walls of the channel is investigated with the help of heat conduction equations. The Math software package was used to solve the nonlinear system of differential equations. The results of this work can be used to solve the cooling problem in various industries.

The work [2] presents the results of nanofluid research in the heat engineering sector. The system of nonlinear differential equations that describes this problem is first converted into a system of ordinary differential equations. Then it is solved using the bvp4c algorithm. Current work has many applications in various fields of engineering, biotechnology, nanotechnology and medicine.

The paper [3] examines the heat exchange of the boundary layer fluid. Similarity transformations are again applied to transform the governing partial differential equations for mass, momentum, and energy into a system of ordinary differential equations. The MATLAB package is used to obtain the solution. The effect of variable thermal conductivity on an exponentially stretched surface is investigated in [4]. The bvp4c algorithm is used to solve this boundary value problem.

The thermal conductivity equation is used to predict the thermo-physiological state of a person. The purpose of the article [5] was to review existing thermophysiological models for the entire human body and its individual segments. The paper provides an overview of the most recognized thermal models, such as Fiala, Berkeley Comfort Model, Tanabe and ThermoSem.

The study [6] proposed clothes with fluid circulation, which can effectively reduce the thermal stress of people in a hot environment. This work used a model of human thermoregulation together with a model of warmth and comfort. Thus, in the article [7], a complete virtual three-dimensional model of human thermophysiological reactions is investigated. The numerical solution of the model was developed using the ANSYS Fluent platform. This platform works due to parallel processes [8].

The purpose and objectives of the research

The purpose of the research is to develop method of effective numerical solution of the system of equations of thermal conductivity. This will make it possible to obtain a numerical value for problems that include the heat conduction equation. To achieve the goal, the following tasks were solved:

 to solve the heat conduction equation, apply the finite difference method together with the parallel Gaussian algorithm;

 to investigate the effectiveness of the method.

As can be seen from the review, the boundary value problems, which include the heat conduction equation, are first analytically translated into the Cauchy problem, and then calculated. This is due to the high computational complexity of boundary value problems. Let's consider one of these.

Mathematical model for predicting the thermophysiological state of human

A multi-compartmental model of human thermoregulation and heat transfer is used for mathematical modeling of human physiological condition prediction [5][6][7][8][9][10][11].

According to the definition given in [9], a compartment is a certain amount of a substance that is released in a biological system and has the property of unity, so in the processes of transport and chemical transformations it can be considered as a whole. Models in which the studied system is presented as a set of compartments are called multi-compartmental.

The class of multicompartmental dynamic models of human heat exchange and thermoregulation describes heat production in organs and tissues, heat transfer by blood flow, conduction, convection, radiation, evaporation from the skin and upper respiratory tract, afferent and efferent thermoregulatory processes.

The structure of the models is based on the compartment. The compartment has an energy source, the ability to receive and transfer heat to neighboring compartments and exchange energy with the environment. The set of compartments and the task of connections between them determines the nature and degree of approximation of the human body and the physiological processes occurring in it. For the sample of a single compartment the cylinder for representation of all parts of a body and sphere -for the head is taken. Each compartment corresponds to a specific organ, body part or body tissue. The models take into account human anatomical and physiological parameters: weight, height, body surface area, biophysical characteristics of tissues and organs, metabolism, oxygen consumption, blood flow, pulmonary ventilation, heat transfer and heat transfer coefficients, cardiac output, sweating, fluid loss, etc. [11]. Approximation of the human body can be performed on any number of geometric parts, for this model the human body was approximated by 13 cylinders and sphere: head, torso, left and right arm (shoulders, forearms, hands), left and right leg (thighs, legs, feet). Each cylinder that approximates a part of the human body has nested cylinders that describe the properties of other tissues of the human body. The heat balance equation is written to each compartment and the result is the following system of differential equations in partial derivatives:

(1)

where cspecific heat, kcal / (kg • ° C); mmass, kg; Ttemperature, ° C; ttime; Mmetabolism, kcal; Qheat flux, kcal; Vvolume, l; Wcardiac output, l / h. Indices: *initial value; bblood; phphysical activity; RSheat loss with respiration; shcold tremor of skeletal muscles; Kconduction; airair; waterwater. To solve this system requires additional conditions: initial and boundary. Initial conditions. At the beginning of time, temperatures, metabolism, blood flow, etc. are set for each compartment: if t = 0: Tij = T 0 ij, for each i,j when 0≤i≤N-1, 0≤j≤L-1. Boundary conditions. To solve this system, boundary conditions of the fourth kind are used. For neighboring compartments, conditions must be met at the boundary between them.

If 𝑥 ∈ 𝐺 𝑥 𝑖𝑗 :

𝜆 𝑖𝑗 𝜕𝑇 𝑖𝑗 𝜕𝑥 = 𝜆 𝑖,𝑗−1 𝜕𝑇 𝑖,𝑗−1 𝜕𝑥 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘,(2)𝜆 𝑏 𝜕𝑇 𝑏 𝜕𝑥 = 𝜆 𝑖,𝑗 𝜕𝑇 𝑖,𝑗

𝜕𝑥 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘.

(3)

If 𝑦 ∈ 𝐺 𝑦 𝑖𝑗 : 𝜆 𝑖𝑗 𝜕𝑇 𝑖𝑗 𝜕𝑦 = 𝜆 𝑖,𝑗−1 𝜕𝑇 𝑖,𝑗−1 𝜕𝑦 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘,(4)𝜆 𝑏 𝜕𝑇 𝑏 𝜕𝑦 = 𝜆 𝑖,𝑗 𝜕𝑇 𝑖,𝑗

𝜕𝑦 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘.

(

)5

If 𝑧 ∈ 𝐺 𝑧 𝑖𝑗 :

𝜆 𝑖𝑗 𝜕𝑇 𝑖𝑗 𝜕𝑧 = 𝜆 𝑖,𝑗−1 𝜕𝑇 𝑖,𝑗−1 𝜕𝑧 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘,(6)𝜆 𝑏 𝜕𝑇 𝑏 𝜕𝑧 = 𝜆 𝑖,𝑗 𝜕𝑇 𝑖,𝑗 𝜕𝑧 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘.(7)

For j = sk, 𝑥 ∈ 𝐺 𝑥 𝑖,𝑠𝑘 :

𝜕𝑇 𝑖,𝑠𝑘 𝜕𝑥 = 𝑎 𝑖𝑗 𝜆 𝑖,𝑠𝑘 (𝑇 𝑖,𝑠𝑘 − 𝑇 𝑎𝑖𝑟 ), ∀𝑡 ≥ 0,(8)

where G ij surface of ij-compartment, Gx, Gy, Gz -the projection of these surfaces on the coordinate axis [12].

Numerical method of finite differences for solving a system of differential parabolic equations

To solve the model (1), the finite difference method was used. This is a method of solving differential equations based on the replacement of differential operators with their approximate values at individual points. When using the finite difference method, the differential problem is replaced by the difference problem.

The equations given in problem (1) are parabolic differential equations in partial derivatives. It is proved that finite difference methods can be used to solve parabolic differential equations. To solve the system of equations, the method of finite differences with an implicit scheme is taken. For parabolic differential equations, the solution of such a scheme is stable under any conditions.

To solve the differential equation by the finite difference method, first the area on which the solution is sought is replaced by a discrete set of points (nodes of the grid). The difference operators corresponding to the differential equation are written in the internal nodes of the grid. The difference operators corresponding to the boundary conditions are written in the boundary nodes of the grid. As a result, we obtain a system of algebraic equations, the number of which is proportional to the number of internal nodes of the grid region. To obtain a numerical solution, we need to solve this system of equations.

Numerical method of finite differences for solving a system of differential parabolic equations

To solve the model (1), the finite difference method was used. This is a method of solving differential equations based on the replacement of differential operators with their approximate values at individual points. When using the finite difference method, the differential problem is replaced by the difference problem. The equations given in problem (1) are parabolic differential equations in partial derivatives. It is proved that finite difference methods can be used to solve parabolic differential equations. To solve the system of equations, the method of finite differences with an implicit scheme is taken. For parabolic differential equations, the solution of such a scheme is stable under any conditions. To solve the differential equation by the finite difference method, first the area on which the solution is sought is replaced by a discrete set of points (nodes of the grid). The difference operators corresponding to the differential equation are written in the internal nodes of the grid. The difference operators corresponding to the boundary conditions are written in the boundary nodes of the grid. As a result, we obtain a system of algebraic equations, the number of which is proportional to the number of internal nodes of the grid region. To obtain a numerical solution, we need to solve this system of equations.

Parallel algorithm of the Gaussian method

To find the temperatures at each step of integration, you need to solve a system of linear algebraic equations. The Gaussian method was used for this purpose.

The Gaussian method is based on the possibility of performing transformations of linear equations, which do not change the solution of the considered system. The Gaussian method involves the sequential execution of two steps. At the first step -the direct course of the Gaussian method -the initial system of linear equations by successive exclusion of unknowns is reduced to the upper triangular form. In the reverse course of the Gaussian method (the second stage of the algorithm) the values of the unknowns are determined. The calculations performed on the elements of the matrix A and the vector b are determined by the following relations:

𝑎 𝑘𝑗 ′ = 𝑎 𝑘𝑗 − 𝑎 𝑘𝑗 𝑎 𝑖𝑖 * 𝑎 𝑖𝑗 ,(14)𝑏 𝑘 ′ = 𝑏 𝑘 − 𝑎 𝑘𝑗 𝑎 𝑖𝑖 * 𝑏 𝑖 ,(15)

where 𝑖 ≤ 𝑗 ≤ 𝑛 − 1, 𝑖 ≤ 𝑘 ≤ 𝑛 − 1, 0 ≤ 𝑖 ≤ 𝑛 − 1.

After bringing the matrix of coefficients to the upper triangular form, it becomes possible to determine the values of the unknown. From the last equation of the transformed system the value of the variable xn-1 can be calculated, then from the penultimate equation it becomes possible to determine the variable xn-2, etc. In general, the calculations of the reverse Gaussian method can be represented by the relations:

𝑥 𝑛−1 = 𝑏 𝑛−1 𝑎 𝑛−1,𝑛−1 ,(16)𝑥 𝑖 = (𝑏 𝑖 − ∑ 𝑎 𝑖𝑗 𝑥 𝑗 𝑛−1 𝑗=𝑖+1 ) 𝑎 𝑖𝑖 ⁄ ,(17)

where i = n-2,…,0. Multicore computers are now commonplace. Multiprocessor machines have almost completely supplanted single-core counterparts and are used everywhere.

Many users do not like working with braking applications. Also annoying is when running a task in an application reduces the performance or sensitivity of another part of the program. Thus, the requirements for improving the efficiency of software are growing every day, which in turn necessitates mechanisms that allow rational use of available resources of modern equipment.

The solution of problem (1) is reduced to the solution of a system of algebraic systems. The dimension of such a system is equal to the number of integration steps in variable directions. The computational complexity of the Gaussian algorithm is of the order of O(n 3 ). In order to speed up the algorithm, it was decided to use a parallel calculation of the Gaussian algorithm.

Since the solution by the Gaussian method is reduced to a sequence of similar computational operations of multiplication and addition over the rows of the matrix, as subtasks can be taken calculations related to the processing of one or more rows of the matrix A and the corresponding element of the vector b. Each iteration associated with the solution of the next subtask begins with the selection of the leading line. The row with the largest absolute value among the elements of the i-th column corresponding to the excluded variable xi is searched.

Because the rows of the matrix A are assigned to different subtasks, to find the maximum value in the subtask column numbered k, the elements must be exchanged for the variable xi, which is excluded. After collecting all these coefficients, it can be determined which of the subtasks contains the leading line and which value is the leading element.To continue the calculations, the leading subtask must send its row to the matrix A, which determines the coefficients of the system of linear equations and the corresponding element of the vector b to all other subtasks with numbers k, xki.

When performing the reverse of the Gaussian method, the subtasks perform the necessary calculations to find the values of the unknowns. As soon as any subtask i, 1≤i≤n, which determines the value of its variable xi, this value must be used by all subtasks with numbers k, k <i: subtasks substitute the obtained value of the new unknown and perform value adjustments for elements of the vector b. Selected basic subtasks are characterized by the same computational complexity. However, the size of the matrix that describes the system of linear equations is significantly larger than the number of threads in the program (p << n), and the basic subtasks can be enlarged by combining several rows of the matrix in one subtask. The use of a sequential data separation scheme for the parallel solution of systems of linear equations will lead to an uneven computational load between the flows: as the exclusion (forward) or determination (reverse) unknown in the Gaussian method for most flows all necessary calculations will be completed and they will stand still. A possible way to solve the problem of balancing calculations may be to use a tape loop scheme for the distribution of data between the enlarged subtasks: in this case, the matrix A is divided into sets (strips) of rows of the form (Figure 2):

where where n is the number of unknowns in the system of linear equations; p is the number of processes that the operating system allows to create a system of linear equations to solve. Comparing the data division scheme and the order of calculations in the Gaussian method, it can be noted that the use of a cyclic method of forming tapes allows to better balance the computational load between subtasks. Therefore, after analyzing the sequential version of the Gaussian algorithm, we can conclude that parallelization is possible for the following computational procedures:

 search for the leading line;  subtraction of the leading line from all lines to be processed;  execution of the return stroke. This approach significantly speeds up data processing than sequential computation.

The study materials and methods

The characteristics of an average person were used for calculations. The work [17] was the source of data. It describes the thermophysical characteristics of a person with a body weight of 74.1 kg and a body surface area of 1.89 m 2 . In order to obtain these average characteristics, a study was conducted on 10 men and 10 women (Table 1). The work [17] also shows that 25 % of the tissues in the human body are skin and fat, 11 % are internal organs, 3 % are nervous tissue, 43 % are muscles, and 18 % are the skeleton. At the same time, the specific heat capacity of the skeleton is 0.5 kcal/g•°С, the heat capacity for fat is 0.6 kcal/g•°С, and 0.9 kcal/g•°С for all other tissues. Thermal conductivity is given in Table 2. Calculations were performed on the language platform Visual C++ using the library MPI. The choice As can be seen from the first diagram (Figure 2) at a dimension of the matrix of 50, the operating time of the parallel Gaussian algorithm is greater than the operating time of the sequential algorithm. This is because the parallelization process is too expensive because the process itself takes time. But with increasing dimension of the matrix (Figure 3 -Figure 4), the parallel algorithm justifies itself and works three times faster, as seen in the last diagram.

According to the table, diagrams are constructed for better visualization (Figure 2 -Figure 4).