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 −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.
- bound water
- latent heat
- thermal energy
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 freezing—thermal 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.
2. Mathematical model of the 2D temperature distribution in logs subjected to refrigeration
2.1 Mechanism of the temperature distribution in logs during their refrigeration
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 :
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
2.2 Mathematical description of the heat sources in logs during their freezing
The volume internal heat source in the logs, , 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 .
That is why as a methodology for the determination of during the freezing of logs, in , a perspective is used, which is already applied for determination of the volume heat source, , 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 :
A numerical approach and an algorithm for the computation of and are given in . The difference in the right-hand part of Eq. (7) reflects the entire mass of free water (in kg), which is contained in 1 m3 of the logs. The wood densities and , 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 in Eq. (8) is determined according to the following equation in relation to the present entirely liquid quantity of nonfrozen water in the wood, , corresponding to the current wood temperature Т < 272.15 K [2, 10]:
2.3 Introducing the mathematical description of heat sources into the mathematical model
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 in Eq. (16) is equal to
where is the effective specific heat capacity of the wood during freezing only of the free water in it, J·kg−1·K−1, and 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.
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 in Eq. (20) is equal to
where is the effective specific heat capacity of the wood during freezing only of the bound water in it, J·kg−1·K−1, and 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.
2.4 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, , can be determined as a relationship of the mass of the ice formed by the bound water in 1 kg wood, , to the sum of the mass of this ice and of the mass of the nonfrozen water in 1 kg wood at and at T < 272.15 K, , according to the following equation:
where is the standardized fiber saturation point at T = 293.15 K, i.e., at t = 20°C.
The sign “minus” in the right parts of Eq. (25) reflects mathematically the circumstances that the icing degree decreases with an increasing temperature Т.
2.5 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 > . 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 to , only a cooling of the logs with full liquid water in them occurs. During the second range from to , 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 .
During the third range from to , 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, , 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 , , , and and also of the thermal conductivity of nonfrozen and frozen wood, , have been suggested by the first co-author earlier [9, 10] using the experimentally determined data in the dissertations by Kanter  and Chudinov  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:
at cylindrical surface of the horizontally situated logs
at frontal surface of the logs
where 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.
2.6 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 m3 of logs, , can be calculated according to the following equation:
where is the energy needed for the phase transition of the bound water from liquid to solid state in 1 m3 of logs subjected to freezing, kWh·m−3, and is the latent thermal energy of the bound water in 1 m3 of logs released in them during its crystallization, kWh·m−3.
It is known that the energy consumption for heating 1 m3 of wood materials, , with an initial mass temperature to a given average mass temperature is determined using the following equation [2, 5, 10]:
The multiplier 3.6·106 in the denominator of Eq. (37) ensures that the values of are obtained in kWh·m−3, instead of in J·m−3.
Based on Eq. (33), it is possible to calculate the energy according to the equation
and the area of ¼ of the longitudinal section of the log subjected to freezing, , is equal to
Based on Eq. (33), it is possible to calculate also the energy according to the equation
3. 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 > . 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 .
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 , 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:/
In Figure 3, the change in the temperature of the processing air medium, , and in its humidity, , and also in the temperature in 4 characteristic points of pine log named below as P1 with u = 0.33 kg·kg−1 and = 470 kg·m−3 and spruce log named as S1 with u = 0.36 kg·kg−1 and = 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.
4. 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 , 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 kg·kg−1 and volume shrinkage S v = 11.8% and also of spruce sapwood with kg·kg−1 and S v = 11.4% were used [8, 10].
4.1 Mathematical description of T in the freezer during logs’ refrigeration
The curvilinear change in the freezing air medium temperature, , 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  by the following equation:
whose coefficients are = 309.7863391, = 0.007125039, = 1.321533597, and = −2.769·10−6 for log P1 and to = 305.6335660, = 0.005833651, = 1.061216339, and = −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.
4.2 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 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 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 , 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, , 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 = 1.67°C for P1 and = 1.54°C for S1. The minimum values of were obtained with the values of = 0.52 for P1 and = 0.48 for S1 in Eqs. (30) and (31).
Figure 4 presents, as an example, the calculated change in , log’s surface temperature , 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.
4.3 Change of the relative icing degree of logs Ψice-bw-avg and its derivative
Figure 5 presents the calculated change of the logs’ icing degree and the derivative 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 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 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 = 0 and also = 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 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 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.
4.4 Change of the thermal energies Q fr-bw and Q Lat-bw
It can be seen that the change of the energies and is happening according to complex curvilinear dependences on the freezing time. The change of , depending on the freezing time, is similar to that of the icing degree , and the change of is similar to that of the derivative .
At the beginning of the logs’ refrigeration process, when = 0, both energies and are also equal to 0. After that 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 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, = 0.487 for P1 and = 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 for S1 is larger than that of P1, but also the maximal value of = 0.431 kWh·m−3 for S1 is larger than the maximal value of = 0.405 kWh·m−3 for P1. After reaching the maximal values, the energy 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 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 (see Figure 6, left). At the end of the 30 h freezing, the energy reaches the following values: 5.992 kWh·m−3 for P1 and 6.657 kWh·m−3 for S1.
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 m3 wood from liquid into solid state, , and latent thermal energy of the bound water, , which is released in 1 m3 of the logs during water crystallization.
Mathematical descriptions and an approach for computing of the energies and 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 and 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 and and their difference, , 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, represents the energy needed for crystallization of that part of the bound water, which amount depends on , , and .
The logs subjected to refrigeration were with the following combinations of their initial temperature t w0, basic density ρb, and moisture content u: t w0 = 25.2°C, ρb = 470 kg·m−3, and u = 0.33 kg·kg−1 for log P1 and t w0 = 23.5°C, ρb = 479 kg·m−3, and u = 0.36 kg·kg−1 for log S1.
It has been determined that the values of the energies ,, and 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: = 6.224 kWh·m−3, = 0.232 kWh·m−3, and = 5.992 kWh·m−3 for P1 and = 6.925 kWh·m−3, = 0.268 kWh·m−3, and = 6.657 kWh·m−3 for S1.
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.).
specific heat capacity, J·kg−1·K−1
latent heat, J·kg−1, or length, m
internal heat source, W·m−3
thermal energy, kWh·m−3
radial coordinate: 0 ≤ r ≤ R, m
shrinkage, %, or aria of ¼ of log’s longitudinal section, m2
moisture content, kg·kg−1 = %/100
longitudinal coordinate: 0 ≤ z ≤ L/2, m
heat transfer coefficients between log’s surfaces and the surrounding air medium, W·m−2·K−1
thermal conductivity, W·m−1·K−1
root square mean error (RSME), oC
step along the coordinates r and z for solving of the model, m
step along the time coordinate for solving of the model, s
relative icing degree of logs or relative degree of solidification of the metal, —
average (for relative icing degree or for root square mean error)
basic (for wood density, based on dry mass divided to green volume)
end of freezing
fiber saturation point
generalized (for specific heat capacity)
ice (for logs’ icing degrees)
medium (for cooling substance)
volume of the metal
parallel to the wood fibers
total (for specific energy needed for freezing of the bound water)
wood effective (for specific heat capacity)
wood with liquid water in it
wood with solid state of water (ice) in it
wood at u = u fsp
wood at u = u nfw
at T = 272.15 K, i.e., at t = −1°C
at T = 293.15 K, i.e., at t = 20°C