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

Abstract

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

1. Introduction

It is known that the duration of the thermal treatment of the frozen logs in the winter aimed at their plasticizing for the production of veneer and the very high energy consumption needed for this treatment depend on the degree of the logs’ icing [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. For example, for the defrosting and plasticizing of beech and oak logs with an initial temperature of −10°C and moisture content of 0.6 kg·kg⁻¹ approximately 68 and 81 kWh·m⁻³ thermal energy respectively [10] is needed.

In the specialized literature there are limited reports about the temperature distribution in subjected to defrosting frozen logs [8, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22] and there is very little information about research of the temperature distribution in logs during their freezing given by the authors only [23, 24, 25]. That is why the modeling and the multi-parameter study of the mutually connected freezing and defrosting processes of logs are of considerable scientific and practical interest.

Considerable contribution to the calculation of the non-stationary distribution of the temperature in frozen and non-frozen logs and to the duration of their heating at conductive boundary conditions has been made by Steinhagen [11, 12] and later one-dimensional and two-dimensional models have been created and solved [13, 14, 15, 16, 17]. The thermal energy, which is needed for the melting of the ice, which has been formed from the freezing of the bound water in the wood, has not been taken into account in the models of cited references.

The models assume that the fiber saturation point is identical for all wood species (i.e., u_fsp = 0.3 kg·kg⁻¹ = const). However, it is known that there are significant differences between the fiber saturation points of the different wood species [1, 2, 3, 4, 5, 6, 7, 8]. The indicated complications and incompleteness in these models have been overcome in the suggested by Deliiski [9, 18, 19] 2-dimensional mathematical model of the transient non-linear heat conduction in frozen logs during their heating and defrosting.

For different engineering calculations it is needed to be able to determine the icing degree of the wood materials depending on the temperature of the influencing on them gas or liquid medium and on the duration of their staying in this medium. Such calculations are carried out using mathematical models, which describe adequately the complex processes of the freezing of both the free and bound water in the wood.

The computer solutions of these models give the non-stationary distribution of the temperature in the materials during their cooling below temperatures, at which a freezing of the whole amount of the free water and a freezing of respective, depending on the temperature, part of the bound water in the wood occurs [24, 25].

This chapter presents the creation, numerical solving and validation of two mutually connected two-dimensional mathematical models of the transient non-linear heat conduction in logs during their freezing and subsequent defrosting at convective boundary conditions.

The model of the freezing process takes into account for the first time the impact of the internal sources of latent heat of both the free and bound water on the temperature distribution. The both models reflect the impact of the temperature on the fiber saturation point of each wood species, with whose participation the current values of the thermo-physical characteristics in each separate volume point of the subjected to freezing and subsequent defrosting logs are computed.

The chapter also presents and visualizes the results from experimental and simulative investigation of the 2D non-stationary temperature distribution in the longitudinal section of beech and poplar logs with a diameter of 0.24 m, length of 0.48 m, and different moisture content during many hours their freezing in a freezer and subsequent defrosting at curvilinear changing temperature of the processing air medium.

2. Mechanism of the temperature distribution in logs during their freezing and subsequent defrosting

2.1 Mathematical model of the 2D temperature distribution in logs subjected to freezing

When the length of the logs, L, is larger than their diameter, D, by not more than 3–4 times, for the calculation of the change in the temperature in the logs’ longitudinal sections (i.e., along the coordinates r and z of these sections) during their freezing in air medium the following 2D mathematical model can be used [24]:

$$\begin{array}{ll}{c}_{\text{we}‐\mathrm{fr}}\cdot {\rho }_{\mathrm{w}}\frac{\partial T\left(r,z,\tau \right)}{\partial \tau }=& {\lambda }_{\mathrm{wr}}\left[\frac{{\partial }^{2}T\left(r,z,\tau \right)}{\partial r^2}+\frac{1}{r}.\frac{\partial T\left(r,z,\tau \right)}{\partial r}\right]+\frac{\partial {\lambda }_{\mathrm{wr}}}{\partial T}{\left[\frac{\partial T\left(r,z,\tau \right)}{\partial r}\right]}^2\\ & + {\lambda }_{\mathrm{wp}}\frac{{\partial }^{2}T\left(r,z,\tau \right)}{\partial z^2}+\frac{\partial {\lambda }_{\mathrm{wp}}}{\partial T}{\left[\frac{\partial T\left(r,z,\tau \right)}{\partial z}\right]}^2+q_{\mathrm{v}}\end{array}$$ E1

with an initial condition

$$T\left(r,z,0\right)=T_{\mathrm{w}0}$$ E2

and boundary conditions for convective heat transfer:

• along the radial coordinate r on the logs’ frontal surface during the freezing process (refer to Figure 3):

  $$\frac{\partial T\left(r,0,\tau \right)}{\partial r}=-\frac{{\alpha }_{\mathrm{wp}‐\mathrm{fr}}\left(r,0,\tau \right)}{{\lambda }_{\mathrm{wp}}\left(r,0,\tau \right)}\left[T\left(r,0,\tau \right)-T_{\mathrm{m}‐\mathrm{fr}}\left(\tau \right)\right],$$ E3

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

$$\frac{\partial T\left(0,z,\tau \right)}{\partial z}=-\frac{{\alpha }_{\mathrm{wr}‐\mathrm{fr}}\left(0,z,\tau \right)}{{\lambda }_{\mathrm{wr}}\left(0,z,\tau \right)}\left[T\left(0,z,\tau \right)-T_{\mathrm{m}‐\mathrm{fr}}\left(\tau \right)\right].$$ E4

Eqs. (1)–(4) represent a common form of a mathematical model of the logs’ freezing process, i.e., of the 2D temperature distribution in logs subjected to freezing.

2.2 Mathematical model of the 2D temperature distribution in logs subjected to defrosting

In cases when the length of the logs, L, is larger than their diameter, D, by not more than 3–4 times, for the calculation of the change in T in the longitudinal sections of the logs during their defrosting in air processing medium the following 2D mathematical model can be used [8, 9]:

$$\begin{array}{ll}{c}_{\text{we}‐\mathrm{dfr}}\cdot {\rho }_{\mathrm{w}}\frac{\partial T\left(r,z,\tau \right)}{\partial \tau }=& {\lambda }_{\mathrm{wr}}\left[\frac{{\partial }^{2}T\left(r,z,\tau \right)}{\partial r^2}+\frac{1}{r}.\frac{\partial T\left(r,z,\tau \right)}{\partial r}\right]+\frac{\partial {\lambda }_{\mathrm{wr}}}{\partial T}{\left[\frac{\partial T\left(r,z,\tau \right)}{\partial r}\right]}^2\\ & + {\lambda }_{\mathrm{wp}}\frac{{\partial }^{2}T\left(r,z,\tau \right)}{\partial z^2}+\frac{\partial {\lambda }_{\mathrm{wp}}}{\partial T}{\left[\frac{\partial T\left(r,z,\tau \right)}{\partial z}\right]}^2\end{array}$$ E5

with an initial condition

$$T\left(r,z,0\right)=T\left(r,z,{\tau }_{\mathrm{fr}}\right)$$ E6

and boundary conditions for convective heat transfer:

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

  $$\frac{\partial T\left(r,0,\tau \right)}{\partial r}=-\frac{{\alpha }_{\mathrm{wp}‐\mathrm{dfr}}\left(r,0,\tau \right)}{{\lambda }_{\mathrm{wp}}\left(r,0,\tau \right)}\left[T\left(r,0,\tau \right)-T_{\mathrm{m}‐\mathrm{dfr}}\left(\tau \right)\right],$$ E7

• along the longitudinal coordinate z on the logs’ cylindrical surface during the defrosting process:

$$\frac{\partial T\left(0,z,\tau \right)}{\partial z}=-\frac{{\alpha }_{\mathrm{wr}‐\mathrm{dfr}}\left(0,z,\tau \right)}{{\lambda }_{\mathrm{wr}}\left(0,z,\tau \right)}\left[T\left(0,z,\tau \right)-T_{\mathrm{m}‐\mathrm{dfr}}\left(\tau \right)\right].$$ E8

Eqs. (5)–(8) represent a mathematical model of the logs’ defrosting process, i.e., of the 2D temperature distribution in logs subjected to defrosting immediately after their freezing.

2.3 Mathematical description of the thermo-physical characteristics of logs

On Figures 1 and 2 the temperature ranges are presented, at which the processes of the logs’ freezing and subsequent defrosting respectively are carried out when u >  $u_{\mathrm{fsp}}$ .

Figure 1.

Temperature ranges of the logs’ freezing process at u > $u_{\mathrm{fsp}}$ and thermo-physical characteristics of the wood and of the frozen free and bound water in it.

Figure 2.

Temperature ranges of the logs’ defrosting process at u > $u_{\mathrm{fsp}}$ and thermo-physical characteristics of the wood and of the frozen free and bound water in it.

There thermo-physical characteristics of the logs and of both the frozen free and bound water in them have also been shown. The information on these characteristics is very important for the solving of the mutually connected mathematical models given above.

During the first range of the logs’ freezing process from $T_{\mathrm{w}0}$ to $T_{\mathrm{fr}-\mathrm{fw}}$ only a cooling of the logs with fully liquid water in them occurs (Figure 1). During the second range from $T_{\mathrm{fr}-\mathrm{fw}}$ to $T_{\mathrm{fr}-\mathrm{bwm}}$ a further cooling of the logs occurs until reaching of the state needed for starting of the crystallization of the free water. During that range also the phase transition of this water into ice is carried out. The second range is absent when u is less than $u_{\mathrm{fsp}}$ . During the third range from $T_{\mathrm{fr}-\mathrm{bwm}}$ to $T_{\mathrm{w}-\mathrm{fre}-\mathrm{avg}}$ a further cooling of the logs is carried out until reaching of the state needed for starting of the crystallization of the bound water. During this range also the phase transition of the bound water into ice is performed.

During the first range of the logs’ defrosting process from $T_{\mathrm{w}-\mathrm{fre}-\mathrm{avg}}$ to $T_{\mathrm{dfr}-\mathrm{bwm}}$ a heating of the frozen logs is carried out until reaching of the state needed for starting and realization of the gradually melting of the frozen bound water in them (Figure 2). During the second range from $T_{\mathrm{dfr}-\mathrm{bwm}}$ to $T_{\mathrm{dfr}-\mathrm{fw}}$ a further heating of the logs occurs until reaching of the state needed for starting and realization of the melting of the frozen free water in them. During the third range from $T_{\mathrm{dfr}-\mathrm{fw}}$ to $T_{\mathrm{w}-\text{dfre}-\mathrm{avg}}$ a further heating of the logs with fully liquid water in them occurs.

The effective specific heat capacities of the logs during the pointed 3 ranges of the freezing process (see Figure 1) above the hygroscopic range, $c_{\text{we}‐\mathrm{fr}}$ , are equal to the below:

$$\mathrm{I}. \text{range}:c_{\text{we}‐\mathrm{fr}1}=c_{\mathrm{w}-\mathrm{nfr}},$$ E9

$$\mathrm{II}. \text{range}:c_{\text{we}‐\mathrm{fr}2}=c_{\mathrm{w}-\mathrm{nfr}}+c_{\mathrm{fw}}-c_{\mathrm{Lat}-\mathrm{fw}},$$ E10

$$\mathrm{III}. \text{range}:c_{\mathrm{wе}‐\mathrm{fr}3}=c_{\mathrm{w}-\mathrm{fr}}+c_{\mathrm{bwm}}-c_{\mathrm{Lat}-\mathrm{bwm}}.$$ E11

The effective specific heat capacities of the logs during the pointed 3 ranges of the defrosting process (see Figure 2) above the hygroscopic range, $c_{\text{we}‐\mathrm{dfr}}$ , are equal to the below:

$$\mathrm{I}. \text{range}:c_{\mathrm{wе}‐\mathrm{dfr}1}=c_{\mathrm{w}-\mathrm{fr}}+c_{\mathrm{bwm}},$$ E12

$$\mathrm{II}. \text{range}:c_{\text{we}‐\mathrm{dfr}2}=c_{\mathrm{w}-\mathrm{nfr}}+c_{\mathrm{fw}},$$ E13

$$\mathrm{III}. \text{Range}:c_{\text{we}‐\mathrm{dfr}3}=c_{\mathrm{w}-\mathrm{nfr}}.$$ E14

Mathematical descriptions of the specific heat capacities $c_{\mathrm{w}-\mathrm{nfr}}$ , $c_{\mathrm{w}-\mathrm{fr}}$ , $c_{\mathrm{fw}}$ , $c_{\mathrm{bwm}}$ , and also of the thermal conductivities of non-frozen, ${\lambda }_{\mathrm{w}‐\mathrm{nfr}}$ , and frozen wood, ${\lambda }_{\mathrm{w}‐\mathrm{fr}}$ , have been suggested by the first co-author earlier [9, 10, 19] using the experimentally determined in the dissertations by Kanter [26] and Chudinov [2] data for their change as a function of t and u. These relations are used in both the European [5, 6, 7, 8, 9, 10] and the American specialized literature [11, 12, 13, 14, 15, 16, 17] when calculating various processes of the wood thermal treatment.

Mathematical descriptions of the specific heat capacities, which are formed by the release of the latent heat of both the free and bound water during their crystallization in the wood, $c_{\mathrm{Lat}-\mathrm{fw}}$ and $c_{\mathrm{Lat}-\mathrm{bwm}}$ respectively (refer to Eqs. (10) and (11)) have been given by the authors in [23]. Mathematical descriptions of the internal heat sources separately for the free and bound water, $q_{\mathrm{v}‐\mathrm{fw}}$ and $q_{\mathrm{v}‐\mathrm{bw}}$ respectively, which participate in the right-hand part of Eq. (1), have been also given.

Our study has shown that for the calculation (in W·m⁻²·K⁻¹) of the radial and longitudinal transfer coefficients of the logs, which participate in the boundary conditions of the model the following equations are most suitable [27]:

• in the radial direction on the cylindrical surface of the horizontally situated logs:

$${\alpha }_{\mathrm{wr}-\mathrm{fr}}=2.56{\left[T\left(0,z,\tau \right)-T_{\mathrm{m}-\mathrm{fr}}\left(\tau \right)\right]}^{E_{\mathrm{fr}}},$$ E15

$${\alpha }_{\mathrm{wr}-\mathrm{dfr}}=2.56{\left[T\left(0,z,\tau \right)-T_{\mathrm{m}-\mathrm{fr}}\left(\tau \right)\right]}^{E_{\mathrm{dfr}}}$$ E16

• in the longitudinal direction on the frontal surface of the logs:

$${\alpha }_{\mathrm{wp}-\mathrm{fr}}=1.123{\left[T\left(r,0,\tau \right)-T_{\mathrm{m}-\mathrm{fr}}\left(\tau \right)\right]}^{E_{\mathrm{fr}}},$$ E17

$${\alpha }_{\mathrm{wp}-\mathrm{dfr}}=1.123{\left[T\left(r,0,\tau \right)-T_{\mathrm{m}-\mathrm{fr}}\left(\tau \right)\right]}^{E_{\mathrm{dfr}}},$$ E18

where Е_fr and Е_dfr are exponents, whose values are determined during the solving and verification of the models through minimization of the root square mean error (RSME) between the calculated by the models and experimentally obtained results about the change of the temperature fields in subjected to freezing and subsequent defrosting logs.

2.4 Transformation of the models to a form suitable for programming

Due to the correct cylindrical shape of the logs, for the solution of the above presented two mutually connected mathematical models an explicit form of the finite-difference method has been used. That form allows for the exclusion of any simplifications of the models and also of the necessity for any linearization of the mathematically described variables in them [8, 9, 18, 19].

The large calculation resources of the contemporary computers eliminate the inconvenience, which creates the limitation for the value of the step along the time coordinate $\mathit{\Delta \tau }$ by using the explicit form [9].

According to the idea of the finite-difference method, the temperature, which is a uninterrupted function of space and time, is presented using a grid vector, and the derivatives $\frac{\partial T}{\partial r}$ , $\frac{\partial T}{\partial z}$ and $\frac{\partial T}{\partial \tau }$ are approximated using the built computational mesh along the spatial and time coordinates through their finite-difference (discrete) analogues.

2.4.1 Transformation of the equations of thermo-conductivity in the models

The transformation of Eqs. (1) and (5) of the models in their suitable for programming discrete analogues has been realized using the presented on Figure 3 coordinate system. That system shows the positioning of the knots of the calculation mesh, in which the non-stationary distribution of the temperature in the longitudinal section of subjected to freezing and subsequent defrosting log has been calculated. The calculation mesh has been built on ¼ of the longitudinal section of the log due to the circumstance that this ¼ is mirror symmetrical towards the remaining 3/4 of the same section.

Figure 3.

Calculation mesh for the solution of the models (left) and positioning of the knots of 2D calculation mesh on 1/4 of longitudinal section of a log subjected to freezing and subsequent defrosting (right).

Taking into consideration the relationships [9, 19]

$${\lambda }_{\mathrm{w}}={\lambda }_{\mathrm{w}0\mathrm{r},\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T-273.15\right)\right],$$ E19

$${\mathrm{\lambda }}_{\mathrm{w}0\mathrm{r}}=K_{\mathrm{wr}}\cdot v\cdot \left[0.165+\left(1.39+3.8u\right)\cdot \left(3.3\cdot {10}^{-7}{\rho }_{\mathrm{b}}^2+1.015\cdot {10}^{-3}{\rho }_{\mathrm{b}}\right)\right],$$ E20

$${\mathrm{\lambda }}_{\mathrm{w}0\mathrm{p}}=K_{\mathrm{wp}}\cdot v\cdot \left[0.165+\left(1.39+3.8u\right)\cdot \left(3.3\cdot {10}^{-7}{\rho }_{\mathrm{b}}^2+1.015\cdot {10}^{-3}{\rho }_{\mathrm{b}}\right)\right],$$ E21

and using the coefficient

$$K_{\mathrm{wp}/\mathrm{wr}}=\frac{K_{\mathrm{wp}}}{K_{\mathrm{wr}}},$$ E22

after applying the explicit form of the finite-differences method to Eq. (1), it is transformed into the following system of equations:

$$\begin{array}{l}{T}_{i,k}^{n+1}=T_{i,k}^n+\frac{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \mathrm{\Delta }\tau }{c_{\text{we}‐\mathrm{fr}1,2,3}^n\cdot {\rho }_{\mathrm{w}}\cdot \mathrm{\Delta }{r}^2}.\\ .\left\{\begin{array}{l}\left[1+\beta \cdot \left(T_{i,k}^n-273.15\right)\right].\left[\begin{array}{l}{T}_{i-1,k}^n+T_{i+1,k}^n+K_{\mathrm{p}/\mathrm{r}}\left(T_{i,k-1}^n+T_{i,k+1}^n\right)-\\ -\left(2+2K_{\mathrm{p}/\mathrm{r}}\right)T_{i,k}^n+\frac{1}{i-1}\left(T_{i-1,k}^n-T_{i,k}^n\right)\end{array}\right]\\ + \beta \cdot \left[{\left(T_{i-1,k}^n-T_{i,k}^n\right)}^2+K_{\mathrm{wp}/\mathrm{wr}}{\left(T_{i,k-1}^n-T_{i,k}^n\right)}^2\right]\end{array}\right\}\end{array},$$ E23

where λ_wor and λ_wop are the wood thermal conductivities at t = 0°C in radial and longitudinal direction respectively, W·m⁻¹·K⁻¹; γ, β, and ν—coefficients for the determination of λ_w, λ_wor and λ_wop. Equation for calculation of γ, β, and ν for non-frozen and frozen wood are given in [9, 19]; K_wr and K_wp—empirically determined coefficients, which take into account the influence of the radial and longitudinal anatomical directions on the wood thermal conductivity, -; Δτ—interval between time levels, i.e., step along the time coordinate, with which the model is solving, s; n—current number of the step Δτ along the time coordinate: n = 0, 1, 2, 3,… ${\tau }_{\mathrm{fr}}$ /Δτ; ${\tau }_{\mathrm{fr}}$ —duration of the freezing process of logs, s; Δr = Δz—step along the coordinates r and z, by which the model is being solved, m; i—current number of the knots of the calculation mesh along the coordinate r: 1 ≤ i ≤ 1 + (R/Δr); k—current number of the knots along the coordinate z: 1 ≤ k ≤ 1 + [(L/2)/Δr]; R—log’s radius, m; L—log’s length, m.

The effective specific heat capacities of the log during the pointed above three ranges of its freezing process (see Figure 1), $c_{\text{we}‐\mathrm{fr}1}$ , $c_{\text{we}‐\mathrm{fr}2}$ , and $c_{\text{we}‐\mathrm{fr}3}$ , which are unitedly represented as $c_{\text{we}‐\mathrm{fr}1,2,3}$ in Eq. (23), are computed according to Eqs. (9)–(11) separately for each knot of the calculation mesh.

The density of the log’s wood above the hygroscopic range, ρ_w, which participates in Eq. (23), is calculated (in kg·m⁻³) according to the following widely accepted in the literature equation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 28]:

$${\rho }_{\mathrm{w}}={\rho }_{\mathrm{b}}\cdot \left(1+u\right),$$ E24

where ${\rho }_{\mathrm{b}}$ is the basic density of the wood, based on dry mass divided to green volume, kg·m⁻³; u—wood moisture content, kg·kg⁻¹.

After applying the explicit form of the finite-differences method to Eq. (5), it is transformed into the following system of equations, which is similar to Eq. (23):

$$\begin{array}{l}{T}_{i,k}^{n+1}=T_{i,k}^n+\frac{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \mathrm{\Delta }\tau }{c_{\text{we}‐\mathrm{dfr}1,2,3}^n\cdot {\rho }_{\mathrm{w}}\cdot \mathrm{\Delta }{r}^2}.\\ .\left\{\begin{array}{l}\left[1+\beta \cdot \left(T_{i,k}^n-273.15\right)\right].\left[\begin{array}{l}{T}_{i-1,k}^n+T_{i+1,k}^n+K_{\mathrm{p}/\mathrm{r}}\left(T_{i,k-1}^n+T_{i,k+1}^n\right)-\\ -\left(2+2K_{\mathrm{p}/\mathrm{r}}\right)T_{i,k}^n+\frac{1}{i-1}\left(T_{i-1,k}^n-T_{i,k}^n\right)\end{array}\right]\\ + \beta \cdot \left[{\left(T_{i-1,k}^n-T_{i,k}^n\right)}^2+K_{\mathrm{wp}/\mathrm{wr}}{\left(T_{i,k-1}^n-T_{i,k}^n\right)}^2\right]\end{array}\right\}\end{array},$$ E25

where the effective specific heat capacities of the log during the pointed above three ranges of its defrosting process (see Figure 2), $c_{\text{we}‐\mathrm{dfr}1}$ , $c_{\text{we}‐\mathrm{dfr}2}$ , and $c_{\text{we}‐\mathrm{dfr}3}$ , which are unitedly represented as $c_{\text{we}‐\mathrm{dfr}1,2,3}$ in Eq. (25), are computed according to Eqs. (12)–(14) separately for each knot of the calculation mesh.

The current number of the step Δτ along the time coordinate, n, in Eq. (25) is equal to: n =  ${\tau }_{\mathrm{fr}}$ /Δτ, ( ${\tau }_{\mathrm{fr}}$ /Δτ) + 1, ( ${\tau }_{\mathrm{fr}}$ /Δτ) + 2, ( ${\tau }_{\mathrm{fr}}$ /Δτ) + 3, …., ( ${\tau }_{\mathrm{fr}}$ + ${\tau }_{\mathrm{dfr}}$ )/Δτ, where ${\tau }_{\mathrm{dfr}}$ is the duration of the log’s defrosting process, s.

2.4.2 Transformation of the equations of initial conditions in the models

The initial condition (2) of the model of logs’ freezing process obtains the following discrete finite-difference form:

$${Т}_{i,k}^0=T_{\mathrm{w}0}.$$ E26

The initial condition (6) of the model of logs’ defrosting process obtains the following discrete finite-difference form:

$${Т}_{i,k}^0=T_{i,k}^{n={\tau }_{\mathrm{fr}}/\mathrm{\Delta }\tau }.$$ E27

where ${\tau }_{\mathrm{fr}}$ is the duration of the freezing process of logs, which precedes the beginning of their subsequent defrosting, s.

2.4.3 Transformation of the equations of boundary conditions in the models along the radial coordinate

The discrete finite difference analogue of the left-hand part of Eq. (3), which is suitable for programming in FORTRAN, has the following form [8, 9]:

$$\frac{\partial T\left(r,0,\tau \right)}{\partial r}\approx \frac{T_{i,1}^{n+1}-T_{i,2}^n}{\mathrm{\Delta }r}.$$ E28

The discrete analogue of the right-hand part of Eq. (3) has the following form:

$$-\frac{{\alpha }_{\mathrm{wp}‐\mathrm{fr}}\left(r,0,\tau \right)}{{\lambda }_{\mathrm{wp}}\left(r,0,\tau \right)}\left[T\left(r,0,\tau \right)-T_{\mathrm{m}‐\mathrm{fr}}\left(\tau \right)\right]\approx -\frac{{\alpha }_{\mathrm{wp}-\mathrm{fr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{i,1}^n-273.15\right)\right]}\left(T_{i,1}^{n+1}-T_{\mathrm{m}-\mathrm{fr}}^{n+1}\right),$$ E29

where according to Eq. (17)

$${\alpha }_{\mathrm{wp}‐\mathrm{fr}}^n=1.123{\left[T_{i,1}^n-T_{\mathrm{m}-\mathrm{fr}}^n\right]}^{E_{\mathrm{fr}}}.$$ E30

After alignment of Eq. (28) with Eq. (29), it is obtained that

$$\frac{{Т}_{i,1}^{n+1}-T_{i,2}^n}{\mathrm{\Delta }r}\approx -\frac{{\alpha }_{\mathrm{wp}-\mathrm{fr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{i,1}^n-273.15\right)\right]}\left(T_{i,1}^{n+1}-T_{\mathrm{m}-\mathrm{fr}}^{n+1}\right).$$ E31

During the solving of the model it is needed to determine the temperature on the log’s frontal surface, $T_{i,1}^{n+1}$ , for each next step n + 1 along the time coordinate. Using Eq. (31), $T_{i,1}^{n+1}$ is equal to

$${Т}_{i,1}^{n+1}\approx T_{i,2}^n-\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wp}-\mathrm{fr}}^n\cdot T_{i,1}^{n+1}}{{\lambda }_{0\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{i,1}^n-273.15\right)\right]}+\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wp}-\mathrm{fr}}^n\cdot T_{\mathrm{m}-\mathrm{fr}}^{n+1}}{{\lambda }_{\mathrm{w}0\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{i,1}^n-273.15\right)\right]}.$$ E32

After designation of

$$G_{i,1-\mathrm{fr}}^n=\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wp}-\mathrm{fr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{i,1}^n-273.15\right)\right]},$$ E33

the boundary condition (3) of the logs’ freezing process obtains the following final form, suitable for programming:

$${Т}_{i,1}^{n+1}=\frac{T_{i,2}^n+G_{i,1-\mathrm{fr}}^n\cdot T_{\mathrm{m}-\mathrm{fr}}^{n+1}}{1+G_{i,1-\mathrm{fr}}^n}.$$ E34

Analogously, the boundary condition (7) of the logs’ defrosting process obtains the following final form, suitable for programming:

$${Т}_{i,1}^{n+1}=\frac{T_{i,2}^n+G_{i,1-\mathrm{dfr}}^n\cdot T_{\mathrm{m}-\mathrm{dfr}}^{n+1}}{1+G_{i,1-\mathrm{dfr}}^n}.$$ E35

The variable $G_{i,1-\mathrm{dfr}}^n$ in Eq. (35) is equal to

$$G_{i,1-\mathrm{dfr}}^n=\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wp}-\mathrm{dfr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{p}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{i,1}^n-273.15\right)\right]},$$ E36

where according to Eq. (18)

$${\alpha }_{\mathrm{wp}‐\mathrm{dfr}}^n=1.123{\left[T_{i,1}^n-T_{\mathrm{m}-\mathrm{fr}}^n\right]}^{E_{\mathrm{dfr}}}.$$ E37

2.4.4 Transformation of the equations of boundary conditions in the models along the longitudinal coordinate

The discrete finite difference analogue of the left-hand part of Eq. (4) at Δz = Δr has the following suitable for programming form [8, 9]:

$$\frac{\partial T\left(0,z,\tau \right)}{\partial z}\approx \frac{T_{1,k}^{n+1}-T_{1,k}^n}{\mathrm{\Delta }r}.$$ E38

The discrete analogue of the right-hand part of Eq. (4) has the following form:

$$-\frac{{\alpha }_{\mathrm{wr}‐\mathrm{fr}}\left(0,z,\tau \right)}{{\lambda }_{\mathrm{wr}}\left(0,z,\tau \right)}\left[T\left(0,z,\tau \right)-T_{\mathrm{m}‐\mathrm{fr}}\left(\tau \right)\right]\approx -\frac{{\alpha }_{\mathrm{wr}-\mathrm{fr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{1,k}^n-273.15\right)\right]}\left(T_{1,k}^{n+1}-T_{\mathrm{m}-\mathrm{fr}}^{n+1}\right).$$ E39

where according to Eq. (15)

$${\alpha }_{\mathrm{wr}‐\mathrm{fr}}^n=2.56{\left[T_{1,k}^n-T_{\mathrm{m}-\mathrm{fr}}^n\right]}^{E_{\mathrm{fr}}}.$$ E40

After alignment of Eq. (38) with Eq. (35), it is obtained that

$$\frac{{Т}_{1,k}^{n+1}-T_{2,k}^n}{\mathrm{\Delta }r}\approx -\frac{{\alpha }_{\mathrm{wr}-\mathrm{fr}}^n}{{\lambda }_{0\mathrm{r}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{1,k}^n-273.15\right)\right]}\left(T_{1,k}^{n+1}-T_{\mathrm{m}-\mathrm{fr}}^{n+1}\right).$$ E41

Using Eq. (41), the temperature $T_{1,k}^{n+1}$ is equal to

$${Т}_{1,k}^{n+1}\approx T_{2,k}^n-\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wr}-\mathrm{fr}}^n\cdot T_{1,k}^{n+1}}{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{1,k}^n-273.15\right)\right]}+\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wr}-\mathrm{fr}}^n\cdot T_{\mathrm{m}-\mathrm{fr}}^{n+1}}{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{1,k}^n-273.15\right)\right]}.$$ E42

After designation of

$$G_{1,k-\mathrm{fr}}^n=\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wr}-\mathrm{fr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{1,k}^n-273.15\right)\right]},$$ E43

the boundary condition (4) of the logs’ freezing process obtains the following final form, suitable for programming:

$${Т}_{1,k}^{n+1}=\frac{T_{2,k}^n+G_{1,k-\mathrm{fr}}^n\cdot T_{\mathrm{m}-\mathrm{fr}}^{n+1}}{1+G_{1,k-\mathrm{fr}}^n}.$$ E44

Analogously, the boundary condition (8) of the logs’ defrosting process obtains the following final form, suitable for programming in FORTRAN:

$${Т}_{1,k}^{n+1}=\frac{T_{2,k}^n+G_{1,k-\mathrm{dfr}}^n\cdot T_{\mathrm{m}-\mathrm{fr}}^{n+1}}{1+G_{1,k-\mathrm{dfr}}^n}.$$ E45

The variable $G_{1,k-\mathrm{dfr}}^n$ in Eq. (45) is equal to

$$G_{1,k-\mathrm{dfr}}^n=\frac{\mathrm{\Delta }r\cdot {\alpha }_{\mathrm{wr}-\mathrm{dfr}}^n}{{\lambda }_{\mathrm{w}0\mathrm{r}}\cdot \gamma \cdot \left[1+\beta \cdot \left(T_{1,k}^n-273.15\right)\right]},$$ E46

where according to Eq. (16)

$${\alpha }_{\mathrm{wr}‐\mathrm{dfr}}^n=2.56{\left[T_{1,k}^n-T_{\mathrm{m}-\mathrm{dfr}}^n\right]}^{E_{\mathrm{dfr}}}.$$ E47

3. Experimental research of the logs’ freezing and defrosting

3.1 Experimental research of the 2D temperature distribution in beech and poplar logs during their freezing and subsequent defrosting

For the validation of the suggested above mathematical models, it is necessary to have experimentally obtained data about the 2D temperature distribution in logs during their freezing and subsequent defrosting. That is why we carried out such experiments.

The logs subjected to freezing and subsequent defrosting in our experimental research were with a diameter of 240 mm, length of 480 mm, and u >  $u_{\mathrm{fsp}}$ . This means, that the logs contained the maximum possible amount of bound water. They were produced from the sap-wood of a freshly felled beech (Fagus sylvatica L.) and poplar (Populus nigra L.) trunks.

Before the experiments, four holes with diameters of 6 mm and different lengths were drilled in each log parallel to its axis until reaching of the characteristic points of the log [23].

The coordinates of the characteristic points of the logs are given on Figure 4. These coordinates of the points allow for the determination of the 2D temperature distribution in logs during their freezing and subsequent defrosting.

Figure 4.

Radial (left) and longitudinal (right) coordinates of four characteristic points for measurement of the temperature in logs subjected to freezing and subsequent defrosting.

For the freezing of the logs according to the suggested methodology by the authors [23], a horizontal freezer was used with adjustable temperature range from −1 to −30°C.

Sensors Pt100 with long metal casings were positioned in the drilled four holes of the logs. The automatic measurement and record of t_m, φ_m, and t in the characteristic points of the logs during the experiments was realized by Data Logger type HygroLog NT3 produced by the Swiss firm ROTRONIC AG (http:/www.rotronic.com). The Data Logger has software HW4 for graphical presentation of the experimentally obtained data. After reaching of about −28°C in the log’s center during the freezing, the freezer was switch off. Then its lid was opened and a defrosting of the log at room temperature was carried out.

In Figures 5 and 6 the change in the temperature of the processing air medium, $t_m$ , and in its humidity, φ_m, and also in the temperature in 4 characteristic points of a beech and of a poplar log respectively during their separately 50 and 70 h is presented. The record of all data was made automatically by Data Logger with intervals of 15 min. The left coordinate axis on the figures is graduated at % of φ_m, and the left one is graduated at °C of t.

Figure 5.

Experimentally determined change in t_m, φ_m, and t in four points of the studied beech log during its 50 h freezing and 50 h subsequent defrosting.

Figure 6.

Experimentally determined change in t_m, φ_m, and t in four points of the studied poplar log during its 50 h freezing and 70 h subsequent defrosting.

The initial temperature, t_w0, basic density, ρ_b, moisture content, u, duration of the freezing, τ_fr, and duration of the subsequent defrosting, τ_dfr, of the logs during the experiments were as follow:

• for the beech log: t_w0 = 22.4°C, ρ_b = 684 kg·kg⁻¹, u = 0.63 kg·kg⁻¹, τ_fr = 50 h, and τ_dfr = 50 h;

• for the poplar log: t_w0 = 19.8°C, ρ_b = 359 kg·kg⁻¹, u = 1.44 kg·kg⁻¹, τ_fr = 50 h, and τ_dfr = 70 h.

3.2 Mathematical description of the air medium temperature during logs’ freezing and subsequent defrosting

The change in the shown on Figures 5 and 6 freezing and defrosting air medium temperatures, ${Т}_{\mathrm{m}-\mathrm{fr}}$ and ${Т}_{\mathrm{m}-\mathrm{dfr}}$ , with very high accuracy (correlation 0.98 and Root Square Mean Error (RSME), σ < 1.0°C) has been approximated with the help of the software package Table Curve 2D (http://www.sigmaplot.co.uk/products/tablecurve2d/tablecurve2d.php) by the following equations:

• during the freezing of the beech log:

  $${Т}_{\mathrm{m}-\mathrm{fr}}=\frac{a_{\mathrm{fr}}+c_{\mathrm{fr}}{\tau }^{0.5}}{1+b_{\mathrm{fr}}{\tau }^{0.5}},$$ E48

  whose coefficients are equal to:

  а_fr = 294.3352069, b_fr = 0.010648218, c_fr = 2.468350514;

• during the freezing of the poplar log:

  $${Т}_{\mathrm{m}-\mathrm{fr}}=\frac{a_{\mathrm{fr}}+c_{\mathrm{fr}}\tau +{е}_{\mathrm{fr}}{\tau }^2+g_{\mathrm{fr}}{\tau }^3}{1+b_{\mathrm{fr}}\tau +d_{\mathrm{fr}}{\tau }^2+f_{\mathrm{fr}}{\tau }^3+h_{\mathrm{fr}}{\tau }^4},$$ E49

  whose coefficients are equal to: а_fr = 301.8210985, b_fr = 0.0004515197, c_fr = 0.111484207, d_fr = −6.6585073·10⁻⁹, e_fr = −1.6653·10⁻⁶, f_fr = 2.52712·10⁻¹⁴, g_fr = 6.46801·10⁻¹², h_fr = 2.94924·10⁻²¹;

• during the defrosting of both the beech and poplar logs:

$${Т}_{\mathrm{m}-\mathrm{dfr}}=\frac{a_{\mathrm{dfr}}+c_{\mathrm{dfr}}{\tau }^{0.5}}{1+b_{\mathrm{dfr}}{\tau }^{0.5}},$$ E50

whose coefficients are equal to as follows:

а_dfr = 297.1420433, b_dfr = −0.00237763, c_dfr = −0.70526837 for the beech log;

а_dfr = 296.3637194, b_dfr = −0.00236425, c_dfr = −0.69281743 for the poplar log.

Eqs. (48) and (49) were used for the solving of Eqs. (3) and (7) and Eq. (50) was used for the solving of Eqs. (4) and (8) of the model.

4. Numerical solution of the mathematical models of the logs’ freezing and subsequent defrosting processes

For the numerical solution of the mutually connected mathematical models, a software package was prepared in Visual FORTRAN Professional developed by Microsoft. Using the package, computations were carried out for the determination of the 2D non-stationary change of t in the characteristic points of ¼ of the longitudinal sections of the studied logs, whose experimentally determined temperature fields are presented on Figures 3 and 4.

The models have been solved with the help of explicit schemes of the finite difference method in a way, analogous to the one used and described in [8, 9, 10, 18, 19]. For the computation of the temperature distribution in ¼ of the longitudinal section of the logs, which is mirror symmetrical towards the remaining 3/4 of the same section, the models were solved with step Δr = Δz = 0.006 m along the coordinates r and z and with the same initial and boundary conditions, as they were during the experimental research.

During the solving of the models, mathematical descriptions of the thermo-physical characteristics of beech wood with $u_{\mathrm{fsp}}^{293.15}=0.31$ kg·kg⁻¹, K_wr = 1.35, and K_wp = 2.40, and also of poplar sapwood with $u_{\mathrm{fsp}}^{293.15}=0.35$ kg·kg⁻¹, K_wr = 1.48, and K_wp = 2.88 [8, 10].

The mathematical models of the logs’ freezing and subsequent defrosting processes have been solved with different values of the exponents ${Е}_{\mathrm{fr}}$ and ${Е}_{\mathrm{dfr}}$ in Eqs. (15)–(18). The calculated by the models change of the temperature in the four characteristic points of the longitudinal logs’ sections with each of the used values of ${Е}_{\mathrm{fr}}$ and ${Е}_{\mathrm{dfr}}$ during the freezing and defrosting has been compared mathematically with the corresponding one experimentally determined change of t in the same points with an interval of 15 min. The aim of this comparison was to find that the values of ${Е}_{\mathrm{fr}}$ and ${Е}_{\mathrm{dfr}}$ , which ensures the best qualitative and quantitative compliance between the calculated and experimentally determined temperature fields in the logs’ longitudinal sections.

As a criterion of the best compliance between the compared values of the temperature total for the four characteristic points the minimum value of RSME, ${\sigma }_{\mathrm{avg}}$ , has been used.

The average value of RSME, ${\sigma }_{\mathrm{avg}}$ , has been calculated according to the following equation:

$${\sigma }_{\mathrm{avg}}=\sqrt{\frac{{\sum }_{n=1}^N{\sum }_{p=1}^P{\left(t_{p,n}^{\text{calc}}-t_{p,n}^{\exp }\right)}^2}{P\cdot \left(N-1\right)}},$$ E51

where $t_{p,n}^{\text{calc}}$ and $t_{p,n}^{\exp }$ are the calculated and experimentally established temperatures in the characteristic points; p—number of the characteristic points in the log’s longitudinal section: р = 1, 2, 3, 4, i.e., Р = 4 is inputted into Eq. (51); n—number of the moments of the freezing and subsequent defrosting processes: n = 1, 2, 3,,….., N = ( ${\tau }_{\mathrm{fr}}$ + ${\tau }_{\mathrm{dfr}}$ )/(150Δτ) because of the circumstance that the comparison of the calculated values of t with experimentally determined values in the same points has been made with an interval of 15 min = 900 s = 150Δτ.

For the determination of ${\sigma }_{\mathrm{avg}}$ software program in the calculation environment of MS Excel has been prepared. At ${\tau }_{\mathrm{fr}}$ + ${\tau }_{\mathrm{dfr}}$ = 100 h for the beech log and at ${\tau }_{\mathrm{fr}}$ + ${\tau }_{\mathrm{dfr}}$ = 120 h for the poplar log, RSME has been calculated with the help of the program simultaneously for total (N + 1)·P = 1604 temperature–time points for the beech log and for total 1924 such points for the poplar log during their freezing and subsequent defrosting. It was determined that the minimum values of RSME overall for the studied 4 characteristic points are equal to σ_avg = 1.29°C for the beech log and to σ_avg = 1.50°C for the poplar log.

The minimum values of ${\sigma }_{\mathrm{avg}}$ were obtained with the following values of the exponents in Eqs. (15)–(18):

• for the beech log: ${Е}_{\mathrm{fr}}$ = 0.52 and ${Е}_{\mathrm{dfr}}$ = 0.32;

• for the poplar log: ${Е}_{\mathrm{fr}}$ = 0.43 and ${Е}_{\mathrm{dfr}}$ = 0.22.

Figures 7 and 8 presents the calculated change in $t_{\mathrm{m}-\mathrm{fr}}$ and $t_{\mathrm{m}-\mathrm{dfr}}$ , which are represented unitedly as $t_{\mathrm{m}}$ , and also in logs’ surface temperature $t_s$ and t of 4 characteristic points of the studied beech and poplar logs respectively.

The comparison to each other of the analogical curves on Figures 5 and 7, and also on Figures 6 and 8 shows good conformity between the calculated and experimentally determined changes in the very complicated temperature fields of the studied logs during their freezing and subsequent defrosting.

Figure 7.

Experimentally determined and calculated change in t_m, φ_m, and t in four points of the beech log during its 50 h freezing and 50 h subsequent defrosting.

Figure 8.

Experimentally determined and calculated change in t_m, φ_m, and t in four points of the poplar log during its 50 h freezing and 70 h subsequent defrosting.

During our wide simulations with the mathematical models, we observed good compliance between computed and experimentally established temperature fields of logs various wood species with different moisture content.

The overall RSME for the studied 4 characteristic points in the logs does not exceed 5% of the temperature ranges between the initial and the end temperatures of the logs subjected to freezing or subjected to subsequent defrosting.

5. Conclusions

This chapter describes the creation, solution, and validation of two mutually connected 2D non-linear mathematical models for the transient heat conduction in subjected to freezing and subsequent defrosting logs with any u ≥ u_fsp. The model of the freezing process takes into account the impact of the internal sources of latent heat of both the free and bound water on the temperature distribution. The both models reflect the impact of the temperature on the fiber saturation point of each wood species, with whose participation the current values of the thermo-physical characteristics in each separate volume point of the subjected to freezing and subsequent defrosting logs are computed.

The mechanism of the temperature distribution in the longitudinal section of the logs during their freezing and subsequent defrosting has been mathematically described by 2D equations of heat conduction. Boundary conditions for convective heat transfer have been implemented in the models. For the transformation of both models into discrete analogues, which are suitable for programming, an explicit form of the finite-difference method has been used.

For the numerical solution of the discrete analogues of the models a software package has been prepared using the programming language FORTRAN, which has been input in the calculation environment of Visual Fortran Professional.

A validation of the models towards own experimentally determined 2D temperature distribution in beech and poplar logs with a diameter of 0.24 m, length of 0.48 m during their separate 50 h freezing in a freezer and many hours subsequent defrosting at room temperature has been carried out. The influence of the curvilinear changing temperature of the air medium in the freezer until reaching of approximately −30°C and also of the air processing medium during the defrosting of logs has been investigated. The following minimum values of the average RSME total for the temperature change in four characteristic points in each of the logs have been obtained:

• σ_avg = 1.29°C for the beech log with t_wo = 22.4°C, ρ_b = 684 kg·m⁻³, and u = 0.63 kg·kg⁻¹;

• σ_avg = 1.50°C for the poplar log with t_wo = 19.8°C, ρ_b = 359 kg·m⁻³, and u = 1.44 kg·kg⁻¹.

During our experimental research it has been determined that in situated on the logs’ inner layers characteristic points the specific practically horizontal sections of retention of the temperature for many hours in the range from 0 to −1°С arise, while in these layers a complete freezing of the free water occurs. Analogous retention of the temperature in the range from −1 to 0°C arises during the logs’ defrosting. The further the point is distanced from the logs’ surfaces and the larger the amount of the free water in the wood is, that much more these sections with temperature retention are extended. Our simulations show that this phenomenon of the freezing and defrosting processes has been correctly reflected in the models (see Figures 7 and 8).

Good adequacy and precision of the models towards the results from wide own experimental studies allow for the carrying out of various calculations with the models, which are connected to the non-stationary temperature distribution in logs from different wood species during their freezing and subsequent defrosting.

The validation of the models with curvilinear change in the temperature of both the freezing and defrosting air mediums will allow us in the future to solve the models with curvilinear changing of the climate temperature [29] over many winter days and nights. The solution will allow for the calculation of the temperature distribution, icing degrees from both the free and bound water, and also different energy characteristics of logs for each desired moment of the their freezing and subsequent defrosting.

The mutually connected models of the freezing and defrosting processes can be applied for the development of scientifically based and energy saving optimized regimes for thermal treatment of frozen logs and also in the software for controllers used for model predictive automatic control [21, 22, 30] of this treatment.

Nomenclature

c

specific heat capacity, J·kg⁻¹·K⁻¹

D

diameter, m

E

exponent, −

L

latent heat, J·kg⁻¹, or length, m

q

internal heat source, W·m⁻³

R

radius, m

r

radial coordinate: 0 ≤ r ≤ R, m

T

temperature, K

t

temperature, °C

u

moisture content, kg·kg⁻¹ = %/100

z

longitudinal coordinate: 0 ≤ z ≤ L/2, m

α

heat transfer coefficients between log’s surfaces and the surrounding air medium, W·m⁻² ·K⁻¹

λ

thermal conductivity, W·m⁻¹·K⁻¹

φ

relative humidity, %

ρ

density, kg·m⁻³

σ

root square mean error (RSME), °C

τ

time, s

Δr

step along the coordinates r and z for solving of the models, m

Δτ

step along the time coordinate for solving of the models, s

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

bwm

maximum possible amount of the bound water in the wood

cr

crystallization

dfre

end of defrosting

fr

freezing

fre

end of freezing

fsp

fiber saturation point

fw

free water

i

knot of the calculation mesh in the direction along the logs’ radius: i = 1, 2, 3,…, 21.

k

knot of the calculation mesh in longitudinal direction of the logs: k = 1, 2, 3, …, 41.

Lat

latent heat

m

medium (for the air during the logs’ freezing and defrosting)

p

parallel to the wood fibers

r

radial direction

v

volume

w

wood

we

wood effective (for specific heat capacity)

w-fr

wood with frozen water in it

w-nfr

wood with fully liquid water in it

0

initial or at 0°C

wop

parallel to the wood fibers at °C

wor

radial direction of wood at °C

Superscripts n

current number of the step Δτ along the time coordinate: n = 0, 1, 2,…

272.15

at 272.15 K, i.e. at −1°C

293.15

at 293.15 K, i.e. at 20°C