Boundary estimation of hyperbolic bio-heat conduction

Boundary estimation of hyperbolic bio-heat conduction

Ching-yu Yang

Department of Mold and Die Engineering, National Kaohsiung University of Applied Sciences, Kaohsiung City 807, Taiwan, ROC

a r t i c l e

i n f o

Article history:

Received 30 August 2010

Received in revised form 7 January 2011 Accepted 7 January 2011

Available online 25 February 2011 Keywords:

Inverse bio-heat conduction problem Hyperbolic bio-heat transfer Biological material

a b s t r a c t

A sequential method is proposed for estimating the boundary condition in hyperbolic bio-heat conduc-tion. The estimated solution is deduced from a numerical approach combined with the concept of future time. The problem with inverse bio-heat conduction is the slow heat-wave propagation speed, resulting in no temperature measurements obtained. Three cases are presented to demonstrate the features and the validity of the proposed method. Comparison between the exact value and the estimated result is made to confirm the validity and accuracy of the proposed method.

Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction

During the past decades, the non-Fourier effect on solving prac-tical engineering problems has attracted more and more attention. The mathematical representation of the non-Fourier law is a hyperbolic form with a finite wave propagation speed. Many meth-ods for solving direct and inverse non-Fourier problems concerning engineering materials have been proposed, such as the non-homogenous solids conduction process, the rapid heating process, and the slow conduction process[1–22]. In contrast, approaches applicable to the domain of biological materials are scarce[23–27]. Owing to the highly heterogeneous nature of living biological tissues, the speed of non-Fourier bio-heat wave propagation is much slower and less predicable than that in engineering materi-als. Despite that, there are still some studies on non-Fourier bio-heat bio-heat conduction. Yang[23]investigated the mechanisms of thermal shock and analyzed the interface thermal resistance. Liu et al. [24]examined skin burn evaluation, which was simulated using a finite difference method. Shih et al.[25]also used the finite difference method to simulate thermal dose distribution during thermal therapy. Liu and Cheng[26,27]used the Laplace transform with a numerical spatial scheme to analyze the non-Fourier tissue problem. The above studies are mainly direct analysis of non-Fourier bio-heat conduction[23–27], while the inverse non-Fou-rier boundary remains largely unknown. In the inverse domain, the boundary condition is estimated from the temperature mea-sured at a location different from that of the boundary applied. The slow speed of heat-wave propagation in the bio-tissue makes it difficult for temperature measurements. In other words, indirect measurement at the sensor location cannot acquire the real-time

boundary input, posing computational difficulty to the inverse analysis. To overcome the problem, this study developed a robust and stable method for estimating the boundary condition in in-verse hyperbolic bio-heat conduction.

The proposed approach is a sequential method combined with the concept of future time. It is also a numerical method that can estimate boundary condition sequentially without the need of sen-sitivity analysis. In the proposed method, a closed-form is derived from a numerical model to represent explicitly the unknown boundary conditions. A finite-element-difference method com-bined with the concept of future time[28]for boundary condition estimation, which is determined step by step along with the tem-poral coordinate.

The rest of the paper is organized as follows. Section 2 intro-duces the current development in non-Fourier conduction and states the features of the proposed method. Section 3 delineates the characteristics of inverse problem and describes the proposed method. In Section 4, a computational algorithm is developed to implement the proposed method. Example problems are presented in Section 5 to demonstrate the application of the proposed meth-od with discussion of the analyzed results. Finally, Section 6 de-scribes the overall contribution and possible applications of this research to inverse bio-heat conduction problems.

2. Problem statement

The inverse hyperbolic bio-heat problem is about finding the boundary condition at one side of the tissue when the temperature measurement is available at a different location. Consider a tissue with length l and finite heat propagation effect. Assume that the value of metabolic heat generation is constant, without any spatial heating source, the value of thermal conductivity is constant and 0017-9310/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved.

doi:10.1016/j.ijheatmasstransfer.2011.02.011

E-mail address:cyyang@cc.kuas.edu.tw

International Journal of Heat and Mass Transfer 54 (2011) 2506–2513

Contents lists available atScienceDirect

International Journal of Heat and Mass Transfer

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j h m t

the initial temperature is distributed all over the tissue medium. The thermal model of bio-heat transfer is shown as follows[24,26]:

q

c

s

@ 2 h @t2þ ð

q

c þ

s

wbcbÞ @h @tþ wbcbh k @2h @x2¼ 0 ð1Þ hðx; 0Þ ¼ 0 0 6 x 6 l ð2Þ @hðx; tÞ @t ¼ 0 at t ¼ 0 and 0 6 x 6 l ð3Þ hðx; tÞ ¼ f ðtÞ at x ¼ 0 and t > 0 ð4Þ qðx; tÞ ¼ 0 at x ¼ l and t > 0 ð5Þ

where

q

, c,

s

and k are density, specific heat, relaxation time and thermal conductivity of tissue, respectively. wband cbare the

perfu-sion and specific heat of blood, respectively. Moreover, the elevation temperature h is defined as h = T  T0, where T is the temperature

response and T0is the initial temperature distribution of the tissue.

Inverse analysis involves estimating the boundary condition f(t) when the temperature is measured at a location other than x = 0. 3. Method for determining boundary conditions

In this section, a numerical formulation is employed to represent Eqs.(1)–(5), in which the undetermined boundary condition is ex-pressed as a vector form in Eq.(4). The temperature distribution can then be represented as a function of the boundary condition. The concept of future time is incorporated into the proposed method. Finally, the inverse solution can be represented as a matrix equation. The proposed method uses a finite-element method with linear elements to discretize the spatial coordinate. By the conventional finite-element procedure with np= 0 grids at t =tj[29], Eqs.(1)–

(5)can be converted into the following discrete form:

½Bf€hjg þ ½Cf _hjg ¼ fRjg  ½Afhjg ð6Þ

where

[A] is the heat matrix of the problem with npdimensions.

[B] and [C] are the transient matrixes of the problem with np

dimensions.

{Rj} is the boundary vector with npcomponents.

{hj} is the temperature vector with npcomponents.

1

D

t2½Bfhj 2hj1þ hj2g þ 1

D

t½Cfhj hj1g ¼ fRjg  ½Afhjg ð7Þ ½Kfhjg ¼ ½Dfhj1g þ ½Efhj2g þ fRjg ð8Þ where ½K ¼ ½A þ 1 Dt2½B þD1t½C. ½D ¼ 2 Dt2½B þD1t½C. ½E ¼  1 Dt2½B. Rjg ¼ fRknownj n o þ Runknownj n o . Rknownj n o

is the vector of the known boundary condition. Runknown

j

n o

is the vector of the unknown boundary.

When t = tj, the temperature distribution {hj} can be derived

from Eq.(8)as follows:

fhjg ¼ ½K1½Dfhj1g þ ½K1½Efhj2g þ ½K1fRjg

¼ ½Ufhj1g þ ½Vfhj2g þ ½WfRjg ð9Þ

where [U] = [K]1_{[D], [V] = [K]}1_{[E] and [W] = [K]}1_.

Similarly, the temperature distribution at t = tm, tm+a, . . . , tm+r1

can be represented as follows: Nomenclature

[A] heat matrix [B], [C] transient matrix [D] 2 Dt2½B þDt1½C [E]  1 Dt2½B [K] ½A þ 1 Dt2½B þ 1 Dt½C

Nt number of the temporal measurements

[U] [K]1_[D] [V] [K]1_[E] [W] [K]1 [fU] intermediate matrix [fV] intermediate matrix [fW] intermediate matrix

[p] number of spatial measurement

q heat flux

r number of the future time

{R} boundary vector at th temporal coordinate

fRknown g vector of the known boundary at th temporal

coordi-nate Runknown



n o

unknown boundary vector at th temporal coordinate t temporal coordinate

T temperature

T0 initial temperature

bu_{c unit row vector with a unit at -component}

{u_{} unit column vectors with a unit at th component}

x spatial coordinate

a

, b,

c

intermediate variables

l

relative error of estimated solution

q

density of tissue

c specific heat of tissue

s

relaxation time of tissue k thermal conductivity of tissue wb perfusion of blood

cb specific heat of blood

Dt increment of temporal domain

k random number

r

standard deviation of measurement error / unknown temperature condition at th grid h evaluation temperature

hmeas measured temperature hsimu simulated temperature

_h derivative of evaluation temperature h €

h second derivative of evaluation temperature h

h vector from previous state and the present known boundary

X intermediate matrix

U vector form of estimated boundary [] matrix

{} column vector bc row vector Subscripts

i, j, k, l, m, n indices

N number of unknown boundary np number of grids at spatial coordinate

Superscripts i, l, n indices i, j, m, u, v indices

fhmg ¼ ½Ufhm1g þ ½Vfhm2g þ ½WfRmg ¼ f0 U   fhm1g þ fV0   fhm2g þ fW0   fRmg fhmþ1g ¼ ½Ufhmg þ ½Vfhm1g þ ½WfRmþ1g ¼ f1 U   fhm1g þ fV1   fhm2g þ fW1   fRmg þ fW0   fRmþ1g fhmþ2g ¼ ½Ufhmþ1g þ ½Vfhmg þ ½WfRmþ2g ¼ fU2   fhm1g þ fV2   fhm2g þ fW2   fRmg þ fW1   fRmþ1g þ fW0   fRmþ2g . . . fhmþr1g ¼ ½Ufhmþr2g þ ½Vfhmþr3g þ ½WfRmþr1g ¼ fr1 U   fhm1g þ fVr1   fhm2g þ fWr1   fRmg þ    þ f2 W   fRmþr3g þ fW1   fRmþr2g þ fW0   fRmþr1g ð10Þ where f0 U   ¼ ½U f0 V   ¼ ½V f0 W   ¼ ½W f1 W   ¼ ½U fW0   fiþ1 U   ¼ ½U fi U   þ fi V   fiþ1 V   ¼ ½U fVi   fiþ2 W   ¼ ½U fiþ1 W   þ ½V fi W   ð11Þ

Therefore, we have the temperature vector at (m + n)-temporal grid

fhmþng ¼ ½Ufhmþn1g þ ½Vfhmþn2g þ ½WfRmþng ¼ fUn   fhm1g þ fVn   fhm2g þ fWn   fRmg þ    þ fW2   fRmþn2g þ fW1   fRmþn1g þ fW0   fRmþng ¼ fUn   fhm1g þ fVn   fhm2g þ Xn k¼0 fnk W   fRmþkg ¼ fUn   fhm1g þ fVn   fhm2g þ Xn k¼0 fnk W   Rknown mþk n o þX n k¼0 fnk W   Runknown_mþk n o ¼ fn U   fhm1g þ fVn   fhm2g þ Xn k¼0 fnk W   Rknown_mþk n o þX n k¼0 fnk W   Runknown_mþk n o ð12Þ

where n is an integer and n = 0, 1, 2, . . . , r  1.

After multiplying a unit row vector bui_{c with both sides of}

Eq.(12), the temperature at i-spatial grid can be calculated as:

hi_mþn¼ bui_c½Ufh mþn1g þ buic½Vfhmþn2g þ buic½WfRmþng ¼ bui_{c f}n U   fhm1g þ buic fVn   fhm2g þ buic Xn k¼0 fnk W   Rknown_mþk n o þ bui_cX n k¼0 fnk W   Runknown mþk n o ¼ bui_{c f}n U   fhm1g þ buic fVn   fhm2g þ buic Xn k¼0 fnk W   Rknown_mþk n o þX N j¼1 Xn k¼0 bui_{c f}nk W   fuj_g/j mþk ð13Þ where Runknownmþk n o ¼PNj¼1fujg/ j

mþk and N is the number of the

esti-mated quantities. Here bui

c is a unit row vector with a unit at i-component; the value of i is the grid number of the measured grid. {uj_{} is the unit}

column vector with a unit at jth, component. Here, the value of /j_mþkis denoted as the unknown boundary and j is the grid number

of the location of the estimated boundary. The value of N is the number of unknown boundary conditions.

The temperatures at N-spatial grid at (m + n)-temporal grid can then be expressed as follows:

hi_mþn¼

a

i mþnþ b i mþnþ XN j¼1 Xn k¼0

c

i;j mþn;k/ j mþk ð14Þ where

a

i mþn¼ bu i_{c f}n U   fhm1g þ buic fVn   fhm2g bi_mþn¼ buicX n k¼0 fnk W   Rknown_mþk n o

c

i;j mþn;k¼ bu i_{c f}nk W   fuj_g

The value of n ranging from 0 to r  1 is substituted into Eq.(14),

him¼

a

i mþ b i mþ XN j¼1

c

i;j m;0/ j m hi_mþ1¼

a

i mþ1þ b i mþ1þ XN j¼1

c

i;j mþ1;0/ j mþ

c

i;j mþ1;1/ j mþ1 hi_mþ2¼

a

i mþ2þ b i mþ2þ XN j¼1

c

i;j mþ2;0/ j mþ

c

i;j mþ2;1/ j mþ1þ

c

i;j mþ2;2/ j mþ2 . . . hi_mþr1¼

a

i mþr1þ b i mþr1þ XN j¼1

c

i;j mþr1;0/ j mþ

c

i;j mþr1;1/ j mþ1 þ

c

i;j mþr1;2/ j mþ2þ    þ

c

i;j mþr1;r1/ j mþr1 ð15Þ

Here, superscripts i and j of

c

represent the measured grid and the estimated grid, respectively; and the subscript of

c

denotes the future time step.

The values of

c

i;j

mþn;k vary with the measured location and the

estimated location. Moreover, it also vary with the number of fu-ture time steps but not with the time step in the global temporal coordinate. In other words, the values of

c

i;j

mþn;k are constants in

each evaluation step and they need to be calculated only once when the location of the measured point and the input boundary are fixed. On the other hand, the coefficients in Eq.(14)

a

i

mþnand

bimþn are derived from the previous states {hm2}, {hm1} and the

present known boundary fRknown

mþk g. Therefore, these coefficients

need to be evaluated at each successive time step.

When t = tm, the estimated condition between t = t1and t = tm1

has been evaluated and the strength of the boundary condition at t = tmis to be estimated. For stabilizing the estimated result in the

in-verse algorithm, several future values of the estimation are temporally assumed to be constant in the subsequent procedure[28]. Then, the jth unknown condition at future time is assumed to be equal, i.e.,

/jm¼ / j mþ1¼    ¼ / j mþr2¼ / j mþr1 ð16Þ

Here r is the number of future time steps.

Substituting Eqs.(15) and (16)into Eq.(14), we get

him¼

a

imþ b i mþ XN j¼1

c

i;j m;0/ j m hi_mþ1¼

a

i mþ1þ b i mþ1þ XN j¼1

c

i;j mþ1;0þ

c

i;j mþ1;1   /jm hi_mþ2¼

a

i mþ2þ b i mþ2þ XN j¼1

c

i;j mþ2;0þ

c

i;j mþ2;1þ

c

i;j mþ2;2   /jm . . . hi_mþr1¼

a

i mþr1þ b i mþr1þ XN j¼1

c

i;j mþr1;0þ

c

i;j mþr1;1  þ

c

i;j mþr1;2þ    þ

c

i;j mþr1;r1  /jm ð17Þ

The estimated results obtained from the measurement ofFigs. 1–3are shown inFigs. 4–6, respectively. As can be seen, an accu-rate result can be approached after the heat-wave arrive at the sen-sor location, even before the heat-wave reaches the sensen-sor location, the estimated boundary can also be captured (seeFigs. 5 and 6). For instance, the heat wave arrives at the sensor location x = 0.0050333 m at t = 103.1526377 s but the boundary can be esti-mated from the time of boundary input (seeFigs. 2 and 5). More-over, the estimated result is accurate after t = 103.1526377 s when the temperature is measured at x = 0.0050333 m, indicating that the proposed inverse algorithm can estimate boundary condi-tion in hyperbolic bio-heat conduccondi-tion.

When measurement error is included in the solutions ofFigs. 1–3, the relative average error of the estimated results are shown inTable 1. The measurement error is set to range between 0.2576 and 0.2576, which implies that the average standard deviation of mea-surements is

r

= 0.1 for a 99% confidence bound. As can be seen, the larger the relative measurement error (i.e., 2:576

r

=hwhere his the mean value of the measurement), the less accurate the estima-tion is. For example, in Case 2, the values of the relative estimated er-rors are 0.0178951 and 0.0357065 when the relative measurement errors are 0.0297138 and 0.0573450, respectively (seeTable 1).

The feature of both direct and inverse analyses is discussed in this section. In direct analysis, Case 1 demonstrates the validity of the proposed method, Case 2 shows the magnitude of temperature re-sponse along the measured locations, and Case 3 reveals that the pro-file of temperature response becomes relatively smooth with increasing distance of the measurement location. In inverse analysis, an inverse result can be obtained even before the heat wave arrives at the sensor location. Moreover, an accurate result is obtained when the heat wave reaches the sensor location (seeFigs. 5 and 6). Inverse estimation with measurement error is also examined. The simulation results show that the larger the relative error, the greater the mea-sured error is obtained (seeTable 1). Our numerical results show that the proposed method is robust and stable. It can be concluded that the proposed method can deal effectively and accurately with the in-verse non-Fourier bio-heat conduction problems.

6. Conclusion

This study proposed an efficient algorithm for determining the boundary condition in inverse non-Fourier bio-heat conduction. The inverse solution is represented as a closed form derived by a finite-element-difference method when the temperature measure-ments are available. The inverse solution at each time grid is solved step by step along the temporal coordinate. The special feature of this method is that no pre-selected functional form for the un-known boundary is required and no sensitivity analysis is needed in the algorithm. Three cases are presented to illustrate the appli-cability of the proposed method. Result obtained shows that the proposed method can estimate accurately the boundary strength even at a slow heat-wave propagation speed. In conclusion, the numerical results show that the proposed method is an accurate and stable inverse technique for solving inverse non-Fourier bio-heat problems. The proposed method is applicable to other kinds of inverse bio-heat problems, such as source strength estimation in the non-Fourier domain.

References

[1] W.-B. Lor, H.-S. Chu, Effect of interface thermal resistance on heat transfer in a composite medium using the thermal wave model, Int. J. Heat Mass Transfer 43 (5) (2000) 653–663.

[2] P.J. Antaki, Importance of non-Fourier heat conduction in solid-phase reactions, Combust. Flame 112 (3) (1998) 329–341.

[3] T. Sanderson, C. Ume, J. Jarzynski, Hyperbolic heat equations in laser-generated ultrasound models, Ultrasonics 33 (6) (1995) 415–421.

[4] L.H. Liu, H.P. Tan, T.W. Tong, Non-Fourier effects on transient temperature response in semitransparent medium caused by laser pulse, Int. J. Heat Mass Transfer 44 (17) (2001) 3335–3344.

[5] A.M. Mullis, Rapid solidification and a finite velocity for the propagation of heat, Mater. Sci. Eng. A 226 (1997) 28–32.

[6] A.M. Mullis, Rapid solidification within the framework of a hyperbolic conduction model, Int. J. Heat Mass Transfer 40 (17) (1997) 4085–4094. [7] R.K. Sahoo, W. Roetzel, Hyperbolic axial dispersion model for heat exchangers,

Int. J. Heat Mass Transfer 45 (6) (2002) 1261–1270.

[8] W. Roetzel, S.K. Das, Hyperbolic axial dispersion model: concept and its application to a plate heat exchanger, Int. J. Heat Mass Transfer 38 (16) (1995) 3065–3076.

[9] W. Roetzel, C.N. Ranong, Consideration of maldistribution in heat exchangers using the hyperbolic dispersion model, Chem. Eng. Process. 38 (4–6) (1999) 675–681.

[10] J.-Y. Lin, The non-Fourier effect on the fin performance under periodic thermal conditions, Appl. Math. Model. 22 (8) (1998) 629–640.

[11] A.-H. Bishri, Modelling non-Fourier heat conduction with periodic thermal oscillation using the finite integral transform, Appl. Math. Model. 23 (12) (1999) 899–914.

[12] D.W. Tang, N. Araki, Non-Fourier heat conduction in a finite medium under periodic surface thermal disturbance, Int. J. Heat Mass Transfer 39 (18) (1996) 1585–1590.

[13] Z.-M. Tan, W.-J. Yang, Propagation of thermal waves in transient heat conduction in a thin film, J. Franklin Inst. 336 (1) (1999) 85–197.

[14] Z.-M. Tan, W.-J. Yang, Heat transfer during asymmetrical collision of thermal waves in a this film, Int. J. Heat Mass Transfer 40 (17) (1997) 3999–4006. [15] C. Weber, Analysis and solution of the ill-posed inverse heat conduction

problem, Int. Heat Mass Transfer 24 (1981) 1783–1791.

[16] N. Al-Khalidy, On the solution of parabolic and hyperbolic inverse heat conduction problems, Int. J. Heat Mass Transfer 41 (1998) 3731–3740. [17] N. Al-Khalidy, Analysis of boundary inverse heat conduction problems using

space marching with Savitzky–Gollay digital filter, Int. Commun. Heat Mass Transfer 26 (2) (1999) 99–208.

[18] H.T. Chen, S.Y. Peng, L.C. Fang, Numerical method for hyperbolic inverse heat conduction problems, Int. Commun. Heat Mass Transfer 28 (6) (2001) 847–856. [19] C.-y. Yang, Direct and inverse solutions of hyperbolic heat conduction

problems, J. Thermophys. Heat Transfer 19 (3) (2005) 217–225.

[20] C.-y. Yang, Estimation of the periodic thermal conditions on the non-Fourier fin problem, Int. J. Heat Mass Transfer 48 (17) (2005) 3506–3515.

[21] C.H. Huang, C.Y. Lin, Inverse hyperbolic conduction problem in estimating two unknown surface heat fluxes simultaneously, J. Thermophys. Heat Transfer 22 (4) (2008) 766–77.

[22] C.-C. Wang, Direct and inverse solutions with non-Fourier effect on the irregular shape, Int. J. Heat Mass Transfer 53 (13–14) (2010) 2685–2693. [23] W.-H. Yang, Thermal (heat) shock biothermomechanical viewpoint, ASME J.

Biomech. Eng. 115 (4b) (1993) 617–621.

[24] J. Liu, X. Chen, L.X. Xu, New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating, IEEE Trans. Biomed. Eng. 46 (4) (1999) 420–428.

[25] T.-C. Shih, H.-S. Kou, C.-T. Liauh, W.-L. Lin, The impact of thermal wave characteristics on thermal dose distribution during thermal therapy: a numerical study, Med. Phys. 32 (9) (2005) 3029–3036.

[26] K.-C. i Liu, P.-J. Cheng, Finite propagation of heat transfer in a multilayer tissue, J. Thermophys. Heat Transfer 22 (4) (2008) 775–782.

[27] K.-C. Liu, Thermal propagation analysis for living tissue with surface heating, Int. J. Thermal Sci. 47 (5) (2008) 507–513.

[28] J.V. Beck, B. Blackwell, C.R.St. Clair, Inverse Heat Conduction: Ill-Posed Problem, first ed., Wiley, New York, 1985.

[29] M.A. Celia, W.G. Gray, Numerical Methods for Differential Equations, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1992.

[30] B. Carnahan, H.A. Luther, J.O. Wilkes, Applied Numerical Methods, first ed., John Wiley & Sons, Inc., New York, 1977.

[31] User’s Manual: Math Library Version 1.0, IMSL Library Edition 10.0, IMSL, Houston, TX, 1987.

Table 1

Relative errors of example problems whenr= 0.1.

Measured location 0.0010070 m 0.0020130 m 0.0030200 m 0.0040266 m 0.0050333 m 0.0060400 m Measurement error 0.0297138 0.0410410 0.0573450 0.0821660 0.1222040 0.1906510 Estimated error – Case 1 0.0139403 0.0134107 0.0181791 0.0226911 0.0267768 0.0372616 Estimated error – Case 2 0.0178951 0.0246515 0.0357065 0.0459909 0.0584220 0.0827630 Estimated error – Case 3 0.0152264 0.0172583 0.0264882 0.0368945 0.0498845 0.0669141