Modeling of the Two-Dimensional Thawing of Logs in an Air Environment

A two-dimensional mathematical model has been created, solved, and verified for the transient nonlinear heat conduction in logs during their thawing in an air environment. For the numerical solution of the model, an explicit form of the finite-difference method in the computing medium of Visual FORTRAN Professional has been used. The chapter presents solutions of the model and its validation towards own experimental studies. During the validation of the model, the inverse task of the heat transfer has been solved for the determination of the logs ’ heat transfer coefficients in radial and longitudinal directions. This task has been solved also in regard to the logs ’ surface temperature, which depends on the mentioned coefficients. The results from the experimental and simulative investigation of 2D nonstationary temperature distribution in the longitudinal section of poplar logs with a diameter of 0.24 m, length of 0.48 m, and an initial temperature of approximately – 30°C during their many hours thawing in an air environment at room temperature are presented, visualized, and analyzed.


Introduction
The duration and the energy consumption of the thermal treatment of frozen logs aimed at their thawing and plasticizing for the production of veneer in winter are very high [1][2][3][4][5][6][7][8][9]. For example, thawing and plasticizing of poplar and pine logs with an initial temperature of -10°C and moisture content of 0.6 kgÁkg À1 about 53 kWhÁm À3 and 64 kWhÁm À3 thermal energy, respectively, are needed [9].
The computation of the temperature field in logs during their thawing in water or steam is carried out using mathematical models, which solve the so-called direct task of the heat transfer. This is the task when all variables in the model are known, and this allows computing the temperature field in the body [23,24].
The computation of the temperature field in logs during their thawing in an air environment requires solving of the so-called inverse task of the heat transfer. This is the task when the model of the studied object and the experimentally obtained temperature field in it are known, but one or more variables in the model need to be determined during the solving and validation of the model [24].
The modeling and the multiparameter study of the thawing process of logs in air environment is of considerable scientific and practical interest. For example, as a result of such a study, it is possible to determine the real initial temperature of logs depending on their dimensions, wood species, moisture content, and the temperature of the air near the logs during their many days staying in an open warehouse before the thermal treatment in the production of veneer. The information about the real value of that immeasurable parameter can be used for scientifically based computing of the optimal, energy saving regimes for thermal treatment of each specific batch of logs.
This chapter presents the creation, numerical solving and validation of a two-dimensional nonlinear mathematical model of the transient heat conduction in frozen logs during their thawing at convective boundary conditions in an air environment. A validation of the models towards own experimentally determined 2D temperature distribution in poplar logs with a diameter of 0.24 m, length of 0.48 m, initial temperature of approximately -30°C, and moisture content above the hygroscopic range during their 70 h thawing at room temperature has been carried out.
During the validation of the model, the inverse task has been solved for the determination of the unknown logs' heat transfer coefficients in radial and longitudinal directions. This task has been solved also in regard to the logs' surface temperature, which depends on the mentioned coefficients.

Mathematical model of the 2D temperature distribution in frozen logs during their thawing in an air environment
In [8] the following common form of a model, which describes the 2D nonstationary temperature distribution subjected to thawing frozen logs in an air environment, has been suggested: with an initial condition and boundary conditions for convective heat transfer: • Along the radial coordinate r on the logs' frontal surface during thawing • Along the longitudinal coordinate z on the logs' cylindrical surface during thawing In [8] solutions of Eq. (1) only at conductive boundary conditions for the case of autoclave steaming of logs aimed at their plasticizing in the production of veneer have been realized and graphically presented.
An approach for solving Eq. (1) at much more complicated convective boundary conditions and verification of model (1) to (4) is considered below.

Mathematical description of the thermophysical characteristics of logs
In Figure 1 the three temperature ranges are presented, at which the process of the logs' thawing above the hygroscopic range is carried out, i.e., when u >u fsp .
There thermophysical characteristics of the logs and of both the frozen and nonfrozen free and bound water in them during the separate temperature ranges are also shown. The information on these characteristics is very important for the solving of the model given above.
The mathematical descriptions of the thermal conductivities of nonfrozen wood, λ w-nfr , and frozen wood, λ w-fr , and also of the specific heat capacities of nonfrozen wood, c wÀnfr , frozen wood, c wÀfr , and frozen free and bound water in the wood, c fw and c bwm , have been suggested in [8,9,18] using the experimentally determined data in the dissertations by Kanter [25] and Chudinov [2] for their change as a function of t and u. According to the suggested in [8,9,18] approach, the wood thermal conductivity during thawing of logs with moisture content u above the hygroscopic range can be calculated with the help of the following equations for λ w T, u, ρ b , u fsp À Á : v ¼ 0:1284 À 0:013u: The coefficients γ and β in Eq. (5) are calculated using the next equations: • For nonfrozen wood at u > u 272:15 fsp and 272:15K < T ≤ 423:15K The fiber saturation points of the wood species, u fsp and u 272:15 fsp , are calculated according to following Eqs. [9]: and consequently is the fiber saturation point at T = 272.15 K (i.e., at t = -1°C), kgÁkg À1 . At t = -1°C, the melting of the frozen bound water in the wood is fully completed, and the melting of the free water in the wood starts [22,26].
The effective specific heat capacities of the logs during the pointed three ranges of the thawing process, c we-1,2,3 , which participate in Eq. (1), are equal to the following: Second range : Third range : c we-3 ¼ c wÀnfr : According to the suggested in [8,9,22] mathematical description, the effective specific heat capacities of the logs during their thawing can be calculated with the help of the following equations for c we-1,2,3 T, u, u fsp À Á : 2862u þ 555 The wood density, ρ w , which participates in Eq. (1), is determined above the hygroscopic range according to the following Equation 27]:

Mathematical description of the heat transfer coefficients of logs
For solving of the 2D mathematical model given above, it is needed to have values for the heat transfer coefficients of the logs in radial and longitudinal directions, α wr and α wp , respectively, which participate in Eqs. (3) and (4).
As it was mentioned in the Introduction, the values of α wr and α wp can be computed by solving the inverse task of the heat transfer between the logs and surrounding air environment.
The calculation of α wr and α wp can be carried out with the help of the following equations of the similarity theory, which are valid for the cases of heating of horizontally situated cylindrical bodies in conditions of free air convection [28]: Nu For the usage of Eqs. (19)- (26), it is needed to have a mathematical description of the thermophysical characteristics of the air, λ, β, w, and а, depending on T and φ. The temperature of the air near the logs subjected to thawing during our experiments described below changes in the range from 243.15 to 303.15 K (i.е., from -30°С to 30°С), and φ changes from 40-100% (see Figures 2 and 3).
For the calculation of λ a , β a , w a , and а a , the temperature of the air, Т a , must be used, but for the calculation of w s and а s in Eq. (26), the temperature of the surface of the logs, T s , has to be used.
In the accessible specialized sources, we did not find suitable mathematical descriptions of λ, β, w, and а of the air, depending on T and φ, which could be applied for the precise determination of α wr and α wp according to Eqs. (19)- (26).
Our further study has shown that for solving the inverse task of the heat transfer between the logs and surrounding air, i.e., for the calculation of the heat transfer coefficients of the logs, which participate in the boundary conditions (3) and (4)  • In the radial direction on the cylindrical surface of the logs • In the longitudinal direction on the frontal surface of the logs where x is an exponent, whose values are determined during the solving and validation of the model through the minimization of the root-square-mean error  (RSME) between the calculated model and experimentally obtained results about the change of the temperature fields subjected to thawing logs.

Experimental research of the 2D temperature distribution in poplar logs during their thawing
For solving the inverse task of the heat transfer aimed at validation of the suggested above mathematical model, it is necessary to have experimentally obtained data about the 2D temperature distribution in logs during their thawing. That is why we realized such experiments using poplar (Populus nigra L.) logs with D = 0.24 m, L = 0.48 m, and u > u fsp [26].
In Figure 4 the coordinates of four representative points of the logs, in which the 2D change in the temperature was measured and registered during the logs' thawing, are shown.
For the freezing of the logs before their thawing, a horizontal freezer was used with adjustable temperature range from -1 to -30°C. Sensors Pt100 with long metal casings were positioned in the drilled four holes of the logs. After 50 h separately freezing each logs, the freezer was switched off. Then its lid was opened, and 70 h thawing of the log at room temperature was carried out.
The automatic measurement and record of t m , φ m , and t in the representative points of the logs during the experiments was accomplished by Data Logger type HygroLog NT3 produced by the Swiss firm ROTRONIC AG.
In Figures 2 and 3, the change in the temperature of the processing air medium, t m , and in its humidity, φ m , and also in the temperature in four representative points of two poplar logs, named below as P1 and P2, respectively, during their separate 70 h thawing is presented.
All curves of the experimentally obtained data on these figures are drawn using the licensed software HW4 of the Data Logger. The left coordinate axis on the figures is graduated at % of φ m , and the right one is graduated at°C of t.

Mathematical description of the air medium temperature during logs' thawing
The change shown in Figures 2 and 3 air medium temperature T m during the logs' thawing with correlation 0.98 and root-square-mean error, σ < 1.5°C, has been approximated with the help of the software package Table Curve 2D by the following equation: whose coefficients are equal to: • For log P1: а = 293.3637194, b = À0.00236425, c = À0.69281743.
and τ is the sum of the time of logs' freezing, equal to 50 h = 180,000 s, and the current time of the subsequent thawing of the logs, s.

Numerical solution of the mathematical model of the logs' thawing process
The mathematical descriptions of the thermophysical characteristics of the logs and also of T m considered above were introduced in the mathematical model (1) to (4). An explicit form of the finite-difference method was used for solving of the model without any simplifications [7,8].

Presentation of Eq. (1) of the model
The presentation of Eq. (1) of the model suitable for programming discrete analogue has been carried out using the given in Figures 5 and 6 coordinate system. These figures show the positioning of the knots of the calculation mesh and four representative points, in which the nonstationary 2D distribution of the temperature in the longitudinal section subjected to thawing log has been calculated. The mesh has been built on ¼ of the longitudinal section of the log due to the fact that this ¼ is mirrored symmetrical towards the remaining ¾ of the same section.
Taking into consideration Eqs. (5) and (6), it can be written that and using the coefficient The discrete finite-difference analogue of the left-hand part of Eq. (1), which is suitable for programming in FORTRAN, has the following form [7,30]: Taking into account Eqs. (30), (31), and (32), the discrete analogue of the righthand part of Eq. (1) has the following form: After alignment of Eq. (34) with Eq. (35) and taking into account Eq. (33), at Δz = Δr, it is obtained that Eq. (1) is transformed into the following system of algebraic equations: we-1,2,3 Á ρ w Á Δr 2 : : The term 1 r in the right-hand part of Eq. (1) is represented as 1 iÀ1 ð ÞÁΔr in Eq. (35). According to the requirements of FORTRAN [7,30], the knots of the calculation mesh, which are situated on the log's surfaces, are denoted by numbers i = 1 and k = 1 along the coordinate axes r and z, respectively.
The temperature in these surface knots is calculated with the help of Eqs. (38) and (42) given below. Since Eq. (36) calculates the temperature in the knots, which are located inside the logs, i.e., in the knots with i ≥ 2 and k ≥ 2, the denominator of the term 1 iÀ1 in this equation is always greater than zero. It can be noted that the effective specific heat capacities of the log during the pointed above three ranges of its thawing process (see Figure 1), c we-1 , c we-2 , and c we-3 , which are unitedly represented as c we-1,2,3 in Eq. (36), are computed according to Eq. (17) separately for each knot of the calculation mesh.

Presentation of the equation of the initial condition in the model
The initial condition (2) of the model of logs' thawing process obtains the following discrete finite-difference form: where T w0-avg is the experimentally determined average mass temperature of the log at the beginning of the thawing process, K.

Presentation of the equation of boundary condition in the model along the radial coordinate
The boundary condition (3) of the logs' thawing process obtains the following final form, suitable for programming in FORTRAN: The variable G n i,1 in Eq. (38) is equal to where according to Eqs. (28) and (29) where

Presentation of the equation of boundary condition in the model along the longitudinal coordinate
Analogously, the boundary condition (4) of the logs' thawing process obtains the following final form, suitable for programming in FORTRAN: The variable G n 1,k in Eq. (42) is equal to where according to Eq. (27) and Т n m is calculated according to Eq. (41).

Input data for solving of the model
The numerical solving and verification of the model (1) to (4) has been realized in the calculation environment of Visual FORTRAN Professional.
Using own software package in that environment, computations were carried out for the determination of the 2D nonstationary change of t in the representative points of the logs P1 and P2, whose experimentally registered temperature fields are presented in Figures 2 and 3, respectively.
As it was mentioned above, the duration of the freezing and duration of the subsequent thawing of the logs were equal to 50 and 70 h, respectively.
The model was solved with step Δr = Δz = 0.006 m along the coordinates r and z, with step Δτ = 6 s [8,23], and with the same initial and boundary conditions, as they were during the experimental research.

Inverse determination of the heat transfer coefficients during solving of the model
The model (1) to (4) was solved with various values of the exponent x in Eqs. (27) and (28). The computed by the model change of t in the four representative points of the logs with each of the tested values of the exponent x during the thawing was compared mathematically with the corresponding one experimentally registered change of t in these points with an interval of 15 min.
The aim of this comparison was to determine the values of x, which ensure the best compliance between the computed and experimentally registered temperature fields in subjected to thawing logs.
As a criterion of the best compliance, the minimum average value of RSME, σ avg , was used, which is equal to where t For the calculation of σ avg , a software program in the calculation environment of MS Excel was prepared. At τ thaw = 70 h = 252,000 s, RSME has been calculated with the help of the program simultaneously for a total of NÁP = 1120 temperature-time points during the thawing of each log.
It was determined that the minimum values of RSME overall for the studied four representative points are equal to σ avg = 1.37°C for log P1 and to σ avg = 1.34°C for log P2. These minimum values of σ avg correspond to the following values of the exponent x in Eqs. (27) and (28), which were obtained during the solving of the inverse task, x = 0.22 for log P1 and x = 0.20 for log P2. Figures 7 and 8 present the calculated change in α wr and α wp during the studied thawing process of the logs P1 and P2, respectively. Figures 9 and 10 present the calculated change in t m and also in the logs' surface temperature t s and t of 4 representative points of the studied logs.
It can be seen that with the decrease of the difference between t m and t s during the logs' thawing, the heat transfer coefficients on Figures 7 and 8 gradually decrease, as follows: • At α wr : from 2.3 to 1.0 WÁm À2 ÁK À1 for P1 and from 1.9 to 1.2 WÁm À2 ÁK À1 for P2.
• At α wp : from 5.1 to 2.2 WÁm À2 ÁK À1 for P1 and from 4.5 to 2.8 WÁm À2 ÁK À1 for P2.  Using the obtained change in the heat transfer coefficients, the change in the logs' surface temperature during the thawing, t s , has been calculated by the model (refer to Figures 9 and 10).
The comparison to each other of the analogical curves in Figures 2 and 9, and also in Figures 3 and 10, shows good conformity between the calculated and experimentally determined changes in the very complicated temperature fields of the studied logs during their thawing.
During our extensive simulations with the model (1) to (4), we established good qualitative and quantitative compliance between computed and experimentally determined temperature fields of logs from numerous wood species with different moisture content above the hygroscopic range [31].
The overall RSME for the studied four representative points in the logs does not exceed 5% of the temperature ranges between the minimal and maximal temperatures of each log during its thawing.

Conclusions
This chapter describes the creation, solving, and validation of a 2D nonlinear mathematical model for the transient heat conduction subjected to thawing frozen logs in an air environment.
The mechanism of the heat distribution in logs during their thawing has been described by a 2D equation of heat conduction at convective boundary conditions. For the numerical solving of the model with the help of explicit form of the finitedifference method, a software package has been prepared in the calculation medium of Visual FORTRAN Professional developed by Microsoft.
A validation of the model towards our own experimentally determined 2D temperature distribution in poplar logs with a diameter of 0.24 m, length of 0.48 m, and initial temperature about -30°C during their 70 h separate thawing at room temperature has been carried out.
During the validation of the model, the inverse problem has been solved for the determination of the logs' heat transfer coefficients in radial and longitudinal directions. This problem has been solved also in regard to the logs' surface temperature, which depends on the mentioned coefficients.
The following minimum values of the average RSME total for the temperature change in four representative points in each of the studied logs have been obtained: • σ avg = 1.37°C for log P1 with ρ b = 359 kgÁm À3 and u = 1.44 kgÁkg À1 .
During the solving of the inverse task, it was determined that the heat transfer coefficients subjected to thawing logs decrease gradually, as follows: • At α wr : from 2.3 to 1.0 WÁm À2 ÁK À1 for P1 and from 1.9 to 1.2 WÁm À2 ÁK À1 for P2.
Good adequacy and precision of the model towards the results from extensive own experimental studies allow for the carrying out of various calculations with it, which are connected to the nonstationary temperature distribution in logs during their thawing in an air environment. For example, as a result of such calculations, it is possible to determine the real initial temperature of logs depending on their dimensions, wood species, moisture content, and the temperature of the air near the logs during their many days staying in an open warehouse before the thermal treatment in the production of veneer.
The information about the real value of that immeasurable parameter is needed for scientifically based computing of the optimal, energy saving regimes for thermal treatment of each specific batch of logs.
The model of the logs' thawing process can be applied also in the software for controllers used for advanced model predictive automatic control [20,21,32] of this treatment. The approach for solving of the inverse task of the heat transfer in this chapter could be further applied in the development and solving of analogous models, for example, for the calculation of the temperature fields during freezing or thawing processes of different wooden and other capillary porous materials.