Oho Lik Chan Aerospace and Mechanical Engineering Department, The University of Arizona, Tucson, AZ 85721
Boundary Element Method Analysis for the Bioheat Transfer Equation In this paper, the boundary element method (BEM) approach is applied to solve the Pennes (1948) bioheat equation. The objective is to develop the BEM formulation and demonstrate its feasibility. The basic BEMformulations for the transient and steady-state cases are first presented. To demonstrate the usefulness of the BEM approach, numerical solutions for 2-D steady-state problems are obtained and compared to analytical solutions. Further, the BEM formulation is applied to model a conjugate problem for an artery imbedded in a perfused heated tissue. Analytical solution is possible when the conduction in the x-direction is negligible. The BEM and analytical results have very good agreement.
Introduction The complexity of the mathematical modeling of the heat the BEM results are more accurate than the FEM. This implies transfer process within a human tissue is enormous. The dif- that for the same desired level of accuracy the BEM requires ficulties arise from the fact that the substance is extremely a coarser mesh than the FEM does. They also found that the heterogeneous (Chen, 1985). The basic heat transfer processes BEM accuracy, with a fixed boundary mesh, does not change in biological systems can be found in a recent review by Chato appreciably with the number of internal sampling points. On (1990). The most commonly used model is the bioheat transfer the other hand, the FEM accuracy is very sensitive to the equation developed by Pennes (1948). He postulated that the internal mesh. effect of the temperature difference between the blood supply In the present work, the BEM formulation for the bioheat and the tissue acts as an energy sink term, Qb = wbcb(T- Ta). transfer equation is first presented. To validate the BEM apThis in essence is very similar to the fin equation where the proach, numerical solutions of a 2-D steady-state problem are sink term represents convective heat loss to the surrounding. presented and compared to an analytical solution. Finally, to The bioheat transfer equation has been applied by many illustrate the capability of the BEM approach, a solution of investigators to hyperthermia. It has been found that the bio- the heat transfer within tissue adjacent to an artery is presented. heat transfer equation provides useful and adequate predic- Analytical solution is possible when the conduction in the xtions (Strohbehn and Roemer, 1984; Roemer, 1990) in many direction is negligible. The BEM and analytical results are cases. A review of numerical solutions of the bioheat transfer presented and compared. equation can be found in (Strohbehn and Roemer, 1984). It suffices to mention that the numerical methods used are mainly Mathematical Formulation the finite difference method and finite element method. The bioheat transfer equation proposed by Pennes (1948) The Boundary Element Method (BEM) is another powerful can be written as, general purpose method (Banerjee and Butterfield, 1981; Mukherjee, 1982; Brebbia et al., 1984). It is far more tolerant of 97"* klV2T*-wbcb(T*-T*) + S*=p,ct — inQ* (1.1) aspect ratio degradation than the Finite Element Method (FEM) and can yield secondary variables as accurately as the primary ones. The internal equations are applied pointwise, hence, Introducing the dimensionless variables, sharp temperature gradients over the domain may be captured (x*,y\z*) „ T'-T* t t*a, S=(1.2) easily. The discretization of the boundary can be tailored to (x,y,z)-, t-- Ll ', 'T=-S L2/k, ' ~ L " S0 0 any irregular shape, thus it is capable of handling complex' geometries. Mukherjee and Morjaria (1982) compared the ef- the governing equation can be expressed in dimensionless form ficiency and accuracy of the BEM and FEM for the Laplace as, equation. They found that for the same level of discretization 3T. V2T-P/T+S = :^ in Q (1.3) dt where Contributed by the Bioengineering Division for publication in the JOURNAL OF BIOMECHANICAI ENGINEERING. Manuscript received by the Bioengineering Division March 11, 1991; revised manuscript received October 26, 1991. Associate Technical Editor: K. R. Diller.
358 / Vol, 114, AUGUST 1992
Pr=
WbCbL2
(1.4)
Transactions of the ASME
Copyright © 1992 by ASME Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/18/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
In addition, appropriate boundary conditions must be specified along the boundary for a well posed problem. These will be presented with specific examples in later sections. Boundary Element Formulation In the Boundary Element Method, the essential steps are (i) derivation of the equivalent boundary integral equation, (ii) determination of the Green's function, and (iii) discretization of the boundary integral. The BEM can be applied most efficiently if an analytical expression for the Green's function of the adjoint equation is available. In this section, we derive the boundary integral equations and the relevant Green's functions for both unsteady and steady problems. We will illustrate the discretization of the boundary integral using the two-dimensional steady-state equation. Transient Formulation. The dimensionless bioheat transfer equation, Eq. (1.3) is rewritten here for sake of convenience, r V2T-PfT+S f
AT
(2.1)
= —- in Q dt
The boundary conditions are T= T0 on IV
(2.2)
dT = dn
(2.3)
where p is a source point and q is the field point. P and Q represent a source point and a field point, respectively, on the boundary. A boundary integral for the determination of any interior point can be derived with the help of the divergence theorem, G(p,Q;t,T)
T(J>
t)r\ J r
dG(p,Q;t,T)
dTQdr + \ \ G(p,q;t,T)S{q,T)dttqdT
G(p,q;tfi)T(qW%
(2.6)
It should be pointed out that Eq. (2.6), the boundary integral suitable for the bioheat equation is identical to that of the transient heat conduction equation. The perfusion term in the bioheat equation is canceled by that in the adjoint equation. However, the fundamental solution G(p,q) to Eq. (2.5) is different from that of the transient heat conduction equation. Equation (2.5) can be solved by using Laplace transform, for 3-D,
o3
8y/it(rBi3(r-r„)onrc;Bi=^ a K,
(2.4)
TT, Tq, and Tc are disjoint and its union forms the entire boundary. To derive the boundary integral for the boundary element formulation, we consider the fundamental solution of the corresponding adjoint equation in an infinite domain, V 2 G - P f G + 5[x,iq) -Xi(p),t-
dtln
G(p,q;t,T)
q0onTq
and dT_
T{Q,t)
dT(Q,t)
T] -•
dG dt
(2.5)
exp
-Pf(j-t)-
4(7-0
H{r-t)
(2.7)
H{r-t)
(2.8)
and for 2-D, G(p,q;t,r) 1 4ir(T-t)
expl-Pf(T-t)-
4(r-0J
where r is the distance between the source point and the field point. A boundary integral equation equivalent to the bioheat equation may now be obtained by taking the limit as p tends to P. This gives,
Nomenclature
matrix coefficients of the final BEM Eq. (2.21) A„ = switched matrix coefficients, see Eq. (2.24) seeEqs. (2.19) and (2.22) B„ matrix coefficients of the final BEM Eq. (2.21) switched matrix coeffiBij cients, see Eq. (2.25) Biot number, see Eq. (2.4) Bi see Eqs. (2.20) and (2.22) dimensionless parameter reQ lated to the cooling length, see Eq. (3.5) heat capacity c shape coefficient c(P) width of the blood vessel d Green's functions G heat transfer coefficient h Heaviside step function H(t) thermal conductivity k characteristic dimension of L the tissue outward normal n dimensionless perfusion Pr coefficient, see Eq. (1.4) A,,
Journal of Biomechanical Engineering
P = source point located on the boundary Peclet number, see Eq. Pe„ (3.7) source point within the domain Q = field point located on the boundary field point with the domain normal temperature gradient on the boundary 2q2)
j=l
dG(P„Q) WiTy + timdTr, dnn
For 3-D, G(p,q) = -
(2.16)
where the boundary of the domain T is divided into N boundary segments. T(P,) represents the temperature at a point P that coincides with the ith node. A suitable shape function i/y must now be chosen for the variation of temperature and flux over the boundary elements ATj. In the present work, both temperature and flux are assumed to be linear over individual boundary elements.
dVr
G(j,,q)S(q)dQq
dG(PhQ) T(Q) dVQ; dnn
dT(Q) q\Q) = dn Q
subject to the appropriate boundary conditions given in Eqs. (2.2) through (2.4). The corresponding fundamental adjoint equation in an infinite domain is V2G-PfG
(.
N
c(P,)nP,)-.
(2.18)
Here r; is a dimensionless local coordinate over individual boundary segments while Tu T2, q\, and q2 are the nodal quantities of they'th boundary segment. Defining the following variables,
4 = o.
Numerical Implementation. In this section, the numerical implementation of the boundary integral equations (2.12 and 2.15) for the steady-state bioheat equation is presented. For simplicity, the two-dimensional versions are used in this work. The procedure can be easily extended to treat three-dimensional and transient formulations. Further, the source term is as3 6 0 / V o l . 114, AUGUST 1992
bf/= \
ypkG{PhQ)dTQ
(2.20)
Ary
The discretized BEM equation can be written as
where Aij^ay-i
+ alj and Bij=bfj-1
+ bjj
(2.22)
Integrals of kernels over elements in Eq. (2.20) must be evaluated with caution. In the present work, the diagonal element of G(P„ 0 , Eq. (2.20), which are log(r) singular may be evaluated numerically through an algorithm for improper integrals (Press et al., 1986). At each location over the entire boundary of the domain, either T, q, or a combination of both is prescribed for a wellposed problem. Equation (2.21) may be rearranged as, N
N
E-W^E^./
y=i
(2.23)
J'=I
Transactions of the ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/18/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
where ii
for q" specified (2.24)
for 7} specified
-B„
q = 0
(1.1:
- Nu Br, for convective condition ! for q" specified for Tj specified
B„
Bij =
(2.25)
- Nu By for convective condition,
y«0_ 1 J ~
(Tj for q" specified
^
qj for Tj specified
>
Uniform boundary mesh 0.1 length, total of 40 elements
(2.26)
\Tj for convective conditionj
->- x
Cqj for q" specified Yjk)=
q = 0
7} for 7} specified
(2.27)
( T„ for convective condition ,
Fig. 1 Computational domain for the BEM formulation of the one-dimensional bioheat equation
This is generally referred as the switching process. It should be pointed out that the switching process for a geometrical corner is more complicated. As defined earlier, the flux q" is equal to the normal temperature gradient. For a geometrical corner, there is an ambiguity in the normal direction. A corner algorithm has been developed by Chan and Chandra (1991b) to resolve this problem. This corner algorithm is also applied to the BEM formulation for the bioheat transfer equation. Equation (2.23) may now be used to solve for the unknown temperature and flux. The temperature at an internal point may be calculated from discretized version of Eq. (2.12).
T(p) = 2^[BW-A?Tj]
(2.28)
j=i
where Apj= af.! + a? and B* = b
Boundary Element Method Analysis for the Bioheat Transfer Equation In this paper, the boundary element method (BEM) approach is applied to solve the Pennes (1948) bioheat equation. The objective is to develop the BEM formulation and demonstrate its feasibility. The basic BEMformulations for the transient and steady-state cases are first presented. To demonstrate the usefulness of the BEM approach, numerical solutions for 2-D steady-state problems are obtained and compared to analytical solutions. Further, the BEM formulation is applied to model a conjugate problem for an artery imbedded in a perfused heated tissue. Analytical solution is possible when the conduction in the x-direction is negligible. The BEM and analytical results have very good agreement.
Introduction The complexity of the mathematical modeling of the heat the BEM results are more accurate than the FEM. This implies transfer process within a human tissue is enormous. The dif- that for the same desired level of accuracy the BEM requires ficulties arise from the fact that the substance is extremely a coarser mesh than the FEM does. They also found that the heterogeneous (Chen, 1985). The basic heat transfer processes BEM accuracy, with a fixed boundary mesh, does not change in biological systems can be found in a recent review by Chato appreciably with the number of internal sampling points. On (1990). The most commonly used model is the bioheat transfer the other hand, the FEM accuracy is very sensitive to the equation developed by Pennes (1948). He postulated that the internal mesh. effect of the temperature difference between the blood supply In the present work, the BEM formulation for the bioheat and the tissue acts as an energy sink term, Qb = wbcb(T- Ta). transfer equation is first presented. To validate the BEM apThis in essence is very similar to the fin equation where the proach, numerical solutions of a 2-D steady-state problem are sink term represents convective heat loss to the surrounding. presented and compared to an analytical solution. Finally, to The bioheat transfer equation has been applied by many illustrate the capability of the BEM approach, a solution of investigators to hyperthermia. It has been found that the bio- the heat transfer within tissue adjacent to an artery is presented. heat transfer equation provides useful and adequate predic- Analytical solution is possible when the conduction in the xtions (Strohbehn and Roemer, 1984; Roemer, 1990) in many direction is negligible. The BEM and analytical results are cases. A review of numerical solutions of the bioheat transfer presented and compared. equation can be found in (Strohbehn and Roemer, 1984). It suffices to mention that the numerical methods used are mainly Mathematical Formulation the finite difference method and finite element method. The bioheat transfer equation proposed by Pennes (1948) The Boundary Element Method (BEM) is another powerful can be written as, general purpose method (Banerjee and Butterfield, 1981; Mukherjee, 1982; Brebbia et al., 1984). It is far more tolerant of 97"* klV2T*-wbcb(T*-T*) + S*=p,ct — inQ* (1.1) aspect ratio degradation than the Finite Element Method (FEM) and can yield secondary variables as accurately as the primary ones. The internal equations are applied pointwise, hence, Introducing the dimensionless variables, sharp temperature gradients over the domain may be captured (x*,y\z*) „ T'-T* t t*a, S=(1.2) easily. The discretization of the boundary can be tailored to (x,y,z)-, t-- Ll ', 'T=-S L2/k, ' ~ L " S0 0 any irregular shape, thus it is capable of handling complex' geometries. Mukherjee and Morjaria (1982) compared the ef- the governing equation can be expressed in dimensionless form ficiency and accuracy of the BEM and FEM for the Laplace as, equation. They found that for the same level of discretization 3T. V2T-P/T+S = :^ in Q (1.3) dt where Contributed by the Bioengineering Division for publication in the JOURNAL OF BIOMECHANICAI ENGINEERING. Manuscript received by the Bioengineering Division March 11, 1991; revised manuscript received October 26, 1991. Associate Technical Editor: K. R. Diller.
358 / Vol, 114, AUGUST 1992
Pr=
WbCbL2
(1.4)
Transactions of the ASME
Copyright © 1992 by ASME Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/18/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
In addition, appropriate boundary conditions must be specified along the boundary for a well posed problem. These will be presented with specific examples in later sections. Boundary Element Formulation In the Boundary Element Method, the essential steps are (i) derivation of the equivalent boundary integral equation, (ii) determination of the Green's function, and (iii) discretization of the boundary integral. The BEM can be applied most efficiently if an analytical expression for the Green's function of the adjoint equation is available. In this section, we derive the boundary integral equations and the relevant Green's functions for both unsteady and steady problems. We will illustrate the discretization of the boundary integral using the two-dimensional steady-state equation. Transient Formulation. The dimensionless bioheat transfer equation, Eq. (1.3) is rewritten here for sake of convenience, r V2T-PfT+S f
AT
(2.1)
= —- in Q dt
The boundary conditions are T= T0 on IV
(2.2)
dT = dn
(2.3)
where p is a source point and q is the field point. P and Q represent a source point and a field point, respectively, on the boundary. A boundary integral for the determination of any interior point can be derived with the help of the divergence theorem, G(p,Q;t,T)
T(J>
t)r\ J r
dG(p,Q;t,T)
dTQdr + \ \ G(p,q;t,T)S{q,T)dttqdT
G(p,q;tfi)T(qW%
(2.6)
It should be pointed out that Eq. (2.6), the boundary integral suitable for the bioheat equation is identical to that of the transient heat conduction equation. The perfusion term in the bioheat equation is canceled by that in the adjoint equation. However, the fundamental solution G(p,q) to Eq. (2.5) is different from that of the transient heat conduction equation. Equation (2.5) can be solved by using Laplace transform, for 3-D,
o3
8y/it(rBi3(r-r„)onrc;Bi=^ a K,
(2.4)
TT, Tq, and Tc are disjoint and its union forms the entire boundary. To derive the boundary integral for the boundary element formulation, we consider the fundamental solution of the corresponding adjoint equation in an infinite domain, V 2 G - P f G + 5[x,iq) -Xi(p),t-
dtln
G(p,q;t,T)
q0onTq
and dT_
T{Q,t)
dT(Q,t)
T] -•
dG dt
(2.5)
exp
-Pf(j-t)-
4(7-0
H{r-t)
(2.7)
H{r-t)
(2.8)
and for 2-D, G(p,q;t,r) 1 4ir(T-t)
expl-Pf(T-t)-
4(r-0J
where r is the distance between the source point and the field point. A boundary integral equation equivalent to the bioheat equation may now be obtained by taking the limit as p tends to P. This gives,
Nomenclature
matrix coefficients of the final BEM Eq. (2.21) A„ = switched matrix coefficients, see Eq. (2.24) seeEqs. (2.19) and (2.22) B„ matrix coefficients of the final BEM Eq. (2.21) switched matrix coeffiBij cients, see Eq. (2.25) Biot number, see Eq. (2.4) Bi see Eqs. (2.20) and (2.22) dimensionless parameter reQ lated to the cooling length, see Eq. (3.5) heat capacity c shape coefficient c(P) width of the blood vessel d Green's functions G heat transfer coefficient h Heaviside step function H(t) thermal conductivity k characteristic dimension of L the tissue outward normal n dimensionless perfusion Pr coefficient, see Eq. (1.4) A,,
Journal of Biomechanical Engineering
P = source point located on the boundary Peclet number, see Eq. Pe„ (3.7) source point within the domain Q = field point located on the boundary field point with the domain normal temperature gradient on the boundary 2q2)
j=l
dG(P„Q) WiTy + timdTr, dnn
For 3-D, G(p,q) = -
(2.16)
where the boundary of the domain T is divided into N boundary segments. T(P,) represents the temperature at a point P that coincides with the ith node. A suitable shape function i/y must now be chosen for the variation of temperature and flux over the boundary elements ATj. In the present work, both temperature and flux are assumed to be linear over individual boundary elements.
dVr
G(j,,q)S(q)dQq
dG(PhQ) T(Q) dVQ; dnn
dT(Q) q\Q) = dn Q
subject to the appropriate boundary conditions given in Eqs. (2.2) through (2.4). The corresponding fundamental adjoint equation in an infinite domain is V2G-PfG
(.
N
c(P,)nP,)-.
(2.18)
Here r; is a dimensionless local coordinate over individual boundary segments while Tu T2, q\, and q2 are the nodal quantities of they'th boundary segment. Defining the following variables,
4 = o.
Numerical Implementation. In this section, the numerical implementation of the boundary integral equations (2.12 and 2.15) for the steady-state bioheat equation is presented. For simplicity, the two-dimensional versions are used in this work. The procedure can be easily extended to treat three-dimensional and transient formulations. Further, the source term is as3 6 0 / V o l . 114, AUGUST 1992
bf/= \
ypkG{PhQ)dTQ
(2.20)
Ary
The discretized BEM equation can be written as
where Aij^ay-i
+ alj and Bij=bfj-1
+ bjj
(2.22)
Integrals of kernels over elements in Eq. (2.20) must be evaluated with caution. In the present work, the diagonal element of G(P„ 0 , Eq. (2.20), which are log(r) singular may be evaluated numerically through an algorithm for improper integrals (Press et al., 1986). At each location over the entire boundary of the domain, either T, q, or a combination of both is prescribed for a wellposed problem. Equation (2.21) may be rearranged as, N
N
E-W^E^./
y=i
(2.23)
J'=I
Transactions of the ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/18/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
where ii
for q" specified (2.24)
for 7} specified
-B„
q = 0
(1.1:
- Nu Br, for convective condition ! for q" specified for Tj specified
B„
Bij =
(2.25)
- Nu By for convective condition,
y«0_ 1 J ~
(Tj for q" specified
^
qj for Tj specified
>
Uniform boundary mesh 0.1 length, total of 40 elements
(2.26)
\Tj for convective conditionj
->- x
Cqj for q" specified Yjk)=
q = 0
7} for 7} specified
(2.27)
( T„ for convective condition ,
Fig. 1 Computational domain for the BEM formulation of the one-dimensional bioheat equation
This is generally referred as the switching process. It should be pointed out that the switching process for a geometrical corner is more complicated. As defined earlier, the flux q" is equal to the normal temperature gradient. For a geometrical corner, there is an ambiguity in the normal direction. A corner algorithm has been developed by Chan and Chandra (1991b) to resolve this problem. This corner algorithm is also applied to the BEM formulation for the bioheat transfer equation. Equation (2.23) may now be used to solve for the unknown temperature and flux. The temperature at an internal point may be calculated from discretized version of Eq. (2.12).
T(p) = 2^[BW-A?Tj]
(2.28)
j=i
where Apj= af.! + a? and B* = b
