Modeling of the Energy for Bound Water Freezing in Logs Subjected to Refrigeration

This study suggests an approach for modeling of the total thermal energy needed for freezing the bound water in logs subjected to refrigeration. The approach maximally considers the physics of the freezing process of the bound water in wood. It considers the nonstationary change in the icing degree of logs formed by the crystallization of the bound water in them, as well as the influence of the fiber saturation point of each wood species on its current amount of nonfrozen bound water depending on temperatures below (cid:1) 1°C. Mathematical descriptions of the thermal energy of the phase transition of bound water in logs and also of the latent thermal energy of the bound water released in logs during their freezing have been executed. These descriptions were introduced in our own 2D nonlinear mathematical model of the freezing process of logs. The model was transformed in a form suitable for programming with the help of explicit schemes of the finite difference method. For the solution of the model and energy simulations with it, a software program was prepared, which was input into the calculation environment of Visual Fortran Professional.


Introduction
It is known that the duration and the energy consumption of the thermal treatment of frozen logs in winter, aimed at their plasticizing for the production of veneer, depend on the degree of the logs' icing [1][2][3][4][5][6][7][8][9][10].
In the accessible specialized literature, there are limited reports about the temperature distribution in frozen logs subjected to defrosting [8,[11][12][13][14][15][16][17][18][19][20][21], and there is very little information about the research of the temperature distribution in logs during their freezing given by the authors only [22][23][24]. That is why the modeling and the multiparameter study of the freezing process of logs are of considerable scientific and practical interest.
For different engineering, technological, and energy calculations, it is necessary to be able to determine the nonstationary temperature field in logs depending on the temperature of the gaseous or liquid medium influencing them 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 freezing both the free and bound water in the wood. In the specialized literature, information about 2D temperature distribution in logs subjected to freezing was given by the authors [22][23][24]. In the available literature for hydrothermal treatment of frozen wood materials, there is no information at all about the quantitative determination of the energy characteristics of both the free and bound water during their freezing in the wood or in other capillary porous materials.
The aim of the present work is to suggest a methodology for mathematical modeling and research of two mutually connected problems: 2D nonstationary temperature distribution in logs subjected to refrigeration and change in two important energy characteristics of the bound water in logs during its freezingthermal energy of the gradual phase transition of bound water from liquid into solid state and latent thermal energy of the bound water released in the logs at temperatures below À1°C.

Mathematical model of the 2D temperature distribution in logs subjected to refrigeration 2.1 Mechanism of the temperature distribution in logs during their refrigeration
The mechanism of the temperature distribution in logs during their refrigeration can be described by the equation heat conduction [2, 5,[8][9][10].
When the length of the logs is less than their diameter by at least 3-4 times, for the calculation of the change in the temperature in the longitudinal sections of the logs (i.e., along the coordinates r and z of these sections) during their refrigeration in the air medium, the following 2D mathematical model can be used [23]: with an initial condition and boundary conditions for convective heat transfer: • along the radial coordinate r on the logs' frontal surface during the cooling process • along the longitudinal coordinate z on the logs' cylindrical surface during the cooling Equations (1)-(4) represent a common form of a mathematical model of the logs' freezing process.

Mathematical description of the heat sources in logs during their freezing
The volume internal heat source in the logs, q v , reflects in Eq. (1) the influence of the latent heat of water in the wood on the logs' freezing process. In the available literature for hydrothermal treatment of frozen wood materials, no information can be found on the approaches for quantitative determination of the heat source q v .
That is why as a methodology for the determination of q v during the freezing of logs, in [23], a perspective is used, which is already applied for determination of the volume heat source, q vM , during the process of solidification of melted metal [25][26][27][28]. According to this methodology, the following equations for determination of the volume internal sources of latent heat separately of the free and the bound water in wood during the logs' freezing process have been obtained [23]: where L crÀice = 3.34Á10 5 JÁkg À1 [2, 5, 22-24] and A numerical approach and an algorithm for the computation of Ψ ice-fw and Ψ ice-bw are given in [23]. The difference ρ w À ρ wUfsp in the right-hand part of Eq. (7) reflects the entire mass of free water (in kg), which is contained in 1 m 3 of the logs. The wood densities ρ w and ρ wUfsp , which participate in Eqs. (7) and (8), are determined above the hygroscopic range according to the equations below [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]29]: The density of the wood ρ wUnfw in Eq. (8) is determined according to the following equation in relation to the present entirely liquid quantity of nonfrozen water in the wood, u nfw , corresponding to the current wood temperature Т < 272.15 K [2, 10]: For the freezing of the free water in the wood, Eq. (14) is transformed in the following form: Equation (15) can be represented, as follows: The effective specific heat capacity c we-fr-fw * in Eq. (16) is equal to where c weÀfrÀfw is the effective specific heat capacity of the wood during freezing only of the free water in it, JÁkg À1 ÁK À1 , and c LatÀfw is the specific heat capacity, which is formed by the release of the latent heat of the free water during its crystallization in the wood, JÁkg À1 ÁK À1 .
After substituting q v in Eq. (1) by the expression for q v-bw from Eq. (6) and of ∂ψ ice-bw ∂τ in this expression by ∂ψ ice-bw ∂τ Á ∂Т ∂Т , it is obtained that For the case of freezing the bound water in the wood, Eq. (18) is transformed in the following form: Equation (19) can be represented, as follows: The effective specific heat capacity c we-fr-bw * in Eq. (20) is equal to where c weÀfrÀbw is the effective specific heat capacity of the wood during freezing only of the bound water in it, JÁkg À1 ÁK À1 , and c LatÀbwm is the specific heat capacity, which is formed by the release of the latent heat of the maximum possible amount of bound water during its crystallization in the wood, JÁkg À1 ÁK À1 .

Mathematical description of the relative icing degree of logs caused from the freezing of the bound water
The relative icing degree of the logs, which is caused from the freezing of only the bound water in them, Ψ iceÀbw , can be determined as a relationship of the mass of the ice formed by the bound water in 1 kg wood, m ice-bw , to the sum of the mass of this ice and of the mass of the nonfrozen water in 1 kg wood at u ¼ u nfw and at T < 272.15 K, m nfw , according to the following equation: Using experimental data of Stamm [30,31], in [10] the following equation was suggested for the calculation of the fiber saturation point of the wood species, depending on T: where u 293:15 fsp is the standardized fiber saturation point at T = 293.15 K, i.e., at t = 20°C.
After substituting u nfw and u fsp in Eq. (22) by the expressions for u nfw from Eq. (13) and for u fsp from Eq. (23), it is obtained that The sign "minus" in the right parts of Eq. (25) reflects mathematically the circumstances that the icing degree Ψ iceÀbw decreases with an increasing temperature Т.

Mathematical description of the thermophysical characteristics of logs
In Figure 1, the temperature ranges are presented, at which the refrigeration process of logs is carried out when u > u fsp . The thermophysical 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 solving the mathematical model given above.
During the first range from T w0 to T frÀfw , only a cooling of the logs with full liquid water in them occurs. During the second range from T frÀfw to T frÀbwm , a further cooling of the logs occurs until reaching the state needed for starting the crystallization of the free water. During this range also, the phase transition of this water into ice is carried out. The second range is absent when the wood moisture content 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 the state needed for starting the crystallization of the bound water. During this range also, the phase transition of the bound water into ice is gradually performed.
The generalized effective specific heat capacities of the logs, c we-frÀgen , during the above pointed three ranges of the freezing process above the hygroscopic range, are equal to the following: Mathematical descriptions of the specific heat capacities c wÀnfr , c wÀfr , c fw , and c bwm and also of the thermal conductivity of nonfrozen and frozen wood, λ w , have been suggested by the first co-author earlier [9,10] using the experimentally determined data in the dissertations by Kanter [32] and Chudinov [2] for their change as a function of t and u.
Using Eqs. (26)-(28), Eqs. (16) and (20) can be united in the following equation, with the help of which, together with the initial condition (2) and the boundary conditions (3), and (4) the 2D change of the temperature in subjected to freezing logs can be calculated: For the convective boundary conditions (3) and (4) of the logs' freezing process, the following relations for the heat transfer coefficients are most suitable [33]: • at cylindrical surface of the horizontally situated logs • at frontal surface of the logs where E fr is determined during the validation of the nonlinear mathematical model of the freezing process through minimization of the root square mean error (RSME) between the computed by the model and experimentally obtained transient temperature fields in logs subjected to freezing.

An approach for modeling of the energy needed for freezing the bound water in logs
The total energy consumption, which is needed for freezing the bound water in 1 m 3 of logs, Q fr-bw-total , can be calculated according to the following equation: where Q fr-bw is the energy needed for the phase transition of the bound water from liquid to solid state in 1 m 3 of logs subjected to freezing, kWhÁm À3 , and Q Lat-bw is the latent thermal energy of the bound water in 1 m 3 of logs released in them during its crystallization, kWhÁm À3 .
It is known that the energy consumption for heating 1 m 3 of wood materials, Q w , with an initial mass temperature T w0 to a given average mass temperature T wÀavg is determined using the following equation [2, 5, 10]: The multiplier 3.6Á10 6 in the denominator of Eq. (37) ensures that the values of Q w are obtained in kWhÁm À3 , instead of in JÁm À3 .
Based on Eq. (33), it is possible to calculate the energy Q fr-bw according to the equation where according to [8,9] c bwÀavg ¼ 1:8938 Á 10 4 u 272:15 fsp À 0:12 Á exp 0:0567 T wÀfrÀavg À 272:15 and the area of ¼ of the longitudinal section of the log subjected to freezing, S w , is equal to Based on Eq. (33), it is possible to calculate also the energy Q Lat-bw according to the equation where the specific heat capacity, which is formed by the release of the latent heat of the bound water, c LatÀbw , is calculated according to Eq. (21) and the density ρ w is calculated according to Eq. (9).

Experimental research of 2D temperature distribution in logs subjected to freezing
For the validation of the suggested above mathematical model, it is necessary to have experimentally obtained data about the temperature distribution in logs during their freezing. The logs subjected to freezing 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 for the separate wood species. They were produced from the sapwood of a freshly felled pine (Pinus sylvestris L.) and spruce (Picea abies 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 the characteristic points of the log [22].
The coordinates of the characteristic points of the logs are given in Figure 2. These coordinates of the points allow the determination of the 2D temperature distribution in logs during their freezing. For refrigeration of the logs according to the suggested methodology by the authors [22], a horizontal freezer was used with adjustable temperature range from À1 to À30°C.
Sensors Pt100 with long metal casings were positioned in the four drilled 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 were realized by data logger-type HygroLog NT3 produced by ROTRONIC AG (http:/www.rotronic. com). The data logger has software HW4 for graphical presentation of the data.
In Figure 3, 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 pine log named below as P1 with u = 0.33 kgÁkg À1 and ρ b = 470 kgÁm À3 and spruce log named as S1 with u = 0.36 kgÁkg À1 and ρ b = 479 kgÁm À3 during their separate 30 h refrigeration is presented. The record of all data was made automatically by the data logger with intervals of 5 min.

Numerical solution of the mathematical model of the logs' freezing process
For numerical solution of the mathematical model, a software package was prepared in Visual Fortran Professional developed by Microsoft. Using the package, computations were carried out for the calculation of the 2D nonstationary change of t in the characteristic points of ¼ of the longitudinal sections of the studied logs, whose experimentally determined temperature fields are presented in Figure 3.
The model has been solved with the help of explicit schemes of the finite difference method in a way analogous to the one used and described in [9,10,33]. For the computation of the temperature distribution in ¼ of the longitudinal section of the logs, which is symmetrical towards the remaining ¾ of the same section, the model was 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 interval between the time levels, Δτ (i.e., the step along the time coordinate), has been determined by the software package according to the condition of stability for explicit schemes of the finite difference method [10], and in our case it was equal to 6 s.
During the solving of the model, the mathematical descriptions of the thermophysical characteristics of pine wood with u 293:15 fsp ¼ 0:30 kgÁkg À1 and volume shrinkage S v = 11.8% and also of spruce sapwood with u 293:15 fsp ¼ 0:32 kgÁkg À1 and S v = 11.4% were used [8,10].

Mathematical description of T in the freezer during logs' refrigeration
The curvilinear change in the freezing air medium temperature, T mÀfr , which is shown in Figure 1, with high accuracy (correlation 0.98 for the both studied logs and root square mean error (RSME) σ = 1.28°C for P1 and σ = 1.22°C for S1) has been approximated with the help of the software package TableCurve 2D [34] by the following equation: whose coefficients are a fr = 309.7863391, b fr = 0.007125039, c fr = 1.321533597, and d fr = À2.769Á10 À6 for log P1 and to a fr = 305.6335660, b fr = 0.005833651, c fr = 1.061216339, and d fr = À2.275Á10 À6 for log S1. Equation.
(39) and its coefficients were introduced in the software for solving Eqs. (3) and (4) of the model.

Computation of 2D temperature distribution in logs during their refrigeration
The mathematical model of the logs' freezing process has been solved with different values of the exponent Е fr in Еqs. (30) and (31). The calculated by the model change of the temperature in four characteristic points of the longitudinal logs' sections with each of the used values of Е fr during the freezing has been compared mathematically with the corresponding experimentally determined change of t in the same points with an interval of 5 min. The aim of this comparison was to find the value of Е fr , 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. For the determination of RSME, a software program in the calculation environment of MS Excel has been prepared. With the help of the program, RSME simultaneously for a total of 1440 temperature-time points during a separate 30 h refrigeration of the logs has been calculated. During the simulations the same initial and boundary conditions have been used as during the experiments. It was determined that the minimum values of RSME overall for the studied four characteristic points are σ avg = 1.67°C for P1 and σ avg = 1.54°C for S1. The minimum values of σ avg were obtained with the values of Е fr = 0.52 for P1 and Е fr = 0.48 for S1 in Eqs. (30) and (31). Figure 4 presents, as an example, the calculated change in t mÀfr , log's surface temperature t s , and t of four characteristic points of the studied pine log P1.
The comparison to each other of the analogical curves in Figure 3-above and Figure 4 shows good conformity between the calculated and experimentally determined changes in the very complicated temperature field of the pine log during its refrigeration. During our wide simulations with the mathematical model, we observed good compliance between computed and experimentally established temperature fields during refrigeration of logs from different wood species with different moisture content. The overall RSME for the studied four characteristic points in the logs does not exceed 5% of the temperature ranging between the initial and the end temperatures of the logs subjected to refrigeration.

Change of the relative icing degree of logs Ψ ice-bw-avg and its derivative
The calculation of the average mass icing degree of the logs caused from the freezing the bound water, Ψ iceÀbwÀavg , is carried out according to the following equation, which was obtained in [23] using Eq. (24): (40) Figure 5 presents the calculated change of the logs' icing degree Ψ iceÀbwÀavg and the derivative Ψ iceÀbw =dT according to Eqs. (40) and (25), respectively, during the 30 h freezing process of the studied pine and spruce logs. The graphs show that the change of these variables is happening according to complex dependences on the freezing time.
It can be seen that the values of Ψ iceÀbwÀavg increase gradually from 0 to 0.487 for P1 and to 0.498 for S1 at the end of the 30 h freezing (Figure 5-left). These values of Ψ iceÀbwÀavg mean that 1-0.487 = 0.513 relative parts (i.e., 51.3%) of the bound water in P1 and 1-0.498 = 0.502 relative parts (i.e., 50.2%) of the bound water in S1 remains in a liquid state in the cell walls of the wood at the end of the 30th h of the logs' freezing process when the calculated average logs' mass temperatures are equal to À27.59°C for P1 and À26.82°C for S1.
It can be pointed that during the first 2.42 h for P1 and 2.75 h for S1 of the freezing process the whole amount of the bound water in the logs is in a liquid state and because of that the icing degree Ψ iceÀbwÀavg = 0 and also Ψ iceÀbw =dT= 0 K À1 . From 2.42th h for P1 and from 2.75th h for S1 the crystallization of the bound water in the peripheral layers of the logs starts. This causes а jump in the change of Ψ iceÀbw =dT from 0 to maximal values, equal to 0.0056 K À1 for P1 and 0.0057 K À1 for S1.
These maximal values remain practically unchanged until reaching 10.00th h for P1 and 10.75th h for S1, when the crystallization of the whole amount of the free water in the logs ends and the freezing of the bound water even in the logs' center starts.
The further decrease of the temperature in the freezer (refer to Figures 3 and 4) causes a gradual crystallization of the bound water in the logs. Because of that the derivatives Ψ iceÀbw =dT decrease slowly, and at the end of the 30 h freezing, they reach the value of 0.0013 K À1 for P1 and 0.0014 K À1 for S1. It can be seen that the change of the energies Q fr-bw and Q Lat-bw is happening according to complex curvilinear dependences on the freezing time. The change of Q fr-bw , depending on the freezing time, is similar to that of the icing degree Ψ iceÀbwÀavg , and the change of Q Lat-bw is similar to that of the derivative Ψ iceÀbw =dT.

Change of the thermal energies Q fr-bw and Q Lat-bw
At the beginning of the logs' refrigeration process, when Ψ iceÀbwÀavg = 0, both energies Q fr-bw and Q Lat-bw are also equal to 0. After that Q fr-bw increases gradually from 0 to 6.224 kWhÁm À3 for P1 and 6.925 kWhÁm À3 for S1 at the end of the 30 h freezing (see Figure 6, left). The larger value of Q fr-bw for S1 is caused from the larger amount of the frozen bound water in S1 in comparison with P1 at the end of the freezing process.
As it was mentioned above, the fiber saturation point of the pine wood is equal to 0.30 kgÁkg À1 , and of the spruce wood, it is 0.32 kgÁkg À1 [8,10]. According to Figure 5, left, the relative icing degree, Ψ iceÀbwÀavg = 0.487 for P1 and Ψ iceÀbwÀavg = 0.498 for S1 at the end of the 30 h freezing. This means that the amount of the crystallized bound water is equal to 0.487 Â 0.30 = 0.146 kgÁkg À1 in P1 and it is equal to 0.498 Â 0.32 = 0.159 kgÁkg À1 in S1. Because of this not only the value of Q fr-bw for S1 is larger than that of P1, but also the maximal value of Q Lat-bw = 0.431 kWhÁm À3 for S1 is larger than the maximal value of Q Lat-bw = 0.405 kWhÁm À3 for P1. After reaching the maximal values, the energy Q Lat-bw decreases gradually and at the end of the 30 h freezing obtains a value of 0.232 kWhÁm À3 for P1 and of 0.268 kWhÁm À3 for S1.  Figure 7 presents the change of the thermal energy Q fr-bw-total during the 30 h refrigeration of the studied logs, which is calculated according to Eq. (32). The change of this energy is similar to that of the energy Q fr-bw (see Figure 6, left). At the end of the 30 h freezing, the energy Q fr-bw-total reaches the following values: 5.992 kWhÁm À3 for P1 and 6.657 kWhÁm À3 for S1.

Discussions
This paper presents a methodology for mathematical modeling and research of two mutually connected problems: 2D nonstationary temperature distribution in logs subjected to refrigeration and change in two important energy characteristics of the bound water in logs during its freezing-thermal energy of the phase transition of the bound water in 1 m 3 wood from liquid into solid state, Q fr-bw , and latent thermal energy of the bound water, Q Lat-bw , which is released in 1 m 3 of the logs during water crystallization.
Mathematical descriptions and an approach for computing of the energies Q fr-bw and Q Lat-bw during the freezing of the bound water at temperatures below À1°C have been carried out. These descriptions are introduced in our own 2D nonlinear mathematical model of the 2D heat distribution in logs during their refrigeration at convective boundary conditions. The model was transformed in a form suitable for programming with the help of explicit schemes of the finite difference method, which excludes the necessity of any simplifications of the model.
A software program for numerical solution of the mathematical model and computation of 2D nonstationary change of the temperature in logs subjected to refrigeration and of the thermal energies Q fr-bw and Q Lat-bw has been prepared in Fortran, which has been input in the calculation environment of Visual Fortran Professional developed by Microsoft.
With the help of the program, computations for the determination of the energies Q fr-bw and Q Lat-bw and their difference, Q fr-bw-total ¼ Q fr-bw À Q Lat-bw , have been completed as an example for the case of a pine log P1 and a spruce log S1 with a diameter of 0.24 m and length of 0.48 m subjected to 30 h refrigeration in a freezer at approximately À30°C. Practically, Q fr-bw-total represents the energy needed for crystallization of that part of the bound water, which amount depends on u fsp , t mÀfr , and τ fr .
It has been determined that the values of the energies Q fr-bw ,Q Lat-bw , and Q fr-bw-total of the studied logs change according to complex relationships depending on the freezing time and after 30 h freezing of the logs reach the following values: Q fr-bw = 6.224 kWhÁm À3 , Q Lat-bw = 0.232 kWhÁm À3 , and Q fr-bw-total = 5.992 kWhÁm À3 for P1 and Q fr-bw = 6.925 kWhÁm À3 , Q Lat-bw = 0.268 kWhÁm À3 , and Q fr-bw-total = 6.657 kWhÁm À3 for S1.

Conclusions
Good adequacy and precision of the model toward the results from wide own experimental studies allow the carrying out of various calculations with the model, which are connected to the nonstationary temperature distribution and energy characteristics of logs from different wood species during their refrigeration. The mathematical model, after its connection with other our model of the logs' defrosting process [9,10], could be input into the software of programmable controllers for optimized model-based automatic control [8,20,21,35] of thermal treatment of frozen logs in the production of veneer.
The approach for the computation of the thermal energies of the bound water in logs during their refrigeration could be used for the creation of analogous models for the computation of the temperature distribution and the energy required for the refrigeration of different capillary porous materials (fruits, vegetables, meet, meet products, etc.).