Mathematical modeling of heat transfer in anisotropic biophysical materials, taking into account the phase transition boundary [0000-0003-4866-2575] [0000-0002-6767-104х] Yaroslav Sokolovskyy1 , Iryna Boretska2 , 3[0000-0002-1008-8427] 4[0000-1111-2222-3333] Bogdana Gayvas , Igor Kroshnyy [0000-0002-7100-3263] 6[0000-0002-8156-8962] Volodymyr Shymanskyi5 and Michal Gregus 1, 2, 4, 5 Ukrainian National Forestry University, Lviv, Ukraine 1 sokolowskyyyar@yahoo.com 2 iryna.boretska@gmail.com 4 kroshny.igor@gmail.com 5 vshymanskiy@gmail.com 3 Center of Mathematical Modeling, Institute of Applied Problems of Mechanics and Mathematics them Y. S. Pidstryhach, Lviv, Ukraine, gaivas.bogdana@gmail.com 6 Comenius University in Bratislava, Bratislava, Slovakia michal.gregus@fm.uniba.sk Abstract. A two-dimensional mathematical model of heat transfer of anisotropic biophysical materials is constructed, taking into account the motion of the boundary of phase transitions. The influence of the main components and the orientation of the main axes of the heat transfer tensor on non-stationary temperature fields in the biophysical material is determined, taking into account the motion of the boundaries of phase transitions. An analytical-numerical method has been developed for determining heat transfer in a biophysical material with a moving boundary of phase transitions, and the moving boundaries of a phase transition in a rectangular region have been established, taking into account the main axes of anisotropy. Algorithms of a nonlinear mathematical model are constructed under variable temperature conditions of the medium. The integrals along the phase transition boundary are determined numerically. All other values included in this equation are calculated according to the physical and thermal characteristics of a particular material. Keywords: phase transitions, mathematical model, heat transfer, analytical- numerical method, anisotropy, biophysical material. 1 Mathematical model of heat transfer in a biophysical material of rectangular cross section The purpose of this study is to construct a two-dimensional mathematical model of heat transfer in biophysical materials of rectangular cross section 2L1 ,2L2  Copyright © 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0) 2019 IDDM Workshops. (L1  x1  L1 ,  L2  x2  L2 , ) taking into account the phase transition boundary. The outer contour of such material in variables x1 , x2 is described by the surface equation   F0 x1 , x2   x12  L12 x22  L22  0 .  (1) In the process of heat exchange of biophysical material with the environment, a dried zone is formed, extending from the outer surface to the depth of the body. We assume that the dried and moisture zones of the biophysical material are separated by a cylindrical surface, the generatrices of which are parallel to the axis of the material [1]. The equation of such a surface is represented in the form Fm x, y,   F0 x1 , x2      , (2) where    is unknown time function. We also assume that in the dried zone the moisture has been removed, and in the remaining volume it has been preserved. The moisture content remaining in the biophysical material is determined by the formula W   L V  Vm  / V  , where V is the body volume; Vm is the dried zone volume;  L is the moisture density. The temperature distribution T x1 , x2 ,  of the biophysical material in the dried zone is described by the equation:     (Cv v  Ca a )  1   Cs s  T    11 T     22 T   12  21  T  F x1, x2 ,  , 2 (3)  x1   x1  x2  x2  x1x2 The indices v, a, s denote the components of vapour, air, and skeleton, and , Cv , Ca , Cs ,  v ,  a ,  s are porosity, heat capacity, the density of vapour, air, skeleton of the material, respectively; ij are components of the thermal conductivity tensor; T is temperature; F x1 , x2 ,  is the summand characterizing the internal source. The main thermal conductivity coefficients are determined through the thermal conductivity coefficients of anisotropic biophysical material and a one-to-one coordinate transformation is established: x1  12 x1 /   22  2   x2 / ; x2  11  1 x1 /   21x2 / . (4)   11  1 2  12 2 ;   22  2 2  221 . The variables x1 , x2  coincide with the main directions of the anisotropy of the thermal conductivity of the plate. If we pass on to variables 1   / 1  x1 ,  2   / 2 1 / 2 x 2 , then in the context of 1/ 2 heat transfer we use the equation   2T T  2T    . (5)    * 2  2 2   1 It is important to obtain the boundary conditions on the surfaces of the anisotropic biophysical material in variables 11 ,  2 : T  H i  T  u t   0, (6)  i l1 where H i*  i /  H i , H i  ~i / i , ~i are heat exchange coefficients. m1l 2  m2 l1 From this surface, temperature reduction process moves inside the body. The boundary conditions of heat exchange (6) are given on the side surfaces. The surface that separates the dry and moist zones will have an oval cylindrical shape, and when completely dried, it is pulled into a line which is the axis of the bar. Given that the volume of the dried-up zone of the bar is a function of time, we present the equation of the boundary of the dried and moist zones in the form:      Fm 1 ,  2 ,   12  21 12  22  22  23  22  24      0 , (7) where   2122 23 24 ,    is so far unknown function of time  . We introduce the following values and variables: T 1 ,  2 ,   Tm / TП  Tm   ,    m cm a / m  a / am ,  *   a ,  m cm / m  1 / am , where TП , Tm is temperature on the contour of the biophysical material and at the phase transition boundary, respectively. Given the conditions of continuity of heat flow between the surfaces F0 and Fm, we obtain an expression for the value  :        12  21 12  22  22  23  22  24     /     , (8) The value of the so far unknown function of time    at the initial moment of the heat transfer process is zero  0  0 . After the process is completed, we have 1  0,  2  0,     1 . Thus, for the accepted notation, we get that   1 at  *     0 and   0 on the phase transition line Fm . For further studies, the equation of the contour of the cross section of the biophysical material and the line separating the dry and moist zones can be written as:   2    2  F0   12  1   22    22   3   42   0 ,  (9)   2    2  Fm   12   1   22    22   3   42       0 ,  (10)         where 1  21  22 / 2,  2  21  22 / 2,  3  23  24 / 2 ,  4  23  24 / 2 . From (10), we define explicitly the equation of the phase transition curve in the biophysical material:     2    3   42   12  21  12  22 .  (11) The signs "+" and "-" before the root refer to cases  2  0 and  2  0 , respectively. The sign "  " under the root corresponds to the case  22   3  0 . 2 Construction of an analytical-numerical method for imple- menting a mathematical model. Identification of moving boundaries of phase transition To construct an analytical-numerical method for implementing the mathematical model of anisotropic biophysical material with allowance for the moisture evaporation zone, the heat balance equation in the region bounded by the outer contour of the biophysical material and the contour of phase transition boundary plays an important role [6, 7]: F0 0 d     ds   dl   dl . (12) Fm 0 d * F0 0 n Fm 0 n If over time  , the volume of the dried zone of the material during heat transfer increases by V , then the amount of heat absorbed due to the phase transition is determined by the dependence Q   m cm T  Tm V Fm , F0  , (13) where V Fm , F0  is the volume of dried zone per unit length of the material, S П is the surface area of the cross-section of the specimen, SФ is the cross-sectiona area of the moist zone surrounded by the contour Fm =0; the m index denotes various previously identified technological characteristics of a biophysical material at the boundary of a phase transition. The given flow Q is defined by the flow on the surface of the phase transition  T  * Q  m   ds  . (14)  F 0 n   m  We pass on to the boundary at  *  0 . Then from expressions (13), (14) we obtain  T  V Fm , F0  m   ds    m c m T  Tm  , (15) n   *  Fm 0  On the basis of (12), (14), (15), we obtain the equation of the heat balance taking into account the moving boundary of the phase transition. F0 0 F0 0 d  V   ds   dl   * ; V Fm , F0     ds   ds   ds (16) d F0 n  * Fm Fm 0 S S This is the main equation for constructing an analytical-numerical method for implementing the mathematical model, taking into account the moving boundary of the moist and dried zones of biophysical material in conditions of heat transfer. To calculate the integrals in formula (16), it is necessary to have explicitly the equations of the line of the phase transition contour, and also to establish the boundaries of the corresponding integrals. In particular, given the equation of the boundary of moist and dried zones of biophysical material (7), we find    1 1,2   1   22    1   22  2122   ,   32   42      2 1,2    3   42            . 1 2 2 1 1 2 2 2 1 Since at the initial point in heat transfer time  0  0 , we get: 1 1,2   1   2 ,  2 1,2    3   4 . The double integrals in the heat balance equation (16) over the surface between the closed contour Fm and the outer contour F0 will be found as the difference between the integral over the surface of the full cross-section of the material and the integral over the surface SФ bounded by the contour Fm . We calculate the integral over the outer contour. So, we have       dl  2  dl   dl , (17) L0 n    L n L3 n  Thus, based on the above mathematical transformations, we obtain an expression for determining the first additive component of the integral over the outer contour (17).   2  l1 m2  l 2 m1 1  2 2  1       dl      4  1 1 2     1  2      L1  l1    1      3      2 4  (18) n l l 1  l      L 1 2 2 2 1  2 3  4 1 l1     1       l2  2   l2 2  4 2             L1  l1 1 1    3  L1  l1 1 1  12  1   22  d1 . 2        Similarly to the above considerations and mathematical transformations, we present the contour integral dl  dx1  in (17). Then we get    2   1 l1 m2  l2 m1 3  2   1   2     dl   2    4 2  2   3      L2  m 2    2  1    2   2 (19) L3 n m1 m2  4    1 m1  2          3  2 m2  1     2  2     4   L2  m2 2 2    1  L2  m2 2   1    m1  m1 1  4 2               22   3   42 d 2 .  2   Thus, based on the performed mathematical calculations, we obtain the expression for finding the second integral over the contour of the biophysical material in formula (17) to calculate the integral over the outer contour. Thus, the complete closed-contour integral L0 is calculated by formula (17) by substituting relations (18) and (19). Now let us pass on to finding the volume of the dried zone of the biophysical material, which is assigned to a unit of length, and which is between the planes F0  0, Fm  0 . This volume is determined by the formulas: F0 0 V F0 , Fm     dx1 dx 2  (20) Fm 0                    1/ 2    3 2 2 2 2 2 2 2 2 2 1 2  1 1 1 2 4 1 1 1 2  4  d1   0    1 12  21 12  22 4   where   1   22    /  32   42  1   22  2122 ,   1    3   42    / 12  21 12  22  .      Let us determine the derivative of the volume of the dried zone of the biophysical material in time, taking into account the time dependence of the magnitude   * and   the time dependence of the upper boundary of the integral. If the integration boundaries are functions of  , then applying the rule of differentiation of a complex function from several variables, we define  V 1  2 d          d  0       *   3  1 2  2  1 1 2  2 2    2 2 4 1  2  1 1 2  2 2   1 (21)  d1     4 1  1 1   2    4 12  21 12  22 2 2 2 2 2            1/ 2    2  2  2  2   2  2  2  2  2      3    1 2 4 1 2  21 22 .   1  4  2  2  2  2 2        1 2 2 2 2 2 1 2  22   21 22   We represent the derivative of the volume of the dried zone in the following way V 1 2 d  J V  J      , (22)  *   d where J   is an integral in (22), A  is the second additive component in (23),     dd   is functional dependence determined by the second additive *   f   * , * * component (22). Through J V – integral which is the first additive component (22). Here, time functions are the limit of integration    12   22  2122 .   The time derivative of   * takes the form 1   d  ' *   2  12   22  2122   2122 / 2  22  2122  * . (23)     d Further, according to (16), it is necessary to determine the integral over the region SФ bounded by the phase transition line. After cumbersome transformations we get   2   ds  4 1   (24) SФ  *    2   5     2 2  2 2  2 2 2 2 2   2    1    3 1  1 1   2   4 1  1 1   2               1   5     1  1 1   2 2 2 2 2 4     3 2           1  3   3 12  21 12  22   42 12  21 12  22    12  21 12  22 2 4 3   1 3  2324   3    12  21 12  22     42 12  21 12  22       12  21 12  22 4 d1 ,  2    where I  1 /  2 /  is transition Jacobian from variables x1 , x2 to 1 ,  2 .  We define the integral  * ds over the outer surface of the contour of the S  biophysical material   * ds  1 /  2 /     2 S 2 2 2 2 2  2 2 2   1  1 1  1  2   3  2   4 d 2 d1  (25)   S П  П  1       2  L1 l1 1 1    1 2              d   1 l 2 2  2    12  21 12  22  2 2 2 2     2 3 2 4 2  L1    1 l1 0   L1 1      L2 m1 1 1  2                 d  1 l1 m2   1 2 2 1 1 2 2 2  2 2 2 3 2 2 2 4 2  1 0 1 λ  λ    λ2  L2  m1 1 ξ1  λ          Δ2 m2     ξ1  Δ1 ξ1  Δ2 2 2 2 2  ξ 2  Δ3 ξ 2  Δ4 dξ 2 dξ1   2 2 2 2 λ L1 1 λ  λ    L1 l1 1 ξ1  λ1 l1 l2 λ2  λ    1 2   2  ( J1  J 2  J 3 )    2 Thus, all the components that are included in the heat balance equation (16) are determined taking into account the moving phase transition boundary. The analysis of the dependences (18) - (25) shows that the integrals along the boundary of the phase transition (21), (25) need to be calculated numerically. To obtain a relation for determining the moving boundary of a phase transition, we introduce some notation and carry out the following transformations    dl  c1 J L ,  dl  c 2 J L3 , (26) L n L3 n l m l m l m l m where c1   1 2 2 1 ; c2   1 2 2 1 .  l1l 2  m1m2 We write the expressions for J   and A  in formula (22) in the following form   J V      (27)   0  3 12  21 12  22   42 12  21 12  22      1  d1 ;   42 12  21   4        12  22 1 2 2 1 1 2 2 2                      1/ 2 2 2 2 2 2 2 2 2     3 1 2 2 4 1 2 A    . 2 2 4                    1 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 2 1 2 2 1 2 According to (25), the integral over the outer surface of the contour is written in the form     ds  I J S П , (28) S П   2 * where J S П  J1  J 2  J 3 . The integral (24) over the region Sф bounded by the phase transition line is represented as:    ds  I  , JФ (29) SФ   2 * where         5/ 2      3 12  21 12  22   42 12  21 12  22      1 JФ         5    2 14 1   1  1   2 2 2 2         2 3                     3/ 2   3  2 2 2 2 2 2 2 2 2 3  1 1 1 2 4 1 1 1 2                      2 2 2 2 14 2 2 2 2 2 2  1 1 1 2 3 4 3 1 1 1 2       d 1/ 2   42 12  21 12  22    34  12  21 12  22  1 Taking into account the above relations (27) - (28), as well as the equation of heat balance (16) taking into account the moving boundary of the phase transition, we obtain F0 0 d V 1 d   ds    I 2 J 1  J 2  J 3   Fm 0 d *  *  d * (30) d d J Ф    *I J V    A  1 I 2  d * d * According to (28), (29), the integral over the outer contour in the balance equation (17) can be written as   l m l m 1 l m l m 1   dl  2  1 2 2 1 JL  1 2 2 1 J L3  . (31) F0 n  l1l 2  m1m2   So, taking into account (16), (28), (29), we get d I 1 J1  J 2  J 3   I 1 2 d* J Ф   2 d *  d (32)    I   d J    A    c J  c J . * d * V 2 1 L 2 L3 Let us assume that J1  J 2  J 3   J S then the equation (32), taking into П account phase transition boundary      1 , takes the form *    I J sП  J    c1 J L  c2 J L3   2 I J V  A . (33)  * Note that the value J Ф  , A , J V   is a function of  . Their boundaries of the phase transition need to be calculated numerically [4, 5]. All other values included in equation (32) are calculated by the given physical and thermal characteristics of a particular material, namely, the cross-sectional dimensions, the main coefficients of thermal conductivity, the main directions, which are calculated by the coefficients of thermal conductivity of a particular anisotropic biophysical material, by transition Jacobian [11, 12] when the equations of thermal conductivity are converted to the canonical form. A two-dimensional mathematical model of heat transfer for non-stationary modes under conditions of heat transfer of anisotropic biophysical material was synthesized and investigated taking into account the moving phase transition boundary [3]. An approximate analytic-numerical solution of a nonlinear problem is constructed for a three-step mode of heat transfer for the case when the solution of a two-dimensional problem T 1 ,  2 ,  is represented as a product of one-dimensional problems T 1 ,  2 ,   T1 1 , T2  2 ,  .The desired solution to the temperature determination problem is represented as: T T1T2  T  BB1   T2 1 1  T2 1 /  H1 T1  u   H1* T1  u1   , (34) n x1 1 x1 p n 1   T1i  i ,    0i  i     ni  i  cos ni2    n  i sin ni2  , i  1,2 (35)  Ti *  i ,    Ani n 1 sin  ni ( i   mi )  ni  exp   ni 2  , (36) where the coefficients Ani are functions of the frequency characteristics  k ,  k associated with the characteristics of the heat transfer process and the magnitude of the phase transition;  ni ,  ni are determined by hyperbolic-trigonometric functions. Formulas were obtained to determine the non-stationary temperature at an arbitrary point of anisotropic biophysical material, depending on the coordinate of the phase transition plane, and the change in ambient temperature. 3 Numerical analysis of the study The given algorithm [9, 10] is tested on a model for the ambient temperature [2, 8] t c  65 оС, relative humidity   0.8 and the speed of the air v  2 m/s, density   581 kg/m3, density of absolutely dry body  0  457 kg/m3, basis density  б  415 kg/m3, porosity П = 0,6, saturated vapour density  n  0.013188 kg/m3,  a0  1.29 kg/m3,  L  1000 kg/m3, moisture exchange coefficient β = 0.000976 m/s, mass exchange coefficient   22.32599 Watt/(m2·degree), thermal conductivity coefficient   0.299993 Watt/(m·degree). Fig. 1. The dependence of the coordinate of phase transition in a material with a density of  0  480 kg/m3 for different values of the temperature of the environment (curve 1 – Tc  25 0С.; 2 – Tc  45 0С.; 3 – Tc  55 0С.; 4 – Tc  65 0С.; 5 – Tc  75 0С.; 6 – Tc  85 0С.). Conclusion A new two-dimensional nonlinear mathematical model of the heat transfer process in anisotropic biophysical materials was constructed taking into account the moving phase transition boundaries. The arbitrary orientation of the principal axes of the thermal conductivity tensor is taken into account, and the influence of the principal components and orientations of the principal axes of the thermal conductivity tensor and non-stationary temperature fields in anisotropic biophysical material are determined. An analytical-numerical method has been developed for determining heat transfer in an anisotropic biophysical material with a moving boundary of phase transitions as well as for establishing moving boundaries of a phase transition in a rectangular region, taking into account the main axes of anisotropy. An algorithm is constructed to determine the moving boundary of the evaporation zone in a biophysical material, depending on the anisotropic characteristics of the material, the parameters of the environment. References 1. Fedotkin, I.M., Burlyay, I.Yu., Rumshin, N.A., Burlyay,Yu.I.: Mathematical modeling of technological processes. Technics, Kyiv, 338 p. (2004). 2. Borovikov, A.N., Ugolev, B.N.: Handbook of wood. Forest industry, Moskow, 296 p. (1989). 3. Kartashov, E.M.: Analytical methods in the theory of thermal conductivity of solids. Higher school, Moscow, 480 p. (1985). 4. Marchuk, G.I.: Methods of splitting. Science, Moscow, 264 p. (1988). 5. Nikitenko, N.I.: Theory of heat and mass transfer. Naukova Dumka, Kyiv, 351 p. (1983). 6. Kafarov, V.V.: Mathematical modeling of the main processes of chemical production. Higher. School, Moscow, 400 p. (1991). 7. Lykov, A.V. Heat and mass transfer: a reference. Energy, Moscow, 560 p. (1971). 8. John F. Sian. Wood: influence of moisture on physical properties. Virginia, 227 p. (1995). 9. Sokolovskyy, Y.; Sinkevych, O.: Software for automatic calculation and construction of chamber drying wood and its components. In: Perspective Technologies and Methods in MEMS Design, MEMSTECH 2016 - Proceedings of 12th International Conference, pp. 209-213. Lviv-Polyana, Ukraine (2016). 10. Sokolowskyi, Y.; Shymanskyi, V.; Levkovych, M.; Yarkun, V.: Mathematical and software providing of research of deformation and relaxation processes in environments with fractal structure. In: Proceedings of the 12th International Scientific and Technical Conference on Computer Sciences and Information Technologies, CSIT 2017, pp. 24-27. Lviv, Ukraine (2016). 11. Sokolovskyy, Y.; Boretska, I.; Gayvas, B.; Kroshnyy, I.: Mathematical modeling of the heat-mass-exchange in anisotropic environments taking into account the boundary of phase transition. In: 13th IEEE International Scientific and Technical Conference on Computer Sciences and Information Technologies, CSIT 2018, Lviv, Ukraine, pp. 147- 150. 12. Sokolovskyy, Y.; Boretska, I.; Kroshnyy, I.; Gayvas, B.: Mathematical models and analysis of the heat-mass-transfer in anisotropic materials taking into account the boundaries of phase transition. In: 15th International Conference on the Experience of Designing and Application of CAD Systems, CADSM 2019 – Proceedings. Polyana, Ukraine (2019).