Open access peer-reviewed chapter

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

Written By

Nencho Deliiski and Natalia Tumbarkova

Submitted: 19 September 2018 Reviewed: 21 December 2018 Published: 25 April 2019

DOI: 10.5772/intechopen.83772

From the Edited Volume

## Low-temperature Technologies

Edited by Tatiana Morosuk and Muhammad Sultan

Chapter metrics overview

View Full Metrics

## 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]:

c we ρ w T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 +   λ wp 2 T r z τ z 2 + λ wp T T r z τ z 2 + q v E1

with an initial condition

T r z 0 = T w 0 E2

and boundary conditions for convective heat transfer:

• along the radial coordinate r on the logs’ frontal surface during the cooling process

T r 0 τ r = α wp fr r 0 τ λ wp r 0 τ T r 0 τ T m fr τ , E3

• along the longitudinal coordinate z on the logs’ cylindrical surface during the cooling

T 0 z τ z = α wr fr 0 z τ λ wr 0 z τ T 0 z τ T m fr τ . E4

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, 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]:

q v fw = K ψ fw ρ w L cr ice ψ ice fw τ , E5
q v bw = K ψ bw ρ w L cr ice ψ ice bw τ , E6

where L cr ice = 3.34·105 J·kg−1 [2, 5, 22, 23, 24] and

K ψ fw = ρ w 1 Ψ ice fw ρ wUfsp Ψ ice fw ρ w , E7
K ψ bw = ρ wUfsp 1 Ψ ice bw ρ wUnfw Ψ ice bw ρ w . E8

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 m3 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]:

ρ w = ρ b 1 + u , E9
ρ wUfsp = ρ b 1 + u fsp . E10

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]:

ρ wUnfw = ρ b 1 + u nfw 1 S v 100 u fsp 272.15 u nfw , E11

where

u fsp 272.15 = u fsp 293.15 + 0.021 , E12
u nfw = 0.12 + u fsp 272.15 0.12 exp 0.0567 T 272.15 @ 213.15 K T 272.15 K . E13

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

After substituting q v in Eq. (1) by the expression for q v fw from Eq. (5) and of ψ ice fw τ in this expression by ψ ice fw τ Т Т , it is obtained that

ρ w c we fr fw T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ wz 2 T r z τ z 2 + λ w z T T r z τ z 2 + K ψ fw ρ w L cr ice ψ ice fw τ T r z τ T . E14

For the freezing of the free water in the wood, Eq. (14) is transformed in the following form:

ρ w c we fr fw K ψ fw L cr ice ψ ice fw T T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ wz 2 T r z τ z 2 + λ w z T T r z τ z 2 . E15

Equation (15) can be represented, as follows:

ρ w c we fr fw T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ wz 2 T r z τ z 2 + λ w z T T r z τ z 2 . E16

The effective specific heat capacity c we fr fw in Eq. (16) is equal to

c we fr fw = c we fr fw K ψ fw L cr ice ψ ice fw T = c we fr fw c Lat fw , E17

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

ρ w c we fr bw T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ wp 2 T r z τ z 2 + λ wp T T r z τ z 2 + K ψ bw ρ w L cr ice ψ ice bw τ T r z τ T . E18

For the case of freezing the bound water in the wood, Eq. (18) is transformed in the following form:

ρ w c we fr bw K ψ bw L cr ice ψ ice bw T T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ wp 2 T r z τ z 2 + λ wp T T r z τ z 2 . E19

Equation (19) can be represented, as follows:

ρ w c we fr bw T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ wp 2 T r z τ z 2 + λ wp T T r z τ z 2 . E20

The effective specific heat capacity c we fr bw in Eq. (20) is equal to

c we fr bw = c we fr bw K ψ bw L cr ice ψ ice bw T = c we fr bw c Lat bwm , E21

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.

### 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, Ψ 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:

Ψ ice bw = m ice bw m ice bw + m nfw = u fsp u nfw u fsp u nfw + u nfw = 1 u nfw u fsp . E22

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:

u fsp = u fsp 293.15 0.001 T 293.15 , E23

where u fsp 293.15 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

Ψ ice bw = 1 0.12 + u fsp 272.15 0.12 exp 0.0567 T 272.15 u fsp 293.15 0.001 T 293.15 . E24

and

Ψ ice bw T = 0.0567 u fsp 272.15 0.12 exp 0.0567 T 272.15 u fsp 293.15 0.001 T 293.15 u fsp 293.15 0.001 T 293.15 2 + 0.001 0.12 + u fsp 272.15 0.12 exp 0.0567 T 272.15 u fsp 293.15 0.001 T 293.15 2 . E25

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 Т.

### 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 > 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 w 0 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:

I . Range : c we fr gen 1 = c w nfr , E26
II . Range : c we fr gen 2 = c w nfr + c fw c Lat fw , E27
III . Range : c fr gen 3 = c w fr + c bwm c Lat bwm . E28

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.

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:

ρ w c we fr gen 1 , 2 , 3 T r z τ τ = λ wr 2 T r z τ r 2 + 1 r . T r z τ r + λ wr T T r z τ r 2 + λ 2 T r z τ z 2 + λ T T r z τ z 2 . E29

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

α wr fr = 2.56 T 0 z τ T m fr τ E fr , E30

• at frontal surface of the logs

α wp fr = 1.123 T r 0 τ T m fr τ E fr , E31

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.

### 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, Q fr bw total , can be calculated according to the following equation:

Q fr bw total = Q fr bw Q Lat bw , E32

where Q fr bw 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 Q Lat bw 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, Q w , with an initial mass temperature T w 0 to a given average mass temperature T w avg is determined using the following equation [2, 5, 10]:

Q w = c w ρ w 3.6 10 6 Т w avg T w 0 . E33

The multiplier 3.6·106 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

Q fr bw = c bw avg ρ w 3.6 10 6 272.15 T w fr avg , E34

where according to [8, 9]

c bw avg = 1.8938 10 4 u fsp 272.15 0.12 exp 0.0567 T w fr avg 272.15 1 + u @ T w fre avg T 272.15 K , E35
T w fr avg = 1 S w S w T r z τ d S w @ T w fre avg T r z τ 272.15 K , E36

and the area of ¼ of the longitudinal section of the log subjected to freezing, S w , is equal to

S w = D L 4 . E37

Based on Eq. (33), it is possible to calculate also the energy Q Lat bw according to the equation

Q Lat bw = c Lat bw ρ w 3.6 10 6 272.15 T w fr avg , E38

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

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

## 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 u fsp 293.15 = 0.30 kg·kg−1 and volume shrinkage S v = 11.8% and also of spruce sapwood with u fsp 293.15 = 0.32 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, 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:

Т m fr = a fr + c fr τ 0.5 1 + b fr τ 0.5 + d fr τ , E39

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.

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

### 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, Ψ ice bw avg , is carried out according to the following equation, which was obtained in [23] using Eq. (24):

Ψ ice bw avg = 1 S w . S w 1 0.12 + u fsp 272.15 0.12 exp 0.0567 [ T r z τ 272.15 ] u fsp 293.15 0.001 T r z τ 293.15 ] d S w @ T w fre avg T r z τ 272.15 K . E40

Figure 5 presents the calculated change of the logs’ icing degree Ψ ice bw avg and the derivative Ψ ice bw / d T 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 / d T = 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 / d T 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 / d T 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 Qfr-bw and QLat-bw

Figure 6 presents according to Eqs. (34) and (38) the calculated change of the thermal energies Q fr bw and Q Lat bw during the 30 h refrigeration of the studied logs.

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 / d T .

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.

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

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

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

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

### Subscripts

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

total

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 = u fsp

wUnfw

wood at u = u nfw

### Superscripts

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

## References

1. 1. Vorreiter L. Holztechnologisches Handbuch. Vienna: Fromm; 1949. 2080 p
2. 2. Chudinov BS. Theoretical research of thermo physical properties and thermal treatment of wood [thesis for DSc.]. Krasnojarsk, USSR: SibLTI; 1966 (in Russian)
3. 3. Kollmann FF. Technologie des Holzes and Holzwerkstoffe. I. Bd., 2 Auflage. Berlin/Göttingen, Heidelberg: Springer; 1951. 2233 p
4. 4. Kollmann FF, Côté WA Jr. Principles of Wood Science and Technology. I. Solid Wood. New York: Springer-Verlag, Berlin, Heidelberg; 1984. 592 p
5. 5. Shubin GS. Drying and Thermal Treatment of Wood. Moscow, USSR: Lesnaya Promyshlennost; 1990. 337 p. (in Russian)
6. 6. Požgaj A, Chovanec D, Kurjatko S, Babiak M. Structure and Properties of Wood. 2nd ed. Bratislava: Priroda a.s; 1997. 485p. (in Slovak)
7. 7. Trebula P, Klement I. Drying and Hydro-Thermal Treatment of Wood. Slovakia: Technical University in Zvolen; 2002. 449 p. (in Slovak)
8. 8. Deliiski N, Dzurenda L. Modelling of the Thermal Processes in the Technologies for Wood Thermal Treatment. Slovakia: Technical University in Zvolen; 2010. 224 p. (in Russian)
9. 9. Deliiski N. Transient heat conduction in capillary porous bodies. In: Ahsan A, editor. Convection and Conduction Heat Transfer. Rijeka, Croatia: InTech Publishing House; 2011. pp. 149-176. DOI: 10.5772/21424
10. 10. Deliiski N. Modelling of the Energy Needed for Heating of Capillary Porous Bodies in Frozen and Non-frozen States. Saarbrücken, Germany: Lambert Academic Publishing, Scholars’ Press; 2013. 116 p. Available from: http://www.scholars-press.com//system/covergenerator/build/1060
11. 11. Steinhagen HP. Computerized finite-difference method to calculate transient heat conduction with thawing. Wood and Fiber Science. 1986;18(3):460-467
12. 12. Steinhagen HP. Heat transfer computation for a long, frozen log heated in agitated water or steam—A practical recipe. Holz als Roh- und Werkstoff. 1991;49(7–8):287-290. DOI: 10.1007/BF02663790
13. 13. Steinhagen HP, Lee HW. Enthalpy method to compute radial heating and thawing of logs. Wood and Fiber Science. 1988;20(4):415-421
14. 14. Khattabi A, Steinhagen HP. Numerical solution to two-dimensional heating of logs. Holz als Roh- und Werkstoff. 1992;50(7-8):308-312. DOI: 10.1007/BF02615359
15. 15. Khattabi A, Steinhagen HP. Analysis of transient non-linear heat conduction in wood using finite-difference solutions. Holz als Roh- und Werkstoff. 1993;51(4):272-278. DOI: 10.1007/BF02629373
16. 16. Khattabi A, Steinhagen HP. Update of numerical solution to two-dimensional heating of logs. Holz als Roh- und Werkstoff. 1995;53(1):93-94. DOI: 10.1007/BF02716399
17. 17. Deliiski N. Modelling and automatic control of heat energy consumption required for thermal treatment of logs. Drvna Industrija. 2004;55(4):181-199
18. 18. Deliiski N. Computation of the 2-dimensional transient temperature distribution and heat energy consumption of frozen and non-frozen logs. Wood Research. 2009;54(3):67-78
19. 19. Deliiski N, Dzurenda L, Tumbarkova N, Angelski D. Computation of the temperature conductivity of frozen wood during its defrosting. Drvna Industrija. 2015;66(2):87-96. DOI: 10.5552/drind.2015.1351
20. 20. Hadjiski M, Deliiski N. Cost oriented suboptimal control of the thermal treatment of wood materials. IFAC-PapersOnLine. 2015;48(24):54-59. DOI: 10.1016/j.ifacol.2015.12.056
21. 21. Hadjiski M, Deliiski N. Advanced control of the wood thermal treatment processing. Cybernetics and Information Technologies, Bulgarian Academy of Sciences. 2016;16(2):179-197. DOI: 10.1515/cait-2016-0029
22. 22. Deliiski N, Tumbarkova N. A methodology for experimental research of the freezing process of logs. Acta Silvatica et Lignaria Hungarica. 2016;12(2):145-156. DOI: 10.1515/aslh-2016-0013
23. 23. Deliiski N, Tumbarkova N. An approach and an algorithm for computation of the unsteady icing degrees of logs subjected to freezing. Acta Facultatis Xylologiae Zvolen. 2017;59(2):91-104
24. 24. Deliiski N, Tumbarkova N. An approach for computing the heat sources in logs subjected to freezing. Acta Silvatica et Lignaria Hungarica. 2018;14(1):35-49. http://aslh.nyme.hu/index.php?id=29435&L=0
25. 25. Salcudean M, Abdullah Z. On the numerical modelling of heat transfer during solidification processes. International Journal for Numerical Methods in Engineering. 1988;25:445-473. DOI: 10.1002/nme.1620250212
26. 26. Dantzig J. Modelling liquid-solid phase changes with melt convection. International Journal for Numerical Methods in Engineering. 1989;28(8):1769-1785. DOI: 10.1002/nme.1620280805
27. 27. Hu H, Argyropoulos S. Mathematical modelling of solidification and melting: A review. Modelling and Simulation in Materials Science and Engineering. 1996;4:371-396. DOI: 10.1088/0965-0393/4/4/004
28. 28. Mihailov E, Petkov V. Cooling parameters and heat quantity of the metal during continuous casting of blooms. International Review of Mechanical Engineering. 2010;4(2):176-184
29. 29. Hrčka R. Model in free water in wood. Wood Research. 2017;62(6):831-837
30. 30. Stamm AJ. Wood and Cellulose Science. New York: Ronald Press; 1964
31. 31. Siau JF. Transport Processes in Wood. NewYork: Springer; 1984
32. 32. Kanter KR. Investigation of the thermal properties of wood [thesis]. Moscow, USSR: MLTI; 1955 (in Russian)
33. 33. Telegin AS, Shvidkiy BS, Yaroshenko UG. Heat- and Mass Transfer. Moscow: Akademkniga; 2002. 456 p. (in Russian)
34. 34. Table Curve 2D. Available from: http://www.sigmaplot.co.uk/products/tablecurve2d/tablecurve2d.php
35. 35. Hadjiski M, Deliiski N, Grancharova A. Spatiotemporal parameter estimation of thermal treatment process via initial condition reconstruction using neural networks. In: Hadjiski M, Atanasov KT, editors. Intuitionistic Fuzziness and Other Intelligent Theories and Their Applications. Cham, Switzerland: Springer International Publishing AG; 2019. pp. 51-80. DOI: 10.1007/978-3-319-78931-6

Written By

Nencho Deliiski and Natalia Tumbarkova

Submitted: 19 September 2018 Reviewed: 21 December 2018 Published: 25 April 2019