Numerical Solution to Two-Dimensional Freezing and Subsequent Defrosting of Logs

Two-dimensional mutually connected mathematical models have been created, solved, and verified for the transient non-linear heat conduction in logs during their freezing and subsequent defrosting. The models reflect the influence of the internal sources of latent heat of both the free and bound water on the logs ’ freezing process and also the impact of the temperature on the fiber saturation point of wood species, with whose participation the current values of the thermo-physical characteristics in each separate volume point of the subjected to freezing and subsequent defrosting logs are computed. The chapter presents solutions of the models with explicit form of the finite-difference method and their validation towards own experimental studies. Results from experimental and simulative investigation of 2D non-stationary temperature distribution in the longitudinal section of beech and pine logs with a diameter of 0.24 m and length of 0.48 m during their many hours freezing in a freezer and subsequent defrosting at room temperature are presented, visualized, and analyzed.


Introduction
It is known that the duration of the thermal treatment of the frozen logs in the winter aimed at their plasticizing for the production of veneer and the very high energy consumption needed for this treatment depend on the degree of the logs' icing [1][2][3][4][5][6][7][8][9][10]. For example, for the defrosting and plasticizing of beech and oak logs with an initial temperature of À10°C and moisture content of 0.6 kgÁkg À1 approximately 68 and 81 kWhÁm À3 thermal energy respectively [10] is needed.
In the specialized literature there are limited reports about the temperature distribution in subjected to defrosting frozen logs [8,[11][12][13][14][15][16][17][18][19][20][21][22] and there is very little information about research of the temperature distribution in logs during their freezing given by the authors only [23][24][25]. That is why the modeling and the multi-parameter study of the mutually connected freezing and defrosting processes of logs are of considerable scientific and practical interest.
Considerable 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 at conductive boundary conditions has been made by Steinhagen [11,12] and later one-dimensional and two-dimensional models have been created and solved [13][14][15][16][17]. The thermal energy, which is needed for the melting of the ice, which has been formed from the freezing of the bound water in the wood, has not been taken into account in the models of cited references.
The models assume that the fiber saturation point is identical for all wood species (i.e., u fsp = 0.3 kgÁkg À1 = const). However, it is known that there are significant differences between the fiber saturation points of the different wood species [1][2][3][4][5][6][7][8]. The indicated complications and incompleteness in these models have been overcome in the suggested by Deliiski [9,18,19] 2-dimensional mathematical model of the transient non-linear heat conduction in frozen logs during their heating and defrosting.
For different engineering calculations it is needed to be able to determine the icing degree of the wood materials depending on the temperature of the influencing on them gas or liquid medium and on the duration of their staying in this medium. Such calculations are carried out using mathematical models, which describe adequately the complex processes of the freezing of both the free and bound water in the wood.
The computer solutions of these models give the non-stationary distribution of the temperature in the materials during their cooling below temperatures, at which a freezing of the whole amount of the free water and a freezing of respective, depending on the temperature, part of the bound water in the wood occurs [24,25].
This chapter presents the creation, numerical solving and validation of two mutually connected two-dimensional mathematical models of the transient nonlinear heat conduction in logs during their freezing and subsequent defrosting at convective boundary conditions.
The model of the freezing process takes into account for the first time the impact of the internal sources of latent heat of both the free and bound water on the temperature distribution. The both models reflect the impact of the temperature on the fiber saturation point of each wood species, with whose participation the current values of the thermo-physical characteristics in each separate volume point of the subjected to freezing and subsequent defrosting logs are computed.
The chapter also presents and visualizes the results from experimental and simulative investigation of the 2D non-stationary temperature distribution in the longitudinal section of beech and poplar logs with a diameter of 0.24 m, length of 0.48 m, and different moisture content during many hours their freezing in a freezer and subsequent defrosting at curvilinear changing temperature of the processing air medium.
2. Mechanism of the temperature distribution in logs during their freezing and subsequent defrosting 2.1 Mathematical model of the 2D temperature distribution in logs subjected to freezing When the length of the logs, L, is larger than their diameter, D, by not more than 3-4 times, for the calculation of the change in the temperature in the logs' longitudinal sections (i.e., along the coordinates r and z of these sections) during their freezing in air medium the following 2D mathematical model can be used [24]: with an initial condition T r; z; 0 ð Þ¼T w0 (2) and boundary conditions for convective heat transfer: • along the radial coordinate r on the logs' frontal surface during the freezing process (refer to Figure 3): • along the longitudinal coordinate z on the logs' cylindrical surface during the freezing: Eqs.
(1)-(4) represent a common form of a mathematical model of the logs' freezing process, i.e., of the 2D temperature distribution in logs subjected to freezing.

Mathematical model of the 2D temperature distribution in logs subjected to defrosting
In cases when the length of the logs, L, is larger than their diameter, D, by not more than 3-4 times, for the calculation of the change in T in the longitudinal sections of the logs during their defrosting in air processing medium the following 2D mathematical model can be used [8,9]: with an initial condition and boundary conditions for convective heat transfer: • along the radial coordinate r on the logs' frontal surface during the defrosting process: • along the longitudinal coordinate z on the logs' cylindrical surface during the defrosting process: Eqs. (5)-(8) represent a mathematical model of the logs' defrosting process, i.e., of the 2D temperature distribution in logs subjected to defrosting immediately after their freezing.

Mathematical description of the thermo-physical characteristics of logs
On Figures 1 and 2 the temperature ranges are presented, at which the processes of the logs' freezing and subsequent defrosting respectively are carried out when u > u fsp .
There thermo-physical characteristics of the logs and of both the frozen free and bound water in them have also been shown. The information on these characteristics is very important for the solving of the mutually connected mathematical models given above.
During the first range of the logs' freezing process from T w0 to T frÀfw only a cooling of the logs with fully liquid water in them occurs ( Figure 1). During the second range from T frÀfw to T frÀbwm a further cooling of the logs occurs until reaching of the state needed for starting of the crystallization of the free water. During that range also the phase transition of this water into ice is carried out. The second range is absent when u is less than u fsp . During the third range from T frÀbwm to T wÀfreÀavg a further cooling of the logs is carried out until reaching of the state needed for starting of the crystallization of the bound water. During this range also the phase transition of the bound water into ice is performed.
During the first range of the logs' defrosting process from T wÀfreÀavg to T dfrÀbwm a heating of the frozen logs is carried out until reaching of the state needed for starting and realization of the gradually melting of the frozen bound water in them  ( Figure 2). During the second range from T dfrÀbwm to T dfrÀfw a further heating of the logs occurs until reaching of the state needed for starting and realization of the melting of the frozen free water in them. During the third range from T dfrÀfw to T wÀdfreÀavg a further heating of the logs with fully liquid water in them occurs.
The effective specific heat capacities of the logs during the pointed 3 ranges of the freezing process (see Figure 1) above the hygroscopic range, c we-fr , are equal to the below: The effective specific heat capacities of the logs during the pointed 3 ranges of the defrosting process (see Figure 2) above the hygroscopic range, c we-dfr , are equal to the below: Mathematical descriptions of the specific heat capacities c wÀnfr , c wÀfr , c fw , c bwm , and also of the thermal conductivities of non-frozen, λ w-nfr , and frozen wood, λ w-fr , have been suggested by the first co-author earlier [9,10,19] using the experimentally determined in the dissertations by Kanter [26] and Chudinov [2] data for their change as a function of t and u. These relations are used in both the European [5][6][7][8][9][10] and the American specialized literature [11][12][13][14][15][16][17] when calculating various processes of the wood thermal treatment.
Mathematical descriptions of the specific heat capacities, which are formed by the release of the latent heat of both the free and bound water during their crystallization in the wood, c LatÀfw and c LatÀbwm respectively (refer to Eqs. (10) and (11)) have been given by the authors in [23]. Mathematical descriptions of the internal heat sources separately for the free and bound water, q v-fw and q v-bw respectively, which participate in the right-hand part of Eq. (1), have been also given.
Our study has shown that for the calculation (in WÁm À2 ÁK À1 ) of the radial and longitudinal transfer coefficients of the logs, which participate in the boundary conditions of the model the following equations are most suitable [27]: • in the radial direction on the cylindrical surface of the horizontally situated logs: • in the longitudinal direction on the frontal surface of the logs: where Е fr and Е dfr are exponents, whose values are determined during the solving and verification of the models through minimization of the root square mean error (RSME) between the calculated by the models and experimentally obtained results about the change of the temperature fields in subjected to freezing and subsequent defrosting logs.

Transformation of the models to a form suitable for programming
Due to the correct cylindrical shape of the logs, for the solution of the above presented two mutually connected mathematical models an explicit form of the finite-difference method has been used. That form allows for the exclusion of any simplifications of the models and also of the necessity for any linearization of the mathematically described variables in them [8,9,18,19].
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 [9].
According to the 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 derivatives ∂T ∂r , ∂T ∂z and ∂T ∂τ are approximated using the built computational mesh along the spatial and time coordinates through their finite-difference (discrete) analogues.

Transformation of the equations of thermo-conductivity in the models
The transformation of Eqs. (1) and (5) of the models in their suitable for programming discrete analogues has been realized using the presented on Figure 3 coordinate system. That system shows the positioning of the knots of the calculation mesh, in which the non-stationary distribution of the temperature in the longitudinal section of subjected to freezing and subsequent defrosting log has been calculated. The calculation mesh has been built on ¼ of the longitudinal section of the log due to the circumstance that this ¼ is mirror symmetrical towards the remaining 3/4 of the same section. Taking into consideration the relationships [9,19] and using the coefficient after applying the explicit form of the finite-differences method to Eq. (1), it is transformed into the following system of equations: where λ wor and λ wop are the wood thermal conductivities at t = 0°C in radial and longitudinal direction respectively, WÁm À1 ÁK À1 ; γ, β, and ν-coefficients for the determination of λ w , λ wor and λ wop . Equation for calculation of γ, β, and ν for nonfrozen and frozen wood are given in [9,19]; K wr and K wp -empirically determined coefficients, which take into account the influence of the radial and longitudinal anatomical directions on the wood thermal conductivity, -; Δτ-interval between time levels, i.e., step along the time coordinate, with which the model is solving, s; n-current number of the step Δτ along the time coordinate: n = 0, 1, 2, 3,… τ fr /Δτ; τ fr -duration of the freezing process of logs, s; Δr = Δz-step along the coordinates r and z, by which the model is being solved, m; i-current number of the knots of the calculation mesh along the coordinate r: 1 ≤ i ≤ 1 + (R/Δr); k-current number of the knots along the coordinate z: The effective specific heat capacities of the log during the pointed above three ranges of its freezing process (see Figure 1), c we-fr1 , c we-fr2 , and c we-fr3 , which are unitedly represented as c we-fr1, 2, 3 in Eq. (23), are computed according to Eqs. (9)-(11) separately for each knot of the calculation mesh.
After applying the explicit form of the finite-differences method to Eq. (5), it is transformed into the following system of equations, which is similar to Eq. (23): we-dfr1, 2, 3 Á ρ w Á Δr 2 : : where the effective specific heat capacities of the log during the pointed above three ranges of its defrosting process (see Figure 2), c we-dfr1 , c we-dfr2 , and c we-dfr3 , which are unitedly represented as c we-dfr1, 2, 3 in Eq. (25), are computed according to Eqs. (12)-(14) separately for each knot of the calculation mesh.

Transformation of the equations of initial conditions in the models
The initial condition (2) of the model of logs' freezing process obtains the following discrete finite-difference form: The initial condition (6) of the model of logs' defrosting process obtains the following discrete finite-difference form: where τ fr is the duration of the freezing process of logs, which precedes the beginning of their subsequent defrosting, s.

Transformation of the equations of boundary conditions in the models along the radial coordinate
The discrete finite difference analogue of the left-hand part of Eq. (3), which is suitable for programming in FORTRAN, has the following form [8,9]: Δr : The discrete analogue of the right-hand part of Eq. (3) has the following form: where according to Eq. (17) After alignment of Eq. (28) with Eq. (29), it is obtained that During the solving of the model it is needed to determine the temperature on the log's frontal surface, T nþ1 i, 1 , for each next step n + 1 along the time coordinate. Using Eq. (31), T nþ1 i, 1 is equal to After designation of the boundary condition (3) of the logs' freezing process obtains the following final form, suitable for programming: Analogously, the boundary condition (7) of the logs' defrosting process obtains the following final form, suitable for programming: The variable G n i, 1Àdfr in Eq. (35) is equal to where according to Eq. (18)

Transformation of the equations of boundary conditions in the models along the longitudinal coordinate
The discrete finite difference analogue of the left-hand part of Eq. (4) at Δz = Δr has the following suitable for programming form [8,9]: The discrete analogue of the right-hand part of Eq. (4) has the following form: where according to Eq. (15) α n wr-fr ¼ 2:56 T n 1, k À T n mÀfr Â Ã E fr : After alignment of Eq. (38) with Eq. (35), it is obtained that Using Eq. (41), the temperature T nþ1 1, k is equal to After designation of the boundary condition (4) of the logs' freezing process obtains the following final form, suitable for programming: Analogously, the boundary condition (8) of the logs' defrosting process obtains the following final form, suitable for programming in FORTRAN: The variable G n 1, kÀdfr in Eq. (45) is equal to where according to Eq. (16) α n wr-dfr ¼ 2:56 T n 1, k À T n mÀdfr Â Ã E dfr : 3. Experimental research of the logs' freezing and defrosting

Experimental research of the 2D temperature distribution in beech and poplar logs during their freezing and subsequent defrosting
For the validation of the suggested above mathematical models, it is necessary to have experimentally obtained data about the 2D temperature distribution in logs during their freezing and subsequent defrosting. That is why we carried out such experiments.
The logs subjected to freezing and subsequent defrosting in our experimental research were with a diameter of 240 mm, length of 480 mm, and u > u fsp . This means, that the logs contained the maximum possible amount of bound water. They were produced from the sap-wood of a freshly felled beech (Fagus sylvatica L.) and poplar (Populus nigra L.) trunks.
Before the experiments, four holes with diameters of 6 mm and different lengths were drilled in each log parallel to its axis until reaching of the characteristic points of the log [23].
The coordinates of the characteristic points of the logs are given on Figure 4. These coordinates of the points allow for the determination of the 2D temperature distribution in logs during their freezing and subsequent defrosting.
For the freezing of the logs according to the suggested methodology by the authors [23], 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. The automatic measurement and record of t m , φ m , and t in the characteristic points of the logs during the experiments was realized by Data Logger type HygroLog NT3 produced by the Swiss firm ROTRONIC AG (http:/www.rotronic.c om). The Data Logger has software HW4 for graphical presentation of the experimentally obtained data. After reaching of about À28°C in the log's center during the freezing, the freezer was switch off. Then its lid was opened and a defrosting of the log at room temperature was carried out.
In Figures 5 and 6 the change in the temperature of the processing air medium, t m , and in its humidity, φ m , and also in the temperature in 4 characteristic points of a beech and of a poplar log respectively during their separately 50 and 70 h is presented. The record of all data was made automatically by Data Logger with  The initial temperature, t w0 , basic density, ρ b , moisture content, u, duration of the freezing, τ fr , and duration of the subsequent defrosting, τ dfr , of the logs during the experiments were as follow: • for the beech log: t w0 = 22.4°C, ρ b = 684 kgÁkg À1 , u = 0.63 kgÁkg À1 , τ fr = 50 h, and τ dfr = 50 h; • for the poplar log: t w0 = 19.8°C, ρ b = 359 kgÁkg À1 , u = 1.44 kgÁkg À1 , τ fr = 50 h, and τ dfr = 70 h.  3.2 Mathematical description of the air medium temperature during logs' freezing and subsequent defrosting The change in the shown on Figures 5 and 6 freezing and defrosting air medium temperatures, Т mÀfr and Т mÀdfr , with very high accuracy (correlation 0.98 and Root Square Mean Error (RSME), σ < 1.0°C) has been approximated with the help of the software package Table Curve 2D (http://www.sigmaplot.co.uk/products/table curve2d/tablecurve2d.php) by the following equations: • during the freezing of the beech log: whose coefficients are equal to: а fr = 294.3352069, b fr = 0.010648218, c fr = 2.468350514; • during the freezing of the poplar log: whose coefficients are equal to: а fr = 301.8210985, b fr = 0.0004515197, c fr = 0.111484207, d fr = À6.6585073Á10 À9 , e fr = À1.6653Á10 À6 , f fr = 2.52712Á10 À14 , g fr = 6.46801Á10 À12 , h fr = 2.94924Á10 À21 ; • during the defrosting of both the beech and poplar logs: whose coefficients are equal to as follows: а dfr = 297.1420433, b dfr = À0.00237763, c dfr = À0.70526837 for the beech log; а dfr = 296.3637194, b dfr = À0.00236425, c dfr = À0.69281743 for the poplar log. Eqs. (48) and (49) were used for the solving of Eqs. (3) and (7) and Eq. (50) was used for the solving of Eqs. (4) and (8) of the model.

Numerical solution of the mathematical models of the logs' freezing and subsequent defrosting processes
For the numerical solution of the mutually connected mathematical models, a software package was prepared in Visual FORTRAN Professional developed by Microsoft. Using the package, computations were carried out for the determination of the 2D non-stationary change of t in the characteristic points of ¼ of the longitudinal sections of the studied logs, whose experimentally determined temperature fields are presented on Figures 3 and 4.
The models have been solved with the help of explicit schemes of the finite difference method in a way, analogous to the one used and described in [8-10, 18, 19]. For the computation of the temperature distribution in ¼ of the longitudinal section of the logs, which is mirror symmetrical towards the remaining 3/4 of the same section, the models were solved with step Δr = Δz = 0.006 m along the coordinates r and z and with the same initial and boundary conditions, as they were during the experimental research.
The mathematical models of the logs' freezing and subsequent defrosting processes have been solved with different values of the exponents Е fr and Е dfr in Eqs. (15)- (18). The calculated by the models change of the temperature in the four characteristic points of the longitudinal logs' sections with each of the used values of Е fr and Е dfr during the freezing and defrosting has been compared mathematically with the corresponding one experimentally determined change of t in the same points with an interval of 15 min. The aim of this comparison was to find that the values of Е fr and Е dfr , which ensures the best qualitative and quantitative compliance between the calculated and experimentally determined temperature fields in the logs' longitudinal sections.
As a criterion of the best compliance between the compared values of the temperature total for the four characteristic points the minimum value of RSME, σ avg , has been used.
The average value of RSME, σ avg , has been calculated according to the following equation: where t calc p, n and t exp p, n are the calculated and experimentally established temperatures in the characteristic points; p-number of the characteristic points in the log's longitudinal section: р = 1, 2, 3, 4, i.e., Р = 4 is inputted into Eq. (51); n-number of the moments of the freezing and subsequent defrosting processes: n = 1, 2, 3,,….., N = (τ fr +τ dfr )/(150Δτ) because of the circumstance that the comparison of the calculated values of t with experimentally determined values in the same points has been made with an interval of 15 min = 900 s = 150Δτ.
For the determination of σ avg software program in the calculation environment of MS Excel has been prepared. At τ fr +τ dfr = 100 h for the beech log and at τ fr +τ dfr = 120 h for the poplar log, RSME has been calculated with the help of the program simultaneously for total (N + 1)ÁP = 1604 temperature-time points for the beech log and for total 1924 such points for the poplar log during their freezing and subsequent defrosting. It was determined that the minimum values of RSME overall for the studied 4 characteristic points are equal to σ avg = 1.29°C for the beech log and to σ avg = 1.50°C for the poplar log.
The minimum values of σ avg were obtained with the following values of the exponents in Eqs. (15)-(18): • for the beech log: Е fr = 0.52 and Е dfr = 0.32; • for the poplar log: Е fr = 0.43 and Е dfr = 0.22. Figures 7 and 8 presents the calculated change in t mÀfr and t mÀdfr , which are represented unitedly as t m , and also in logs' surface temperature t s and t of 4 characteristic points of the studied beech and poplar logs respectively.
The comparison to each other of the analogical curves on Figures 5 and 7, and also on Figures 6 and 8 shows good conformity between the calculated and experimentally determined changes in the very complicated temperature fields of the studied logs during their freezing and subsequent defrosting.
During our wide simulations with the mathematical models, we observed good compliance between computed and experimentally established temperature fields of logs various wood species with different moisture content.
The overall RSME for the studied 4 characteristic points in the logs does not exceed 5% of the temperature ranges between the initial and the end temperatures of the logs subjected to freezing or subjected to subsequent defrosting.

Conclusions
This chapter describes the creation, solution, and validation of two mutually connected 2D non-linear mathematical models for the transient heat conduction in  subjected to freezing and subsequent defrosting logs with any u ≥ u fsp . The model of the freezing process takes into account the impact of the internal sources of latent heat of both the free and bound water on the temperature distribution. The both models reflect the impact of the temperature on the fiber saturation point of each wood species, with whose participation the current values of the thermo-physical characteristics in each separate volume point of the subjected to freezing and subsequent defrosting logs are computed.
The mechanism of the temperature distribution in the longitudinal section of the logs during their freezing and subsequent defrosting has been mathematically described by 2D equations of heat conduction. Boundary conditions for convective heat transfer have been implemented in the models. For the transformation of both models into discrete analogues, which are suitable for programming, an explicit form of the finite-difference method has been used.
For the numerical solution of the discrete analogues of the models a software package has been prepared using the programming language FORTRAN, which has been input in the calculation environment of Visual Fortran Professional.
A validation of the models towards own experimentally determined 2D temperature distribution in beech and poplar logs with a diameter of 0.24 m, length of 0.48 m during their separate 50 h freezing in a freezer and many hours subsequent defrosting at room temperature has been carried out. The influence of the curvilinear changing temperature of the air medium in the freezer until reaching of approximately À30°C and also of the air processing medium during the defrosting of logs has been investigated. The following minimum values of the average RSME total for the temperature change in four characteristic points in each of the logs have been obtained: • σ avg = 1.29°C for the beech log with t wo = 22.4°C, ρ b = 684 kgÁm À3 , and u = 0.63 kgÁkg À1 ; • σ avg = 1.50°C for the poplar log with t wo = 19.8°C, ρ b = 359 kgÁm À3 , and u = 1.44 kgÁkg À1 .
During our experimental research it has been determined that in situated on the logs' inner layers characteristic points the specific practically horizontal sections of retention of the temperature for many hours in the range from 0 to À1°С arise, while in these layers a complete freezing of the free water occurs. Analogous retention of the temperature in the range from À1 to 0°C arises during the logs' defrosting. The further the point is distanced from the logs' surfaces and the larger the amount of the free water in the wood is, that much more these sections with temperature retention are extended. Our simulations show that this phenomenon of the freezing and defrosting processes has been correctly reflected in the models (see Figures 7 and 8).
Good adequacy and precision of the models towards the results from wide own experimental studies allow for the carrying out of various calculations with the models, which are connected to the non-stationary temperature distribution in logs from different wood species during their freezing and subsequent defrosting.
The validation of the models with curvilinear change in the temperature of both the freezing and defrosting air mediums will allow us in the future to solve the models with curvilinear changing of the climate temperature [29] over many winter days and nights. The solution will allow for the calculation of the temperature distribution, icing degrees from both the free and bound water, and also different energy characteristics of logs for each desired moment of the their freezing and subsequent defrosting.
The mutually connected models of the freezing and defrosting processes can be applied for the development of scientifically based and energy saving optimized regimes for thermal treatment of frozen logs and also in the software for controllers used for model predictive automatic control [21,22,30]