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