Nonlinear Behavior of Thermal Lagging in Laser-Irradiated Layered Tissue

Abstract

This work theoretically studies the temperature rise behavior in laser-irradiated layered tissue, which was stratified into skin, fat, and muscle. Taking the effect of microstructural interaction into account, a modified non-Fourier equation of bioheat transfer was developed based on the nonlinear form of the dual-phase-lag model. The broad beam irradiation method was adopted to calculate the laser energy deposition in tissues. The nonlinear high order partial differential equation and the complicated boundary conditions at the interface between two adjacent layers cause some mathematical difficulties. It is not easy to deal with such a problem. Therefore, a hybrid numerical scheme has been proposed to solve the present problem. The comparison between the linear and nonlinear effects of phase lag times on bioheat transfer has been made.

1. Introduction

From the viewpoint of the thermal treatment, if well controlled, the high-intensity pulsed heating can produce accurate and appropriate heating amount. The entire treatment time can be shortened, and the treatment quality will be enhanced. This characteristic is particularly important to the bulky objects. Therefore, in the last decade, the applications of pulse heating in the skin treatment are booming, and a revolutionary impact is on the skin beauty. For example, in the past surgical excision, dermabrasion and other methods were used to remove tattoos, which can be easily removed with a laser. In addition, a variety of skin vascular lesions, pigmented birthmarks, moles, spots, scars, even hair, and wrinkles are removed with a laser.

The Pennes bioheat transfer model is commonly used to simulate thermal behavior in biological bodies due to simplicity and validity. The Pennes bioheat equation describes the thermal behavior based on the classical Fourier's law which depicts an infinitely fast propagation of thermal signal. In reality, the living tissues are highly nonhomogenous and need a relaxation time to accumulate enough energy to transfer to the nearest element. As a result, to solve the paradox which occurred in the Pennes model, the thermal wave model of bioheat transfer was proposed for the investigation of physical mechanisms and the behaviors in thermal wave propagation in living tissues [1].

The properties of hyperbolic diffusion in bioheat transfer have attracted the relevant researchers' attention. Liu [2] and Özen et al. [3], respectively, studied the thermal wave propagation bioheat transfer in homogenous and multilayer tissues. Shih et al. [4] explored the impact of thermal wave characteristics on thermal dose distribution during thermal therapy. However, in order to consider the effect of microstructural interactions in the fast transient process of heat transport, a phase lag for temperature gradient absent in the thermal wave model was introduced [5]. The corresponding model is called the dual-phase-lag (DPL) model. Antaki [6] used the DPL model to interpret heat conduction in processed meat that was interpreted with the thermal wave model. Antaki [6] estimated the phase lag time of heat flux to be 14–16 s and the phase lag time of temperature gradient to be 0.043–0.056 s for processed meat. Xu et al. [7] presented a system discussion on the application of the DPL model in the biothermomechanical behavior of skin tissue. A more rational prediction of temperature distribution is always needed in the development of hyperthermia.

For more realistic predictions, this work would employ the DPL model of bioheat transfer to analyze the thermal behavior in tissue, which was stratified into skin, fat, and muscle. The DPL equation of bioheat transfer was always developed with the first-order Taylor series expansion of DPL model. For a more general form, this paper develops a modified DPL equation of bioheat transfer based on the second-order Taylor expansion. For convenience of statement, this paper calls the former first-order DPL equation and another second-order DPL equation. The second-order DPL equation is a fourth-order partial differential equation. Due to the difference in physiological and thermal properties, the boundary conditions at the interface between two adjacent layers become complicated. There are mathematical difficulties in dealing with such a problem. The hybrid numerical scheme [8] based on the Laplace transform and the modified discretization technique is extended to solve the present problem. The deviations of the results from the bioheat transfer equations based on Pennes' model, thermal wave model, and dual-phase-lag model are presented and discussed.

2. Problem Formulation

In order to solve the paradox which occurred in the classical heat flux model and to consider the effect of microstructural interactions, the DPL model was suggested [5] with

\overset{⇀}{q} (t + τ_{q}) = - k \nabla T (t + τ_{T}),

(1)

where T is the temperature, k is the heat conductivity, q is the heat flux, and t is the time. τ_q and τ_T can be interpreted as phase lags arising from “thermal inertia” and “microstructural interaction,” respectively [6]. τ_q means the phase lag of the heat flux and τ_T means the phase lag of the temperature gradient. The heat flux precedes the temperature gradient for τ_q < τ_T. The temperature gradient precedes the heat flux for τ_q > τ_T. Equation (1) reduces to the thermal wave model by setting τ_T = 0 and reduces to Fourier's heat equation by setting τ_q = τ_T.

Equation (1) is, usually, developed in the first-order Taylor series expansion. As [7] did, this paper would develop (1) in the second-order Taylor series expansion for a more general form. Thus, it is rewritten as

(1 + τ_{q} \frac{\partial}{\partial t} + \frac{τ_{q}^{2}}{2} \frac{\partial^{2}}{\partial t^{2}}) \overset{⇀}{q} = - (1 + τ_{T} \frac{\partial}{\partial t} + \frac{τ_{T}^{2}}{2} \frac{\partial^{2}}{\partial t^{2}}) k \nabla T .

(2)

In a local energy balance, the energy conservation equation of bioheat transfer is described as

ρ C \frac{\partial T}{\partial t} = - \nabla \cdot \overset{⇀}{q} + W_{b} C_{b} (T_{b} - T) + q_{m} + q_{r},

(3)

where t is time. ρ, c, and T denote density, specific heat, and temperature of tissue. c_b and w_b are, respectively, the specific heat and perfusion rate of blood. q_m is the metabolic heat generation and q_r is the heat source for spatial heating. T_b is the arterial temperature.

Substituting (2) into the energy conservation equation (3) leads to the second-order DPL equation of bioheat transfer with constant physiological parameters as follows:

(1 + τ_{T} \frac{\partial}{\partial t} + \frac{τ_{T}^{2}}{2} \frac{\partial^{2}}{\partial t^{2}}) k \nabla^{2} T = (1 + τ_{q} \frac{\partial}{\partial t} + \frac{τ_{q}^{2}}{2} \frac{\partial^{2}}{\partial t^{2}}) \times [ρ c \frac{\partial T}{\partial t} - w_{b} c_{b} (T_{b} - T) - q_{m} - q_{r}] .

(4)

As τ_T = 0, (4) undergoes the hyperbolic heat transfer equation with negating second-order Taylor series expansion for τ_q. Deleting the second-order Taylor series expansions for τ_q and τ_T leads to the first-order DPL equation of bioheat transfer. On the other hand, (4) would reduce to the Pennes equation for τ_q = τ_T = 0.

Consider a broad laser beam with uniform irradiance (φ_in) is applied normally to the human body at time t = 0⁺. When the spot size of the broad laser beam is much larger than the thickness of the thermally affected zone for the time period of interest, a 1D model would be sufficient for analyzing the thermal response of the heated medium [9, 10]. The tissue of the human body was often described as a three-layer model [10]. Figure 1 shows a one-dimensional analytical human body model that consists of skin, fat, and muscle.

Figure 1:

Three-layer human body model.

The 1D form of (4) with constant thermal parameters is written as

(1 + τ_{T} \frac{\partial}{\partial t} + \frac{τ_{T}^{2}}{2} \frac{\partial^{2}}{\partial t^{2}}) k \frac{\partial^{2} T}{\partial x^{2}} = (1 + τ_{q} \frac{\partial}{\partial t} + \frac{τ_{q}^{2}}{2} \frac{\partial^{2}}{\partial t^{2}}) \times [ρ c \frac{\partial T}{\partial t} - w_{b} c_{b} (T_{b} - T) - q_{m} - q_{r}] .

(5)

When laser irradiation is highly absorbed in the tissue as for some UV and IR wavelengths, scattering may be neglected. In this situation, the laser light can be absorbed within a very small depth from the irradiated boundary surface, laser heating can be reasonably considered as a boundary surface heat flux, and the terms involving q_r in (4) and (5) are neglected. Laser heating serves as a spatially varied body heat source that heats up the tissue for scattering. Under broad beam irradiance, the light distribution in scattering tissues is described with two exponentials as [11]

φ (x) = φ_{in} [C_{1} \exp (- \frac{d_{1} x}{δ}) - C_{2} \exp (- \frac{d_{2} x}{δ})],

(6)

where the coefficients C₁, d₁, C₂, and d₂ depend on the diffuse reflectance and can be obtained from the Monte Carlo solutions. δ is the effective optical penetration depth, which is defined based on the diffusion theory as

δ = \frac{1}{\sqrt{3 μ_{a} [μ_{a} + μ_{s} (1 - g)]}},

(7)

where μ_a is the absorption coefficient and μ_s is the scattering coefficient. The reduced scattering coefficient μ_s′ is introduced as μ_s′ = μ_s(1 – g), in which g is the anisotropy factor.

This study will treat the body tissue in two ways: highly absorbing tissue and strongly scattering tissue. The forced surface cooling is prosecuted in heating process. For circulation effect of blood, the temperature deep in tissue is often regarded to keep at the core temperature. For highly absorbing tissue, the boundary conditions at the ends for t > 0 are described as

\begin{matrix} q (0, t) = h (T_{f} - T) + φ_{in} (1 - R_{d}), \\ T (L_{3}, t) = T_{b}, \end{matrix}

(8)

where h is the heat convection coefficient between the cooling medium and the skin surface, T_f is the temperature of the cooling medium, and the location x = L₃ is defined as the body core. R_d is the diffuse reflectance of light at the irradiated surface.

For strongly scattering tissue, the boundary condition at x = 0 becomes

q (0, t) = h (T_{f} - T)

(9)

and the laser volumetric heat source q_r is described as

q_{r} = μ_{a} φ (x) .

(10)

This study assumes the laser irradiance φ_in is deposited as shown in Figure 2.

Figure 2:

Laser pulse sequences.

The boundary conditions at the interfaces of two adjacent layers are obtained from the assumption that temperature and heat flux are continuous.

The initial conditions are regarded as

T (x, 0) = T_{b}, \frac{\partial T (x, 0)}{\partial t} = \frac{\partial^{2} T (x, 0)}{\partial t^{2}} = 0 .

(11)

3. Analytical Method

First, the Laplace transform method is employed to map the transient problem into steady one. Equation (5) and the boundary conditions can be written in the transform domain based on the initial conditions (11) as

\begin{matrix} \frac{d^{2} {\tilde{T}}_{j}}{d x^{2}} - λ_{j}^{2} {\tilde{T}}_{j} = - g_{j} for j = 1, 2, 3, \\ \frac{d \tilde{T} (0)}{d x} = \frac{h}{K_{1}} (\frac{T_{f}}{s} - \tilde{T}), \\ \tilde{T} (L_{3}) = T_{b} \frac{1}{s} \end{matrix}

(12)

and at the interfaces of two adjacent layers

\begin{matrix} {\tilde{T}}_{j} = {\tilde{T}}_{j + 1} for j = 1, 2, \\ \frac{d {\tilde{T}}_{j}}{d x} = \frac{K_{j + 1}}{K_{j}} \frac{d {\tilde{T}}_{j + 1}}{d x} for j = 1, 2, \end{matrix}

(13)

where

\begin{matrix} λ_{j}^{2} = \frac{1}{K_{j}} [ρ_{j} c_{j} s + w_{b} c_{b}] for j = 1, 2, 3, \\ K_{j} = \frac{k_{j} (1 + τ_{T, j} s + τ_{T, j}^{2} s^{2} / 2)}{1 + τ_{q, j} s + τ_{q, j}^{2} s^{2} / 2} for j = 1, 2, 3, \\ g_{j} = \frac{1}{k_{j} (1 + τ_{T, j} s + τ_{T, j}^{2} s^{2} / 2)} (\frac{q_{m}}{s} + {\tilde{q}}_{r}) for j = 1, 2, 3 . \end{matrix}

(14)

The function T is written as $\tilde{T}$ in the Laplace domain. s is the Laplace transform parameter with respect to t.

A modified discretization scheme is applied to numerically solve the present problem. Please refer to [8] for the details of the present numerical scheme.

4. Results and Discussion

The optical and thermal physical properties used to perform all present computations in the present paper are shown in Tables 1 and 2, respectively. The relaxation times and the arterial blood temperature are specified as τ_q = 15.5 s, τ_T = 10 s, and T_b = 37°C. The incident laser irradiance is considered as q₀ = 2 W/cm². For highly absorbing tissues, the diffuse reflectance of light at the irradiated surface R_d is 0.05 [9].

Table 1:

The optical properties of tissues [12].

Parameters	Skin	Fat	Muscle

Absorption coefficient (mm⁻¹)	0.033	0.013	0.121
Reduced scattering coefficient (mm⁻¹)	2.73	1.26	0.89
Effective optical penetration depth (mm⁻¹)	1.91	4.49	1.65

Table 2:

The thermal physical properties of tissues.

Parameters	Blood [7]	Skin [7]	Fat [7]	Muscle [9]

Specific heat (J/kg °C)	3770.0	3300.0	2700.0	4187. 0
Density (kg/m³)	1060.0	1116.0	971.0	1000.0
Thermal conductivity (W/m °C)	—	0.445	0.185	0.628
Metabolic heat generation (W/m³)	—	368.1	368.3	1190.0
Thickness (m)	—	0.0012	0.0048	0.054
Perfuse rate (m³/s/m³)		0.05	0.0	0.00187

Figure 3 shows the variations of temperature predicted with various equations for h = 0 and h = 150 W/m² s at the irradiated surface of the highly absorbing tissue. It is found that the thermal response is almost simultaneous for various equations during heating period. In the thermal wave model, the effect of τ_q = 15.5 s induced the higher temperature at the early times. And then, the reflection of thermal wave created at the interfaces between two adjacent layers makes several subthermal waves appear as shown in Figure 3. Correspondingly, in the other models, the effects of thermal diffusion make the temperatures lower and depress the reflection phenomenon happening at the interfaces. As h = 150 W/m² s, the temperatures are reduced for the surface cooling function.

Figure 3:

Variations of temperature at the irradiated surface of the highly absorbing tissue.

The comparison of temperature variation at the interface of fat and muscle (x = 6 mm) in highly absorbing tissue and strongly scattering tissue is offered in Figure 4. It is found that the starting times of thermal response are almost simultaneous for various heat conduction models in the strongly scattering tissue, as shown in Figure 4(b), although the Pennes equation depicts an infinitely fast propagation of thermal signal. Because the heating source was regarded as the volumetric heat source, the distance from the heating source to the thermal response location is zero. As a result, no thermal transport delay is created. However, in the highly absorbing tissue, the heating source is the surface heat flux. Due to the difference in the velocity of thermal transport, the starting times of thermal response are obviously different at the location x = 6 mm for various heat conduction models, as shown in Figure 4(a).

Figure 4:

Comparison of temperature variation at the interface of fat and muscle (x = 6 mm) in (a) highly absorbing tissue and (b) strongly scattering tissue.

The effect of the heat convection on skin surface with h = 150 W/m² s on the variations of temperature at the location x = 6 mm is shown in Figure 5. The problem would become a nonlinear problem as h ≠ 0. Solving this problem becomes more difficult. It is observed from Figure 4(b) that the temperatures at the location x = 6 mm are slightly reduced due to the convection effect of h = 150 W/m² s. The convection cooling function of skin surface can work at h = 150 W/m² s for such laser heating. In other words, the quality of treatment can be controlled using surface cooling.

Figure 5:

Variations of temperature predicted with various equations at the location x = 6 mm for h = 150 W/m² s in the strongly scattering tissue.

Figure 6 shows the temperature distribution in the strongly scattering tissue at t = 25 s. Two vertical lines are used to show the position of the layer interfaces. From Table 1, it is known that the absorption coefficient of muscle is higher, so the muscle layer is able to efficiently absorb the laser energy. As a result, in the muscle layer, the affected region is larger, and the temperatures are also higher. Because the non-Fourier effect has been weakened at t = 25 s, the temperature distributions are nearly coincident.

Figure 6:

Temperature distribution in the strongly scattering tissue for h = 0 at t = 25 s.

5. Conclusions

This study treats the body tissue in two ways: highly absorbing tissue and strongly scattering tissue. The comparisons of the results from the bioheat transfer equations based on Pennes' model, thermal wave model, and dual-phase-lag model were done. It is found that the starting times of thermal response are almost simultaneous for various heat conduction models in the strongly scattering tissue, although the Pennes equation depicts an infinitely fast propagation of thermal signal. Because the heating source was regarded as the volumetric heat source, the distance from the heating source to the thermal response location is zero. As a result, no thermal transport delay is created. However, in the highly absorbing tissue, the heating source is the surface heat flux. Due to the difference in the velocity of thermal transport, the starting times of thermal response would be different at the location inner medium for various heat conduction models. The higher temperature is induced in the thermal wave model at the early times for the effect of τ_q. The comparison between the predicted results with the first-order and second-order DPL equations further depicts that the effect of microstructural interactions is obvious only in a short time period. The convection cooling function of skin surface can work for such continued pulse heating.

Footnotes

Acknowledgments

This work is funded by the National Science Council under Grants nos. NSC-100-2221-E-269-015 and NSC-102-2221-E-269-013.

References

Liu

Chen

, and Xu

L. X.

, “New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 4, pp. 420–428, 1999.

Liu

K.-C.

, “Thermal propagation analysis for living tissue with surface heating,” International Journal of Thermal Sciences, vol. 47, no. 5, pp. 507–513, 2008.

Özen

Ş.

Helhel

, and Cerezci

, “Heat analysis of biological tissue exposed to microwave by using thermal wave model of bio-heat transfer (TWMBT),” Burns, vol. 34, no. 1, pp. 45–49, 2008.

Shih

T.-C.

Kou

H.-S.

Liauh

C.-T.

, and Lin

W.-L.

, “The impact of thermal wave characteristics on thermal dose distribution during thermal therapy: a numerical study,” Medical Physics, vol. 32, no. 9, pp. 3029–3036, 2005.

Tzou

D. Y.

, Macro- to Microscale Heat Transfer: The Lagging Behavior, Taylor & Francis, Washington, DC, USA, 1996.

Antaki

P. J.

, “New interpretation of non-Fourier heat conduction in processed meat,” Journal of Heat Transfer, vol. 127, no. 2, pp. 189–193, 2005.

Seffen

K. A.

, and Lu

T. J.

, “Non-Fourier analysis of skin biothermomechanics,” International Journal of Heat and Mass Transfer, vol. 51, no. 9–10, pp. 2237–2259, 2008.

Liu

K.-C.

and Cheng

P-J.

, “Finite propagation of heat transfer in a multilayer tissue,” Journal of Thermophysics and Heat Transfer, vol. 22, no. 4, pp. 775–782, 2008.

Zhou

Chen

J. K.

, and Zhang

, “Dual-phase lag effects on thermal damage to biological tissues caused by laser irradiation,” Computers in Biology and Medicine, vol. 39, no. 3, pp. 286–293, 2009.

10.

Kanezaki

Hirata

Watanabe

, and Shirai

, “Parameter variation effects on temperature elevation in a steady-state, onedimensional thermal model for millimeter wave exposure of one-and three-layer human tissue,” Physics in Medicine and Biology, vol. 55, no. 16, pp. 4647–4659, 2010.

11.

Gardner

C. M.

Jacques

S. L.

, and Welch

A. J.

, “Light transport in tissue: accurate, heuristic expressions for one dimensional fluence rate and escape function based upon monte carlo simulations,” Lasers in Surgery and Medicine, pp. 129–138, 1996.

12.

Simpsonyz

C. R.

Kohlyx

Essenpreisk

, and Copey

, “Near-infrared optical properties of ex Pvivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Physics in Medicine and Biology, vol. 43, pp. 2465–2478, 1998.