Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

52, 4, pp. 927-935, Warsaw 2014

APPLICATION OF THE CONTROL VOLUME METHOD USING THE

VORONOI POLYGONS FOR NUMERICAL MODELING OF BIO-HEAT

TRANSFER PROCESSES

Mariusz Ciesielski

Czestochowa University of Technology, Częstochowa, Poland; e-mail: mariusz.ciesielski@icis.pcz.pl

Bohdan Mochnacki

Higher School of Labour Safety Management, Katowice, Poland; e-mail: bmochnacki@wszop.edu.pl

In the paper, the problem concerning the numerical modeling of thermal processes in the
domain of a biological tissue being in thermal contact with the environment is discussed.
The changing ambient temperature causes that the non-steady heat transfer process is con-
sidered. The cross-section of the forearm (2D problem) is treated as a non-homogeneous
domain in which the sub-domains of skin tissue, fat, muscle and bone are distinguished.
From themathematical point of view, the boundary-initial problemdescribedby the system
of energy equations (the Pennes equations), boundary conditions on the external surface of
the system, boundary conditions on the surfaces limiting the successive sub-domains and the
initial condition is analyzed. At the stage of numerical computations, the Control Volume
Method using the Voronoi polygons is applied. In the final part of the paper, examples of
computations are shown.

Keywords:bioheat transfer,numericalmodeling,ControlVolumeMethod,Voronoi’sdiagram

1. Introduction

The thermal processes proceeding in the domain of a biological tissue are, as a rule, described
by the well known Pennes equations (e.g. Mochnacki and Majchrzak, 2003; Majchrzak, 2013).
The typical heat diffusion equation is supplemented by components determining the capacities
of internal heat sources connected with the blood perfusion andmetabolism. Themathematical
form of perfusion heat source results from the assumption that the tissue is supplied by a large
number of blood capillaries uniformly distributed in the area under consideration. In the case
of tissue freezingmodeling, the right hand side of the Pennes equation should be supplemented
by the third internal heat source controlling the evolution of latent heat (e.g. Majchrzak et al.,
2011). The Pennes equation constitutes a basis of the so-called tissue models (Mochnacki and
Majchrzak, 2003; Majchrzak, 2013; Majchrzak et al., 2011; Mochnacki and Piasecka-Belkhayat,
2013). Atkin et al. (1994) talk over different bio-heat transfer tissue models and conclude that
the Pennes equation is the best approach to the modeling of bio-heat transfer because of its
simplicity.

In literature, one can also find the so-called vascular models (Majchrzak, 2013; Zhu and
Weinbaum, 1995; Wang et al., 2007). The vascular models are applied because of the need to
include the thermally significant blood vessels of considerable size. In this paper, the presence of
blood vessels (arteries and veins) in the domain considered is taken into account, but the blood
temperature is optionally established on the basis of literature data (it is not calculated), and
the model presented belongs rather to the tissue ones.
There have recently been attempts to use descriptions of thermal processes in the biological

tissue on the basis of porous media theory. The domain considered is divided into two parts,
that this the vascular and extravascular region (cells and interstitial tissue). The blood vessels



928 M. Ciesielski, B. Mochnacki

are treated as voids, while the remaining area as thematrix (solid body). The equations concer-
ning the vascular and extravascular regions contain the parameter ε (porosity) (Zhang, 2009;
Nakayama andKuwahara, 2008; Khaled andVafai, 2003). Assuming certain simplifications, the
above mentioned equations can be substituted by a single equation close to the dual phase lag
model (Khaled and Vafai, 2003).

In this paper, the cross section of the forearm (middle part) shown in Fig. 1 (Schuenke et
al., 2010) is considered. The domain is non-homogeneous and constitutes a composition of skin,
fat, muscle, bone and blood vessels.

Fig. 1. Forearm cross section and a simplified 2D geometrical model

The mathematical model of thermal processes proceeding in the area considered subjected
to the time-dependent external heat source is presented in Section 2.

At the stage of numerical modeling, the Control VolumeMethod using the Voronoi tessella-
tion is used (Domanski et al., 2010). The Voronoi polygons are characterized by the geometric
properties of a well meeting the requirements for the shape of control volumes. The details of
the approach proposed are discussed in Section 3. In the final part of the paper, examples of
computations are presented.

2. Governing equations

Heat transfer processes proceeding in the tissue domain are described by the system of the
Pennes equations of the following form

ce(T)
∂Te(x,t)

∂t
=∇[λe(T)∇Te(x,t)]+Qpere(T)+Qmete(T) e=1, . . . ,4 (2.1)

where e = 1, . . . ,4 identifies the tissue sub-domains (skin, fat, muscle and bone respectively),
ce is the volumetric specific heat, λe is the thermal conductivity, Qper and Qmet are the capa-
cities of volumetric internal heat sources connected with the blood perfusion and metabolism,
[W/m3], T , x = {x1,x2}, t denotes temperature, spatial co-ordinates and time, respectively.
The perfusion heat source is given by the formula

Qpere(T)= cbGbe(T)[Tb−Te(x,t)] Tb =
1

2
(Tbartery +Tbvein) (2.2)

where Gbe is the blood perfusion [m
3blood/(sm3 tissue)], cb is the blood volumetric specific

heat and Tbartery and Tbvein are the arterial and vein blood temperatures. Themetabolic heat
source Qmet can be treated both as a constant value and a temperature-dependent function.



Application of the Control Volume Method using the Voronoi polygons... 929

On the contact surface between the tissue sub-domains, the IV type of boundary conditions
are assumed

x∈Γk−l :











−λk
∂Tk(x,t)

∂n
=−λl

∂Tl(x,t)

∂n

Tk(x,t)=Tl(x,t)

(k,l)∈{(1,2),(2,3),(3,4)} (2.3)

where ∂/∂n denotes the normal derivative.

On the outer surface of the skin (e = 1), the III type of boundary condition is taken into
account

x∈Γout : −λ1
∂T1(x,t)

∂n
=−αout[Tamb(t)−T1(x,t)] (2.4)

where αout is the heat transfer coefficient, Tamb is the ambient temperature. The same type of
boundary conditions is given on the contact surfaces between the blood vessels and soft tissue
sub-domains, in particular

x∈Γartery : −λ3
∂T3(x,t)

∂n
=−αartery[Tbartery −T3(x,t)] (2.5)

and

x∈Γvein : −λe
∂Te(x,t)

∂n
=−αvein[Tbvein−Te(x,t)] e= {2,3} (2.6)

The initial conditions are also given

t=0 : Te(x,t)=Tsteady(x) e=1, . . . ,4 (2.7)

where Tsteady is the temperature distribution corresponding to the steady state conditions in
the domain considered for the given ambient temperature and the initial external heat transfer
coefficient.

3. Numerical algorithm

For the purpose of heat transfer modeling in the domain considered, the tissue sub-domains
are divided into small cells (control volumes) known as the Voronoi polygons (also called the
Thiessen or Dirichlet cells in two dimensions) (Okabe et al., 2000). The polygon that contains
the point xi (central point) is denoted by ∆Vi, see Fig. 2. All of theVoronoi regions are convex
polygons, and each polygon is defined by lines that bisect the sectors between the central point
and its neighbouring points. The bisecting lines and the connection lines are perpendicular to
each other (it is very convenient at the stage of CVMequations construction).Whenwe use this
rule for every point in the area, the area will be completely covered by adjacent polygons.Many
algorithms to construct the Voronoi polygons can be found in literature, see e.g. the Delaunay
triangulation (Watson, 1981).

In Fig. 3, an example of the control volumemesh in the domain of the forearm cross-section
is shown. The domain is divided into 2966 control volumes. The positions of CV central points
close to the contact surface between the tissue sub-domains are analytically determined in order
to achieve a better approximation of the tissue shape.

The control volumemethod (CVM)constitutes an effective tool for numerical computation of
heat transfer processes. The domain analyzed is divided into N volumes. The CVM algorithm
allows one to find the transient temperature field at the set of nodes corresponding to the



930 M. Ciesielski, B. Mochnacki

Fig. 2. The Voronoi polygons for a set of arbitrarily distributed points

Fig. 3. Control volumemesh: cross-section of the forearm

Fig. 4. Control volume ∆Vi

central points of the control volumes. The nodal temperatures can be found on the basis of
energy balances for the successive volumes.

Let us consider the control volume ∆Vi with the central node xi. It is assumed here that
the thermal capacities and capacities of the internal heat sources are concentrated at the nodes
representing the elements, while the thermal resistances are concentrated on the sectors joining
thenodes.The energybalances corresponding to theheat exchange between the analyzed control
volume ∆Vi and the adjoining control volumes results from the integration of energy equation
(2.1) with respect to time and volume CVi. Let us consider the interval of time ∆t= t

f+1−tf.



Application of the Control Volume Method using the Voronoi polygons... 931

Then

tf+1
∫

tf

∫

CVi

ce(T)
∂Te(x,t)

∂t
dV dt=

tf+1
∫

tf

∫

CVi

∇[λe(T)∇Te(x,t)] dV dt

+

tf+1
∫

tf

∫

CVi

[Qpere(T)+Qmete(T)] dV dt

(3.1)

Using Gauss-Ostrogradsky’s theorem, one obtains

tf+1
∫

tf

∫

CVi

ce(T)
∂Te(x,t)

∂t
dV dt=

tf+1
∫

tf

∫

Ai

n · [λe(T)∇Te(x,t)] dAdt

+

tf+1
∫

tf

∫

CVi

[Qpere(T)+Qmete(T)] dV dt

(3.2)

where Ai is the surface (perimeter) limiting the CVi. The approximation of the left-hand side
of equation (3.2) can be taken in form

tf+1
∫

tf

∫

CVi

ce(T)
∂Te(x,t)

∂t
dV dt∼= c

f
i (T
f+1
i −T

f
i )∆Vi (3.3)

where c
f
i is an integralmean of thermal capacity, and this value is approximated by the volume-

tric specific heat corresponding to the temperature T
f
i . In a similar way, one can approximate

the last component in equation (3.2), namely

tf+1
∫

tf

∫

CVi

(Qpere(T)+Qmete(T)) dV dt∼= [(Qper)
f
i +(Qmet)

f
i ]∆Vi∆t (3.4)

The term determining the heat conduction between ∆Vi and its neighbourhoods can bewritten
in form

tf+1
∫

tf

∫

Ai

n · [λe(T)∇Te(x,t)] dAdt=

tf+1
∫

tf

(

ni
∑

j=1

∫

Ai(j)

ni(j) · [λ(T)∇T(x,t)]i(j) dAi(j)

)

dt

=

tf+1
∫

tf

(

ni
∑

j=1

ni(j) · [λ(T)∇T(x,t)]i(j)Ai(j)

)

dt∼=

tf+1
∫

tf

(

ni
∑

j=1

λij
Ti(j)−Ti

hi(j)
Ai(j)

)

dt

=
ni
∑

j=1

λ
f
ij

T
f
i(j)
−T

f
i

hi(j)
Ai(j)∆t

(3.5)

where λij is the mean thermal conductivity between the nodes i and i(j), in particular

λij =
2λiλi(j)

λi+λi(j)
(3.6)



932 M. Ciesielski, B. Mochnacki

Such an assumption causes that in the denominators of conductional terms, the well known
thermal resistances appear. The energy balance written in convention of ‘explicit’ scheme takes
the form

c
f
i (T
f+1
i −T

f
i )∆Vi =

ni
∑

j=1

λ
f
ij

T
f
i(j)
−T

f
i

hi(j)
Ai(j)∆t+

[

(Qper)
f
i +(Qmet)

f
i

]

∆Vi∆t (3.7)

fromwhich (introducing (Qper)
f
i =(Gb)

f
i cb(Tb−T

f
i ))

T
f+1
i =T

f
i +

∆t

c
f
i∆Vi

ni
∑

j=1

λ
f
ij

T
f
i(j)
−T

f
i

hi(j)
Ai(j)+

∆t

c
f
i

[

(Gb)
f
i cb(Tb−T

f
i )+(Qmet)

f
i

]

(3.8)

or

T
f+1
i =T

f
i +

ni
∑

j=1

Wi(j)(T
f
i(j)
−T

f
i )+

∆t

c
f
i

[

(Gb)
f
i cb(Tb−T

f
i )+(Qmet)

f
i

]

(3.9)

where

Wi(j) =
λ
f
ijAi(j)∆t

c
f
ihi(j)∆Vi

(3.10)

The stability condition

1−
ni
∑

j=1

Wi(j)−
∆t(Gb)

f
i cb

c
f
i

> 0 (3.11)

(for all i) allows one to determine the critical time step.
In the case of external CV (the boundary of CVi between the nodes i and i(j) lies on the

surface Γout, Γartery and Γvein), the following approximation of conductional term in equation
(3.5) is used

ni(j) · [λ(T)∇T(X,t)]i(j)Ai(j)
∼=
Ta−Ti
hi(j)
2λi
+ 1
α

Ai(j) (3.12)

where Ta = {Tamb,Tbartery,Tbvein} and α= {αout,αartery,αvein}, respectively. Then, the coef-
ficient Wi(j) in formula (3.10) takes the form

Wi(j) =
Ai(j)∆t

c
f
i

(

hi(j)
2λi
+ 1
α

)

∆Vi
(3.13)

and simultaneously T
f
i(j)
=Ta. One can see that the denominator in formula (3.13) corresponds

to the thermal resistance between the central point of CVi and its environment in the i(j) di-
rection.

4. Example of computations

The forearm domain stays in the thermal contact with the environment whose temperature is
equal to Tamb = 20

◦C, while the heat transfer coefficient is αout = 3.7W/(m
2K). The initial

temperature distribution is found using theGaussmethod (simple iterationmethod). The ther-
mophysical parameters of the successive sub-domains are taken from Fiala et al. (1999), see
Table 1.



Application of the Control Volume Method using the Voronoi polygons... 933

Table 1.Thermal parameters of the tissues

λ [W/mK] ρ [kg/m3] c [J/(kgK)] Gb [1/s] Qmet [W/m
3]

Bone 0.75 1357 1700 0.0000 ·10−3 0

Muscle 0.42 1085 3768 0.5380 ·10−3 684

Fat 0.16 850 2300 0.0036 ·10−3 58

Skin 0.47 1085 3680 1.1000 ·10−3 368

Blood 1069 3650

Additionally, the followingblood temperatures are assumed: Tbartery =36
◦C,Tbvein =35

◦C,
while αartery =αvein =5000W/(m

2K).At themoment t=0theambient temperature increases
to 40◦C (the heat transfer coefficient is the same). Such a situationmay correspond to the input
from a room outside in a summer hot day. The numerical simulation concerns the process of
tissue heating.
In Fig. 5, the initial temperature distribution is shown. The next figures present temporary

solutions for times 15min. and 300min. (this result practically corresponds to the steady state
conditions).

Fig. 5. The initial temperature distribution

Fig. 6. Temporary solutions for 15min



934 M. Ciesielski, B. Mochnacki

Fig. 7. Temporary solutions for 300min

5. Final remarks

In the paper, a problem fromthe scope of bio-heat transfermodeling is discussed.Themathema-
tical model belongs to the group of tissuemodels basing on the Pennes equation. The numerical
solution is obtained using the developed by the authors version of the Control VolumeMethod
using the Voronoi tessellation. According to our opinion, the Voronoi polygons characterized
by geometric properties perfectly meet the conditions required by the 2D control volumes. The
mesh generation and the computations of the control volumes geometrical parameters are a dif-
ficult task, but the authors have developed procedures to carry out such investigations. In the
case discussed, the consideration of internal heat sources has been necessary.
The numerical simulation of forearm heating in conditions of natural convection is only an

example of possible application of the method proposed in the scope of bio-heat transfer. The
algorithm presented can be, among others, used for numerical modeling of thermal processes
proceeding in the domain of a living tissue subjected to the strong external thermal interactions
(e.g. modeling of frostbites or burns). The shape and location of tissue sub-domains can be
accurately reproduced (if the assumption of 2D approximation is acceptable). The application
of the method for numerical analysis of the vascular models seems also be possible.

Acknowledgement

The paper has been a part of Project PB3/2013 sponsored by WSZOPKatowice.

References

1. ArkinH., Xu L.X,HolmesK.R., 1994,Recent development inmodeling heat transfer in blood
perfused tissues, IEEE Transactions on Biomedical Engineering, 41, 97-107

2. Domanski Z., Ciesielski M., Mochnacki B., 2010, Applications of control volume method
using theVoronoi tessellation in numericalmodeling of solidification process, [In:]Current Themes
in Engineering Science 2009, KorsunskyA. (Ed.), BookSeries:AIPConferenceProceedings,1220,
17-26

3. Fiala D., LomasK.J., Stohrer M., 1999,A computermodel of human thermoregulation for a
wide range of environmental conditions: the passive system, Journal of Applied Physiology, 87, 5,
1957-1972

4. Khaled R.A., Vafai K., 2003, The role of porous media in modeling flow and heat transfer in
biological tissues, International Journal of Heat and Mass Transfer, 46, 4989-5003



Application of the Control Volume Method using the Voronoi polygons... 935

5. Majchrzak E., 2013, Application of different variants of the BEM in numerical modeling of
bio-heat transfer processes,MCB: Molecular and Cellular Biomechanics, 10, 3, 201-232

6. Majchrzak E., Mochnacki B., Dziewonski M., Jasinski M., 2011, Numerical modeling of
hyperthermia and hypothermia processes,Advanced Materials Research, 268-270, 257-262

7. Mochnacki B., Majchrzak E., 2003, Sensitivity of the skin tissue on the activity of external
heat sources,CMES – Computer Modeling in Engineering and Sciences, 4, 3/4, 431-438

8. MochnackiB.,Piasecka-BelkhayatA., 2013,Numericalmodeling of skin tissue heating using
the interval finite difference method,MCB: Molecular and Cellular Biomechanics, 10, 3, 233-244

9. Nakayama A., Kuwahara F., 2008, A general bioheat transfer model based on the theory of
porous media, International Journal of Heat and Mass Transfer, 51, 3190-3199

10. Okabe A., Boots B., Sugihara K., Chiu S.N., 2000, Spatial Tessellations: Concepts and
Applications of Voronoi Diagrams, Second Edition,Wiley

11. Schuenke M., Schulte E., Schumacher U. Ross L., Lamperti E., 2010,General Anatomy
and Musculoskeletal System (Atlas of Anatomy), Thieme, Stuttgart-NewYork

12. WangH., DaiW., BejanA., 2007,Optimal temperature distribution in a 3D triple-layered skin
structure embedded with artery and vein vasculature and induced by electromagnetic radiation,
International Journal of Heat and Mass Transfer, 50, 1843-1854

13. Watson D.F., 1981, Computing the n-dimensional Delaunay tessellation with application to Vo-
ronoi polytopes,The Computer Journal, 24, 2, 167-172

14. Zhang Y., 2009, Generalized dual-phase lag bioheat equations based on nonequilibrium he-
at transfer in living biological tissues, International Journal of Heat and Mass Transfer, 52,
4829-4834

15. Zhu L., Weinbaum S., 1995, A model for heat transfer from embedded blood vessels in two
dimensional tissue preparations, Journal of Biomechanical Engineering, 117, 64-73

Manuscript received April 1, 2014; accepted for print April 29, 2014