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{split} \rho\_{\text{we-fr}} \cdot \rho\_{\text{w}} \frac{\partial T(r, z, \tau)}{\partial \tau} = \lambda\_{\text{wr}} \left[ \frac{\partial^2 T(r, z, \tau)}{\partial r^2} + \frac{1}{r} \frac{\partial T(r, z, \tau)}{\partial r} \right] + \frac{\partial \lambda\_{\text{wr}}}{\partial T} \left[ \frac{\partial T(r, z, \tau)}{\partial r} \right]^2 \\ + \lambda\_{\text{wp}} \frac{\partial^2 T(r, z, \tau)}{\partial z^2} + \frac{\partial \lambda\_{\text{wp}}}{\partial T} \left[ \frac{\partial T(r, z, \tau)}{\partial z} \right]^2 + q\_{\text{v}} \end{split} \tag{1}$$

Numerical Solution to Two-Dimensional Freezing and Subsequent Defrosting of Logs DOI: http://dx.doi.org/10.5772/intechopen.84706

with an initial condition

$$T(r, z, \mathbf{0}) = T\_{\mathbf{w}\mathbf{0}} \tag{2}$$

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(r,\mathbf{0},\boldsymbol{\tau})}{\partial \boldsymbol{r}} = -\frac{a\_{\mathrm{wp-fr}}(r,\mathbf{0},\boldsymbol{\tau})}{\lambda\_{\mathrm{wp}}(r,\mathbf{0},\boldsymbol{\tau})} [T(r,\mathbf{0},\boldsymbol{\tau}) - T\_{\mathrm{m-fr}}(\boldsymbol{\tau})],\tag{3}$$

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

$$\frac{\partial T(\mathbf{0}, z, \tau)}{\partial \mathbf{z}} = -\frac{a\_{\text{wr-fr}}(\mathbf{0}, z, \tau)}{\lambda\_{\text{wr}}(\mathbf{0}, z, \tau)} [T(\mathbf{0}, z, \tau) - T\_{\text{m-fr}}(\tau)].\tag{4}$$

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{split} c\_{\text{we-dfr}} \cdot \rho\_{\text{w}} \frac{\partial T(r, \mathbf{z}, \tau)}{\partial \tau} = \lambda\_{\text{wr}} \left[ \frac{\partial^2 T(r, \mathbf{z}, \tau)}{\partial r^2} + \frac{1}{r} \frac{\partial T(r, \mathbf{z}, \tau)}{\partial r} \right] + \frac{\partial \lambda\_{\text{wr}}}{\partial T} \left[ \frac{\partial T(r, \mathbf{z}, \tau)}{\partial r} \right]^2 \\ + \lambda\_{\text{wp}} \frac{\partial^2 T(r, \mathbf{z}, \tau)}{\partial \mathbf{z}^2} + \frac{\partial \lambda\_{\text{wp}}}{\partial T} \left[ \frac{\partial T(r, \mathbf{z}, \tau)}{\partial \mathbf{z}} \right]^2 \end{split} \tag{5}$$

with an initial condition

$$T(r, z, \mathbf{0}) = T(r, z, \tau\_{\text{fr}}) \tag{6}$$

and boundary conditions for convective heat transfer:

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

$$\frac{\partial T(r,\mathbf{0},\boldsymbol{\tau})}{\partial r} = -\frac{a\_{\mathrm{wp-dfr}}(r,\mathbf{0},\boldsymbol{\tau})}{\lambda\_{\mathrm{wp}}(r,\mathbf{0},\boldsymbol{\tau})}[T(r,\mathbf{0},\boldsymbol{\tau}) - T\_{\mathrm{m-dfr}}(\boldsymbol{\tau})],\tag{7}$$

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

$$\frac{\partial T(\mathbf{0}, z, \tau)}{\partial \mathbf{z}} = -\frac{a\_{\text{wr-dfr}}(\mathbf{0}, z, \tau)}{\lambda\_{\text{wr}}(\mathbf{0}, z, \tau)} [T(\mathbf{0}, z, \tau) - T\_{\text{m-dfr}}(\tau)].\tag{8}$$

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

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 Tw0 to Tfr˜fw only a cooling of the logs with fully liquid water in them occurs (Figure 1). During the second range from Tfr˜fw to Tfr˜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 ufsp. During the third range from Tfr˜bwm to T<sup>w</sup>˜fre˜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<sup>w</sup>˜fre˜avg to Tdfr˜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 1.

Temperature ranges of the logs' freezing process at u > ufsp 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 > ufsp and thermo-physical characteristics of the wood and of the frozen free and bound water in it.

Numerical Solution to Two-Dimensional Freezing and Subsequent Defrosting of Logs DOI: http://dx.doi.org/10.5772/intechopen.84706

(Figure 2). During the second range from Tdfr�bwm to Tdfr�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 Tdfr�fw to Tw�dfre�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 <sup>1</sup>) above the hygroscopic range, cwe‐fr, are equal to the below:

$$\text{I. range}: \mathcal{c}\_{\text{we-fr1}} = \mathcal{c}\_{\text{w-nfrv}} \tag{9}$$

$$\text{II. range}: c\_{\text{we-fr2}} = c\_{\text{w-nfr}} + c\_{\text{fw}} - c\_{\text{Lat-fw}} \tag{10}$$

$$\text{III. range}: \mathcal{c}\_{\text{we-fr}3} = \mathcal{c}\_{\text{w-fr}} + \mathcal{c}\_{\text{bwm}} - \mathcal{c}\_{\text{Lat-bwm}}.\tag{11}$$

The effective specific heat capacities of the logs during the pointed 3 ranges of the defrosting process (see Figure <sup>2</sup>) above the hygroscopic range, cwe‐dfr, are equal to the below:

$$\text{I. range}: c\_{wc\text{-}df\text{r}1} = c\_{\text{w-fr}} + c\_{\text{bwm}}.\tag{12}$$

$$\text{III. range}: \mathcal{c}\_{\text{we-dfr2}} = \mathcal{c}\_{\text{w-nfr}} + \mathcal{c}\_{\text{fw}} \tag{13}$$

$$\text{III. Range } \colon \mathcal{c}\_{\text{we-dfr3}} = \mathcal{c}\_{\text{w-nfr}}.\tag{14}$$

Mathematical descriptions of the specific heat capacities c<sup>w</sup>�nfr, c<sup>w</sup>�fr, cfw, cbwm, and also of the thermal conductivities of non-frozen, λ<sup>w</sup>‐nfr, and frozen wood, λ<sup>w</sup>‐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–10] and the American specialized literature [11–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, cLat�fw and cLat�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<sup>v</sup>‐fw and <sup>q</sup><sup>v</sup>‐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�<sup>2</sup> �K�<sup>1</sup> ) 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:

$$a\_{\rm wr-fr} = 2.56[T(0, z, \tau) - T\_{\rm m-fr}(\tau)]^{E\_{\rm fr}},\tag{15}$$

$$a\_{\rm wr-dfr} = 2.56[T(\mathbf{0}, z, \tau) - T\_{\rm m-fr}(\tau)]^{E\_{\rm dfr}} \tag{16}$$

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

$$a\_{\rm wp-fr} = \mathbf{1.123} [T(r, \mathbf{0}, \tau) - T\_{\rm m-fr}(\tau)]^{E\_{\rm fr}},\tag{17}$$

$$a\_{\rm wp-dfr} = \mathbf{1.123} [T(r, \mathbf{0}, \tau) - T\_{\rm m-fr}(\tau)]^{E\_{\rm dfr}},\tag{18}$$

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 Δτ 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 <sup>∂</sup><sup>T</sup> derivatives <sup>∂</sup><sup>T</sup>, and <sup>∂</sup><sup>T</sup> are approximated using the built computational mesh along <sup>∂</sup><sup>r</sup> <sup>∂</sup><sup>z</sup> <sup>∂</sup><sup>τ</sup> 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).

Numerical Solution to Two-Dimensional Freezing and Subsequent Defrosting of Logs DOI: http://dx.doi.org/10.5772/intechopen.84706

Taking into consideration the relationships [9, 19]

$$
\lambda\_{\rm w} = \lambda\_{\rm w0r,p} \cdot \boldsymbol{\gamma} \cdot [\mathbf{1} + \boldsymbol{\beta} \cdot (T - 27 \mathbf{3.15})],\tag{19}
$$

$$\lambda\_{\rm w0r} = K\_{\rm wr} \cdot v \cdot \left[ 0.165 + (1.39 + 3.8u) \cdot \left( 3.3 \cdot 10^{-7} \rho\_{\rm b}^2 + 1.015 \cdot 10^{-3} \rho\_{\rm b} \right) \right], \tag{20}$$

$$\lambda\_{\rm w0p} = K\_{\rm wp} \cdot v \cdot \left[ \mathbf{0.165} + (\mathbf{1.39} + \mathbf{3.8u}) \cdot \left( \mathbf{3.3 \cdot 10^{-7}} \rho\_{\rm b}^2 + \mathbf{1.015 \cdot 10^{-3}} \rho\_{\rm b} \right) \right], \tag{21}$$

and using the coefficient

$$K\_{\rm wp/wrr} = \frac{K\_{\rm wp}}{K\_{\rm wr}},\tag{22}$$

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

$$\begin{aligned} T\_{i,k}^{u+1} &= T\_{i,k}^{u} + \frac{\lambda\_{\text{wbr}} \cdot \chi \cdot \Delta \tau}{c\_{\text{we-fr1},2,3}^{u} \cdot \rho\_{\text{w}} \cdot \Delta r^{2}}, \\\\ \left\{ \left[ 1 + \beta \cdot \left( T\_{i,k}^{u} - 273.15 \right) \right] \cdot \begin{bmatrix} T\_{i-1,k}^{u} + T\_{i+1,k}^{u} + K\_{\text{p}/r} \left( T\_{i,k-1}^{u} + T\_{i,k+1}^{u} \right) - \\\\ - \left( 2 + 2K\_{\text{p}/r} \right) T\_{i,k}^{u} + \frac{1}{i-1} \left( T\_{i-1,k}^{u} - T\_{i,k}^{u} \right) \end{bmatrix} \right\}, \\\\ &+ \beta \cdot \left[ \left( T\_{i-1,k}^{u} - T\_{i,k}^{u} \right)^{2} + K\_{\text{wp}/\text{wr}} \left( T\_{i,k-1}^{u} - T\_{i,k}^{u} \right)^{2} \right] \end{aligned} \tag{23}$$

where λwor and λwop are the wood thermal conductivities at t = 0°C in radial and longitudinal direction respectively, W�m�<sup>1</sup> �K�<sup>1</sup> ; γ, β, and ν—coefficients for the determination of λw, λwor and λwop. Equation for calculation of γ, β, and ν for nonfrozen and frozen wood are given in [9, 19]; Kwr and Kwp —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,… τfr/Δτ; τ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 <sup>1</sup>), cwe‐fr1, cwe‐fr2, and cwe‐fr3, which are unitedly represented as cwe‐fr1,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�<sup>3</sup> ) according to the following widely accepted in the literature equation [1–17, 28]:

$$
\rho\_{\mathbf{w}} = \rho\_{\mathbf{b}} \cdot (\mathbf{1} + \boldsymbol{\mu}),
\tag{24}
$$

where ρb is the basic density of the wood, based on dry mass divided to green volume, kg�m�<sup>3</sup> ; u—wood moisture content, kg�kg�<sup>1</sup> .

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{aligned} T\_{i,k}^{u+1} &= T\_{i,k}^{u} + \frac{\lambda\_{\text{w0r}} \cdot \boldsymbol{\chi} \cdot \Delta \boldsymbol{\pi}}{c\_{\text{we-diffr},2,3}^{u} \cdot \rho\_{\text{w}} \cdot \Delta \boldsymbol{\pi}^{2}}. \\\\ &\cdot \left\{ \left[ \mathbbm{1} + \boldsymbol{\beta} \cdot \left( T\_{i,k}^{u} - 273.15 \right) \right] \cdot \begin{bmatrix} T\_{i-1,k}^{u} + T\_{i+1,k}^{u} + K\_{\text{p}/r} \left( T\_{i,k-1}^{u} + T\_{i,k+1}^{u} \right) - \\ - \left( 2 + 2K\_{\text{p}/r} \right) T\_{i,k}^{u} + \frac{1}{i-1} \left( T\_{i-1,k}^{u} - T\_{i,k}^{u} \right) \end{bmatrix} \right\}, \\\ &\cdot \boldsymbol{\beta} \cdot \left[ \left( T\_{i-1,k}^{u} - T\_{i,k}^{u} \right)^{2} + K\_{\text{wp}/\text{wr}} \left( T\_{i,k-1}^{u} - T\_{i,k}^{u} \right)^{2} \right] \end{aligned} \tag{25}$$

where the effective specific heat capacities of the log during the pointed above three ranges of its defrosting process (see Figure <sup>2</sup>), cwe‐dfr1, cwe‐dfr2, and cwe‐dfr3, which are unitedly represented as cwe‐dfr1, <sup>2</sup>, 3in 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 = τfr/Δτ, (τfr/Δτ) + 1, (τfr/Δτ) + 2, (τfr/Δτ) + 3, …., (τfr+τdfr)/Δτ, where τ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:

$$T\_{i,k}^{0} = T\_{\mathbf{w}0}.\tag{26}$$

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

$$T\_{i,k}^{0} = T\_{i,k}^{n=\mathbf{r}\_{\text{fr}}/\Delta \mathbf{r}}.\tag{27}$$

where τ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(r, \mathbf{0}, \tau)}{\partial r} \approx \frac{T\_{i,1}^{n+1} - T\_{i,2}^{n}}{\Delta r}. \tag{28}$$

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

$$-\frac{a\_{\rm wp\,fr}(r,\mathbf{0},\tau)}{\lambda\_{\rm wp}(r,\mathbf{0},\tau)}[T(r,\mathbf{0},\tau)-T\_{\rm m\,fr}(\tau)] \approx -\frac{a\_{\rm wp\,-fr}^{\rm u}}{\lambda\_{\rm w0p}\cdot\chi\cdot\left[\mathbf{1}+\beta\cdot\left(T\_{i,1}^{\rm u}-273.15\right)\right]}\left(T\_{i,1}^{\rm u+1}-T\_{\rm m-fr}^{\rm u+1}\right),\tag{29}$$

Numerical Solution to Two-Dimensional Freezing and Subsequent Defrosting of Logs DOI: http://dx.doi.org/10.5772/intechopen.84706

where according to Eq. (17)

$$a\_{\rm wp-fr}^{\rm n} = \mathbf{1}.\mathbf{123} \left[ T\_{i,1}^{\rm n} - T\_{\rm m-fr}^{\rm n} \right]^{E\_{\rm fr}}.\tag{30}$$

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

$$\frac{T\_{i,1}^{n+1} - T\_{i,2}^{n}}{\Delta r} \approx -\frac{a\_{\text{wp}-\text{fr}}^{n}}{\lambda\_{\text{w0p}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{i,1}^{n} - 273.15\right)\right]} \left(T\_{i,1}^{n+1} - T\_{\text{m-fr}}^{n+1}\right). \tag{31}$$

During the solving of the model it is needed to determine the temperature on the log's frontal surface, T<sup>n</sup> i, þ 1 1 , for each next step n + 1 along the time coordinate. Using Eq. (31), T<sup>n</sup> i, þ 1 1 is equal to

$$T\_{i,1}^{n+1} \approx T\_{i,2}^{n} - \frac{\Delta r \cdot d\_{\rm wp-fr}^{n} \cdot T\_{i,1}^{n+1}}{\lambda\_{\rm 0p} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{i,1}^{n} - 273.15\right)\right]} + \frac{\Delta r \cdot d\_{\rm wp-fr}^{n} \cdot T\_{\rm m-fr}^{n+1}}{\lambda\_{\rm w0p} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{i,1}^{n} - 273.15\right)\right]}.\tag{32}$$

After designation of

$$G\_{i, \text{1-fr}}^{\text{n}} = \frac{\Delta r \cdot a\_{\text{wp}-\text{fr}}^{\text{n}}}{\lambda\_{\text{w0p}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{i, \text{1}}^{\text{n}} - 273.15\right)\right]},\tag{33}$$

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

$$T\_{i,1}^{n+1} = \frac{T\_{i,2}^{n} + G\_{i,1-\text{fr}}^{n} \cdot T\_{\text{m-fr}}^{n+1}}{1 + G\_{i,1-\text{fr}}^{n}}.\tag{34}$$

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

$$T\_{i,1}^{n+1} = \frac{T\_{i,2}^n + G\_{i,1-\text{dfr}}^n \cdot T\_{m-\text{dfr}}^{n+1}}{1 + G\_{i,1-\text{dfr}}^n} . \tag{35}$$

The variable G<sup>n</sup> i, <sup>1</sup>�dfr in Eq. (35) is equal to

$$G\_{i,1-\text{dfr}}^{\text{n}} = \frac{\Delta r \cdot a\_{\text{wp}-\text{dfr}}^{\text{n}}}{\lambda\_{\text{w0p}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{i,1}^{\text{n}} - 273.15\right)\right]},\tag{36}$$

where according to Eq. (18)

$$a\_{\rm wp-dfr}^{\rm n} = \mathbf{1.123} \left[ T\_{i,1}^{\rm n} - T\_{\rm m-fr}^{\rm n} \right]^{E\_{\rm dfr}}.\tag{37}$$

#### 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(0, z, \tau)}{\partial z} \approx \frac{T\_{1,k}^{n+1} - T\_{1,k}^{n}}{\Delta r}. \tag{38}$$

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

$$-\frac{a\_{\text{wr-fr}}(0,z,\tau)}{\lambda\_{\text{wr}}(0,z,\tau)}[T(0,z,\tau)-T\_{\text{m-fr}}(\tau)] \approx -\frac{a\_{\text{wr-fr}}^{\text{n}}}{\lambda\_{\text{wr0r}}\cdot\chi\cdot\left[1+\beta\cdot\left(T\_{1,k}^{\text{n}}-273.15\right)\right]}\left(T\_{1,k}^{\text{n+1}}-T\_{\text{m-fr}}^{\text{n+1}}\right). \tag{39}$$

where according to Eq. (15)

$$a\_{\rm wr-fr}^{\rm n} = 2.56 \left[ T\_{1,k}^{\rm u} - T\_{\rm m-fr}^{\rm u} \right]^{E\_{\rm fr}}.\tag{40}$$

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

$$\frac{T\_{1,k}^{n+1} - T\_{2,k}^{n}}{\Delta r} \approx -\frac{a\_{\text{wr}-\text{fr}}^{n}}{\lambda\_{0\text{r}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{1,k}^{n} - 273.15\right)\right]} \left(T\_{1,k}^{n+1} - T\_{\text{m-fr}}^{n+1}\right). \tag{41}$$

Using Eq. (41), the temperature T<sup>n</sup> 1, þ k 1 is equal to

$$T\_{1,k}^{n+1} \approx T\_{2,k}^{n} - \frac{\Delta r \cdot \alpha\_{\text{wr}-\text{fr}}^{n} \cdot T\_{1,k}^{n+1}}{\lambda\_{\text{w0r}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{1,k}^{n} - 273.15\right)\right]} + \frac{\Delta r \cdot \alpha\_{\text{wr}-\text{fr}}^{n} \cdot T\_{\text{m}-\text{fr}}^{n+1}}{\lambda\_{\text{w0r}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{1,k}^{n} - 273.15\right)\right]}.\tag{42}$$

After designation of

$$G\_{1,k-\text{fr}}^{\text{n}} = \frac{\Delta r \cdot a\_{\text{wr}-\text{fr}}^{\text{n}}}{\lambda\_{\text{w0r}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{1,k}^{\text{n}} - 273.15\right)\right]},\tag{43}$$

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

$$T\_{1,k}^{n+1} = \frac{T\_{2,k}^n + G\_{1,k-\text{fr}}^n \cdot T\_{m-\text{fr}}^{n+1}}{1 + G\_{1,k-\text{fr}}^n}.\tag{44}$$

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

$$T\_{1,k}^{n+1} = \frac{T\_{2,k}^{n} + G\_{1,k-\text{dfr}}^{n} \cdot T\_{m-\text{fr}}^{n+1}}{1 + G\_{1,k-\text{dfr}}^{n}}.\tag{45}$$

The variable G<sup>n</sup> <sup>1</sup>,k�dfr in Eq. (45) is equal to

$$G\_{1,k-\text{dfr}}^{\text{n}} = \frac{\Delta r \cdot a\_{\text{wr}-\text{dfr}}^{\text{n}}}{\lambda\_{\text{w0r}} \cdot \chi \cdot \left[1 + \beta \cdot \left(T\_{1,k}^{\text{n}} - 273.15\right)\right]},\tag{46}$$

where according to Eq. (16)

$$a\_{\rm wr-dfr}^{\rm n} = 2.56 \left[ T\_{1,k}^{\rm n} - T\_{\rm m-dfr}^{\rm n} \right]^{E\_{\rm dfr}}.\tag{47}$$
