Transient Heat Conduction in Capillary Porous Bodies

Capillarity is a well known phenomenon in physics and engineering. Porous materials such as soil, sand, rocks, mineral building elements (cement stone or concrete, gypsum stone or plasterboards, bricks, mortar, etc.), biological products (wood, grains, fruit, etc.) have microscopic capillaries and pores which cause a mixture of transfer mechanisms to occur simultaneously when subjected to heating or cooling. In the most general case each capillary porous material is a peculiar system characterized by the extremely close contact of three intermixed phases: gas (air), liquid (water) and solid. Water may appear in them as physically bounded water and capillary water (Chudinov, 1968, Twardowski, Richinski & Traple, 2006). Both the bounded water and the capillary water can be found in liquid or hard aggregate condition. Physically bounded water co-operates with the surface of a solid phase of the materials and has different properties than the free water. The maximum amount of bounded water in porous materials corresponds to the maximal hygroscopicity, i.e. moisture absorbed by the material at the 100% relative vapour pressure. The maximum hygroscopicity of the biological capillary porous bodies is known as fibre saturation point. Capillary water fills the capillary tube vessels, small pores or sharp, narrow indentions of bigger pores. It is not bound physically and is called free water. Free water is not in the same thermodynamic state as liquid water: energy is required to overcome the capillary forces, which arise between the free water and the solid phase of the materials. For the optimization of the heating and/or cooling processes in the capillary porous bodies, it is required that the distribution of the temperature and moisture fields in the bodies and the consumed energy for their heating at every moment of the process are known. The intensity of heating or cooling and the consumption of energy depend on the dimensions and the initial temperature and moisture content of the bodies, on the texture and microstructural features of the porous materials, on their anisotropy and on the content and aggregate condition of the water in them, on the law of change and the values of the temperature and humidity of the heating or cooling medium, etc. (Deliiski, 2004, 2009). The correct and effective control of the heating and cooling processes is possible only when its physics and the weight of the influence of each of the mentioned above as well as of many other specific factors for the concrete capillary porous body are well understood. The summary of the influence of a few dozen factors on the heating or cooling processes of the


Introduction
Capillarity is a well known phenomenon in physics and engineering. Porous materials such as soil, sand, rocks, mineral building elements (cement stone or concrete, gypsum stone or plasterboards, bricks, mortar, etc.), biological products (wood, grains, fruit, etc.) have microscopic capillaries and pores which cause a mixture of transfer mechanisms to occur simultaneously when subjected to heating or cooling. In the most general case each capillary porous material is a peculiar system characterized by the extremely close contact of three intermixed phases: gas (air), liquid (water) and solid. Water may appear in them as physically bounded water and capillary water (Chudinov, 1968, Twardowski, Richinski & Traple, 2006. Both the bounded water and the capillary water can be found in liquid or hard aggregate condition. Physically bounded water co-operates with the surface of a solid phase of the materials and has different properties than the free water. The maximum amount of bounded water in porous materials corresponds to the maximal hygroscopicity, i.e. moisture absorbed by the material at the 100% relative vapour pressure. The maximum hygroscopicity of the biological capillary porous bodies is known as fibre saturation point. Capillary water fills the capillary tube vessels, small pores or sharp, narrow indentions of bigger pores. It is not bound physically and is called free water. Free water is not in the same thermodynamic state as liquid water: energy is required to overcome the capillary forces, which arise between the free water and the solid phase of the materials. For the optimization of the heating and/or cooling processes in the capillary porous bodies, it is required that the distribution of the temperature and moisture fields in the bodies and the consumed energy for their heating at every moment of the process are known. The intensity of heating or cooling and the consumption of energy depend on the dimensions and the initial temperature and moisture content of the bodies, on the texture and microstructural features of the porous materials, on their anisotropy and on the content and aggregate condition of the water in them, on the law of change and the values of the temperature and humidity of the heating or cooling medium, etc. (Deliiski, 2004(Deliiski, , 2009). The correct and effective control of the heating and cooling processes is possible only when its physics and the weight of the influence of each of the mentioned above as well as of many other specific factors for the concrete capillary porous body are well understood. The summary of the influence of a few dozen factors on the heating or cooling processes of the capillary porous bodies is a difficult task and its solution is possible only with the assistance of adequate for these processes mathematical models. There are many publications dedicated to the modelling and computation of distribution of the temperature and moisture content in the subjected to drying capillary porous bodies at different initial and boundary conditions. Drying is fundamentally a problem of simultaneous heat and mass transfer under transient conditions. Luikov (1968) and later Whitaker (1977) defined a coupled system of non-linear partial differential equations for heat and mass transfer in porous bodies. Practically all drying models of capillary porous bodies are based on these equations and include a description of the specific initial and boundary conditions, as well as of thermo-and mass physical characteristics of the subjected to drying bodies (Ben Nasrallah & Perre, 1988, Doe, Oliver & Booker, 1994, Ferguson & Lewis, 1991, Kulasiri & Woodhead 2005, Murugesan et al., 2001, Zhang, Yang & Liu 1999. As a rule the non-defective drying is a very long continuous process, which depending on the dimensions of the bodies and their initial moisture content can last many hours, days or even months. However there are a lot of cases in the practice where the capillary porous bodies are subjected to a relatively short heating and/or cooling, as a result of which a significant change in their temperature and a relatively small change in their moisture content occurs. A typical example for this is the change in temperature in building elements under the influence of the changing surrounding temperature during the day and night or during fire. Widely used during the production of veneer, plywood or furniture parts technological processes of thermal treatment of wood materials (logs, lumber, etc.) with the aim of plasticizing or ennoblement of the wood are characterized by a controlled change in the temperature in the volume of the processed bodies, without it being accompanied by significant changes in their moisture content. In these as well as in other analogous cases of heating and/or cooling of capillary porous bodies the calculation of the non-stationary change in the temperature field in the bodies can be carried out with the assistance of models, which do not take into account the change in moisture content in the bodies. The number of such published models is very limited, especially in comparison to the existing large variety of drying models. Axenenko (1995) presents and uses a 1D model for the computation of the temperature change in the exposed to fire gypsum plasterboards. According to the author, after obtaining the temperature fields inside the plasterboard, changes in the material properties and temperature deformations can be calculated and used as initial data for the study of the structural behaviour of the entire plasterboard assembly. Considerate contribution to the calculation of the non-stationary distribution of the temperature in frozen and non-frozen logs and to the duration of their heating has been made by H. P. Steinhagen. For this purpose, he, alone, (Steinhagen, 1986(Steinhagen, , 1991 or with coauthoring (Steinhagen, Lee & Loehnertz, 1987), (Steinhagen & Lee, 1988) has created and solved a 1-dimensional, and later a 2-dimensional (Khattabi & Steinhagen, 1992, 1993 mathematical model, whose application is limited only for u ≥ 0,3 kg.kg -1 . The development of these models is dominated by the usage of the method of enthalpy, which is rather more complicated than its competing temperature method. The models contain two systems of equations, one of which is used for the calculation of the change in temperature at the axis of the log, and the other -for the calculation of the temperature distribution in the remaining points of its volume. The heat energy, which is needed for the melting of the ice, which has been formed from the freezing of the hygroscopically bounded water in the wood, although the specific heat capacity of that ice is comparable by value to the capacity of the frozen wood itself (Chudinov, 1966), has not been taken into account. These models assume that the fibre saturation point is identical for all wood species and that the melting of the ice, formed by the free water in the wood, occurs at 0ºC. However, it is known that there are significant differences between the fibre saturation point of the separate wood species and that the dependent on this point quantity of ice, formed from the free water in the wood, thaws at a temperature in the range between -2°C and -1°C (Chudinov, 1968). This paper presents the creation and numerical solutions of the 3D, 2D, and 1D mathematical models for the transient non-linear heat conduction in anisotropic frozen and non-frozen prismatic and cylindrical capillary porous bodies, where the physics of the processes of heating and cooling of bodies is taken into account to a maximum degree and the indicated complications and incompleteness in existing analogous models have been overcome. The solutions include the non-stationary temperature distribution in the volume of the bodies for each u ≥ 0 kg.kg -1 at every moment of their heating and cooling at prescribed surface temperature, equal to the temperature of the processing medium or during the time of convective thermal processing.

Mechanism of heat distribution in capillary porous bodies
During the heating or cooling of the capillary porous materials along with the purely thermal processes, a mass-exchange occurs between the processing medium and the materials. The values of the mass diffusion in these materials are usually hundreds of times smaller than the values of their temperature conductivity. These facts determine a not so big change in defunding mass in the materials, which lags significantly from the distribution of heat in them during the heating or cooling. This allows to disregarding the exchange of mass between the materials and the processing medium and the change in temperature in them to be viewed as a result of a purely thermo-exchange process, where the heat in them is distributed only through thermo-conductivity. Because of this the mechanism of heat distribution in capillary porous bodies can be described by the equation of heat conduction (also known as the equation of Fourier-Kirchhoff). Its most compact form is as follows: www.intechopen.com (1) This form holds for each coordinate system and for each processing medium -both for immobile and mobile.
As it was described in the introduction, the water contained in these bodies can be found in liquid or hard aggregate condition. It is known that the specific heat capacity of the liquid (non-frozen) water at 0°С is equal to 4237 J.kg -1 .K -1 , and the specific heat capacity of ice is 2261 J.kg -1 .K -1 , i.e. almost two times smaller (Chudinov, 1966(Chudinov, , 1984. Because of this, the frozen water in capillary porous bodies causes smaller values of c in comparison to the case, when the water in them is completely liquid. The ice in capillary porous bodies can be formed from the freezing of higroscopically bounded water or of the free water in them. It is widely accepted that the phase transition of water into ice and vice versa to be expressed with the help of the so-called "latent heat" in the ice of the frozen body. When solving problems, connected to transient heat conduction in frozen bodies, it makes sense to include the latent heat in the so-called effective specific heat capacity c e (Chudinov, 1966), which is equal to the sum of the own specific heat capacity of the body с and the specific heat capacity of the ice, formed in them from the freezing of the hygroscopically bounded water and of the free water, i.e. eb w f w cc c c = ++. (2) When solving the problems, it must be taken into consideration that the formation and thawing of both types of ice in these bodies takes place at different temperature ranges. Because of this for each of the diapasons in equation (2) the sum of с with bw c and/or fw c (shown below as an example) participates. When modeling processes of transient heat conduction in anisotropic capillary porous bodies it is also necessary to take into consideration that the thermal conductivity of these bodies λ apart from Т and u depends additionally on the direction of the influencing heat flux towards the anatomic directions of the body -radial, tangential and longitudinal to the fibers.

Mathematical models for transient heat conduction in prismatic bodies
If it is assumed, that the anatomical directions of a prismatic capillary porous body coincide with the coordinate axes, in the absence of an internal source of heat q in equation (1), the following form of this equation in the Cartesian coordinate system is obtained: Txyz After the differentiation of the right side of equation (3) on the spatial coordinates x, y, and z, excluding the arguments in the brackets for shortening of the record, the following mathematical model of the process of non-stationary heating or cooling (further calling thermal processing) of the capillary porous bodies with prismatic form is obtained: and boundary conditions: • during the time of thermal processing of the prisms at their prescribed surface temperature, equal to the temperature of the processing medium: [ ] The system of equations (4) ÷ (9) presents a 3D mathematical model, which describes the change in temperature in the volume of capillary porous bodies with prismatic form during the time of their thermal processing at corresponding initial and boundary conditions. When the length of the subjected to thermal processing body exceeds its thickness by at least ( ) p rt 2 22 , 5 λ λ λ ÷ + times, then the heat transfer through the frontal sides of the body can be neglected, because it does not influence the change in temperature in the cross-section, which is equally distant from the frontal sides. In these cases for the calculation of the change in T in this section (i.e. only along the coordinates x and y) the following 2D model can be used: with an initial condition and boundary conditions:

155
• for thermal processing of the prisms at their prescribed surface temperature: • for convective thermal processing of the prisms: When the thickness of the body is smaller than its width by at least 2 ÷ 3 times, and than its length by at least ( ) p r 22 , 5 λ λ ÷ times, then the non-stationary change in T, for example along the radial direction x of the body, coinciding with its thickness in the section, equally distant from the frontal sides, can be calculated using the following 1D model: with an initial condition and boundary conditions: • for thermal processing of the prisms at their prescribed surface temperature: • for convective thermal processing of the prisms:

Mathematical models for transient heat conduction in cylindrical bodies
The mechanism of the heat distribution in the volume of cylindrical capillary porous bodies during their thermal processing can be described by the following non-linear differential equation with partial derivatives, which is obtained from the equation (3) after passing in it from rectangular to cylindrical coordinates (Deliiski, 1979) er  22  2   22  2  2  p  r  p  22   ,, , Trz  Trz  Trz  Trz  cT u  T u  T u  rr  rr   Tu  Tu  Trz  Trz  Trz  Trz  Tu  Tr  T  z If heating or cooling cylindrical bodies of material, which is homogenous in their cross section, the distribution of T in their volume does not depend on φ, but only depends on r and z. Consequently, when excluding the participants in the equation (19), containing φ and when also omitting the arguments in the brackets for the shortening of the record, the following 2D mathematical model is obtained, which describes the change of temperature in the volume of capillary porous bodies with cylindrical form: with an initial condition and boundary conditions: • for thermal processing of the bodies at their prescribed surface temperature: • for convective thermal processing of the bodies: When the length of the body exceeds its diameter by at least ( ) p r 22 , 5 λ λ ÷ times, then the heat transfer through the frontal sides of the body can be neglected, because it does not influence the change in temperature of its cross section, which is equally distant from the frontal sides. In such cases, for the calculation of the change in T only along the coordinate r of this section, the following 1D model can be used: with an initial condition and with boundary conditions, which are identical to the ones in equations (17) and (18), but with derivative of T along r instead of along x in (18).

Transformation of the models for transient heat conduction in suitable form for programming
Analytical solution of mathematical models, which contain non-linear differential equations with partial derivatives in the form of (4) and (20), is practically impossible without significant simplifications of these equations and of their boundary conditions. For the numerical solution of models with such equations the methods of finite differences or of the finite elements can be used. When the bodies have a correct shape -prismatic or cylindrical, the method of finite differences is preferred, because its implementation requires less computational resources from the computers.
For the numerical solution of the above presented models for transient heat conduction it makes sense to use the explicit form of the finite-difference method, which allows for the exclusion of any simplifications. The large calculation resources of the contemporary computers eliminate the inconvenience, which creates the limitation for the value of the step along the time coordinate τ Δ by using the explicit form (refer to equation (47)). According to the main idea of the finite-difference method, the temperature, which is a uninterrupted function of space and time, is presented using a grid vector, and the are approximated using the built computational mesh along the spatial and time coordinates through their finite-difference (discrete) analogues. For this purpose the subjected to thermal processing body, or 1/8 of it in the presence of mirror symmetry towards the other 7/8, is "pierced" by a system of mutually perpendicular lines, which are parallel to the three spatial coordinate axes. The distances between the lines, also called as the step of differentiation, is constant for each coordinate direction (Fig. 1). The knots from 0 to 6 with temperatures accordingly from 0 T to 6 T are centers for the presented and neighboring it volume elements. The length of the sides of the volume element x Δ , y Δ and z Δ are steps of differentiation, the size of which determines the distance between the separate knots of the mesh. The entire time of thermal processing of the body is also separated into definite number n intervals (steps) with equal duration τ Δ . A volume element of a subjected to thermal processing body used for the solution of equation (4) is shown on Fig. 1, together with its belonging part from the rectangular calculation mesh.  (4) using the explicit form of the finite-difference method As a result of the described procedures, the process for the solution of the non-linear differential equation with partial derivatives (4) is transformed to the solution of the equivalent to it system of linear finite-difference equations, requiring the carrying out of numerous one-type, but not complicated algebraic operations. In the built in the body 3D rectangular mesh, the heat transfer is taken into consideration only along its gradient lines. For each cross point of three mesh lines, called a knot, a finite-difference equation is derived, using which the temperature of this knot is calculated. Finite-difference equations are obtained by substituting derivatives in the differential equation (4), with their approximate expressions, which are differences between the values of the function in selected knots of the calculation mesh. The temperature in knot 0 with coordinates ( i x , j y , k z , n τ ) on Fig. 1 is designated as , and the temperature in neighboring knots accordingly as 11 , , T in the preceding moment n τ Δ need to be known.

Discrete analogues of models for transient heat conduction in prismatic bodies
The transformation of the non-linear differential equation with partial derivatives (4) in its discrete analogue with the help of the explicit form of the finite-difference method is carried out using the shown on Fig. 2 coordinate system for the positioning of the knots of the calculation mesh, in which the distribution of the temperature in a subjected to thermal processing capillary porous body with prismatic form is computed.
Then after carrying out the differentiation of λ along Т in equation (4) and substituting in it with finite differences of the first derivative along the time and the first and second derivatives along the spatial coordinates, the following system of equations is obtained: , , 22 2( ) 1( 2 7 3 , 1 5 ) After transformation of the equation (28) it is obtained, that the value of the temperature in whichever knot of the built in the volume of a subjected to thermal processing body 3D mesh at the moment (1 ) n τ + Δ is determined depending on the already computed values for the temperature in the preceding moment n τ Δ using the following system of finite differences equations: ,, 1 ,, 1 ,, , ,, The initial condition (5) in the 3D mathematical model is presented using the following finite differences equation: The boundary conditions (6) ÷ (9) get the following easy for programming form: • for thermal processing of the prisms at their prescribed surface temperature: • for convective thermal processing of the prisms: In the boundary conditions (31) ÷ (34) is reflected the requirement for the used for the solution of the models computation environment of FORTRAN, that the knots of the mesh, which are positioned on the corresponding surface of the subjected to thermal processing prism, to be designated with number 1, i.e. 1 ≤ i ≤ M; 1 ≤ j ≤ N; 1 ≤ k ≤ KD (Dorn & McCracken 1972). Using the system of equations (29) ÷ (34) the distribution of the temperature in subjected to thermal processing prismatic materials can be calculated, whose radial, tangential and parallel to the fibers direction coincide with its coordinate axes х, у and z.
Since the subjected to thermal processing in the practice prismatic materials usually do not have a clear radial or clear tangential sides, but are with partially radial or partially tangential, then in equation (29) instead of the coefficients λ 0 in the observed two anatomical directions their average arithmetic value can be substituted, which determines the thermal conductivity at 0°С cross sectional to the body's fibers λ 0cr : 0r 0t 0cr 2 λ λ λ Also the thermal conductivity at 0°С in the direction parallel to the fibers λ 0p can be expressed through λ 0cr using the equation where the coefficient 0p p/cr 0cr K λ λ = depends on the type of the capillary porous material.
For the unification of the calculations it makes sense to use one such step of the calculation mesh along the spatial coordinates ∆x = ∆y = ∆z (refer to Fig. 2). Taking into consideration this condition, and also of equations (35) and (36), the system of equations (29) becomes 1 ,, , In this case in equations (32) and (33) instead of r α and t α , cr α must be used, and instead of 0r λ and 0t λ , 0cr λ must be used.
• for convective thermal processing of the prisms: ( ) The discrete analogue of the 1D model, in which equations (15) with an initial condition and boundary conditions: • for thermal processing of the prisms at their prescribed surface temperature: • for convective thermal processing of the prisms: Using the system of equations (37) , www.intechopen.com in which min T and max T are correspondingly the smallest and biggest of all values of the temperatures, encountered in the initial and boundary conditions of the heat transfer when solving the mathematical model, and d K is a coefficient, reflecting the dimensioning of the heat flux when calculating the step τ Δ : for a 3D heat flux dp / c r 6 KK = ; for a 2D heat flux d 4 K = and for a 1D heat flux d 2 K = .

Discrete analogues of models for transient heat conduction in cylindrical bodies
For the obtaining of a discrete analogues of equations (20) ÷ (24) the explicit form of the finite-difference method has been used in a way comparable to the reviewed above case of 2D thermal processing of prismatic bodies. The transformation of the non-linear differential equation with partial derivatives (20) in its discrete analogue is carried out using the shown on Fig. 3 coordinate system for the positioning of the knots from the calculation mesh, in which the distribution of temperature in the longitudinal section of subjected to thermal processing cylindrical body is calculated. The calculation mesh for the solution of the 2D model with the help of the finite-differences method is built on ¼ of the longitudinal section of the cylindrical body due to the circumstance that this ¼ is mirror symmetrical towards the remaining ¾ of the same section.
• for convective heating or cooling of the bodies in the processing medium: ( ) The discrete analogue of the 1D model, in which equations (25), (26), (17) and (18) T  TT  cr  TT  T  TT  TT with an initial condition presented through equation (44) and boundary conditions, presented through equations (45) and (46), where in (46) instead of cr α and 0cr λ , r α and 0r λ are used. The value of the τ Δ , which guarantees the obtaining of robust solutions for the presented models, is determined by the condition for stability, described by equation (47). When solving the 2D model dp / r 4 KK = , and when solving the 1D model d 2 K = in (47).

Mathematical description of the thermo-physical characteristics of the capillary porous bodies
The mathematical description of the thermo-physical characteristics of the capillary porous bodies consists of the deduction of summarized dependencies as a function of the factors which have an influence on them, which with maximum precision correspond to the their experimentally determined values in the interesting for the practice sufficiently wide ranges for the change in the factors. It is necessary to have such a description when solving the above presented mathematical models of the transient heat conduction in capillary porous bodies. The following thermo-physical characteristics are present in these models: specific heat capacity, thermal conductivity and density of the capillary porous bodies. The mathematical description of the above mentioned characteristics of the wood, which is a typical representative of the studied bodies, is shown as an example below.

Mathematical description of the specific heat capacity of the wood
The mathematical description of e c is done according to equation (2) and is based on an analysis of the physics of the wood frosting and defrosting process. The description reveals the facts determined by Chudinov (1966Chudinov ( , 1984, which state that the thawing of the ice from the situated in the wood free water occurs in the temperature range between 271,15 and 272,15 K, and the thawing of the ice from the situated in the wood hygroscopically bounded water depends on the temperature, which is Т < 271,15 К. Besides this, at Т < 271,15, a certain portion, nfw u , from the hygroscopically bounded water is found in a non-frozen state, where the value of nfw u is calculated by equation (58). The own specific heat capacity of the wood c is described mathematically using the shown on Fig. 4 experimentally determined in the dissertations of Kanter (1955) and Chudinov (1966) data for its change as a function of t and u. In the description apart from t and u, the independent parameter fsp u has been input, which is different for the separate wood species and reflects the influence of the anatomic characteristics of the wood on c. From the analysis on Fig. 4 it can be seen that at t = -2°C and all values for the content of water 0,3 u > kg.kg -1 during the thawing of wood a jump in the values of c is observed. This jump is explained by the phase transition of the frozen free water in the wood at the given values for t and u, when the influence on c of a significant difference in the specific heat capacity of the water in a liquid and hard aggregate condition is observed. On Fig. 4 the jump of c for 0,2 u = kg.kg -1 occurs at a temperature around -17°С. It is characterized by the phase transition of the frozen part of the hygroscopically bounded water in the wood at these values for t and u. Both at 0 u = , and at 0,1 u = kg.kg -1 there is a lack of a jump in the values for c, because even at -40°С there is no ice in the wood. The following equations for fsp (,, ) cT uu have been derived, which are valid for all wood species for the intervals 0 ≤ u ≤ 1,2 kg.kg -1 and 223 ≤ Т ≤ 373 К (Deliiski 1990(Deliiski , 2003b(Deliiski , 2004 The calculated according to equation (58)  Taking into consideration the latent heat in the phase transition (crystallization) of the water, equal to 5 3, 34.10 J.kg -1 (Chudinov, 1966) the following equation has been obtained for the calculation of fw c (Deliiski, 2003b, 2004, 2009, Dzurenda & Deliiski, 2010 www.intechopen.com The change in the calculated according to equation (59) values for fw c depending on u and fsp u of the wood is shown on Fig. 6. It can be seen that the decrease in fsp u at a given u causes a proportional increase of the capacity fw c due to the increase in the numerator in equation (59). The increase of u causes an increase of fw c , which is caused by the increase of the quantity of ice in the more humid wood, which was formed from the free water. For the calculation of bw c the following equation (Deliiski 2003b(Deliiski , 2009 The change in the calculated according to equation (60) values for bw c depending on T and u is shown on Fig. 7 for the comparatively more often subjected to thermal processing beech wood with fsp 0,31 u = kg.kg -1 . When Т and u are outside of the indicated in equations (59) and (60) ranges, the values of fw c and bw c have been assumed to be equal to zero during the solution of the models.

Mathematical description of the thermal conductivity of the wood
A mathematical description of the thermal conductivity of the wood λ has been done using the shown on Fig. 8 experimentally determined in the dissertation by Kanter (1955) and Chudinov (1966) Fig. 4 find a wide use in both the European (Shubin, 1990, Trebula, 1996, Videlov, 2003 and the American specialized literature (Steinhagen, 1986(Steinhagen, , 1991 when calculating various processes of thermal processing of wood. From the analysis on Fig. 8 it can be seen that at fsp uu > the temperature t = -2°C for the defrosting of the wood is critical, because at this temperature a jump in the values of cr λ is observed. This jump increases with the increase in u. It is explained by a significant difference in the thermal conductivity of the water in liquid and in hard aggregate condition. In the description of cr λ , as well as in that of e c , fsp u has been introduced as an independent variable for different wood species, instead of the ordinarily used averaged value of fsp 0,30 u = kg.kg -1 for all wood species. The following equations for the calculation of the coefficients of λ in equation (27) for the ranges 01 , 2 u ≤ ≤ kg.kg -1 and 223 373 T ≤ ≤ К have been derived (Deliiski, 1994, 2003b, Dzurenda & Deliiski, 2010: 0,165 1, 39 3,8 3, 3.10 1,015.10 Ku The equations, which have been suggested by Chudinov (1966Chudinov ( , 1968) and shown in (Deliiski, 1977) can be used for the determining of the values of the coefficient a K in equation (61), which takes into account the influence on 0 λ of the heat flux towards the separate anatomic directions of the wood, i.e. for the determining of r λ , t λ and p λ . The coefficients and β in equation (27)     Equation (70), together with equation (27), make the calculation of λ applicable for all temperatures larger than -50°С (i.e. 223,15 T = К), because only at 223,15 T = К and the minimum possible for the wood species fsp 0,2 u = kg.kg -1 , the value of nfw u is nfw 0,125 u = kg.kg -1 (Fig. 5).

Mathematical description of the wood density
The mathematical description of the wood density ρ depending on the moisture content u and basic wood density b ρ which influence it has been carried out using the shown on Fig. 9 experimental values set by (Sergovski, 1975) and shown in (Deliiski, 1977(Deliiski, , 2003b. The following equations, which are applicable to all wood species, have been derived: where according to (Siau 1984) 20 fsp The values for 20 fsp u are given in the specialized literature (Siau, 1984, Videlov, 2003.

Computer solution of the models for transient heat conduction in the bodies
For the numerical solution of the above presented models for transient heat conduction in prismatic and cylindrical capillary porous bodies a software package has been prepared in www.intechopen.com FORTRAN (Dorn & McCracken, 1972), which has been input in the developed by Microsoft calculation environment of VISUAL FORTRAN PROFESSIONAL (Deliiski 2003b). The software package can be used for the calculation and colour visualization (either as animation of the whole process or as 3D, 2D, 1D graphs of each desired moment of the process) of the non-stationary distribution of the temperature fields in the materials containing or not containing ice during their thermal processing. The computation of the change in the temperature field in the volume of materials containing ice in the beginning of their thermal processing is interconnected for the periods of the melting of the ice and after that, taking into account the flexible spatial boundary of the melting ice. The computation of the temperature fields is done interconnectedly and for the processes of heating and consequent cooling of the materials, i.e. the calculation of the non-stationary change in temperatures in the volume of the materials during the time of their cooling begins from the already reached during the time of calculations distribution of temperature in the end of the heating. Based on the calculations it can be determined when the moment of reaching in the entire volume of the heated wood has occurred for the necessary optimal temperatures needed for bending of the parts or for cutting the veneer, as well as the stage of the ennoblement of the wood desired by the clients. kg.kg -1 , u = 0,6 kg.kg -1 and 20 fsp 0,29 u = kg.kg -1 . The heating of the prisms continues until the reaching of the minimally required for cutting of veneer temperature in their centre, equal to 0 c 50 C t = . During the time of cooling of the heated prisms a redistribution and equalization of t in their cross section takes place, which is especially appropriate for the obtaining of quality veneer. The change in m t and t is shown on Fig. 11 in 5 characteristic points from the cross section of the prisms with coordinates, which are given in the legend of the graphs. has been also calculated with the given above parameters during the time of a 3-stage high temperature thermal processing in autoclave and during the time of the consequent cooling with surface convection at 0 m 20 C t = (Fig. 13). Using 3D graphs and 2D diagrams a part of the results is shown on Fig. 14 from the simulation studies on the heat transfer in the radial and longitudinal direction of the frozen beech logs with t 0 = -2°C, whose temperature field is shown on the left side on Fig. 12. The non-stationary temperature distribution during specific time intervals of the thermal processing is clearly observed from the 3D graphs (left columns on Fig. 14). The 2D diagrams which show in more detail the results from the simulations can be used rather for qualitative than quantitative analysis of the thermal processing of the materials (right columns on Fig. 14). On the left parts of Fig. 10, Fig. 11, Fig. 12 and Fig. 13 the characteristic non-linear parts can be seen well, which show a slowing down in the change in t in the range from -2°С to -1°С, in which the melting of the ice takes place, which was formed in the wood from the freezing of the free water in it. This signifies the good quality and quantity adequacy of the mathematical models towards the real process of heating of ice-containing wood materials. The calculated with the help of the models results correspond with high accuracy to wide experimental data for the non-stationary change in t in the volume of the containing and not containing ice wood logs, which have been derived in the publications by (Schteinhagen, 1986(Schteinhagen, , 1991 and (Khattabi & Steinhagen, 1992, 1993. The results presented on the figures show that the procedures for calculation of nonstationary change in t in prepared software package, realizing the solution of the mathematical models according to the finite-differences method, functions well for the cases of heating and cooling both for frozen and non-frozen materials at various initial and boundary conditions of the heat transfer during the thermal processing of the materials. The good adequacy and precision of the models towards the results from numerous own and foreign experimental studies allows for the carrying out of various calculations with the models, which are connected to the non-stationary distribution of t in prismatic and cylindrical materials from various wood species and also to the heat energy consumption by the wood at random encountered in the practice conditions for thermal processing.

Conclusion
This paper describes the creation and solution of non-linear mathematical models for the transient heat conduction in anisotropic frozen and non-frozen capillary porous bodies with prismatic and cylindrical shape and at any u ≥ 0 kg.kg -1 . The mechanism of the heat distribution in the entire volume of the bodies is described only by one partial differential equation of heat conduction. For the first time the own specific heat capacity of the bodies and the specific heat capacity of the ice, formed in them from the freezing of the hygroscopically bounded water and of the free water are taken into account in the models. The models take into account the physics of the described processes and allow the 3D, 2D and 1D calculation of the temperature distribution in the volume of subjected to heating and/or cooling anisotropic or isotropic bodies in the cases, when the change in their moisture content during the thermal processing is relatively small. For the solution of the models an explicit form of the finite-difference method is used, which allows for the exclusion of any simplifications in the models. For the usage of the models it is required to have the knowledge and mathematical description of several properties of the subjected to thermal processing frozen and nonfrozen capillary porous materials. In this paper the approaches for mathematical description of thermo-physical characteristics of materials from different wood species, which are typical representatives of anisotropic capillary porous bodies, widely subjected to thermal treatment in the practice are shown as examples.
For the numerical solution of the models a software package has been prepared in FORTRAN, which has been input in the developed by Microsoft calculation environment of Visual Fortran Professional. The software allows for the computations to be done for heating and cooling of the bodies at prescribed surface temperature, equal to the temperature of the processing medium or during the time of convective thermal processing. The computation of the change in the temperature field in the volume of materials containing ice in the beginning of their thermal processing is interconnected for the periods of the melting of the ice and after that, taking into account the flexible boundary of the melting ice. The computation of the temperature fields is done interconnectedly and for the processes of heating and consequent cooling of the materials, i.e. the calculation of the change in temperatures in the volume of the materials during the time of their cooling begins from the already reached during the time of calculations distribution of temperature in the end of the heating. It is shown how based on the calculations it can be determined when the moment of reaching in the entire volume of the heated and after that cooling body has occurred for the necessary optimal temperatures needed, for example, for bending of wood parts or for cutting the veneer from plasticised wooden prisms or logs. The models can be used for the calculation and colour visualization (either as animation of the whole process or as 3D, 2D, 1D graphs of each desired moment of the process) of the distribution of the temperature fields in the bodies during their thermal processing. The development of the models and algorithms and software for their solution is consistent with the possibility for their usage in automatic systems with a model based (Deliiski 2003a(Deliiski , 2003b(Deliiski , 2009 or model predicting control of different processes for thermal treatment.

Acknowledgement
This work was supported by the Scientific Research Sector of the University of Forestry, Sofia, Bulgaria.