## Abstract

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.

### Keywords

- logs
- refrigeration
- bound water
- freezing
- latent heat
- thermal energy

## 1. 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 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

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.

### 2.2 Mathematical description of the heat sources in logs during their freezing

The volume internal heat source in the logs,

That is why as a methodology for the determination of

where ^{5} J·kg^{−1} [2, 5, 22, 23, 24] and

A numerical approach and an algorithm for the computation of ^{3} of the logs. The wood densities

The density of the wood *Т* < 272.15 K [2, 10]:

where

### 2.3 Introducing the mathematical description of heat sources into the mathematical model

After substituting

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

where ^{−1}·K^{−1}, and ^{−1}·K^{−1}.

After substituting

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

where ^{−1}·K^{−1}, and ^{−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, *T* < 272.15 K,

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 *T* = 293.15 K, i.e., at *t* = 20°C.

After substituting

and

The sign “minus” in the right parts of Eq. (25) reflects mathematically the circumstances that the icing degree *Т*.

### 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* >

During the first range from

During the third range from

The generalized effective specific heat capacities of the logs,

Mathematical descriptions of the specific heat capacities *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] when calculating various processes of the wood thermal treatment.

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 *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 m^{3} of logs,

where ^{3} of logs subjected to freezing, kWh·m^{−3}, and ^{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,

The multiplier 3.6·10^{6} in the denominator of Eq. (37) ensures that the values of ^{−3}, instead of in J·m^{−3}.

Based on Eq. (33), it is possible to calculate the energy

and the area of ¼ of the longitudinal section of the log subjected to freezing,

Based on Eq. (33), it is possible to calculate also the energy

where the specific heat capacity, which is formed by the release of the latent heat of the bound water,

## 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* > *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 (

In Figure 3, the change in the temperature of the processing air medium, *u* = 0.33 kg·kg^{−1} and ^{−3} and spruce log named as S1 with *u* = 0.36 kg·kg^{−1} and ^{−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 [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 ^{−1} and volume shrinkage *S*_{v} = 11.8% and also of spruce sapwood with ^{−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,

whose coefficients are ^{−6} for log P1 and to ^{−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 *t* in the same points with an interval of 5 min. The aim of this comparison was to find the value of

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,

Figure 4 presents, as an example, the calculated change in *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

The calculation of the average mass icing degree of the logs caused from the freezing the bound water,

Figure 5 presents the calculated change of the logs’ icing degree

It can be seen that the values of

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 ^{−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 ^{−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 ^{−1} for P1 and 0.0014 K^{−1} for S1.

### 4.4 Change of the thermal energies *Q*_{fr-bw} and *Q*_{Lat-bw}

Figure 6 presents according to Eqs. (34) and (38) the calculated change of the thermal energies

It can be seen that the change of the energies

At the beginning of the logs’ refrigeration process, when ^{−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

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, ^{−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 ^{−3} for S1 is larger than the maximal value of ^{−3} for P1. After reaching the maximal values, the energy ^{−3} for P1 and of 0.268 kWh·m^{−3} for S1.

Figure 7 presents the change of the thermal energy ^{−3} for P1 and 6.657 kWh·m^{−3} for S1.

## 5. 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, ^{3} of the logs during water crystallization.

Mathematical descriptions and an approach for computing of the energies

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

With the help of the program, computations for the determination of the energies

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 ^{−3}, ^{−3}, and ^{−3} for P1 and ^{−3}, ^{−3}, and ^{−3} for S1.

## 6. 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.).

## Nomenclature

c | specific heat capacity, J·kg−1·K−1 |

D | diameter, m |

E | exponent, — |

L | latent heat, J·kg−1, or length, m |

q | internal heat source, W·m−3 |

Q | thermal energy, kWh·m−3 |

R | radius, m |

r | radial coordinate: 0 ≤ r ≤ R, m |

S | shrinkage, %, or aria of ¼ of log’s longitudinal section, m2 |

T | temperature, K |

t | temperature, °C |

u | moisture content, kg·kg−1 = %/100 |

z | 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 |

ρ | density, kg·m−3 |

σ | root square mean error (RSME), oC |

τ | time, s |

Δr | 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, — |

@ | at |

avg | average (for relative icing degree or for root square mean error) |

b | basic (for wood density, based on dry mass divided to green volume) |

bw | bound water |

cr | crystallization |

fr | freezing |

fre | end of freezing |

fsp | fiber saturation point |

fw | free water |

gen | generalized (for specific heat capacity) |

ice | ice (for logs’ icing degrees) |

Lat | latent heat |

m | medium (for cooling substance) |

vM | volume of the metal |

nfw | nonfrozen water |

0 | initial |

p | parallel to the wood fibers |

r | radial direction |

total | total (for specific energy needed for freezing of the bound water) |

v | volume |

w | wood |

we | wood effective (for specific heat capacity) |

wL | wood with liquid water in it |

wS | wood with solid state of water (ice) in it |

wUfsp | wood at u = ufsp |

wUnfw | wood at u = unfw |

272.15 | at T = 272.15 K, i.e., at t = −1°C |

293.15 | at T = 293.15 K, i.e., at t = 20°C |