An immersed interface method for Pennes bioheat transfer equation

We consider an immersed finite element method for solving one dimensional Pennes bioheat transfer equation with discontinuous coefficients and nonhomogenous flux jump condition. Convergence properties of the semidiscrete and fully discrete schemes are investigated in the $L^{2}$ and energy norms. By using the computed solution from the immerse finite element method, an inexpensive and effective flux recovery technique is employed to approximate flux over the whole domain. Optimal order convergence is proved for the immersed finite element approximation and its flux. Results of the simulation confirm the convergence analysis.

The heat transfer analysis of biological bodies during surface heating or cooling has been studied by many researchers [5,7,12,13]. In many diagnostic and therapeutic applications of heat transfer, evaluation of transient and spatial distribution of temperature is required both on the skin surface and inside the living tissues. Investigation of transient temperature distribution requires solving a bioheat transfer equation with relevant boundary and initial conditions. Several numerical methods including finite element (FEM) and finite difference (FDM), have been studied to simulate Pennes bioheat transfer equation [3,15]. Finite difference schemes were more developed to simulate Pennes equation in 1-D and 3-D in a triple layered skin structure allowing variable thermal properties across the layers. We refer the reader to [3,4] and the references therein for more details.
Obtaining the analytical solution of Pennes bioheat transfer equation in a multiple layered structure is difficult due to the variable thermal parameters and the interface conditions. In order to successfully approximate the transient temperature distribution in the tissues and at the skin surface of the biological body, it is necessary to consider the variation of thermal properties that appear in the bioheat transfer equation. Therefore, one of the goals of this article is to simulate Pennes bioheat transfer equation using the immersed finite element method (IFEM), which is suitable for handling variable thermal properties when heat transfer through the different layers of the biological body.
An IFEM requires construction of special basis functions for the interface elements while keeping the standard basis functions for non-interface elements. The interface basis functions are constructed to satisfy the interface jump conditions. Detailed explanation of this IFEM can be found in [8,9,10,11], and references therein. In this paper due to the one dimensional nature of the problem, we used a special function to reduce the nonhomogeneous jump condition to a homogeneous one, making it easier to prove optimal convergence. In higher dimensions it is much harder to find such functions. See He et al. [6] for related issues.
In many applications it is also important to accurately approximate the heat flux u = −βT x as well as the solution T of the initial boundary value problem. Even though the finite element method provides a symmetric positive definite system to solve for the variable T , it does not approximate flux automatically. Chou and Tang [1] introduced a technique to approximate flux of the solution T , when numerical solution is approximated accurately by the standard finite element method. Chou [2] extended this method to IFE spaces (both cases are for elliptic problems). It has the special merit of reproducing the exact fluxes at the nodes and interface. In Section 3, we extend this methodology to parabolic problems and show that the approximate flux has optimal order convergence.
The rest of the paper is organized as follows. The error analysis for the semidiscrete problem and the fully discrete problem is given in Section 2. Section 3 is devoted to the convergence of the flux approximations. Finally, numerical results are presented in Section 4.
2. IFE method. In this section we study an immersed finite element method for the bioheat interface problem (1.1). In the error analysis, we shall use the standard notation for the Sobolev spaces H k (I) endowed with Sobolev norm · k,I . H k 0 (I) is the usual space with homogenous boundary conditions. We use the shorthand H k for H k (I) when I = (a, b). Let α be the interface point in (a, b). For the ensuing analysis we need the function spaces equipped with the norm u 2 2,α := u 2 2,(a,α) + u 2 2,(α,b) .
To derive a natural weak formulation of the problem (1.1), we integrate it against q ∈ H 1 0 and use integration by parts on each subinterval. The resulting formulation is where a(T, q) = (βT , q ) + (T, q), Rewriting T as T = p +p, p ∈ H 1 0 , we see that it is sufficient to find p : J → H 1 0 such that , where the right hand side is a bounded linear functional on H 1 0 , using the Sobolev imbedding theorem on the last term. Note that by using smooth testing functions it can be shown that It is also easy to see that the bilinear form a(·, ·) is bounded and coercive on H 1 0 . To approximate the weak problem (2.2), we consider the one dimensional linear IFE space introduced in [8,10]. Let a = x 0 < x 1 , . . . < x N = b be a partition of I = [a, b] and let the interface point α be located in (x κ , x κ+1 ) for some κ. We term [x κ , x κ+1 ] an interface element and the remaining elements non-interface elements. Let h i = x i − x i−1 , (i = 1, 2, . . . , N ). In the IFE space, for non-interface elements the local shape function is the usual linear function. For the interface element a local shape function takes the form It is easy to see that functionφ is uniquely determined by the above conditions [8]. The linear IFE space is defined as S h = span{φ i } N −1 i=1 , whereφ i (x k ) = δ ik andφ i is piecewise linear and glued together by the local shape basis functions. It is clear that S h ⊂ H 1 0 . Moreover, Lin et al. [11] proved that the interpolation error in this space has optimal order. That is, for . The semidiscrete immersed interface finite element problem based on the above weak formulation is: with the initial condition p h (0) = π h p 0 . Of course, the approximate T h to T can be obtained through T h = p h +p. In the following we will concentrate on the error In this section we study the convergence properties of the semidiscrete IFE method. We shall show that the error p − p h between the solutions of the continuous and spatially discrete problems are of optimal order. Let us first define the elliptic projection of Note that R h w is well defined by the coercivity of a onŜ h due to the conformingness In fact, their proofs are similar to the standard linear elements, provided that we make sure the constants C are independent of the interface parameter α. For example, in the proof of the first inequality we have to use the coercivity constant min{min{β − , β + }, 1}, the norm of a, and the constant in (2.4), which are all independent of α. The second inequality can be proved using the first one and the usual duality argument while checking that none of the constants involved depend on α. Unlike in 2-D IFE (nonconforming) methods, here we are having a conforming method and proving the above type of two inequalities is much simpler. Now we study the parabolic equation. As in the study of standard parabolic equations [14], we introduce the elliptic projection p * = R h p. Moreover, substituting w by p and differentiating (2.6) with respect to t leads to (2.7) Now the total error is written as and the elliptic projection error satisfies α . The next theorem shows that the error between the semidiscrete solution and the continuous solution is of optimal order in the L 2 norm. Theorem 2.1. Let p h and p be the solutions of the semidiscrete problem (2.5) and the parabolic problem (2.2), respectively. Then there exists a positive constant C independent of h such that Proof. We only have to find error bound for θ. It follows from (2.2), (2.5) and (2.6) that Using the fact a(θ, θ) ≥ 0, we get where the 2 term was added since γ 1/2 θ 2 0 may not differentiable when θ = 0. Since the second term is less than the fourth term, we have Canceling the common factor and then integrating both sides with respect to t and let → 0, we get Using the equivalence of the norms we have Letting p h (0) =π h p 0 we can see that, Therefore, substituting the estimates for θ(0) 0 and ρ t 0 in (2.10) gives Finally the total error This completes the proof.
In a similar fashion, we can obtain the error estimate for the semidiscrete solution in the energy norm.

CHAMPIKE ATTANAYAKE AND SO-HSIANG CHOU
Proof. By settingφ = θ t in (2.9) we obtain or, explicitly Let γ max = max{γ − , γ + } and γ min = min{γ − , γ + }. Then Taking small enough to absorb the θ t term into the left side e.g., ≤ γ min /γ max , we conclude 1 2 . Integrating with respect to t and then using the equivalence of the norms, it follows that where C depends on β and γ. For p h (0) =πp 0 , θ(0) 1 ≤ Ch p 0 2,α . Substituting the estimates for ρ t 0 and θ(0) 1 in the above inequality we have which completes the proof.

2.2.
Completely discrete error. In this section we analyze a fully discrete scheme based on the backward Euler finite difference approximation. Optimal order convergence in the L 2 norm is established. We first partition the time interval [0,t] into M subintervals by 0 = t 0 < t 1 < . . . < t M =t points, with t n = nk and k =t/M , the uniform time stepsize. Then the completely discrete backward Euler Galerkin approximation p n N ∈ S h to p n = p(t n ) of problem (1.1) is defined as the solution of where l(q) = (T a , q) − a(p, q) − Qq(α) and Under appropriate regularity assumptions for p, the following theorem shows the fully discrete error has convergence rate of order O(h 2 + k). Proof. As in the previous theorem, we split the error into two terms as , the elliptic projection defined in (2.6). Since the estimates for ρ j are known, we only have to find a bound for θ j . It follows from (2.12), that Now define [14] ω (2.14) For notational convenience, let us defineθ = γ 1/2 θ,ω = γ 1/2 ω and so on. Takinĝ φ = θ j , (∂θ j ,θ j ) + a(θ j , θ j ) = −(ω j ,θ j ).
Summing over j from 1 to n leads to Using the norm equivalences we have As for the ω 2 term we have Proof. By settingφ =∂θ j in (2.13) we obtain, Note that for any sequence η j , ∂ η j 2 0 ≥ 0 implies Therefore, just as in the proof of Theorem 2.2, we can use the -inequality ab ≤ a 2 + 1 4 b 2 , a, b > 0 for > 0 small enough to absorb the ∂ θ j 0 term on the right side into the first term on the left side, when estimating (2.15). Thus using this technique, we again derivē ∂ β 1/2 θ j x 2 0 +∂ θ j 2 0 ≤ C γω j 2 0 , where ω j is as in the proof of the last theorem. Summing the above inequality over j from j = 1 to n and use equivalence of norms, we have The remaining part is just like in the last theorem.
3. Flux recovery. In this section we will construct a method to approximate the exact flux −βT x in (1.1). Since we have concentrated on computing p h to approximate p in T = p +p. We will first transform (1.1) accordingly. The transformed problem is Thus it suffices to approximate u = −βp x , which will be done by a flux recovery process in which once p N ≈ p has been computed, the approximate flux u N is obtained by using a simple formula. The u N is globally continuous and linear elementwise except on the interface element. On the interface element, u N is piecewise linear.
We now describe how to construct the approximate flux on non-interface elements first. Multiply (3.1) byφ j and integrate by parts over Do the same over I j+1 = [x j , x j+1 ] to get Since p N is a good approximate of p for a fixed time t n , it is natural to define the discrete fluxes u N (x − j , t n ) and u N (x + j , t n ) on I j and I j+1 , respectively as, f (x, t n )φ j dx.

Thus, from (3.3) and (3.4) we see that
In a similar fashion by makingφ =φ j+1 in (2.12), it is easy to see that . Therefore, following the flux approximation method in [1], for a given t n ∈ J and for any x ∈ [x j , x j+1 ], j = k, we define the approximate flux u N (x, t n ) as (3.5) To approximate flux on the interface element, I κ , we multiply (3.1) by the interface basis functionφ κ , integrate by parts over [x κ , α] and [α, x κ+1 ], and use the interface conditions to get which has the same form as for a noninterface element. For the adjacent noninterface element Accordingly we define Again the fact that u N (x − κ , t n ) = u N (x + κ , t n ) can be checked by settingφ =φ κ in (2.12). It is then not hard to see that all we need to check is Similarly, it is easy to see that u N (x − κ+1 , t n ) = u N (x + κ+1 , t n ) = u N (x κ+1 , t n ). Now we need to define u N (α, t n ), the approximate flux at the interface point. For this we take clue from (3.1) again. Integrating it over [x κ , α] and using u = −βp x , we get (3.10) Similarly, integrating over [α, Thus accordingly we define Now we claim u N (α − , t n ) = u N (α + , t n ). Recall (3.8) and Substituting these relations into (3.11)-(3.12) and usingφ κ +φ κ+1 = 1 give the desired result. So we can define u N (α, t n ) = u N (α ± , t n ). Finally we define the global approximate flux by interpolation at the nodes and the interface point. In particular, u N (x, t n ) on the interface element [x κ , x κ+1 ] has the formula (3.13) where λ κ , λ α − , and λ κ+1 , λ α + are the linear Lagrange nodal basis functions defined in [x κ , α] and [α, x κ+1 ], respectively. We summarize the above findings in the following theorem.

Numerical examples.
In this section we give two numerical examples to confirm our theory. (4.1) Here,  Table 2. Error analysis of the solution and the flux for the case with β − = 100 and β + = 1.

Example 2.
To illustrate our method applied to Pennes bioheat equation we consider a double layered tissue structure in the domain of 1 mm cross section. The first and the second layers have a thickness of of 0.33 mm and 0.67 mm, respectively. More specifically, consider the following initial boundary value problem in [0, 1] with an interface at 0.13: 1]. The elevation of the temperature at the left end of the tissue is kept constant, and the heat flux approaches 0 at the other end [3]. This is realistic for a situation such as a skin burn in a biological body. Thus, the initial and boundary conditions are T (0, t) = 12 • C for t > 0, T x (1, t) = 0 for t > 0. At the interface, we assume the temperature T satisfies the jump conditions [T ] α = 0, [βT x ] α = Q = 1.
The thermal properties for the two layer are given in Table 3 [3]. Temperature distribution and flux distribution for our example are shown in Figure 1 and 2, respectively.