Finite Elements in Analysis and Design 135 (2017) 11–21

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design

journal homepage: www.elsevier.com/locate/finel

Sensitivity analysis for radiofrequency induced thermal therapies using the MARK
complex finite element method
⁎
Juan F. Monsalvoa, Manuel J. Garcíaa,b, , Harry Millwaterb, Yusheng Fengb
a
Applied Mechanics Research Group, Universidad EAFIT, Medellín, Colombia
b
Mechanical Engineering Department, University of Texas at San Antonio, TX, USA

A R T I C L E I N F O A BS T RAC T

Keywords: In radiofrequency induced thermal procedures for cancer treatment, the temperature of the cancerous tissue is
Sensitivity analysis raised over therapeutic values while maintaining the temperature of the surrounding tissue at normal levels. In
Complex Taylor series expansion order to control these temperature levels during a thermal therapy, it is important to predict the temperature
Bioheat transfer equation distribution over the region of interest and analyze how the variations of the different parameters can affect the
Radio frequency ablation
temperature in the healthy and damaged tissue. This paper proposes a sensitivity analysis of the radiofrequency
Liver cancer treatment
induced thermal procedures using the complex Taylor series expansion (CTSE) finite element method (ZFEM),
which is more accurate and robust compared to the finite difference method. The radiofrequency induced
thermal procedure is modeled by the bioheat and the Joule heating equations. Both equations are coupled and
solved using complex-variable finite element analysis. As a result, the temperature sensitivity with respect to any
material property or boundary condition involved in the process can be calculated using CTSE.
Two thermal therapeutical examples, hyperthermia and ablation induced by radio frequency, are presented
to illustrate the capabilities and accuracy of the method. Relative sensitivities of the temperature were computed
for a broad range of parameters involved in the radiofrequency induced thermal process using ZFEM. The major
feature of the method is that it enables a comprehensive evaluation of the problem sensitivities, including both
model parameters and boundary conditions. The accuracy and efficiency of the method was shown to be
superior to the finite difference method. The computing time of a complex finite element analysis is about 1.6
times the computing time of real finite element analysis; significantly lower than the 2 times of forward/
backward finite differencing or 3 times of central differencing. It was found that the radiofrequency
hyperthermia procedure is very sensitive to the electric field and temperature boundary conditions. In the
case of the radiofrequency ablation procedure, the cooling temperature of the electrodes has the highest liver/
tumor temperature sensitivity. Also, thermal and electrical conductivities of the healthy tissue were the
properties with the highest temperature sensitivities. The result of the sensitive analysis can be used to design
very robust and safe medical procedures as well as to plan specific patient procedures.

1. Introduction The goal of a radiofrequency induced thermal treatment is to raise

the temperature of the tumor to a value in the range from (42 − 46)°C
A radiofrequency induced thermal procedure is a minimally in- for hyperthermia therapy, or to a value in the range from (50 − 100)°C
vasive clinical method for the treatment of hepatocellular carcinoma for radiofrequency ablation therapy. The temperature field is applied
that consists of the heating of biological tissues through the emission of for certain time so that the tumor tissues are destroyed while keeping
an external electric field. The radiofrequency induced thermal proce- the temperature in the healthy tissue region below critical temperature
dure is used when the tumor or other dysfunctional tissues are not to avoid unwanted damage [1–3]. The success of either hyperthermia
surgically resectable due to the location or the potential of significant and ablation induced by radio frequency depend on: (i) how well the
collateral damage in the surgical process. Radio frequency ablation can temperature field is controlled and optimized to ensure the tempera-
be used to treat tumors in the liver as well as those in other organ sites ture is high enough to destroy the cancer cells in the tumor region (ii)
such as lung, kidney, pancreas and bone. how to minimize the damage in the healthy tissues surrounding the

⁎
Corresponding author at: Applied Mechanics Research Group, Universidad EAFIT, Medellín, Colombia.
E-mail addresses: jmonsa13@eafit.edu.co (J.F. Monsalvo), manuel.garciaruiz@utsa.edu (M.J. García), harry.millwater@utsa.edu (H. Millwater),
yusheng.feng@utsa.edu (Y. Feng).

http://dx.doi.org/10.1016/j.finel.2017.07.001
Received 14 February 2017; Received in revised form 1 June 2017; Accepted 2 July 2017
Available online 13 July 2017
0168-874X/ © 2017 Elsevier B.V. All rights reserved.
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

tumor with reasonable margin prescribed by an oncologist or surgeon.

The thermal modeling of the radiofrequency induced thermal
treatment has been studied by several researchers with Pennes [4]
the first to introduce a mathematical equation explaining the heat
transfer present in a biological tissue. This equation is known as the
bioheat transfer equation. Tunç et al. [5] applied the bioheat transfer
equation in the analysis of hyperthermia treatments. Chang and
Nguyen developed a model to predict the size of the thermal lesions
created by the radiofrequency induced thermal procedure using the
Arrhenius equation which relates the temperature and the exposure
time, using a first order kinetics relationship [6]. Fuentes et al.
compared the thermal damage region, recreated in a computer
simulation, with the radiofrequency induced thermal procedure in an Fig. 1. A simplified 2D model of a radiofrequency induced thermal treatment
in vitro perfused bovine liver model, finding a good correlation (hyperthermia) with an embedded tumor.
between the two damage regions [1]. In addition, the bioheat transfer
in some organs may need a hyperbolic term to account for finite time ∂Ω = Γ1 + Γ2 + Γ3 + Γ4 . The electric field is generated from two electro-
dispersion of heat energy [7]. des located on the top Γ1 and bottom Γ2 of the healthy tissue,
Despite significant research efforts, the prediction of the tempera- generating a voltage gradient. The top and bottom of the healthy tissue
ture distribution is still an important issue in the reliability and Γ3 are covered with a cooling pad that preserves a constant tempera-
effectiveness of the treatment. Since the human body is a complex ture at the surface. The left and right side Γ4 are assumed insulated
blend of tissue with heterogeneous electric and thermal characteristics, surfaces. A detail description of the geometry is presented in Section 4.
it is a difficult task to predict and control the temperature distribution The radiofrequency induced thermal treatment can be described as
and to determine the optimum heating configuration during therapy. the coupling of two physical phenomena. The first is the heat induced
Recent advances in this direction include the development of a patient- by the electric field generated from the electrodes; the second is the
specific image-based 3D model for thermal ablation simulation to heat transfer in the tissue and considers the metabolic heat, perfusion
optimize treatment efficacy [8]. Also, Schumann et al. proposed an rate and the heat produced by the electric field. The combination of
access path determination method, for radiofrequency ablation plan- these two models allows one to predict the temperature field distribu-
ning, based on image processing and numerical optimization [9]. tion on the domain (organ). A quasi-static electric field approximation
Therefore, it is important to develop a tool capable of predicting the was used due to the fact that the wave length of the radiofrequency is
thermal results given the variation of the material properties and the much greater compared to the depth of the computational domain,
radiofrequency induced thermal variables. Hence, a thorough sensitiv- therefore, the free charge continuity equation takes the following form
ity analysis of each variable involved in the process will allow one to [28],
identify those variables with the largest impact on the tumor tempera-
∇·J (x) = 0, (1)
ture distribution. The Complex Taylor Series Expansion (CTSE)
method was used here to perform this analysis. It was first proposed where the electric current density J = J (x)[A/m2] is related to the
by [10] and [11]. [12] were the first to obtain a method to approximate electric field through the Ohm's law
derivatives of real functions using complex variables. The CTSE
∇·σ (x) E (x) = 0. (2)
method calculates the first derivative with respect a parameter by
perturbing the parameter along the imaginary axes. There is no Here, σ = σ (x) [s/m ] is the tissue electrical conductivity and E = E (x)
difference operation involved, and therefore it avoids the cancellation [V/ m ] denotes the electric field. The electric field can be defined as the
error presented in the Finite Difference Method (FDM). This allows the gradient of the electric potential ϕ (x) [V]
method to be step size independent, easy to implement and highly E (x) = ∇ϕ (x). (3)
accurate. The method has proven to be very attractive in applications
where a sensitivity analysis is required, and it has been applied in Substituting Eq. (3) into Eq. (2) yields the governing electric potential
several fields of engineering as: the computation of sensitivity deriva- equation,
tives of NavierStokes equations [13], shape sensitivity analysis [14], ∇·σ (x)∇ϕ (x) = 0. (4)
fatigue sensitivity analysis [15], fracture mechanics [16–19], nonlinear
problems in plasticity [20] and creep [21], aerodynamics [22,23], When the electric conductivity σ is constant over the domain, Eq. (4) is
structural dynamics [24,25], and nonlinear structural analysis [26], transformed into the Laplace equation σ ∇2 ϕ = 0 . However, this is not
among others. However, the CTSE has not been exploited in bioengi- the case in the present study as the healthy tissue and tumor tissue
neering or medical research. have different physical properties. The boundary conditions necessary
In this paper, a simplified 2D geometry of the radiofrequency to solve Eq. (4) are as follows: a Dirichlet boundary condition with the
induced thermal process was modelled using the bioheat transfer value of the electric potential relative to ground is specified on each
∂ϕ
equation coupled with the Joule heating equation. The sensitivity was electrode Γ1, Γ2. A Newman boundary condition equal to zero, ∂n = 0
calculated using CTSE, and the results were compared with the representing an ideal electrical isolation, is specified on the remaining
sensitivity found using the finite difference method in order to external boundaries Γ3, Γ4, see Fig. 1. The differential in electric
demonstrate the superior step size independence of the CTSE method. potential produces an electric current through the tissue and releases
heat. This process is called Joule heating and is also referred to as
2. Radiofrequency induced thermal process Ohmic or resistive heating. The amount of heat released is proportional
to the square of the current, or in terms of the electric field
A simplified 2D model of the radiofrequency ablation, adopted and Q J (x ) = σ (x) E (x)·E (x) or
modified from [27], is shown in Fig. 1. The domain Ω is composed of
Q J (x ) = σ (x)∥∇ϕ (x)∥2 . (5)
two subdomains, the healthy tissue domain Ω1, and the tumor domain
Ω2. A point x = (x, y) ∈ Ω , will have different properties depending on On the other hand, the temperature distribution inside the healthy and
whether it is located in the tumor, x ∈ Ω2 , or in the healthy tissue, the tumor tissue is described by the bioheat transfer equation also
x ∈ Ω1. The boundary of the domain ∂Ω is composed of four parts known as the Pennes equation [4]. This equation has an added source

12
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

term Q J (x) [W/m ] which is a product of the Joule heating. The steady 3. Evaluation of sensitivities
state version of the bioheat transfer equation is
A sensitivity analysis consists of determining the impact of a
λ∇2 T (x) + ρb Cb Gb [Tb − T (x)] + Qm (x) + Q J (x) = 0, (6) parameter on the output of a system. These parameters can be tissue
in which T [K ] denotes the temperature, λ (x) [W/mK ] is the thermal properties, initial and boundary conditions. The selection of the
conductivity, ρb [kg/m3] denotes the blood density, Cb [J/kg K ] is the method for computing sensitivities is not an easy task and key issues
specific heat capacity of the blood, Gb [1/s ] defines the perfusion rate, to be considered are the accuracy and the computational cost. It is also
Tb [K] is the temperature of the blood, Qm (x) [W/m3] denotes the important that the differentiation scheme be easy to implement.
metabolic heat source and Q J (x) [W/m3] is the added heat source There are four main methods for computing the sensitivities. First,
product of the radiofrequency induced thermal treatment. the direct differentiation method where the partial derivatives are
In summary, a differential electric potential is applied to the derived analytically, the results from this method are exact. It is
electrodes which produces an electric potential in the tissue that is efficient only when the number of design variables is low (Choi and
solved using Eq. (4). Once the potential is computed, the heating Kim [29]). Second, the adjoint method that can be thought of as the
source due to the Joule heating effect is computed using Eq. (5). reverse of direct differentiation, where the derivative is calculated from
Finally, the temperature distribution is computed from Eq. (6) using the outputs rather than the inputs using the calculus of variations. This
the Joule heating source previously calculated. Fig. 1 shows the method is a feasible solution for problems with many design variables
boundary conditions of the thermal model: a constant temperature and few outputs where the direct differentiation becomes difficult to
along the cooling pad and an adiabatic condition on the remaining implement (Choi and Kim [29]). Third, automatic differentiation or
external boundaries. computational differentiation that consist of the differentiation of the
computer code itself. Automatic differentiation is based on the concept
of the chain rule and the fact that a high-level function is composed of
2.1. FEM formulation multiple smaller functions each having its own partial derivative and
applying the chain rule of differentiation [30]. Fourth, the Finite
The radiofrequency thermal therapy problem characterized by Eqs. Difference method (FDM) is a straightforward method for obtaining
(4)–(6), are solved using the Finite Element Method (FEM). This sensitivities as there is no need of code changes. The calculation of the
method is expressed first in real variables and then followed by a sensitivities using the FDM consists of perturbing the variable of
complex variable representation. interest and recomputing the output, in order to determine the ratio
The discrete form of Eq. (4) that describes the electric potential ϕ of the change. The CTSE method is similar to the finite difference
using the Galerkin method is given by method as it uses a Taylor series expansion to approximate the
Hϕ = q . (7) derivative. However, CTSE method uses complex variable and the
perturbation is made in the imaginary axis. It is easy to code, robust,
where ϕ is a vector of the electric potential at discrete points, q is a accurate and does not suffer from subtraction error.
vector of loads (electric potential applied), and H is the system matrix The Taylor series expansion of f (x + ih ) is defined as
defined by
f ″(x )
f (x + ih ) = f (x ) + ihf ′(x ) − h2 +  (h3).
∂N ∂Nj 2 (15)
Hij = ∫ σ i
Ω ∂xr ∂xr
dx.
(8) Taking the imaginary part of the f (x + ih ), that is, Im ( f (x + ih ))
where Ni (x) is the standard Lagrange basis defined over the nodes of neglecting the higher order terms and dividing by the step size (h)
the element. The r index stands for the coordinate axes (x1, x2 ) = (x, y), yields
and summation is assumed over repeated indices. In a similar way the Im( f (x + ih ))
discrete form of the bioheat transfer equation is given by f ′ (x ) = +  (h2 ),
h (16)
(K + M) T = qb + qm + q J, (9) which is a second-order approximation of the first derivative. In other
words, the information of the first derivative is found in the imaginary
where K is the matrix due to diffusion, M and qb the matrix and vector
part of the function. As can be seen, Eq. (16) does not have a difference
due to the blood the perfusion rate, qm the metabolic heat source and
operation between two numbers. The step size should be extremely
q J the Joule heating source These terms are defined as follows
small in order to recover the exact numerical derivative. Additionally,
∂Nj (x) the method is straightforward to implement in an existing FEM code.
Kij = ∫Ω λ (x) ∂N∂ix(rx) ∂xr
dx ,
(10) This can be accomplished by either using the native complex data type
in a user-written code or by incorporating user-defined finite elements
into a commercial code as it was done in [17]. Eq. (16) yields accurate
Mij = ∫Ω Bb Ni (x) Nj (x) dx, (11) numerical derivatives when h is very small, e.g., h = 10−10 .
Computing the derivative in this way requires that the variables are
q bi = ∫Ω Bb Tb Ni (x) dx, (12)
complex consisting of real and imaginary parts. This implies an
increase in the memory and the computation time. Despite these
apparent drawbacks, the accuracy, robustness and ease of implementa-
q mi = ∫Ω Qm (x) Ni (x) dx, (13) tion of the CTSE method make it a great option for calculating first
order derivatives.
q Ji = ∫Ω QJ (x) Ni (x) dx, (14)
3.1. Implementation of the CTSE method into a finite element code
where Bb = ρb Cb Gb. The energy that the tissue receives via the electric
field is dissipated by the effect of the electrical resistivity of the The conversion of a real-valued FEM solver into a complex valued
material. Once the electric field is found from Eq. (7), the heat can FEM solver involves the transformation of each real number into
be calculated by computing the gradient of the electric field, using Eq. complex. A complex number, a + ib , can be represented as a pair (a,b)
(5), and then the temperature computed by solving the bioheat and numerous languages support operations among complex numbers.
equation, Eq. (9). Also, there are linear algebra modules with solvers that are capable of

13
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

handling complex variables. However, with the aim to support multi- Table 2
complex analysis in the future, this implementation was made using Thermal properties of the human blood. (Data source from reference [36]).
the Cauchy-Riemann matrix representation for a complex number
Temperature Specific Heat Density
[31]. If z is a complex number and z = a + ib , then, Tb [K] Cb [J/kg K] ρb [kg/m3]
⎡ Re(z ) − Im(z )⎤ ⎡ a − b ⎤
z ≔⎢ ⎥=⎢ ⎥. 310.15 4180 1000
⎣ Im(z ) Re(z ) ⎦ ⎣ b a ⎦ (17)

In matrix form the real part is located along the diagonal of the matrix.
electric potential were as follows
Operations between complex numbers are transformed into standard
matrix operations. The CTSE implementation was based on a standard
finite element code using linear triangular elements. The real variables • ϕ1 = 10 [V] on Γ1, the upper electrode.
• ϕ1 = −10 [V] on Γ2, the lower electrode.
were turned into complex variables using the Cauchy-Riemann matrix
notation. For example, a real valued element matrix Ke for a linear • ∂ϕ /∂n = 0 on Γ3 + Γ4 , ideal insulation on the remaining external
boundaries.
element is a 3×3 real matrix; for a complex matrix form it will be:
⎡ Re(K11) Re(K12 ) Re(K13) − Im(K11) − Im(K12 ) − Im(K13) ⎤ For the bioheat transfer the following boundary conditions were
⎢ ⎥ implemented
⎢ Re(K21) Re(K22 ) Re(K23) − Im(K21) − Im(K22 ) − Im(K23) ⎥
⎢ Re(K31) Re(K32 ) Re(K33) − Im(K31) − Im(K32 ) − Im(K33)⎥
• T = 32.5°C on Γ , Γ and Γ ; the cooling pad.
e
K =⎢
Im(K11) Im(K12 ) Im(K13) Re(K11) Re(K12 ) Re(K13) ⎥ 1 1 2 3
⎢
⎢ Im(K21) Im(K22 ) Im(K23) Re(K21) Re(K22 ) Re(K23) ⎥
⎥
• ∂T /∂n = f = 0 on Γ , zero heat flux on the
4 remaining external
⎢⎣ Im(K31) Im(K32 ) Im(K33) Re(K31) Re(K32 ) Re(K33) ⎥⎦ boundaries.
(18)
The computation of sensitivities was accomplished using CTSE with a
Finally, the solutions of the different partial differential equations step size of h = 10−10 for all parameters.
will result in a vector as
⎡ Re(u1) ⎤ 4.1. Methodology
⎡ Re(u1) ⎤ ⎢ ⎥
⎢ ⎥ ⎢ Re(u2 )⎥
⎢ Re( u2 ) ⎥ ⎢ Re(u3) ⎥ Solving the radiofrequency induced thermal process involves three
⎢ Re(u3) ⎥ ⎢ ∂u1 ⎥ steps. (a) Given an electrical potential at the electrodes, solve for the
u=⎢ = h ,
Im(u1) ⎥ ⎢ ∂k ⎥ electrical potential distribution in the whole domain, Eq. (4). (b)
⎢ ⎥ ⎢ ∂u2 ⎥
⎢ Im(u2 )⎥ ⎢ ∂k h ⎥ Determine the gradient of the electrical potential to compute the
⎢⎣ Im(u3) ⎥⎦ ⎢ ∂u3 ⎥ Joule heating source Eq. (5) (c) Given the Joule heating source, solve
⎢⎣ ∂k h ⎥⎦ (19) for the temperature distribution, Eq. (6). These two steady state
where the real part corresponds to the normal output of the system, equations were solved in a sequential fashion using the complex
whereas, the imaginary part divided by the step size h, will give us the variable implementation of the FEM code.
sensitivity with respect to the variable that was perturbed by h in the The problem was solved using four different meshes in order to
imaginary direction. This result is general and can be applied to any examine the convergence of the solution in the sensitivities. The
type of finite element of two- or three-dimensional domains. In those coarsest mesh contained 496 elements, the next mesh contained
cases, the number of degrees of freedom will change accordingly. 1984 elements, the third mesh contained 7936 elements, and the
finest mesh consisted of 13,206 elements. The solution for the
temperature at the center of the tumor and the relativity sensitivity
4. Numerical example of the temperature with respect to λ1 at the center of the tumor was
computed with each mesh and the results are compared in Fig. 2. It can
FEM and ZFEM implementations of the Joule heating and the be observed the asymptotic convergence of both solutions. The
bioheat transfer equation were developed in the Matlab ® language temperature result of the 1984 node mesh differs by 0.3% with respect
[32]. A simple geometric domain was selected to illustrate the to the finest mesh and the sensitivity result differs by 0.99%. The
capabilities of the method. This example, with the same geometry results for the 7936 node mesh differ by 0.04% for the temperature,
and boundary conditions, was used by Majchrzak et al. [27] and is
presented in Fig. 1. Its dimensions are 0.08[m] in the x direction and
0.04[m] in the y direction. The electrodes are located at the upper and
lower boundary from x = 0.032 to x = 0.048[m]. The tumor is repre-
sented as a square domain with size of 0.015[m] and its left inferior
vertex has a coordinate of x = 0.032[m], y = 0.016[m]. The values of the
properties used for the healthy and the tumor tissue are presented in
the Table 1. The values of the properties for the human blood are listed
in the Table 2.
For the purpose of this analysis, these properties were considered
independent of the temperature. The boundary conditions for the

Table 1
Thermal and electrical properties of a healthy (Ω1) and a carcinoma (Ω2) human liver at a
frequency of 100[kHz]. (Data source from references [33,34] and [35]).

Tissue λ [W/m K] Gb [1/s] Qm [W/m3] σ [s/m]

Ω1 0.469 0.0005 4200 0.179

Ω2 0.600 0.002 42000 0.461 Fig. 2. Convergence analysis for the temperature and sensitivity with respect to λ
evaluated at the center of the tumor tissue.

14
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

Fig. 3. Electric potential and temperature field of the radiofrequency induced thermal process. (a) Electric potential field ϕ [V], (b) Temperature field distribution T [°C].

N
and by 0% for the sensitivity, with respect to the results of the finest tumor domain were calculated. That is, max(SkT (xj )) and ∑ j SkT (xj )/ N
mesh. for all xj ∈ Ω , with N the total number of nodes in Ω, and max(SkT (xj ))
Fig. 3 shows the electrical potential and the temperature distribu- N
and ∑ j 2 SkT (xj )/ N2 for all xj ∈ Ω2 , and N2 the total number of nodes in
tion for the example problem. It is observed that using a voltage of Ω2. The results of the maximum and averaged relative sensitivities of
−10 [V] and 10 [V] as boundary conditions, a maximum temperature of the temperature with respect to different variables (electrical conduc-
42.9°C was obtained and was located in the tumor region of the tivity, perfusion rate, metabolic heat and the thermal conductivity for
domain. This temperature value is just below the recommended the healthy and carcinoma tissue) were calculated and the results are
therapeutic value of 43°C . summarized in Table 3. The average sensitivity of the parameters on
At this point, in terms of the desired treatment, it is necessary to the tumor domain is plotted in Fig. 5. It can be observed that the
increase the temperature value at the tumor domain only. The electrical and the thermal conductivity of the healthy tissue are the two
procedure is to increase the difference in the electrical potential of properties with the greatest impact on the temperature, with an
the two electrodes. Furthermore, due to the delicate nature of the average relative sensitivity value of σ1 = 1.246 × 10−2 and
procedure, it is important to determine what is the effect of small λ1 = −1.242 × 10−2 respectively. Furthermore, the thermal conductivity
variations of the different parameters involved in the heating of the of the tumor, with SλT2 = −0.073 × 10−2 , and the electrical conductivity
tissue. Identifying the most critical factor and understanding the effect of the tumor, with SσT2 = 0.088 × 10−2 , seem not to have a significant
of the different variables in the radiofrequency induced thermal process influence on the temperature distribution even over the tumor region.
will lead to better treatments using this technique. Boundary conditions model the effect of the surroundings on the
object of interest. Their values are determined based on physical
4.2. Sensitivity analysis considerations or experimental measurements. Obviously, these data
may have measurement error or uncertainty. Therefore, it is important
The goal of the research is to determine the effect that the to characterize how sensitive are the results to the boundary condi-
parameters and the tissues properties have on the temperature tions. In order to compute these sensitivities, the following perturba-
distribution. This information can be used to obtain the best treatment tions were applied to the boundary conditions:
with a minimal collateral damage. Therefore, a sensitivity analysis for
each parameter and property was carried out using CTSE. Since the
sensitivity found by the CTSE method for each property is an absolute
• Electric field Dirichlet boundary condition: ϕ (x) = 10 + i h at the
top electrode, x ∈ Γ1.
sensitivity which can not be compared, a relative sensitivity is used. It
is defined as [37]
• ∂T (x)
Temperature field Newman boundary condition: ∂n = f = 0 + i h
at the adiabatic boundary, x ∈ Γ4 .

SkT (x) =
∂T (x)/ T (x)
=
∂ ln T (x)
=
k ∂T (x)
,
• Temperature field Dirichlet boundary condition: T (x) = 305.65 + i h
∂k / k ∂ ln k T (x) ∂k (20) at the cooling pad, x ∈ Γ1 + Γ3.

where k represents an arbitrary parameter of the model. Eq. (20) The first analysis will provide the temperature sensitivity with
provides a dimensionless value that can be used to compare the respect to the applied voltage, SϕTΓ . The second one will provide the
1
sensitivities obtained for all parameters. Fig. 4 presents the relative temperature sensitivity with respect to the heat flux through the
temperature sensitivities with respect to the thermal conductivity of the T
boundary S fΓ 4 . The third one will provide the temperature sensitivity
healthy (see Fig. 4a) and the tumor (see Fig. 4b) tissues; and with with respect to the assumed, or measured, temperature at the cooling
respect to the electrical conductivity of the healthy (see Fig. 4c) and the pad STTCP Table 4. The temperature sensitivities over the domain are
tumor (see Fig. 4d) tissues. It can be observed that these two properties shown in Fig. 6. Also, the maximum and average sensitivity over the
have opposite effects on the temperature. Increasing the thermal total domain Ω and the tumor domain Ω2 were calculated and are
conductivity will decrease the temperature and increasing the electric shown in Table 5. It can be observed that the temperature sensitivity
conductivity will increase the temperature. Also, the temperature is less with respect to the temperature boundary condition is 100% near to the
sensitive to the thermal conductivity of the carcinoma than to the cooling pad and decreases away from the boundary. This is an expected
thermal conductivity of the healthy tissue, and its effect seems to be result as this boundary condition sets the temperature and therefore
limited within the tumor region. The thermal conductivity of the has a one-to-one effect on the temperature of the healthy tissue near to
healthy tissue acts as a diffuser, spreading the temperature across the the cooling pad. Over the tumor region, the average sensitivity is
domain. On the other hand, an increase in the electrical conductivity of 36.87% and the maximum is 53.56%. These are large values when
the healthy tissue would lead to an increase in the temperature field on compared to the sensitivities with respect to the other parameters. The
the tumor region as the Fig. 4c suggests. However, it will also increase value of this boundary condition comes from experimental measure-
the temperature in the healthy tissue and its maximum value will be in ment. Therefore, very accurate measurements have to be made in order
the healthy region. to obtain accurate temperature values in the healthy and tumor tissues.
The relative sensitivity is computed at each point in the domain The second most influential boundary condition is the electric potential
SkT (xj ) with xj ∈ Ω representing a finite element node. In order to boundary condition. This value can be controlled in the radiofrequency
compare the sensitivities with respect to the different parameters its induced thermal procedure. An increase in the electric potential of the
maximum value and average value over the total domain and over the

15
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

Fig. 4. Relative sensitivity of temperature with respect to: (a) thermal conductivity of the healthy tissue SλT1 , (b) thermal conductivity of the tumor SλT2 , (c) electrical conductivity of the
healthy tissue SσT1 , (b) electrical conductivity of the tumor SσT2 .

Table 3
Average sensitivity of the temperature with respect to different properties. The values are multiplied by 100.

Thermal conduct. Perfusion rate Metabolic heat Electric conduct.

Properties λ1 λ2 Gb1 Gb2 Qm1 Qm2 σ1 σ2

Maximum sensitivity over Ω −1.77 −0.203 −0.275 −0.807 0.320 0.785 1.500 0.432
Average sensitivity over Ω −0.78 −0.002 −0.056 −0.180 0.170 0.18 0.624 0.050
Maximum sensitivity over Ω2 −1.51 −0.203 −0.240 −0.807 0.200 0.785 1.370 0.260
Average sensitivity over Ω2 −1.242 −0.073 −0.153 −0.646 0.138 0.642 1.246 0.088

Table 4
Average and maximum sensitivity of the temperature with respect to the boundary
conditions. Values are multiplied by 100.

Electric Adiabatic Cooling Temp.

potential boundary

Maximum sensitivity 1.920 3.040 100.00

over Ω
Average sensitivity over 0.670 0.490 61.74
Ω
Maximum sensitivity 1.630 0.110 53.56
over Ω2
Average sensitivity over 1.330 0.075 36.87
Ω2

expansion method for the radiofrequency induced thermal procedure,

the sensitivity of temperature with respect to two parameters (the
Fig. 5. Average tumor sensitivity of temperature with respect to the properties. The thermal conductivity and the electrical conductivity of the healthy
values are multiplied by 100. tissue) were calculated using the Finite Difference Method (FDM). The
comparison was made for different step sizes (h), from h = 1 to
electrodes will increase the temperature in the whole domain (0.67% h = 1e−20 .
on average) and its effect is larger inside the tumor (1.33% on average).
│f ′ − f ′ref │
Finally, the sensitivity with respect to the ideal adiabatic condition is ε= ,
high in the region near the boundary and where the temperature is low. │f ′ref │ (21)
However, at the tumor region, it has an average of 0.0755%. This value
where f = maxx (SλT1 (x )),
The reference derivative was computed using
shows that the adiabatic condition is well imposed.
the CTSE method with a step size of h = 10−50 . Fig. 7 and 8 plot the
results of this comparison. It can be observed that in both cases, the
4.3. A comparison of CTSE against the FDM sensitivity computed using the finite difference method fails at h values
smaller than 10-12. This is due to round-off errors involved in the
In order to verify the accuracy of the complex Taylor series calculation of the sensitivity caused by the subtraction of terms. On the

16
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

Fig. 6. Relative sensitivity of temperature with respect to the boundary conditions: (a) electrical Dirichlet boundary condition SϕTΓ , (b) ideal adiabatic boundary condition S TfΓ 4 , (c)
1
temperature at the cooling pad STTCP .

Table 5
Average and maximum sensitivity of the temperature with respect to the boundary
conditions. Values are multiplied by 100.

Electric Adiabatic Cooling

potential boundary Temperature.

Maximum sensitivity over Ω 1.920 3.040 100.00

Average sensitivity over Ω 0.670 0.490 61.74
Maximum sensitivity over Ω2 1.630 0.110 53.56
Average sensitivity over Ω2 1.330 0.075 36.87

Fig. 8. Normalized error in the electrical conductivity sensitivity estimates ∂T /∂σ1, given
by finite-difference and the complex Taylor series expansion method for different step
sizes.

the solution time. The current implementation uses sparse matrices

where only the elements different from zero are stored. The matrix size
in a sparse matrix is the number of entries different from zero.
Therefore, the memory requirements for the complex-value code will
depend mainly on how the complex perturbation affects the matrix. For
example, the sensitivity with respect to a boundary condition will only
affect the system vector so the size of the sparse matrix is only
Fig. 7. Normalized error in the thermal conductivity sensitivity estimates ∂T /∂λ1, given increased by a factor of two (imaginary part equal zero in Eq. (18).
by finite-difference and the complex Taylor series expansion method for different step On the other hand, if the sensitivity involves a perturbation in the full
sizes.
stiffness matrix then the size of the sparse matrix is increased by a
factor of four. In order to assess the computational cost of a complex
other hand, the sensitivity computed using the CTSE method remains perturbation, these two extreme cases using eight different mesh
constant throughout the entire range of the step size and oscillations densities were evaluated. Table 6 presents the results for the two types
are not present. As opposed to the finite difference method, the CTSE of sensitivities. a) Sensitivity with respect to a boundary condition (e.g.
method still converges when very low step size values are used. temperature boundary condition) and b) Sensitivities with respect to
Therefore, the sensitivity computed with the CTSE method is step size the material properties (e.g thermal conductivity). For each case, the
independent. CPU time and the memory used by the program were computed. It can
be observed case a) that the memory required was doubled and the
4.4. Computational efficiency total computing time was increased by a factor of 1.61 on average when
compared to a single FEM evaluation. In case b) the memory
When using a complex-valued code the number of degrees of requirements were multiplied by a factor of 4. However, the solution
freedom are duplicated, causing an increase in the memory used and

17
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

Table 6
Computing time and memory requirements for complex-coded FEM for: a) Boundary condition perturbation, b) Material property perturbation. Real-coded FEM requirements are
presented for reference but they do not produce sensitivity information alone.

a) Dirichlet boundary condition complex perturbation

Real Complex

Num. Elem. DOF Num. Pos. sparse Total Time (s) DOF Num. Pos. sparse Total Time (s) Time ratio

24 20 106 0.44 40 212 0.72 1.62

96 63 379 1.75 126 758 2.83 1.62
384 221 1,429 6.95 442 2,858 11.28 1.62
1,536 825 5,545 27.74 1,650 11,090 44.85 1.62
6,144 3,185 21,841 110.72 6,370 43,682 178.27 1.61
24,576 12,513 86,689 441.17 25,026 173,378 713.76 1.62
98,304 49,601 345,409 1,778.88 99,202 690,818 2,846.58 1.60
1,572,864 788,225 5,510,401 28,420.28 1,576,450 11,020,802 45,611.86 1.60

b) System matrix complex perturbation

Real Complex

Num. Elem. DOF Num. Pos. sparse Total Time (s) DOF Num. Pos. sparse Total Time (s) Time ratio

24 20 106 0.46 40 424 0.75 1.62

96 63 379 1.82 126 1,516 2.94 1.62
384 221 1,429 7.21 442 5,716 11.68 1.62
1,536 825 5,545 28.76 1,650 22,180 46.52 1.62
6,144 3,185 21,841 114.32 6,370 87,364 185.05 1.62
24,576 12,513 86,689 457.20 25,026 346,756 741.62 1.62
98,304 49,601 345,409 1,778.88 99,202 1,381,636 2,978.00 1.67
1,572,864 788,225 5,510,401 28,420.28 1,576,450 22,041,604 48,147.97 1.69

time increased only by a factor of 1.64 on average. These results determine the influence of the parameters in the resulting tumor/tissue
represent an improvement over the run-time increase reported in the temperature. This example is based on the experimental and numerical
literature for full matrix storage, e.g. Martins et al. [38] and Whitfield study developed by Fuentes et al. [1]. The RFA procedure used
et al. [39] reported increases in run-time by a factor of 3. In addition, consisted of a pulsed voltage emanating through a tri-cooled tip
there are other potential approaches to improve the efficiency of the electrode. The electrode was run in impedance mode for 12 minutes
complex method as the so-called complex variable semi-analytical per ablation and cooled with chilled water. The liver domain was
method [40]. subdivided into two domains the healthy tissue domain Ω1 and the
tumor domain Ω2. The boundary was subdivided into the external
4.4.1. Finite difference comparison boundary Γ2 and the electrodes Γ1 inside the tumor tissue. A 2D section
When forward or backward finite difference is used to compute a of the liver was used for the analysis.
sensitivity one needs to solve the problem at two different points. The
solution time is a factor of two and the truncation error is  (h ). When 5.1. Liver properties
the central finite difference is used to improve the truncation error to
 (h2 ), the problem needs to be solved at three different points. The The liver properties were selected as follows. The electric conduc-
increase in the solution time is a factor of three. In contrast, CTSE has a tivity of the healthy and tumor tissue at a frequency of 460 kHz were
truncation error of  (h2 ) and the overhead time is only 1.6 times. obtained from [35], the heat transfer conductivities from [1] and the
Perfusion rates were from [41].
5. Liver ablation treatment sensitivities With regard of the metabolic heat transfer generation, the authors
found no agreement in the literature. For example, a value of
Fig. 9 shows the geometry of a human liver that was used to Qm = 1091 W/m3 for the healthy tissue is used by several authors like
illustrate the capabilities of the method. A sensitivity analysis of the [42,43]. This value is based on the paper by Mitchell [44]. Elia [45]
radiofrequency ablation (RFA) process was undertaken in order to identified the specific resting metabolic rates of major organs and

Fig. 9. 3D Liver geometry and a 2D section showing the resulting temperature of an radiofrequency ablation process. No cooling at the electrodes was considered in this case.

18
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

tissues in young adults with normal weight. A value of 9.7 W/m3 was Table 8
given for the liver. Fiala et al. [46] presented a computer model of Sensitivity of the temperature with respect to the metabolic heat for different range of
values. Qm (Ω1) is the metabolic heat for the healthy tissue and Qm (Ω2 ) for the tumor
human thermo-regulation for a wide range of environmental condi-
tissue.
tions. A metabolic rate of 4100 W/m3 was given to the human viscera.
Kummitha et al. [47] computed the Energy Expenditure (EE), for Qm (Ω1) − Qm (Ω2 ) [W/m3]
different organs/tissues, using information from mouse and a mathe-
matical model. They proposed a metabolic heat rate for the liver of 0−0 10–100 420–4200 450–29,000
10.8 W/m3
maximum temp at tumor [°C] 54.65 54.65 54.70 55.04
No information was found relative to the liver tumor metabolic maximum temp at healthy [°C] 45.26 45.26 45.30 45.50
heat. However, Gautherie [48] used an electric analogy model of heat maximum sensitivity at tumor [%] 0.00 0.00 0.02 0.13
transfer in cancerous breast to determine the metabolic heat genera- maximum sensitivity at healthy [%] 0.00 0.00 0.01 0.08
tion. He proposed a value of 450 W/m3 for normal tissue and
29000 W/m3 for cancer tissue. Also, Liu and Xu [49] recognized the
lack of experimental measurements of this parameter and assumed a sensitivity of the temperature to the metabolic heat generation. The
value of 420 W/m3 for healthy and 4200 W/m3 for tumor tissues under- voltage was 15 V and no cooling was considered. Table 8 shows the
neath the skin. In [50], these values are multiplied by a factor of ten (no results obtained for the different ranges of metabolic heat. It can be
reason given). In spite of the uncertainty of these values, a number of observed that the maximum temperature changes only for very high
articles use this source for the metabolic heat generation. [51,52] and values of Qm and the change is about 1.3%. However, these high values
[27] are examples of that. have not been reported in liver cancer. It is also observed that local
Finally, some authors line Chang [6] assumed that the metabolic sensitivities obtained by CTSE are very small and therefore metabolic
heat source was an insignificant influence in the ablation simulation heat generation has a minor role in the ablation process.
and was not taken into account when solving the bioheat equation. Fig. 9 shows the temperature distribution for a voltage of 15 V and
Summarizing, the metabolic heat generation possible values are in a no cooling at the electrodes. In this case the maximum temperature at
range from zero to 29000 W/m3. Therefore, sensitivities of this para- the tumor was 54.65 °C and the maximum temperature at the healthy
meter along the whole range need to be computed. tissue was 45.26 °C. It can be observed that the higher temperatures
A summary of the properties of the healthy and tumor tissue are were concentrated in the tumor area. In contrast, Fig. 10 shows the
shown in Table 7 and the properties of the blood are shown in Table 2. temperature distribution for the same applied voltage but with cooling
at the electrode tips. It can be observed that there is a region in the
tumor in between the electrodes where the temperature does not reach
5.2. Boundary conditions
values higher than 35 °C and therefore that part of the tumor will not be
effectively treated. Also, the sensitivity map shows that the highest
The electric potential as varied in a range in order to obtain the
sensitive areas are precisely around the electrodes.
ablation temperatures. An electric potential of 15V, 20V, 30V was used
Based on the previous results, only the sensitivities with respect to
at the electrodes Γ1, and a zero potential at the external boundaries of
the cooling temperature STT∞ , convective cooling coefficient STh, electric
the liver Γ2.
potential boundary condition STV, and electric conductivity of the
The cooling at the electrodes was modelled as a thin cylinder with healthy tissue SσT1 were calculated. Their maximum and averaged values
water inside. As the diameter is very small, 1 mm, the fluid can be were calculated over the complete domain and over the healthy and
consider laminar and the convective coefficient approximated as tumor domains. The results are summarized in Table 9. It can be
h ≈ 2000 W/m2. K [53]. The cooling temperature was set to 30 °C and observed that the parameter with the largest influence in the tempera-
the temperature at the external boundary Γ2 was set to 37 °C . ture is the cooling temperature as it directly prescribe the values of the
temperature of the tissue in the vicinity of the electrodes. The second
5.3. Mesh convergence largest influence in the temperature is the electric potential applied. A
variation of the potential will affect the temperature over the whole
The problem was solved using three different meshes in order to domain. Another important parameter is the electric conductivity of
examine the convergence of the solution. The coarsest mesh contained the healthy tissue σ1. The relative sensitivity is about 5%. As this
3009 elements, the second mesh contained 5505 elements, and the variable is not under control of the process, special care must be taken
finest mesh consisted of 10,510 elements. Due to the unstructured when designing an adequate radiofrequency induced thermal process
nature of the mesh, the maximum norm of the temperature sensitivity as the main interest of the process is to raise the temperature of the
with respect to λ was used to test the solution convergence. The results tumor only. To accomplish this it will be desirable to increase the
were 0.4145 using the mesh with 3009 elements; 0.4139 with 5505 electrical conductivity value of the tumor. Some research in this
elements; and 0.4137 with 10,510 elements. It is observed that the last direction includes the introduction of nanoparticles [51,2,1,54,55]
two results differ only by 0.05%, which represents the relative error in which are capable of modifying the general properties into the tumor
the solution. Based on these results, the mesh with 5505 elements was region and concentrating the energy.
selected for the present study.

6. Conclusions
5.4. Results
Sensitivity analysis of the radiofrequency induced thermal process
The first set of computational experiments was designed to test the
was accomplished using the complex Taylor series expansion finite
element method. A finite element code was written to support complex
Table 7
Physical properties of the liver healthy tissue and liver tumor tissue. Metabolic heat
variables using the Cauchy-Riemann matrix representation.
transfer is evaluated in a range. Data source is presented in Section 5.1. Derivatives with respect to a parameter were obtained by perturbing
the variable of interest by a small amount h along the imaginary axis
Tissue λ [W/m°K] Gb [1/s] Qm [W/m3] σ [s/m] and then computing the imaginary part of the result. This constitutes a
major advantage when designing a medical treatment where the
Ω1 0.5 0.018 0 to 4200 0.26
Ω2 0.6 0.022 0 to 29000 0.504 sensitivities with respect to all the variables have to be considered in
order to design a robust and safe procedure.

19
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

Fig. 10. Temperature distribution and temperature sensitivity for a 15 V RFA procedure and colling temperature of 30 °C at the electrodes.

Table 9 Appendix A. Supplementary data

Maximum and average sensitivities computed over the liver Ω1, tumor Ω1, and healthy
tissue Ω2. The values are multiplied by 100.
Supplementary data associated with this article can be found in the
STh STV online version at http://dx.doi.org/10.1016/j.finel.2017.07.001.
STT∞ SσT1

Average sensitivity over Ω 16.29 −0.18 3.39 1.38 References

Average sensitivity over Ω2 39.45 −0.45 7.36 3.21
Average sensitivity over Ω1 1.21 −0.01 0.88 0.20 [1] D. Fuentes, R. Cardan, R.J. Stafford, J. Yung, G.D. Dodd, Y. Feng, High-fidelity
Maximum Sensitivity over Ω 90.48 0.00 10.75 4.85 computer models for prospective treatment planning of radiofrequency ablation
Maximum sensitivity over Ω2 90.48 −0.04 10.75 4.85 with in vitro experimental correlation, J. Vasc. Interv. Radiol. 21 (11) (2010)
Maximum sensitivity over Ω1 17.05 0.00 5.19 2.26 1725–1732. http://dx.doi.org/10.1016/j.jvir.2010.07.022.
[2] Y. Feng, D. Fuentes, Model-based planning and real-time predictive control for
laser-induced thermal therapy, Int. J. Hyperth. 27 (8) (2011) 751–761.
[3] M.N. Rylander, Y. Feng, K. Zimmermann, K.R. Diller, Measurement and mathe-
The electric field and bioheat equations were solved using the CTSE matical modeling of thermally induced injury and heat shock protein expression
finite element code to automatically obtain sensitivity results. Relative kinetics in normal and cancerous prostate cells, Int. J. Hyperth. 26 (8) (2010)
sensitivities were computed and compared to establish the most 748–764.
[4] H.H. Pennes, Analysis of tissue and arterial blood temperatures in the resting
influential properties in the radiofrequency induced thermal process.
human forearm, J. Appl. Physiol. 1 (2) (1948) 93–122.
The parameters included the physical properties thermal, electrical [5] M. Tunç, U. Çamdali, S. Cikrikci, The bio-heat transfer equation and its applica-
conductivity, perfusion rate, metabolic heat, and heat capacity. The tions in hyperthermia treatments, Eng. Comput. 23 (4) (2006) 451–463.
sensitives of the boundary conditions and cooling temperature were [6] I. Chang, U.D. Nguyen, Thermal modeling of lesion growth with radiofrequency
ablation devices, Biomed. Eng. Online 3 (1) (2004) 27. http://dx.doi.org/10.1186/
also computed. 1475-925X-3-27.
It was found the thermal and electrical conductivity of the healthy [7] M. Zhang, Z. Zhou, S. Wu, L. Lin, H. Gao, Y. Feng, Simulation of temperature field
tissue were the most influential properties in the final temperature. for temperature-controlled radio frequency ablation using a hyperbolic bioheat
equation and temperature-varied voltage calibration: a liver-mimicking phantom
This means that if the thermal energy is not carefully applied to the study, Phys. Med. Biol. 60 (24) (2015) 9455.
right cancerous region healthy tissues can be damaged. It also supports [8] Z. Wang, I. Aarya, M. Gueorguieva, D. Liu, H. Luo, L. Manfredi, L. Wang,
the idea of locally and artificially changing these properties at the D. McLean, S. Coleman, S. Brown, A. Cuschieri, Image-based 3d modeling and
validation of radiofrequency interstitial tumor ablation using a tissue-mimicking
location of the tumor. breast phantom, Int. J. Comput. Assist. Radiol. Surg. 7 (6) (2012) 941–948. http://
In the case of th RFA process, it was found that the highest relative dx.doi.org/10.1007/s11548-012-0769-3 (URL 〈http://dx.doi.org/10.1007/
sensitivity was the cooling temperature with values as high as 90%. The s11548-012-0769-3〉).
[9] C. Schumann, C. Rieder, S. Haase, K. Teichert, P. S¨uss, P. Isfort, P. Bruners,
highest sensitivities were located in the area between the electrodes. As
T. Preusser, Interactive multi-criteria planning for radiofrequency ablation, Int. J.
a consequence, when the low cooling temperatures are maintained, it Comput. Assist. Radiol. Surg. 10 (6) (2015) 879–889. http://dx.doi.org/10.1007/
will result in tumor areas not exposed to ablation temperature. The s11548-015-1201-6 (URL 〈http://dx.doi.org/10.1007/s11548-015-1201-6〉).
[10] J.N. Lyness, C.B. Moler, Numerical differentiation of analytic functions, SIAM J.
second largest temperature sensitivity was found to be with respect to
Numer. Anal. 4 (2) (1967) 202–210. http://dx.doi.org/10.1137/0704019.
the electric potential applied and a variation of the potential will affect [11] J.N. Lyness, Numerical algorithms based on the theory of complex variable, in:
the temperature over the whole domain. The metabolic heat transfers 1967 Proceedings of the 22nd National Conference, ACM '67, ACM, New York, NY,
values for the liver and tumor are not well known. A sensitivity analysis USA, 1967, pp. 125–133. http://dx.doi.org/10.1145/800196.805983.
[12] W. Squire, G. Trapp, Using complex variables to estimate derivatives of real
in a different range of possible values revealed that there is no functions, SIAM Rev. 40 (1) (1998) 110–112. http://dx.doi.org/10.1137/
significant influence of the metabolic heat in the final ablation S003614459631241X.
temperature. [13] V. Vatsa, Computation of sensitivity derivatives of navier-stokes equations using
complex variables, Adv. Eng. Softw. 31 (8–9) (2000) 655–659. http://dx.doi.org/
The sensitivity of the temperature with respect to the thermal and 10.1016/S0965-9978(00)00025-9.
electrical conductivity was compared with the finite difference method. [14] A. Voorhees, H. Millwater, R. Bagley, Complex variable methods for shape
The results showed that as opposed to the finite difference method, sensitivity of finite element models, Finite Elem. Anal. Des. 47 (10) (2011)
1146–1156.
CTSE is step size independent. In contrast to the FDM, the CTSE only [15] A. Voorhees, H. Millwater, R. Bagley, P. Golden, Fatigue sensitivity analysis using
needs one run for each sensitivity of interest. The increase in run time complex variable methods, Int. J. Fatigue 40 (2012) 61–73.
for a complex FEA as compared to a real FEA was only 1.6; significantly [16] D. Wagner, H. Millwater, 2D weight function development using a complex taylor
series expansion method, Eng. Fract. Mech. 86 (2012) 23–37. http://dx.doi.org/
lower than the 2x for forward/backward finite differencing or 3x for
10.1016/j.engfracmech.2012.02.006 (URL 〈http://www.sciencedirect.com/
central differencing. Furthermore, the current CTSE implementation is science/article/pii/S0013794412000689〉).
efficient in memory and uses less computational time that finite [17] H. Millwater, D. Wagner, A. Baines, K. Lovelady, Improved WCTSE method for the
generation of 2D weight functions through implementation into a commercial finite
difference method. CTSE has significant advantages in the sensitivity
element code, Eng. Fract. Mech. 109 (2013) 302–309. http://dx.doi.org/10.1016/
analysis computation. Once the finite element code is written with j.engfracmech.2013.07.012 (URL 〈http://www.sciencedirect.com/science/article/
complex variables, it is easy to compute the sensitive of the variables of pii/S0013794413002579〉).
interest in a precise and robust manner. For this reason, the method is [18] Z. Jing, X. Wu, Wide-range weight functions and stress intensity factors for
arbitrarily shaped crack geometries using complex taylor series expansion method,
an excellent alternative to traditional sensitivity methods. Eng. Fract. Mech. 138 (2015) 215–232. http://dx.doi.org/10.1016/j.engfrac-
mech.2015.03.006 (URL 〈http://www.sciencedirect.com/science/article/pii/

20
J.F. Monsalvo et al. Finite Elements in Analysis and Design 135 (2017) 11–21

S0013794415000788〉). [37] T.F. Edgar, D.M. Himmelblau, Optimization of Chemical Processes, 2nd ed.,
[19] X.C. Zhao, X.R. Wu, D.H. Tong, Weight functions and stress intensity factors for McGraw-Hill, New York, 2001.
pin-loaded single-edge notch bend specimen, Fatigue Fract. Eng. Mater. Struct. 38 [38] J.R.R.A. Martins, P. Sturdza, J.J. Alonso, The complex-step derivative approx-
(12) (2015) 1519–1528. http://dx.doi.org/10.1111/ffe.12343 (URL 〈http://dx.doi. imation, ACM Trans. Math. Softw. 29 (3) (2003) 245–262.
org/10.1111/ffe.12343〉). [39] D. Whitfield, J. Newman, W. Anderson, Step-size independent approach for
[20] A. Montoya, R. Fielder, A. Gomez-Farias, H. Millwater, Finite-element sensitivity multidisciplinary sensitivity analysis, Journal of aircraft 40 (3).
for plasticity using complex variable methods, J. Eng. Mech. 141 (2) (2015) [40] W. Jin, B.H. Dennis, B.P. Wang, Improved sensitivity analysis using a complex
04014118. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000837 (URL variable semi-analytical method, Struct. Multidiscip. Optim. 41 (3) (2010)
〈http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000837〉. 433–439.
[21] A. Gomez-Farias, A. Montoya, H. Millwater, Complex finite element sensitivity [41] B.E. Van Beers, I. Leconte, R. Materne, A.M. Smith, J. Jamart, Y. Horsmans,
method for creep analysis, Int. J. Press. Vessel. Pip. 132–133 (2015) 27–42. http:// Hepatic perfusion parameters in chronic liver disease: dynamic ct measurements
dx.doi.org/10.1016/j.ijpvp.2015.05.006 (URL 〈http://www.sciencedirect.com/ correlated with disease severity, Am. J. Roentgenol. 176 (3) (2001) 667–673.
science/article/pii/S0308016115000587〉). [42] K.N. Rai, S.K. Rai, Effect of metabolic heat generation and blood perfusion on the
[22] W.K. Anderson, J. Newman, D.L. Whitfield, Sensitivity analysis for navier-stokes heat transfer in the tissues with a blood vessel, Heat. Mass Transf. 35 (1) (1999)
equations on unstructured meshes using complex variables, AIAA J. 39 (1) (2001) 75–79. http://dx.doi.org/10.1007/s002310050300 (URL 〈http://dx.doi.org/10.
56–63. 1007/s002310050300〉).
[23] J.C. Newman, E.J. Whitfield, W.K. Anderson, Step-size independent approach for [43] P.K. Gupta, J. Singh, K. Rai, Numerical simulation for heat transfer in tissues
multidisciplinary sensitivity analysis, J. Aircr. 40 (3) (2003) 566–573 (URL during thermal therapy, J. Therm. Biol. 35 (6) (2010) 295–301. http://dx.doi.org/
〈http://arc.aiaa.org/doi/pdf/10.2514/2.3131〉). 10.1016/j.jtherbio.2010.06.007 (URL 〈http://www.sciencedirect.com/science/
[24] B. Wang, A. Apte, Complex variable method for eigensolution sensitivity analysis, article/pii/S030645651000063X〉).
AIAA J. 44 (12) (2006) 2958–2961. [44] J.W. Mitchell, T.L. Galvez, J. Hengle, G.E. Myers, K.L. Siebecker, Thermal response
[25] J. Garza, H. Millwater, Multicomplex newmark-beta time integration method for of human legs during cooling, J. Appl. Physiol. 29 (6) (1970) 859–865.
sensitivity analysis in structural dynamics, AIAA J. 53 (5) (2015) 1188–1198. [45] M. Elia, Organ and tissue contribution to metabolic rate, Energy metabolism: tissue
http://dx.doi.org/10.2514/1.J053282 (URL 〈http://dx.doi.org/10.2514/1. determinants and cellular corollaries, 1992 (1992) 19–60.
J053282〉). [46] D. Fiala, K.J. Lomas, M. Stohrer, A computer model of human thermoregulation for
[26] S. Kim, J. Ryu, M. Cho, Numerically generated tangent stiffness matrices using the a wide range of environmental conditions: the passive system, J. Appl. Physiol. 87
complex variable derivative method for nonlinear structural analysis, Comput. (5) (1999) 1957–1972 (URL 〈http://jap.physiology.org/content/87/5/1957〉).
Methods Appl. Mech. Eng. 200 (1–4) (2011) 403–413. http://dx.doi.org/10.1016/ [47] C.M. Kummitha, S.C. Kalhan, G.M. Saidel, N. Lai, Relating tissue/organ energy
j.cma.2010.09.004 (URL 〈http://www.sciencedirect.com/science/article/pii/ expenditure to metabolic fluxes in mouse and human: experimental data integrated
S0045782510002549〉). with mathematical modeling, Physiol. Rep. 2 (9) (2014) e12159. http://dx.doi.org/
[27] E. Majchrzak, G. Dziatkiewicz, M. Paruch, The modelling of heating a tissue 10.14814/phy2.12159 (URL 〈http://www.ncbi.nlm.nih.gov/pmc/articles/
subjected to external electromagnetic field, Acta Bioeng. Biomech. 10 (2) (2008) PMC4270223/〉).
29–37. [48] M. Gautherie, Thermopathology of breast cancer: measurement and analysis of in
[28] D.J. Griffiths, Introduction to Electrodynamics, 4th ed., Pearson, Boston, 2013. vivo temperature and blood flow, Ann. N. Y. Acad. Sci. 335 (1) (1980) 383–415.
[29] K.K. Choi, N.-H. Kim, Structural Sensitivity Analysis and Optimization 1: Linear [49] J. Liu, L.X. Xu, Boundary information based diagnostics on the thermal states of
Systems, Springer, 2005. http://dx.doi.org/10.1007/b138709. biological bodies, Int. J. Heat. Mass Transf. 43 (16) (2000) 2827–2839.
[30] L.B. Rall, Automatic Differentiation: Techniques and Applications, Springer, Berlin, [50] Z.-S. Deng, J. Liu, Monte carlo method to solve multidimensional bioheat transfer
1981. problem, Numer. Heat. Transf.: Part B: Fundam. 42 (6) (2002) 543–567.
[31] G.B. Price, An Introduction to Multicomplex Spaces and Functions, CRC Press, [51] Y.-G. Lv, Z.-S. Deng, J. Liu, 3-d numerical study on the induced heating effects of
New York, 1990. embedded micro/nanoparticles on human body subject to external medical
[32] The MathWorks Inc., MATLAB and Statistics Toolbox Release 2015a, Natick, electromagnetic field, IEEE Trans. Nanobiosci. 4 (4) (2005) 284–294.
Massachusetts, United States, 2015. [52] Z.-S. Deng, J. Liu, Mathematical modeling of temperature mapping over skin
[33] T.E. Cooper, G.J. Trezek, A probe technique for determining the thermal surface and its implementation in thermal disease diagnostics, Comput. Biol. Med.
conductivity of tissue, J. Heat. Transf. 94 (2) (1972) 133–140. 34 (6) (2004) 495–521.
[34] J. Valvano, J. Cochran, K. Diller, Thermal conductivity and diffusivity of bioma- [53] T.L. Bergman, F.P. Incropera, Introduction to Heat Transfer, John Wiley & Sons,
terials measured with self-heated thermistors, Int. J. Thermophys. 6 (3) (1985) New Jersey, 2011.
301–311. http://dx.doi.org/10.1007/BF00522151. [54] Y. Feng, D. Fuentes, A. Hawkins, J. Bass, M.N. Rylander, A. Elliott, A. Shetty,
[35] D. Haemmerich, D.J. Schutt, A.W. Wright, J.G. Webster, D.M. Mahvi, Electrical R.J. Stafford, J.T. Oden, Nanoshell-mediated laser surgery simulation for prostate
conductivity measurement of excised human metastatic liver tumours before and cancer treatment, Eng. Comput. 25 (1) (2009) 3–13.
after thermal ablation, Physiol. Meas. 30 (5) (2009) 459–466. http://dx.doi.org/ [55] Y. Feng, M. Rylander, J. Bass, J. Oden, K. Diller, Optimal design of laser surgery for
10.1088/0967-3334/30/5/003. cancer treatment through nanoparticle-mediated hyperthermia therapy, in: NSTI-
[36] S. Tungjitkusolmun, Three-dimensional finite-element analyses for radio-frequency Nanotech, vol. 1, 2005, pp. 39–42.
hepatic tumor ablation, IEEE Trans. Bio-Med. Eng. 49 (1) (2002) 3–9.

21