The Method of Effective Numerical Solution of the System of Equations of Thermal Conductivity Oleg Kurchenko, Kseniia Dukhnovska, Oksana Kovtun, Anastasiia Nikolaienko and Iryna Yurchuk Taras Shevchenko National University of Kyiv, Bohdan Hawrylyshyn str. 24, Kyiv, UA-04116, Ukraine Abstract The article considers the method of effective numerical solution of the system of equations of thermal conductivity, which is used for numerical modeling of the forecast of thermophysiological state of man, and its parallel software implementation on a multiprocessor system. The method is suitable for building a three-dimensional model for predicting the thermophysiological state of man. The results of numerical experiments and parallel calculations show the good suitability of the proposed algorithms for parallel implementation and their significant potential for creating on their basis applied computer systems for modeling the prediction of human thermophysiological state. Keywords 1 System of equations of thermal conductivity, parallel algorithm of the Gaussian method, the method of finite differences 1. 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. 2. 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 Information Technology and Implementation (IT&I-2023), November 20-21, 2023, Kyiv, Ukraine EMAIL: oleg.kurchenko@knu.ua (A.1); kseniia.dukhnovska@knu.ua (A.2); kovok@ukr.net (A.3); n_nastja@ukr.net (A.4); i.a.yurchuk@gmail.com (A.5) ORCID: 0000-0002-3507-2392 (O. Kurchenko); 0000-0002-4539-159X (K. Dukhnovska); 0000-0003-0871-5097 (O. Kovtun); 0000-0002- 2402-2947 (A. Nikolaenko); 0000-0001-8206-3395 (I. Yurchuk); ©️ 2023 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) CEUR ceur-ws.org Workshop ISSN 1613-0073 Proceedings 206 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]. 3. 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. 4. 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-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 207 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 c – specific heat, kcal / (kg · ° C); m – mass, kg; T – temperature, ° C; t – time; M – metabolism, kcal; Q – heat flux, kcal; V – volume, l; W – cardiac output, l / h. Indices: * – initial value; b – blood; ph – physical activity; RS – heat loss with respiration; sh – cold tremor of skeletal muscles; K – conduction; air – air; water – water. 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 = T0ij, 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 (2) 𝜆 =𝜆 𝑖𝑗 𝜕𝑥 𝑖,𝑗−1 𝜕𝑥 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘, 𝜕𝑇𝑏 𝜕𝑇𝑖,𝑗 (3) 𝜆𝑏 𝜕𝑥 = 𝜆𝑖,𝑗 𝜕𝑥 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘. 𝑖𝑗 If 𝑦 ∈ 𝐺𝑦 : 𝜕𝑇𝑖𝑗 𝜕𝑇𝑖,𝑗−1 (4) 𝜆𝑖𝑗 𝜕𝑦 = 𝜆𝑖,𝑗−1 𝜕𝑦 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘, 𝜕𝑇 𝜕𝑇𝑖,𝑗 (5) 𝜆𝑏 𝜕𝑦𝑏 = 𝜆𝑖,𝑗 𝜕𝑦 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘. 𝑖𝑗 If 𝑧 ∈ 𝐺𝑧 : 𝜕𝑇𝑖𝑗 𝜕𝑇𝑖,𝑗−1 (6) 𝜆𝑖𝑗 𝜕𝑧 = 𝜆𝑖,𝑗−1 𝜕𝑧 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘, 𝜕𝑇 𝜕𝑇𝑖,𝑗 (7) 𝜆𝑏 𝜕𝑧𝑏 = 𝜆𝑖,𝑗 𝜕𝑧 , ∀𝑡 ≥ 0, 𝑗 ≠ 𝑠𝑘. For j = sk, 𝑥 ∈ 𝐺𝑥 𝑖,𝑠𝑘 : 𝜕𝑇𝑖,𝑠𝑘 (8) 𝜕𝑥 = 𝑎𝑖𝑗 𝜆𝑖,𝑠𝑘 (𝑇𝑖,𝑠𝑘 − 𝑇 𝑎𝑖𝑟 ), ∀𝑡 ≥ 0, 𝜕𝑇 𝜕𝑇𝑖,𝑠𝑘 (9) 𝜆𝑏 𝜕𝑥𝑏 = 𝜆𝑖,𝑠𝑘 𝜕𝑥 , ∀𝑡 ≥ 0. For j = sk, 𝑧 ∈ 𝐺𝑦 𝑖,𝑠𝑘 : 𝜕𝑇𝑖,𝑠𝑘 (10) 𝜕𝑦 = 𝑎𝑖𝑗 𝜆𝑖,𝑠𝑘 (𝑇𝑖,𝑠𝑘 − 𝑇 𝑎𝑖𝑟 ), ∀𝑡 ≥ 0, 208 𝜕𝑇𝑏 𝜕𝑇𝑖,𝑠𝑘 (11) 𝜆𝑏 = 𝜆𝑖,𝑠𝑘 , ∀𝑡 ≥ 0. 𝜕𝑦 𝜕𝑦 for j = sk, 𝑥 ∈ 𝐺𝑧𝑖,𝑠𝑘 : 𝜕𝑇𝑖,𝑠𝑘 (12) = 𝑎𝑖𝑗 𝜆𝑖,𝑠𝑘 (𝑇𝑖,𝑠𝑘 − 𝑇 𝑎𝑖𝑟 ), ∀𝑡 ≥ 0, 𝜕𝑧 𝜕𝑇 𝜕𝑇 (13) 𝜆𝑏 𝑏 = 𝜆𝑖,𝑠𝑘 𝑖,𝑠𝑘 , ∀𝑡 ≥ 0, 𝜕𝑧 𝜕𝑧 where Gij surface of ij-compartment, Gx, Gy, Gz - the projection of these surfaces on the coordinate axis [12]. 5. 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. 6. 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. 7. 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 209 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 (16) 𝑥𝑛−1 = , 𝑎𝑛−1,𝑛−1 𝑛−1 (17) 𝑥𝑖 = (𝑏𝑖 − ∑ 𝑎𝑖𝑗 𝑥𝑗 )⁄𝑎𝑖𝑖 , 𝑗=𝑖+1 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