............... Deliiski: 3D Modeling and Visualization of Non-Stationary Temperature ...
DRVNA INDUSTRIJA 64 (4) 293-303 (2013) 293
Nencho Deliiski
1
3D Modeling and
Visualization of Non-
Stationary Temperature
Distribution during Heating
of Frozen Wood
3D modeliranje i vizualizacija nestacionarne
distribucije temperature tijekom zagrijavanja
smrznutog drva
Received – prispjelo: 27. 1. 2013.
Accepted – prihvaćeno: 6. 11. 2013.
UDK: 630*812.141; 630*847
doi:10.5552/drind.2013.1306
AbSTRAcTA 3-dimensional mathematical model has been developed, solved, and veried for the transient
non-linear heat conduction in frozen and non-frozen wood with prismatic shape at arbitrary initial and boundary
conditions encountered in practice. The model takes into account for the rst time the ber saturation point of
each wood species, u
fsp
, and the impact of the temperature on u
fsp
of frozen and non-frozen wood, which are then
used to compute the current values of the thermal and physical characteristics in each separate volume point of
the material subjected to defrosting.
This paper presents solutions of the model with the explicit form of the nite-difference method. Results of simu-
lation investigation of the impact of frozen bound water, as well as of bound and free water, on 3D temperature
distribution in the volume of beech and oak prisms with dimensions 0.4 x 0.4 x 0.8 m during their defrosting at the
temperature of the processing medium of 80 °C are presented, analyzed and visualized through color contour plots.
Keywords: 3D mathematical model, frozen wood, nite difference method, temperature distribution, contour plots
Kreiran je i riješen 3D matematički model te provjeren za nelinearno provođenje topline u smrz-
nutome i nesmrznutom drvu prizmatičnog oblika pri proizvoljnim početnim i rubnim uvjetima koji se susreću u
praksi. Prvi put model uzima u obzir točku zasićenosti vlakanaca za svaku vrstu drva (u
fsp
) i utjecaj temperature na
u
fsp
smrznutoga i nesmrznutog drva, koji se primjenjuju pri izračunavanju trenutačne vrijednosti termo-zikalnih
svojstava u svakoj posebno deniranoj točki volumena materijala koji se odmrzava.
Rad prikazuje rješenja modela s eksplicitnim oblikom metode konačnih razlika. Rezultati simulacijskih istraživanja
o utjecaju zamrznute vezane vode te vezane i slobodne vode na 3D raspodjelu temperature u volumenu bukovih
i hrastovih prizmi dimenzija 0,4 x 0,4 x 0,8 m tijekom odmrzavanja pri temperaturi procesnog medija od 80
о
C
prezentirani su i analizirani te vizualizirani crtežima u boji.
Ključne riječi : 3D matematički model, smrznuto drvo, metoda konačnih razlika, raspodjela temperature , kon-
turni crteži
1
Author is professor at Faculty of Forest Industry, University of Forestry, Soa, Bulgaria.
1
Autor je profesor Fakulteta šumske industrije Šumarskog sveučilišta, Soja, Bugarska.
Deliiski: 3D Modeling and Visualization of Non-Stationary Temperature ... ...............
294 DRVNA INDUSTRIJA 64 (4) 293-303 (2013)
1 INTRODUCTION
1. UVOD
For the optimization of the control of the heating
process of wood in veneer and plywood mills, it is nec-
essary to know the temperature distribution at every
moment of the process (Shubin, 1990; Trebula and
Klement, 2003; Pervan, 2009). Considerable contribu-
tion was made to the calculation of non-stationary dis-
tribution of temperature in frozen and non-frozen logs,
and to the duration of their heating (Steinhagen, 1986,
1991). Later on, 1-dimensional and 2-dimensional
models were developed and solved (Steinhagen et al.,
1987; Steinhagen and Lee, 1988; Khattabi and Steinha-
gen, 1992, 1993, 1995), whose applications are limited
only to wood with moisture content above ber satura-
tion point.
The heat energy, required for melting the ice,
formed from bound water in the wood, has not been
taken into account in these models. The models assume
that the ber saturation point is identical for all wood
species (i.e. u
fsp
= 0.3 kg·kg
-1
= const) and that the melt-
ing of the ice, formed from free water in the wood, oc-
curs at 0 ºC.
However, it is known that there are signicant
differences between the ber saturation point of differ-
ent wood species (Požgai et. al., 1997; Videlov, 2003)
and that, depending on this point, the quantity of the
ice formed from free water in the wood melts at a tem-
perature in the range between -2 ºC and -1 ºC (Chudi-
nov, 1968, 1984). The complications and deciencies
indicated in these models have been overcome by a
2-dimensional mathematical model of the transient
non-linear heat conduction in frozen and non-frozen
logs suggested by Deliiski (2004, 2011).
This paper presents the development, verication
and solutions of an analog 3-dimensional mathematical
model of the transient non-linear heat conduction in
frozen and non-frozen wood with prismatic shape at
arbitrary initial and boundary conditions encountered
in practice. The model takes into account for the rst
time the ber saturation point of each wood species,
u
fsp
, and the impact of the temperature on u
fsp
of frozen
and non-frozen wood, which are then used to compute
the current values of the thermal and physical charac-
teristics in each separate volume point of the material
subjected to defrosting.
This paper also presents the results of simulation
investigation of the impact of the frozen bound water
and free water on 3D temperature distribution in the
volume of beech and oak prisms with dimensions 0.4 x
0.4 x 0.8 m during their defrosting at the temperature
of the processing medium of 80 °C.
2 MATHERIAL AND METHODS
2. MATERIJAL I METODE
2.1. 3D mathematical model of the defrosting
process of prismatic wood materials
2.1. 3D matematički model procesa odmrzavanja
prizmatičnoga drvnog materijala
The defrosting process of prismatic wood materi-
als during their thermal treatment can be described by
a non-linear differential equation of the thermal-con-
ductivity, using the Cartesian coordinates (Deliiski,
2003a):
( )
( )
( )
( )
( )
( )
( )
( )
∂
t∂
l
∂
∂
+
∂
t∂
l
∂
∂
+
+
∂
t∂
l
∂
∂
=
t∂
t∂
r
z
zyxT
uT
zy
zyxT
uT
y
x
zyxT
uT
x
zyxT
uTuTc
,,,
,
,,,
,
,,,
,
,,,
),(,
pt
re
(1)
After the differentiation of the right side of equa-
tion (1) on the spatial coordinates x, y, and z, excluding
the arguments in the brackets for shortening of the re-
cord, the following mathematical model is obtained of
the non-stationary defrosting of wood materials with
prismatic shape subjected to heating:
2
p
2
2
p
2
t
2
2
t
2
r
2
2
re
∂
∂
∂
l∂
+
∂
∂
l+
∂
∂
∂
l∂
+
∂
∂
l+
+
∂
∂
∂
l∂
+
∂
∂
l=
t∂
∂
r
z
T
T
z
T
y
T
T
y
T
x
T
T
x
TT
c
(2)
with an initial condition
( )
0
0,,, TzyxT =
, (3)
and a boundary condition
( ) ( ) ( ) ( )
t=t=t=t
m
,0,,,,0,,,,0 TyxTzxTzyT
. (4)
For the solution of the system of equations (2) to
(4), it is necessary to make a mathematical description
of thermal and physical characteristics of the wood, c
e
,
l
r
, l
t
, l
p
, and of its density, r. Equations in (Deliiski,
2003a, 2011) and (Deliiski and Dzurenda, 2010) pre-
sent a mathematical description of the effective spe-
cic heat capacity coefcient, c
e
, of the frozen wood as
a sum of the capacities of the wood itself, c, and the ice
produced by freezing of the free water, c
fw
, and of the
hygroscopically bound water, c
bw
. Other equations
quoted by the above authors present mathematical de-
scriptions of wood density,
r
, and of its thermal con-
ductivity, l, in different anatomical directions.
The given mathematical descriptions of
e
c
,
r
l
,
t
l
, and
p
l
(Deliiski, 2011), which are part of the
model (2) to (4), have now been updated by taking into
account, for the rst time, the inuence of the ber
saturation point of wood species on the values of ther-
mal and physical characteristics during wood defrost-
ing, and the inuence of the temperature on ber satu-
ration point of frozen and non-frozen wood. This has
been done using the method presented by Deliiski
(2013) during the update of the mathematical descrip-
tion of l.
2.2. Transformation of 3D model to a form suitable
for programming
2.2. Transformacija 3D modela u odgovarajući oblik
za programiranje
The following system of equations (Equation 5)
has been derived by passing to nal increases in equa-
tion (2) with the usage of the same, as well as by the
explicit form of the nite-difference method described
by Deliiski (2003a, 2011) and taking into account the
............... Deliiski: 3D Modeling and Visualization of Non-Stationary Temperature ...
DRVNA INDUSTRIJA 64 (4) 293-303 (2013) 295
mathematical description of the thermal conductivity,
l, in different anatomical directions.
Since in practice prismatic materials subjected to
thermal treatment usually do not have a clear radial or
clear tangential orientation, and are partially radially or
partially tangentially oriented, then in equation (5) in-
stead of the coefcients
0
l
in the observed two ana-
tomical directions, their average arithmetic value can
be used, as it determines the thermal conductivity at 0
°С perpendicular to the wood bers (Equation 6):
(6)
Also, the thermal conductivity at 0 °С in the di-
rection parallel to the bers λ
0p
can be expressed
through
0cr
l
using the equation
(7)
where the coefcient
depends on the
wood species (Deliiski, 2003a).
(5)
For uniformity of the calculations, it is reasona-
ble to use one step of the calculation mesh along the
spatial coordinates
xD
=
yD
=
zD
(see Fig. 1). Taking
into consideration this condition and equations (6) and
(7), the system of equations (5) becomes equation (8).
The initial condition (3) in the model is presented
using the following nite differences equation:
0
0
,,
TT
kji
=
. (9)
(8)
The boundary conditions (4) acquire the follow-
ing form suitatable for programming:
1
m
1
1,,
1
,1,
1
,,1
++++
===
nn
ji
n
ki
n
kj
TTTT . (10)
The presentation of a non-linear differential
equation (2) from the mathematical model through its
discrete analogue (8) corresponds to the setting of the
coordinate system and positioning of the nodes in the
mesh shown in Fig. 1, in which a non-stationary 3D
temperature distribution in prismatic wood materials
during their defrosting is calculated. The calculation
mesh for the solution of the model through the nite-
difference method is built on a 1/8 part of the prism
volume, because of its mirror symmetry with the re-
maining 7/8 parts of the prism volume.
The setting of the coordinate system, shown in
Fig. 1 allows, with the help of only one system of equa-
tions (8), to calculate the change in the temperature in
any mesh node of the volume of the prism subjected to
defrosting at the moment (n + 1)Dt using the already
calculated values of T at the preceding moment nDt.
Wide experimental studies have been performed
for the determination of a 1-, 2- and 3-dimensional
temperature distribution in the volume of frozen and
non-frozen oak, beech, poplar and pine
»
¼
º
«
¬
ª
w
Ww
O
w
w
»
¼
º
«
¬
ª
w
Ww
O
w
w
»
¼
º
«
¬
ª
w
Ww
O
w
w
Ww
Ww
U
z
zyxT
uT
zy
zyxT
uT
y
x
zyxT
uT
x
zyxT
uTuTc
,,,
,
,,,
,
,,,
,
,,,
),(,
pt
re
2
p
2
2
p
2
t
2
2
t
2
r
2
2
re
¸
¹
·
¨
©
§
w
w
w
Ow
w
w
O
¸
¸
¹
·
¨
¨
©
§
w
w
w
Ow
w
w
O
¸
¹
·
¨
©
§
w
w
w
Ow
w
w
O
Ww
w
U
z
T
T
z
T
y
T
T
y
T
x
T
T
x
TT
c
0
0,,, TzyxT
,
W W W W
m
,0,,,,0,,,,0 TyxTzxTzyT
.
>@
>@
>@
°
°
°
°
°
¿
°
°
°
°
°
¾
½
°
°
°
°
°
¯
°
°
°
°
°
®
»
¼
º
«
¬
ª
EE
'
O
»
»
¼
º
«
«
¬
ª
E
E
'
O
»
»
¼
º
«
«
¬
ª
E
E
'
O
U
W'J
2
1,,,,,,1,,1,,,,
2
0p
2
,1,,,,,,1,,1,
,,
2
0t
2
,,1,,,,,,1,,1
,,
2
0r
e
,,
1
,,
215,2731
2
15,2731
2
15,2731
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
TTTTTT
z
TTTTT
T
y
TTTTT
T
x
c
TT
.
2
0t0r
0cr
OO
O
.
0crp/cr0p
O O K ,
0cr
0p
p/cr
O
O
K
»
¼
º
«
¬
ª
w
Ww
O
w
w
»
¼
º
«
¬
ª
w
Ww
O
w
w
»
¼
º
«
¬
ª
w
Ww
O
w
w
Ww
Ww
U
z
zyxT
uT
zy
zyxT
uT
y
x
zyxT
uT
x
zyxT
uTuTc
,,,
,
,,,
,
,,,
,
,,,
),(,
pt
re
2
p
2
2
p
2
t
2
2
t
2
r
2
2
re
¸
¹
·
¨
©
§
w
w
w
Ow
w
w
O
¸
¸
¹
·
¨
¨
©
§
w
w
w
Ow
w
w
O
¸
¹
·
¨
©
§
w
w
w
Ow
w
w
O
Ww
w
U
z
T
T
z
T
y
T
T
y
T
x
T
T
x
TT
c
0
0,,, TzyxT
,
W W W W
m
,0,,,,0,,,,0 TyxTzxTzyT
.
>@
>@
>@
°
°
°
°
°
¿
°
°
°
°
°
¾
½
°
°
°
°
°
¯
°
°
°
°
°
®
»
¼
º
«
¬
ª
EE
'
O
»
»
¼
º
«
«
¬
ª
E
E
'
O
»
»
¼
º
«
«
¬
ª
E
E
'
O
U
W'J
2
1,,,,,,1,,1,,,,
2
0p
2
,1,,,,,,1,,1,
,,
2
0t
2
,,1,
,,,,,1,,1
,,
2
0r
e
,,
1
,,
215,2731
2
15,2731
2
15,2731
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
TTTTTT
z
TTTTT
T
y
TTTTT
T
x
c
TT
.
2
0t0r
0cr
OO
O
.
0crp/cr0p
O O K
,
0cr
0p
p/cr
O
O
K
»
¼
º
«
¬
ª
w
Ww
O
w
w
»
¼
º
«
¬
ª
w
Ww
O
w
w
»
¼
º
«
¬
ª
w
Ww
O
w
w
Ww
Ww
U
z
zyxT
uT
zy
zyxT
uT
y
x
zyxT
uT
x
zyxT
uTuTc
,,,
,
,,,
,
,,,
,
,,,
),(,
pt
re
2
p
2
2
p
2
t
2
2
t
2
r
2
2
re
¸
¹
·
¨
©
§
w
w
w
Ow
w
w
O
¸
¸
¹
·
¨
¨
©
§
w
w
w
Ow
w
w
O
¸
¹
·
¨
©
§
w
w
w
Ow
w
w
O
Ww
w
U
z
T
T
z
T
y
T
T
y
T
x
T
T
x
TT
c
0
0,,, TzyxT
,
W W W W
m
,0,,,,0,,,,0 TyxTzxTzyT
.
>@
>@
>@
°
°
°
°
°
¿
°
°
°
°
°
¾
½
°
°
°
°
°
¯
°
°
°
°
°
®
»
¼
º
«
¬
ª
EE
'
O
»
»
¼
º
«
«
¬
ª
E
E
'
O
»
»
¼
º
«
«
¬
ª
E
E
'
O
U
W'J
2
1,,,,,,1,,1,,,,
2
0p
2
,1,,,,,,1,,1,
,,
2
0t
2
,,1,
,,,,,1,,1
,,
2
0r
e
,,
1
,,
215,2731
2
15,2731
2
15,2731
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
TTTTTT
z
TTTTT
T
y
TTTTT
T
x
c
TT
.
2
0t0r
0cr
OO
O
.
0crp/cr0p
O O K ,
0cr
0p
p/cr
O
O
K
2
>@
°
°
°
°
°
¿
°
°
°
°
°
¾
½
°
°
°
°
°
¯
°
°
°
°
°
®
»
»
¼
º
«
«
¬
ª
E
»
»
¼
º
«
«
¬
ª
E
'U
W'JO
2
1,,,,p/cr
2
,1,,,
2
,,1,,
,,p/cr1,,1,,p/cr
,1,,1,,,1,,1
,,
2
e
0cr
,,
1
,,
)(
)()(
)24()(
)15,273(1
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
TTK
TTTT
TKTTK
TTTT
T
xc
TT
.
0
0
,,
TT
kji
.
1
m
1
1,,
1
,1,
1
,,1
nn
ji
n
ki
n
kj
TTTT .
Deliiski: 3D Modeling and Visualization of Non-Stationary Temperature ... ...............
296 DRVNA INDUSTRIJA 64 (4) 293-303 (2013)
Figure 1 Positioning of nodes in a 3D calculation mesh of a discretized wooden prism
Slika 1. Pozicioniranje čvorova u 3D računskoj mreži u diskretiziranoj drvenoj prizmi
prismatic materials during their thermal treatment. The
values of the coefcient
in equation (8) have
been determined through the solution of the model
with the same initial and boundary conditions in order
to achieve maximum conformity between the calculat-
ed and experimental results.
It has been determined that the coefcient
p/cr
K
has the following values: for oak K
p/cr
=1.76, for beech
K
p/cr
=1.88, for poplar K
p/cr
=2.03, and for pine K
p/cr
=2.26
(Deliiski, 2003a, 2011).
3 RESULTS AND DISCUSSION
3. REZULTATI I RASPRAVA
3.1. Computation of 3D temperature distribution in
frozen wood during its defrosting
3.1. Izračun 3D raspodjele temperature u smrznutome
drvnom materijalu tijekom njegova odmrzavanja
For the numerical solution of the above presented
mathematical model, a software package has been de-
veloped in FORTRAN and integrated in the calculation
environment of Visual Fortran Professional developed
by Microsoft, as a part of the Windows Ofce software
(Deliiski, 2011).
With the help of this software package, 3D tem-
perature changes of beechwood (Fagus Silvatica L.)
and oakwood (Quercus petraea Liebl.) prisms with di-
mensions d = 0.4 m, b = 0.4 m, L = 0.8 m, initial tem-
perature of t
0
= -40 °C and two values of wood mois-
ture content u = 0.3 kg·kg
-1
and u = 0.6 kg·kg
-1
have
been studied during their 20 h heating with the inter-
mediate stage of melting at the heating temperature of
t
m
= 80 °C. The prisms with u = 0.3 kg·kg
-1
contain the
maximum possible quantity of frozen bound water in
beech and oak wood and contain no ice in the cell lu-
mens (i.e. contain no ice from free water). The prisms
with u = 0.6 kg·kg
-1
not only contain frozen bound wa-
ter but also contain a signicant quantity of frozen free
water.
The heating medium temperature, t
m
, increases
exponentially from t
m0
= t
0
to t
m
= 80 °C = const with
the time constant of 1800 s. This increasing of t
m
at the
beginning of the heating process of prisms can be seen
in Fig. 4 and 5. The values of d, b, L, t
m
, and u have
been selected so as to correspond to cases often en-
countered in practice.
The duration of 20 h of the prism heating at t
m
=
80 °C has been proven suitable for complete melting of
the ice in the studied prisms. The calculations have
been done with average values of r
b
= 560 kg·m
-3
and
= 0.31 kg·kg
-1
of the beech wood and of r
b
= 670
kg·m
-3
and
= 0.29 kg·kg
-1
of the oak wood (Vide-
lov, 2003; Deliiski and Dzurenda, 2010).
The computations have been carried out in a step
on the spatial coordinates
xD
= 0.001 m = 10 mm, i.e.
with the nodes M = 1 + [d/(2
xD
)] = 21 and N = 1 + [b/
(2
xD
)] = 21 along the x and y coordinates, respective-
ly, and KD = 1 + [L/(2
xD
)] = 41 along the z coordi-
nate. This means that the calculation meshes in the vol-
ume of the prisms consist of 21 x 21 х 41 = 18 081
nodes in total.
The step on time coordinate,
tD
, which is deter-
mined by the software package that keeps the stability
condition (Deliiski, 2011) of 3D solution of the explic-
it form of the nite-difference method and takes into
account the maximum values of λ and с
е
during wood
defrosting process, is as follows:
▪ for beech wood: Δτ = 30 at u = 0.3 kg·kg
-1
and Δτ =
25 at u = 0.6 kg·kg
-1
;
▪ for oak wood: Δτ = 40 at u = 0.3 kg·kg
-1
and Δτ = 30
at u = 0.6 kg·kg
-1
.
It takes 30 to 45 s to compute the temperature
distribution in the volume of each of the studied prisms
during a 20 h thermal treatment using the above values
2
>@
°
°
°
°
°
¿
°
°
°
°
°
¾
½
°
°
°
°
°
¯
°
°
°
°
°
®
»
»
¼
º
«
«
¬
ª
E
»
»
¼
º
«
«
¬
ª
E
'U
W'JO
2
1,,,,p/cr
2
,1,,,
2
,,1,,
,,p/cr1,,1,,p/cr
,1,,1,,,1,,1
,,
2
e
0cr
,,
1
,,
)(
)()(
)24()(
)15,273(1
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
n
kji
TTK
TTTT
TKTTK
TTTT
T
xc
TT
.
0
0
,,
TT
kji
.
1
m
1
1,,
1
,1,
1
,,1
nn
ji
n
ki
n
kj
TTTT .
............... Deliiski: 3D Modeling and Visualization of Non-Stationary Temperature ...
DRVNA INDUSTRIJA 64 (4) 293-303 (2013) 297
for
tD
with the help of Intel Pentium (4) CPU 3.0
GHz processor. Using the input data for solving the
model, the value for the interval (INT) is given in sec-
onds. After completing each INT from the beginning of
the process, the calculated temperature distribution in
the prism volume is recorded on computer hard-drive.
The records can be consequently seen on a monitor,
graphically processed, and/or printed. Besides taking
into account the stability condition for solving the 3D
model, the value of the step
tD
is calculated so as to
be divisible by the input value of INT, using the soft-
ware package.
Fig. 2 and 3 show the tables with the computed
temperature distribution in 121 nodes of the calcula-
tion mesh in the central cross-section of the beech
prisms at every 5 h of the defrosting process.
Fig. 4 and 5 shows the temperature change of the
surface of beech and oak prisms subjected to defrost-
ing, which is equal to t
m
, as well as of t in 6 character-
istic points of their volume.
The rst three characteristic points with coordi-
nates (d/4, b/8, L/8), (d/4, b/4, L/4), and (d/4, b/4, L/2)
allow for the tracking of the inuence on the defrosting
process of the gap from the prisms base (see Fig. 1 –
Figure 2 Change in t in the nodes of the calculation mesh, situated in the central cross section of a beech prism with
dimensions 0.4 x 0.4 x 0.8 m and u = 0.3 kg·kg
-1
during every 5 h of defrosting at t
m
= 80°C
Slika 2. Promjene temperature u čvorovima računske mreže smještenima na središnjemu poprečnom presjeku bukove prizme
dimenzija 0,4 x 0,4 x 0,8 m i u = 0,3 kg·kg
-1
tijekom svakih 5 h odmrzavanja pri temperaturi t
m
= 80°C