Analytical Solutions to 3-D Bioheat Transfer Problems with or without Phase Change

1. Introduction

Theoretical analysis on the bioheat transfer process has been an extremely important issue in a wide variety of bioengineering situations such as cancer hyperthermia, burn injury evaluation, brain hypothermia, disease diagnostics, thermal comfort analysis, cryosurgery and cryopreservation etc. In this chapter, the theoretical strategies towards exactly solving the three-dimensional (3-D) bioheat transfer problems for both cases with and without phase change were systematically illustrated based on the authors’ previous works. Typical closed form analytical solutions to the hyperthermia bioheat transfer problems with space or transient heating on skin surface or inside biological bodies were summarized. In addition, exact solutions to the 3-D temperature transients of tissues under various phase change processes such as cryopreservation of biomaterials or cryosurgery of living tissues subject to freezing by a single or multiple cryoprobes were also outlined. Such solution is comprehensive enough by taking full account of many different factors such as generalized initial and boundary conditions, blood perfusion heat transfer, volumetric heating of hyperthermia apparatus or heat sink of cryoprobes etc. For illustrating the applications of the present methods, part of the solutions were adopted to analyze the selected bioheat transfer problems. The versatility of these theoretical approaches to tackle more complex issues was also discussed. The obtained solutions are expected to serve as the basic foundation for theoretically analyzing bioheat transfer problems.

2. Motivations of analytical solutions to bioheat transfer problem

Analytical solutions to bioheat transfer problems are very important in a wide variety of biomedical applications [1]. Especially, understanding the heat transfer in biological tissues involving either raising or lowering of temperature is a necessity for many clinical practices such as tumor hyperthermia [2], burn injury evaluation [3, 4], brain hypothermia resuscitation [5], disease thermal diagnostics [6], thermal comfort analysis [7], cryosurgery planning [8, 9], and cryopreservation programming [10]. The bioheat transfer problems involved in the above applications can generally be divided into two categories: with and without phase change. In this chapter, the phase change especially denotes the solid-liquid phase transition of biological hydrated tissues. The cases without phase change usually include tumor hyperthermia, burn injury evaluation, brain hypothermia resuscitation, disease diagnostics, and thermal comfort analysis, while the cases with phase change include cryosurgery and cryopreservation.

To guarantee optimal clinical outputs for such applications, it is essential to predict in advance the transient temperature distribution of the target tissues. For example, in a tumor hyperthermia process, the primary objective is to raise the temperature of the diseased tissue to a therapeutic value, typically above 43ºC, and then thermally destroy it [11]. Temperature prediction would be used to find an optimum way either to induce or prevent such thermal damage to the target tissues. In contrast to the principle of hyperthermia, cryosurgery realizes its clinical purpose of controlled tissue destruction through deep freezing and thawing [12]. Applications of this treatment are quite wide in clinics owning to its outstanding virtues such as quick, clean, relatively painless, good homeostasis, and minimal scaring. An accurate understanding of the extent of the irregular shape of the frozen region, the direction of ice growth, and the temperature distribution within the ice balls during the freezing process is a basic requirement for the successful operation of a cryosurgery. Therefore, solving the bioheat transfer problems involved is very important for both hyperthermia and cryosurgery. Moreover, in thermal diagnostics, thermal comfort analysis, brain hypothermia resuscitation, and burn injury evaluation, similar bioheat transfer problems are also often encountered [13].

It is commonly accepted that mathematical model is the basis for solving many practical problems. Because modeling bioheat transfer is of the utmost importance in many biomedical applications such as proper device or heating/cooling protocol design, a number of bioheat transfer equations for living tissue have been proposed since the landmark work by Pennes published in 1948 [14], in which the perfusion heat source/sink was introduced. Until now, the classical Pennes equation is also commonly accepted as the best practical approach for modeling bioheat transfer in view of its simplicity and excellent validity [15]. This is because most of the other models either still lack sound experimental grounding or just appear too complex for mathematical solution. Although the real anatomical geometry of a biological body can be incorporated, the Pennes equation remains the most useful model for characterizing the heat transport process in most biomedical applications. For brevity, here only cases for space-dependent thermal properties will be mainly discussed. Then a generalized form of the Pennes equation for this purpose can be written as:

$$\begin{array}{l}\rho c\frac{\partial T(X,t)}{\partial t}=\nabla \cdot k(X)\nabla \left[T(X,t)\right]+{\rho }_b{c}_{\text{b}}{\omega }_{\text{b}}(X)\left[T_{\text{a}}-T(X,t)\right]\\ +Q_{\text{m}}(X,t)+Q_{\text{r}}(X,t)\begin{array}{cc},& X\in \Omega \end{array}\end{array}$$E1

where, $$\rho$$, c are the density and the specific heat of tissue, respectively; $${\rho }_b$$and $$c_{\text{b}}$$ denote the density and the specific heat of blood, respectively; $X$contains the Cartesian coordinates$x$, $y$and$z$; $\Omega$denotes the analyzed spatial domain; $$k(X)$$is the space dependent thermal conductivity; and $${\omega }_{\text{b}}(X)$$ is the space dependent blood perfusion. The value of blood perfusion represents the blood flow rate per unit tissue volume and is mainly from microcirculation including the capillary network plus small arterioles and venules. $$T_{\text{a}}$$is the blood temperature in the arteries supplying the tissue and is often treated as a constant at 37°C; $T(X,t)$is the tissue temperature; $Q_{\text{m}}(X,t)$is the metabolic heat generation; and $Q_{\text{r}}(X,t)$ the distributed volumetric heat source due to externally applied spatial heating.

From the historical viewpoint, we can find that the development of the bioheat transfer’s art and science can be termed as one to modify and improve the Pennes model [15]. Among the many efforts, the blood perfusion term in the Pennes equation has been substantially studied which led to several conceptually innovative bioheat transfer models such as Wulff’s continuum model [16], Chen-Holmes model addressing both the flow and perfusion properties of blood [17], and the Weinbaum-Jiji three-layer model to characterize the heat transfer in the peripheral tissues [18]. The bioheat transfer equation and its extended forms can be directly used to characterize the thermal process of the biological bodies subject to various external or interior factors such as convective interaction with a heated or cooled fluid, radiation by fire or laser, contact with a heating or freezing apparatus, electromagnetic effect, or a combination among them. Such issues can be treated using different boundary conditions as well as spatial heating or freezing patterns. Generally, the geometric shape, dimensions, thermal properties and physiological characteristics for tissues, as well as the arterial blood temperature, can be used as the input to the Pennes equation for a parametric study. According to a specific need in clinics, the bioheat transfer model can even be modified by taking more factors into concern [19]. Traditionally, for solving bioheat transfer problems, people relied too heavy on numerical approaches such as finite difference method, finite element method, and boundary element method etc. Numerical simulation is necessary when the analytical solutions are not available. But if both analytical and numerical solutions can be obtained for the same issue, the analytical one is often preferred. Except for its simplicity being used to compile computer codes, the analytical solution is very attractive since its efficiency depends weakly on the dimensions of the problem, in contrast to the numerical methods. For analytical method, solution at a desired point can be performed independently from that of the other points within the domain, which can be an asset when temperatures are needed at only some isolated sites or times. But for most of the conventional numerical methods (except Monte Carlo simulation), the temperatures at all mesh points must be simultaneously computed even when only the temperatures at a single point are needed [20]. In this sense, the analytical solution will save computational time greatly, which is valuable in clinical practices.

Based on the above considerations, we aimed in this chapter to present several typical closed form analytical solutions to bioheat transfer problems with or without phase change, in which relatively complex boundary or heating/cooling conditions, and existence of discrete large blood vessel were included. Derivation of the solutions was mainly based on the Green’s function method, which is beneficial for dealing with the non-homogeneous problems with spatial or transient heating source and initial temperature distribution, as well as complex cooling or boundary conditions. For generalized and practical purpose, complex bioheat transfer problems encountered in several typical clinical applications as well as basic studies such as tumor hyperthermia, cryosurgery, cryopreservation, and interpretation of physiological phenomena etc. will be especially addressed.

3. Bioheat transfer problems without phase change

3.1. Generalized analytical solutions to 3-D bioheat transfer problems

Derivation of the solutions was based on the Green’s function method, since the Green’s function obtained for the differential equation is independent of the source term. Therefore it can be flexibly used to calculate the temperature distribution for various spatial or temporal source profiles. Furthermore, the Green’s function method is capable of dealing with the transient or space-dependent boundary conditions. Up to now, quite a few studies have applied the Green’s function method to solve the bioheat transfer problems [21-25]. However, in most of the existing analytical studies, the available solutions to the bioheat transfer problem are for the cases with one dimensional geometry, steady state, infinite domain, constant heating, or heat conduction equations not considering blood perfusion, which may not be practical for some real bio-thermal situations. In this section, the generalized analytical solutions, which have incorporated relatively complex situations such as the 3-D tissue domain, the transient or space-dependent boundary conditions, and volumetric heating, were especially addressed. Such solutions are expected to be very useful in a variety of bio-thermal practices. The 3-D computational domain with widths s₁ = s₂ = 0.08m and height L = 0.03m was depicted as the shadowed region in Fig. 1, where s₁ and s₂ were widths of the tissue domain to be analyzed in y and z directions, respectively; the skin surface was defined at x=0 while the body core at x=L.

Figure 1.

Calculation geometry for 3-D case [13]

For brief, only 3-D case with constant thermal parameters will be particularly studied, which is a good approximation when no phase change occurred in tissue. The corresponding 3-D Pennes equation can be derived from Equation (1) as:

$$\rho c\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial x^2}+k\frac{{\partial }^{2}T}{\partial y^2}+k\frac{{\partial }^{2}T}{\partial z^2}+{\omega }_b{\rho }_b{c}_b\left(T_a-T\right)+Q_m+Q_r\left(x,y,z;t\right)$$E2

The generalized boundary conditions (BCs) often encountered in a practical clinical situation can be written as:

$$-k\frac{\partial T}{\partial x}=f_1\left(y,z;t\right),x=0$$E3

or

$-k\frac{\partial T}{\partial x}=h_f\left[f_2(y,z;t)-T\right],x=0$E4

where, $$f_1\left(y,z;t\right)$$is the time-dependent surface heat flux, $$f_2\left(y,z;t\right)$$is the time-dependent temperature of the cooling medium, and $h_f$ is the heat convection coefficient between the medium and the skin surface. In this chapter, Equation (3) was named the second BC and Equation (4) the third BC.

The body core temperature was regarded as a constant ($T_c$) on considering that the biological body tends to keep its core temperature to be stable, i.e.

$$T=T_c\begin{array}{cc},& x=L\end{array}$$E5

The BCs at y and z directions can be expressed as

$$-k\frac{\partial T}{\partial y}=0,y=0$$E6

$$-k\frac{\partial T}{\partial y}=0,y=s_1$$E7

$$-k\frac{\partial T}{\partial z}=0,z=0$$E8

$$-k\frac{\partial T}{\partial z}=0,z=s_2$$E9

The reason for adopting the adiabatic conditions in the two ends of the $y$ and $z$ directions is from the consideration that at the positions far from the beam center of the heat deposition apparatus, the temperatures there were almost not affected by the external heating, which generally has a strong decay in $y$ and $z$ directions. However, it should be mentioned that more generalized boundary conditions in $y$ and $z$ directions can also be dealt with by the present approach, but they were not listed here for brevity.

The initial temperature is

$$T\left(x,y,z;0\right)=T_0\left(x,y,z\right),t=0$$E10

where, $$T_0\left(x,y,z\right)$$can be approximated by the 1-D solution, representing the initial temperature field for the basal state of biological bodies, which can be obtained through solving the following equation sets:

$\{\begin{array}{l}k\frac{d^2{T}_0(x)}{d{x}^2}+{\omega }_b{\rho }_b{c}_b\left[T_a-T_0(x)\right]+Q_m=0\\ T_0(x)=T_c\begin{array}{ccc},& & x=L\end{array}\\ -k\frac{d{T}_0(x)}{dx}=h_0\left[T_f-T_0(x)\right]\begin{array}{ccc},& & x=0\end{array}\end{array}$E11

where, $h_0$is the apparent heat convection coefficient between the skin surface and the surrounding air under physiologically basal state and is an overall contribution from natural convection and radiation, and $T_f$ the surrounding air temperature.

The solution to Equation (11) is:

$$\begin{array}{l}{T}_0(x)=T_a+\frac{Q_m}{{\omega }_b{\rho }_b{c}_b}+\frac{\left(T_c-T_a-\frac{Q_m}{{\omega }_b{\rho }_b{c}_b}\right)\cdot \left[\sqrt{A}ch\left(\sqrt{A}x\right)+\frac{h_0}{k}sh\left(\sqrt{A}x\right)\right]}{\sqrt{A}ch\left(\sqrt{A}L\right)+\frac{h_0}{k}sh\left(\sqrt{A}L\right)}\\ \begin{array}{ccc}& & \end{array}+\frac{\frac{h_0}{k}\left(T_f-T_a-\frac{Q_m}{{\omega }_b{\rho }_b{c}_b}\right)\cdot sh\left[\sqrt{A}\left(L-x\right)\right]}{\sqrt{A}ch\left(\sqrt{A}L\right)+\frac{h_0}{k}sh\left(\sqrt{A}L\right)}\end{array}$$E12

where,$A={\omega }_b{\rho }_b{c}_b/k$. Through using the following transformation [13]

$$T\left(x,y,z;t\right)=T_0\left(x,y,z\right)+W\left(x,y,z;t\right)\cdot \exp \left(-\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}t\right)$$E13

Equation (2) was transformed to the following form:

$$\frac{\partial W}{\partial t}=\alpha \frac{{\partial }^{2}W}{\partial x^2}+\alpha \frac{{\partial }^{2}W}{\partial y^2}+\alpha \frac{{\partial }^{2}W}{\partial z^2}+\frac{Q_r\left(x,y,z;t\right)}{\rho c}\cdot \exp \left(\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}t\right)$$E14

where, $$\alpha =k/\rho c$$is the thermal diffusivity of tissue. The corresponding boundary and initial conditions are:

$$-k\frac{\partial W}{\partial x}=g_1\left(y,z;t\right),x=0$$E15

$-k\frac{\partial W}{\partial x}=h_f\left[g_2(y,z;t)-W\right],x=0$E16

$$W=0,x=L$$E17

$$-k\frac{\partial W}{\partial y}=0,y=0$$E18

$$-k\frac{\partial W}{\partial y}=0,y=s_1$$E19

$$-k\frac{\partial W}{\partial z}=0,z=0$$E20

$$-k\frac{\partial W}{\partial z}=0,z=s_2$$E21

$$W\left(x,y,z;0\right)=0,t=0$$E22

where,

$$g_1\left(y,z;t\right)=\left[{k\frac{\partial T_0}{\partial x}|}_{x=0}+f_1\left(y,z;t\right)\right]\cdot \exp \left(\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}t\right)$$E23

$$g_2\left(y,z;t\right)=\left[{\frac{k}{h_f}\frac{\partial T_0}{\partial x}|}_{x=0}-{T_0|}_{x=0}+f_2\left(y,z;t\right)\right]\cdot \exp \left(\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}t\right)$$E24

Using Green function method, $$W\left(x,y,z;t\right)$$can be solved from the combined Equations (14-24). The Green’s functions $G_1$ and $G_2$ to the second and third BCs can finally be obtained:

$\begin{array}{l}{G}_1(x,y,z,t;\xi ,\vartheta ,\varsigma ,\tau )=\\ \frac{2}{L{s}_1{s}_2}{\displaystyle \sum _{l=1}^{\infty }{\displaystyle \sum _{m=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{R}_1{R}_2{e}^{-\alpha ({\beta }_l^2+{\gamma }_m^2+{\psi }_n^2)(t-\tau )}\cos {\beta }_{l}x\cos {\beta }_l\xi \cos {\gamma }_{m}y\cos {\gamma }_m\vartheta \cos {\psi }_{n}z\cos {\psi }_n\varsigma }}}\end{array}$E25

$\begin{array}{l}{G}_2(x,y,z,t;\xi ,\vartheta ,\varsigma ,\tau )=\\ {\displaystyle \sum _{p=1}^{\infty }{\displaystyle \sum _{m=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{c}_{pmn}{e}^{-\alpha ({\beta }_p^2+{\gamma }_m^2+{\psi }_n^2)(t-\tau )}}}}\sin {\beta }_p(L-x)\sin {\beta }_p(L-\xi )\cos {\gamma }_{m}y\cos {\gamma }_m\vartheta \cos {\psi }_{n}z\cos {\psi }_n\varsigma \end{array}$E26

where,

$${\beta }_l=\frac{2l-1}{2L}\pi ,l=1,2,3,\mathrm{...}$$E27

$${\gamma }_m=\frac{m\pi }{s_1},m=0,1,2,\mathrm{...}$$E28

$${\psi }_n=\frac{n\pi }{s_2},n=0,1,2,\mathrm{...}$$E29

$R_1=\{\begin{array}{l}1\begin{array}{cc},& m=0\end{array}\\ 2\begin{array}{cc},& m=1,2,3,\mathrm{...}\end{array}\end{array}$E30

$R_2=\{\begin{array}{l}1\begin{array}{cc},& n=0\end{array}\\ 2\begin{array}{cc},& n=1,2,3,\mathrm{...}\end{array}\end{array}$E31

$c_{pmn}=\frac{2\left[{\beta }_p^2+{\left(h_f/k\right)}^2\right]R_1{R}_2}{\left\{L\left[{\beta }_p^2+{\left(h_f/k\right)}^2\right]+h_f/k\right\}{s}_1{s}_2}$E32

The Eigen-values ${\beta }_p$ are the positive roots of the following equation

$${\beta }_p\cdot \cot ({\beta }_{p}L)=-h_f/k$$E33

Then, the solution of Equation (14) can be easily obtained. For the second BC at the skin surface, one has

$$\begin{array}{l}W\left(x,y,z;t\right)=\frac{\alpha }{k}{\displaystyle {\int }_0^{t}d\tau {\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}{G}_1{\left(x,y,z,t;\xi ,\vartheta ,\varsigma ,\tau \right)|}_{\xi =0}{g}_1\left(\vartheta ,\varsigma ;\tau \right)d\vartheta d\varsigma }}+}\\ {\displaystyle {\int }_0^{t}d\tau {\displaystyle {\int }_0^L{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}{G}_1\left(x,y,z,t;\xi ,\vartheta ,\varsigma ,\tau \right)}}}}\cdot \frac{Q_r\left(\xi ,\vartheta ,\varsigma ;\tau \right)}{\rho c}\exp \left(\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}\tau \right)d\xi d\vartheta d\varsigma \end{array}$$E34

For the third BC, the solution is

$$\begin{array}{l}W\left(x,y,z;t\right)=\frac{\alpha }{k}{\displaystyle {\int }_0^{t}d\tau {\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}{G}_2{\left(x,y,z,t;\xi ,\vartheta ,\varsigma ,\tau \right)|}_{\xi =0}{h}_f\cdot g_2\left(\vartheta ,\varsigma ;\tau \right)d\vartheta d\varsigma }}+}\\ {\displaystyle {\int }_0^{t}d\tau {\displaystyle {\int }_0^L{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}{G}_2\left(x,y,z,t;\xi ,\vartheta ,\varsigma ,\tau \right)}}}}\cdot \frac{Q_r\left(\xi ,\vartheta ,\varsigma ;\tau \right)}{\rho c}\exp \left(\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}\tau \right)d\xi d\vartheta d\varsigma \end{array}$$E35

Then the tissue temperature field can be constructed as:

$$T\left(x,y,z;t\right)=T_0\left(x,y,z\right)+W\left(x,y,z;t\right)\cdot \exp \left(-\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}t\right)$$E36

Clearly, the above method can also be extended to solve some other three-dimensional problems such as in spherical and cylindrical coordinates. But they will not be listed here for brevity. To illustrate the application of the above analytical solutions, a selective 3-D hyperthermia problem with point heating sources was particularly studied as an example. Accordingly, the temperature distribution of tissue subject to the point heating in volume was analytically solved. Practical examples for the point heating can be found in clinics where heat was deposited though inserting a conducting heating probe in the deep tumor site. Previously, such problems received relatively few attentions in compared with other heating patterns. Here, the point-heating source to be studied can be expressed as:

$Q_r(x,y,z,t)=P_1(t)\delta (x-x_0)\delta (y-y_0)\delta (z-z_0)$E37

where, $P_1(t)$is the point-heating power, $\delta$is the Dirac function, $(x_0,y_0,z_0)$is the position of the point-heating source.

The results were given in Fig. 2, which represent the temperature distribution in biological bodies heated by one and two-point sources, respectively. In calculations, the typical tissue properties were applied as given in Table 1. In Fig. 2(a), the single heating source was fixed at position (0.021m, 0.04m, 0.04m); in Fig. 2(b), the two point-heating sources were at (0.021m, 0.032m, 0.04m) and (0.021m, 0.048m, 0.04m), respectively. It makes clear that the maximum temperatures of the tissues occur at the positions of the point-heating sources. Further, one can still observe that the temperature for the tissues surrounding the point-heating sources can fairly be kept at a lower temperature on the whole. This is very beneficial for the hyperthermia operation since one can then selectively control the temperature level at the diseased tissue sites while the healthy tissues at the surrounding area will just stay below the safe threshold. This may be one of the most attractive features why the invasive heating probes are frequently used to thermally kill the tumor in the deep tissue, although they may cause mechanical injury. The above solutions are expected to be valuable for such hyperthermia treatment planning.

Figure 2.

Temperature distribution at cross-section $z=0.04m$ after 1200s’ heating

	Unit	Value
Air temperature (Tf)	°C	25
Artery blood temperature (Ta)	°C	37
Blood perfusion of tissue (ωb)	ml/s/ml	0.0005
Body core temperature (Tc)	°C	37
Density of tissue (ρ)	Kg/m^3	1000
Density of blood (ρb)	Kg/m^3	1000
Heat convection coefficient (h0)	W/m^2· °C	10
Heat convection coefficient (hf)	W/m^2· °C	100
Metabolic heat generation of tissue (Qm)	W/m^3	33800
Specific heat of tissue (c)	J/Kg· °C	4200
Specific heat of blood (cb)	J/Kg· °C	4200
Temperature of cooling medium (f2)	°C	15
Thermal conductivity of tissue (k)	W/m· °C	0.5

Table 1.

Typical thermophysical properties of soft biological tissues [13].

3.2. Analytical solutions to 3-D bioheat transfer involved in hyperthermia for prostate

Localized transurethral thermal therapy has been widely used as a non-surgical modality for treatment of benign prostatic hyperplasia [26]. One of the critical issues in clinical application is to effectively heat and cause coagulation necrosis in target tissue while simultaneously preserving the surrounding healthy tissue, especially the prostatic urethra and rectum. This requires administration of an optimal thermal dose which can induce the desired three dimensional tissue temperature distributions in the prostate during the therapy. In this section, the analytical approach to solving the transient 3-D temperature field was illustrated, which can be used to predict point-by-point tissue temperature mapping during the heating.

The transurethral microwave catheter (T3 catheter) was used as the heating apparatus in this section. Geometric presentation of the prostate with the inserted T3 catheter was shown in Fig. 3. It was modeled as a cylinder of 3.4cm in diameter and 3cm in length with constant temperature $$T_{\infty }$$ at the surface which is near the body core temperature. The catheter was represented by the inner cylinder. The induced volumetric heating in tissue is [26]:

$$Q_r\left(r,\theta ,z\right)=C_{t}Q\frac{\left[2\epsilon \left(r-s\cos \theta \right)+\left(N-2\right)\right]e^{-2\epsilon \left(r-s\cos \theta \right)}}{{\left(r-s\cos \theta \right)}^N}{e}^{-{\left(z-L/2\right)}^2/z_0^2}$$E38

where Q is the applied microwave power, C_t is a proportional constant, ε is the microwave attenuation constant in tissue, and z₀ is the critical axial decay length along the catheter while L is the total length of the prostate.

Practically, the microwave antenna is located with an offsets from the geometrical center to produce an asymmetric microwave field, which can prevent overheating the rectum. The chilled water at a given temperature flows between the antenna and the inner catheter wall.

The Pennes equation for the 3-D temperature field in the prostate can then be applied as:

$$k\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)+k\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\theta }^2}+k\frac{{\partial }^{2}T}{\partial z^2}+{\omega }_b{\rho }_b{c}_b\left(T_a-T\right)+Q_m+Q_r=\rho c\frac{\partial T}{\partial t}$$E39

Figure 3.

D configuration of the prostate under microwave heating [26]

To obtain an analytical solution, all these parameters were assumed to be uniform throughout the prostate and remained constant except for ω_b which varied with respect to the heat power. The cooling effect from the chilled water running inside the catheter was modeled by an overall convection coefficient h. The external boundary at the capsule was prescribed as the body core temperature$$T_{\infty }$$. Therefore, one has:

$$k\frac{\partial T}{\partial r}=h\left(T-T_f\right),r=R_0$$E40

$$T=T_{\infty },r=R_2$$E41

where, R₀ and R₂ are the urethra and prostate radius, respectively; T_f is the coolant temperature. Considering the Gaussian distribution of microwave power deposition along the z direction, adiabatic conditions can be used at two ends of the prostate:

$$k\frac{\partial T}{\partial z}=0,z=0,L$$E42

The initial temperature is

$$T\left(r,\theta ,z,0\right)=T_0\left(r,\theta ,z\right),t=0$$E43

Using transformation:

$$T=\Delta +T_{\infty }+(T_{\infty }-T_f)\frac{\ln (r/R_2)}{\ln (R_2/R_0)+k/(h{R}_0)}$$E44

One can rewrite the above equations (Equations (39-44) as:

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \Delta }{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\Delta }{\partial {\theta }_{}^2}+\frac{{\partial }^2\Delta }{\partial z^2}-\frac{{\omega }_b{\rho }_b{C}_b\Delta }{k}+\frac{Q*\left(r,\theta ,z\right)}{k}=\frac{1}{\alpha }\frac{\partial \Delta }{\partial t}$$E45

$$k\frac{\partial \Delta }{\partial r}=h\Delta ,r=R_0$$E46

$$\Delta =0,r=R_2$$E47

$$k\frac{\partial \Delta }{\partial z}=0,z=0,L$$E48

$$\Delta \left(r,\theta ,z,0\right)=F\left(r,\theta ,z\right),t=0$$E49

where,

$$Q*\left(r,\theta ,z\right)=Q_r\left(r,\theta ,z\right)+Q_m+{\omega }_b{\rho }_b{C}_b\left[T_a-T_{\infty }-\left(T_{\infty }-T_f\right)\frac{\ln \left(r/R_2\right)}{\ln \left(R_2/R_0\right)+k/\left(h{R}_0\right)}\right]$$E50

$$F\left(r,\theta ,z\right)=T_0\left(r,\theta ,z\right)-T_{\infty }-\left(T_{\infty }-T_f\right)\frac{\ln \left(r/R_2\right)}{\ln \left(R_2/R_0\right)+k/\left(h{R}_0\right)}$$E51

The Green’s function for the above equation sets can be obtained as:

$$\begin{array}{l}G\left(r,\theta ,z,t;\xi ,\eta ,\lambda ,\tau \right)=\\ -{\displaystyle \sum _{m=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \sum _{l=0}^{\infty }\left[J_m\left({\mu }_n^{*}r\right)+A_n{N}_m\left({\mu }_n^{*}r\right)\right]\cos \left[m\left(\theta -\eta \right)\right]\cos \left({\gamma }_{l}z\right)e^{-\alpha {\beta }_l^2\left(t-\tau \right)}\frac{\kappa \alpha F_{mn}\left(\xi \right)\cos \left({\gamma }_l\lambda \right)H\left(t-\tau \right)}{L}}}}\end{array}$$E52

where

$$\kappa =\{\begin{array}{l}1,l=0\\ 2,l=1,2,\mathrm{3...}\end{array}$$E53

$$A_n=-J_m\left({\mu }_n^{*}{R}_2\right)/N_m\left({\mu }_n^{*}{R}_2\right)$$E54

$${\gamma }_l=\frac{l\pi }{L},l=0,1,2,\mathrm{3...}$$E55

$${\beta }_l$$can be found by

$${\beta }_l^2={\gamma }_l^2+{\mu }_n^2,l=0,1,2,\mathrm{3...}$$E56

$${\mu }_n$$can be calculated from

$${\mu }_n=\sqrt{{\mu }_n^{*}{}^2+\frac{{\omega }_b{\rho }_b{C}_b}{k}}$$E57

$${\mu }_n^{*}$$are the positive roots of the following equation:

$$k\left[J_m\text{'}\left({\mu }_n^{*}{R}_0\right)+A_n{N}_m\text{'}\left({\mu }_n^{*}{R}_0\right)\right]=h\left[J_m\left({\mu }_n^{*}{R}_0\right)+A_n{N}_m\left({\mu }_n^{*}{R}_0\right)\right]$$E58

Here J_m and N_m are the Bessel functions of the first kind and the second kind, respectively.

$$F_{mn}=\{\begin{array}{l}\frac{\xi \left[J_0\left({\mu }_n^{*}\xi \right)+A_n{N}_0\left({\mu }_n^{*}\xi \right)\right]}{2\pi }/{\displaystyle {\int }_{R_0}^{R_2}r\left[J_0\left({\mu }_n^{*}r\right)+A_n{N}_0\left({\mu }_n^{*}r\right)\right]{}^{2}dr},m=0\\ \frac{\xi \left[J_m\left({\mu }_n^{*}\xi \right)+A_n{N}_m\left({\mu }_n^{*}\xi \right)\right]}{\pi }/{\displaystyle {\int }_{R_0}^{R_2}r\left[J_m\left({\mu }_n^{*}r\right)+A_n{N}_m\left({\mu }_n^{*}r\right)\right]{}^{2}dr},m=1,2,3,\mathrm{...},\infty \end{array}$$E59

H(t-τ) is a heavy-side unit step function which has the following properties:

$$\frac{dH\left(t\right)}{dt}=\delta \left(t\right)$$E60

$$H\left(t\right)=\{\begin{array}{l}1fort>0\\ 0fort\le 0\end{array}$$E61

Finally, the temperature field was constructed as:

$$\begin{array}{l}\Delta \left(r,\theta ,z,t\right)={\displaystyle {\int }_0^{t}d\tau }{\displaystyle {\int }_{R_0}^{R_2}{\displaystyle {\int }_{{}_{-\pi }}^{\pi }{\displaystyle {\int }_0^L\frac{-Q*\left(\xi ,\eta ,\lambda \right)}{k}\cdot G\left(r,\theta ,z,t;\xi ,\eta ,\lambda ,\tau \right)d\xi d\eta d\lambda }}}\\ +{\displaystyle {\int }_{R_0}^{R_2}{\displaystyle {\int }_{{}_{-\pi }}^{\pi }{\displaystyle {\int }_0^L\frac{-G\left(r,\theta ,z,t;\xi ,\eta ,\lambda ,0\right)\cdot F\left(\xi ,\eta ,\lambda \right)d\xi d\eta d\lambda }{\alpha }}}}\end{array}$$E62

where, $$Q*\left(r,\theta ,z\right)$$and $$F\left(r,\theta ,z\right)$$ were given in Equations (51-52), respectively.

This analytical solution has been applied to perform parametric studies on the bioheat transfer problems involved in prostate hyperthermia [26].

3.3. Analytical solutions to bioheat transfer with temperature fluctuation

Contributed from microcirculation including the capillary network plus small arterioles and venules of less than 100μm in diameter, blood perfusion plays an important role in the transport of oxygen, nutrients, pharmaceuticals, and heat through the body [27]. Although generally treated as a constant, blood perfusion is in fact a transient value even under physiological basal state. This is due to external perturbation and the self-regulation of biological body. The pulsative blood flow behavior also makes blood perfusion a fluctuating quantity. It is well accepted that blood perfusion is a fluctuating quantity around a mean value. Corresponding to the pulsative blood flow and very irregular distribution of blood vessels with various sizes, temperatures in intravital living tissues also appear fluctuating. To address this issue, we have obtained before an analytical model to characterize the temperature fluctuation in living tissues based on the Pennes equation [27]. It provides a theoretical foundation to better understanding the temperature fluctuation phenomena in living tissues. The Pennes equation used for the analysis can be rewritten as

$$\rho C\frac{\partial T}{\partial t}=K{\nabla }^{2}T+W_b{C}_b\left(T_a-T\right)+Q_m$$E63

Considering that small perturbations of arterial blood temperature, blood perfusion, and metabolic heat generation will result in tissue temperature fluctuation, each of these parameters can be expressed as the sum of a mean and a fluctuation value, i.e.

$T=\overline{T}+T^{\prime }$E64

$T_a=\overline{T_a}+{T^{\prime }}_a$E65

$W_b=\overline{W_b}+{W^{\prime }}_b$E66

$$Q_m=\overline{Q_m}+{Q^{\prime }}_m$$E67

where, symbol “ — ” represents the mean value, and “ ' ” the fluctuation one. Here, temporal averaging is adopted and defined as:

$\overline{A\left(x,y,z;t\right)}=\frac{1}{\Gamma }{\displaystyle {\int }_{t-\Gamma /2}^{t+\Gamma /2}A\left(x,y,z;\tau \right)}d\tau$E68

where, $A\left(x,y,z;t\right)$denotes the transient physical quantity at the vicinity of point ($x,y,z$), and $\Gamma$ the temporal averaging period.

Compared with the mean value, the fluctuation value is generally a small quantity. Then one has the following statistical relation:

$\overline{T^{\prime }}=0\begin{array}{cc},& \overline{{T^{\prime }}_a}=0\begin{array}{cc},& \overline{{W^{\prime }}_b}=0\end{array}\end{array}\begin{array}{cc},& \overline{{Q^{\prime }}_m}\end{array}=0$E69

Substituting Equations (65-68) into Equation (64) leads to:

$$\rho C\frac{\partial (\overline{T}+T^{\prime })}{\partial t}=K{\nabla }^2(\overline{T}+T^{\prime })+(\overline{W_b}+{W^{\prime }}_b)C_b\left[(\overline{T_a}+{T^{\prime }}_a)-(\overline{T}+T^{\prime })\right]+\overline{Q_m}+{Q^{\prime }}_m$$E70

Further,

$$\overline{\rho C\frac{\partial (\overline{T}+T^{\prime })}{\partial t}}=\overline{K{\nabla }^2(\overline{T}+T^{\prime })}+\overline{(\overline{W_b}+{W^{\prime }}_b)C_b\left[(\overline{T_a}+{T^{\prime }}_a)-(\overline{T}+T^{\prime })\right]}+\overline{\overline{Q_m}+{Q^{\prime }}_m}$$E71

Using Equation (70), Equation (72) was simplified as:

$$\rho C\frac{\partial \overline{T}}{\partial t}=K{\nabla }^2\overline{T}+\overline{W_b}{C}_b(\overline{T_a}-\overline{T})+C_b\overline{{W^{\prime }}_b{T^{\prime }}_a}-C_b\overline{{W^{\prime }}_b{T}^{\prime }}+\overline{Q_m}$$E72

Subtracting Equation (73) from Equation (71) leads to:

$$\begin{array}{l}\rho C\frac{\partial T^{\prime }}{\partial t}=K{\nabla }^2{T}^{\prime }+{W^{\prime }}_b{C}_b(\overline{T_a}-\overline{T})-C_b\overline{{W^{\prime }}_b{T^{\prime }}_a}+C_b\overline{{W^{\prime }}_b{T}^{\prime }}\\ \begin{array}{ccc}& & \end{array}-{W^{\prime }}_b{C}_b{T}^{\prime }-\overline{W_b}{C}_b{T}^{\prime }+\overline{W_b}{C}_b{T^{\prime }}_a+{W^{\prime }}_b{C}_b{T^{\prime }}_a+{Q^{\prime }}_m\end{array}$$E73

Equations (73) and (74) consist of the theoretical models for characterizing the temperature fluctuation in living tissues. Derivation of the perturbation Equation (74) is similar to that of the well known Reynolds equation in fluid mechanics. Compared with the Pennes equation, there are two additional terms appearing in Equation (73) both of which have explicit physical meaning: $$C_b\overline{{W^{\prime }}_b{T^{\prime }}_a}$$and $$-C_b\overline{{W^{\prime }}_b{T}^{\prime }}$$ represent the mean transferred energy due to perfusion perturbations and temperature fluctuations in tissue and arterial blood, respectively. It is from Equation (74) that the temperature fluctuation ($T^{\prime }$) and the pulsative blood perfusion, arterial blood temperature and metabolic heat generation ($W^{\prime }$, ${T^{\prime }}_a$, and${Q^{\prime }}_m$) were correlated. Clearly, since the mean tissue temperature $\overline{T}$ is a space dependent value, $T^{\prime }$is expected to be different at various tissue positions. To solve for Equation (73) and (74), dimension analysis was performed to simplify the model. One can express the orders of magnitude of the mean and pulsative physical quantities as:

$\begin{array}{l}\overline{T}~O(1)\begin{array}{cc}\begin{array}{cc}& \begin{array}{cc}\overline{T_a}~O(1)& \end{array}\end{array}& \overline{W_b}~\end{array}O(1)\begin{array}{cc}& \overline{Q_m}~O(1)\end{array}\\ T^{\prime }~O(\delta )\begin{array}{cc}\begin{array}{ccc}& {T^{\prime }}_a~O(\delta )& \end{array}& {W^{\prime }}_b~\end{array}O(\delta )\begin{array}{cc}& {Q^{\prime }}_m~O(\delta )\end{array}\end{array}$E74

where $O(\delta )$ stands for the value far less than$O(1)$, since the pulsative physical quantities investigated in this study is a small value.

Omitting those terms less than $O(1)$ in Equation (73) and those less than $O(\delta )$ in Equation (74), the Equations (73) and (74) can be respectively simplified as:

$$\rho C\frac{\partial \overline{T}}{\partial t}=K{\nabla }^2\overline{T}+\overline{W_b}{C}_b(\overline{T_a}-\overline{T})+\overline{Q_m}$$E75

$$\rho C\frac{\partial T^{\prime }}{\partial t}=K{\nabla }^2{T}^{\prime }+{W^{\prime }}_b{C}_b(\overline{T_a}-\overline{T})-\overline{W_b}{C}_b{T}^{\prime }+\overline{W_b}{C}_b{T^{\prime }}_a+{Q^{\prime }}_m$$E76

For the interpretation of temperature fluctuation in living tissues, it is reasonable to apply the 1-D degenerated forms of Equations (76) and (77). The boundary condition at the skin surface can be chosen as convective case which is often encountered in reality, i.e.

$-K\frac{\partial T}{\partial x}=h_f(T_f-T),x=0$E77

At the body core, a symmetrical or adiabatic condition can be used, namely

$-K\frac{\partial T}{\partial x}=0,x=L$E78

Then using the relations in Equations (65-68), the boundary and initial conditions of Equations (76) and (77) can be respectively obtained as:

$-K\frac{\partial \overline{T}}{\partial x}=h_f(T_f-\overline{T})\begin{array}{cc},& x=0\end{array}$E79

$-K\frac{\partial \overline{T}}{\partial x}=0\begin{array}{cc},& x=L\end{array}$E80

$\overline{T}=T_0(x)\begin{array}{cc},& t=0\end{array}$E81

$K\frac{\partial T^{\prime }}{\partial x}=h_f{T}^{\prime }\begin{array}{cc},& x=0\end{array}$E82

$-K\frac{\partial T^{\prime }}{\partial x}=0\begin{array}{cc},& x=L\end{array}$E83

$T^{\prime }=0\begin{array}{cc},& t=0\end{array}$E84

where, $x=0$is defined as the skin surface while $x=L$ the body core; $h_f$is the apparent heat convection coefficient between the skin and the environment which is the contribution from the natural convection and radiation; $T_f$the environment temperature; and $T_0(x)$ the initial temperature distribution. The following transformation was introduced [27]

$\overline{T(x,t)}=\overline{R(x,t)}\cdot \exp \left(-\frac{\overline{W_b}{C}_b}{\rho C}t\right)$E85

Then Equations (76) and (80-82) were respectively rewritten as:

$$\frac{\partial \overline{R}}{\partial t}=\alpha \frac{{\partial }^2\overline{R}}{\partial x^2}+\frac{\overline{W_b}{C}_b\overline{T_a}+\overline{Q_m}}{\rho C}\cdot \exp \left(\frac{\overline{W_b}{C}_b}{\rho C}t\right)$$E86

$-K\frac{\partial \overline{R}}{\partial x}=h_f\left[T_f\cdot \exp \left(\frac{\overline{W_b}{C}_b}{\rho C}t\right)\cdot H(t)-\overline{R}\right]\begin{array}{cc},& x=0\end{array}$E87

$-K\frac{\partial \overline{R}}{\partial x}=0\begin{array}{cc},& x=L\end{array}$E88

$\overline{R}=T_0(x)\begin{array}{cc},& t=0\end{array}$E89

where, $\alpha =K/\rho C$is the diffusivity of tissue, and $H(t)=\{\begin{array}{l}0\begin{array}{cc},& t<0\end{array}\\ 1\begin{array}{cc},& t>0\end{array}\end{array}$ the Heaviside function.

If the Green’s function for the above equation system is obtained, its transient solution can thus be constructed [13]. Through introducing an auxiliary problem corresponding to Equations (87-90), the Green’s function G can finally be obtained as:

$$G_{\overline{R}}\left(x,t;\xi ,\tau \right)={\displaystyle \sum _{n=1}^{\infty }{e}^{-\alpha {\beta }_n^2(t-\tau )}\cdot \frac{2\left[{\beta }_n^2+{\left(h_f/K\right)}^2\right]H\left(t-\tau \right)}{L\left[{\beta }_n^2+{\left(h_f/K\right)}^2\right]+h_f/K}\cdot \cos {\beta }_n(L-x)}\cdot \cos {\beta }_n(L-\xi )$$E90

where, the Eigen-values ${\beta }_n$ are the positive roots of the following equation

$${\beta }_n\cdot \tan ({\beta }_{n}L)=h_f/K$$E91

Then, the solution to Equation (87) can be obtained as

$$\begin{array}{l}\overline{R\left(x,t\right)}=\frac{\alpha }{k}{\displaystyle {\int }_0^t{G}_{\overline{R}}{\left(x,t;\xi ,\tau \right)|}_{\xi =0}\cdot h_f{T}_f\exp \left(\frac{\overline{W_b}{C}_b}{\rho C}\tau \right)\cdot H(\tau )d\tau +}\\ \begin{array}{ccc}& & \end{array}{\displaystyle {\int }_0^{t}d\tau {\displaystyle {\int }_0^L{G}_{\overline{R}}\left(x,t;\xi ,\tau \right)\frac{\overline{W_b}{C}_b\overline{T_a}+\overline{Q_m}}{\rho C}\cdot \exp \left(\frac{\overline{W_b}{C}_b}{\rho C}\tau \right)d\xi }}+\\ \begin{array}{ccc}& & \end{array}{\displaystyle {\int }_0^L{G}_{\overline{R}}{\left(x,t;\xi ,\tau \right)|}_{\tau =0}{T}_0(\xi )d\xi }\end{array}$$E92

Substituting Equation (93) into Equation (86) leads to the mean temperature$\overline{T(x,t)}$. The above solution is applicable to any transient environmental temperature $T_f(t)$ and space-dependent metabolic heat generation$\overline{Q_m(x,t)}$. For the fluctuation temperature, the solving process is as follows. Using the similar transformation as Equation (86)

$T^{\prime }(x,t)=R^{\prime }(x,t)\cdot \exp \left(-\frac{\overline{W_b}{C}_b}{\rho C}t\right)$E93

Equations (77) and (83-85) can be respectively converted to

$$\frac{\partial R^{\prime }}{\partial t}=\alpha \frac{{\partial }^2{R}^{\prime }}{\partial x^2}+\frac{{W^{\prime }}_b{C}_b(\overline{T_a}-\overline{T})+\overline{W_b}{C}_b{T^{\prime }}_a+{Q^{\prime }}_m}{\rho C}\cdot \exp \left(\frac{\overline{W_b}{C}_b}{\rho C}t\right)$$E94

$K\frac{\partial R^{\prime }}{\partial x}=h_f{R}^{\prime }\begin{array}{cc},& x=0\end{array}$E95

$-K\frac{\partial R^{\prime }}{\partial x}=0\begin{array}{cc},& x=L\end{array}$E96

$R^{\prime }=0\begin{array}{cc},& t=0\end{array}$E97

Following the same procedure described above, the Green’s function of this equation set is:

$$G_{R^{\prime }}\left(x,t;\xi ,\tau \right)=G_{\overline{R}}\left(x,t;\xi ,\tau \right)$$E98

where, $$G_{\overline{R}}\left(x,t;\xi ,\tau \right)$$was given in Equation (91).

Consequently, the fluctuation variable $R^{\prime }$ in Equation (95) can be obtained as

$$R^{\prime }\left(x,t\right)={\displaystyle {\int }_0^{t}d\tau {\displaystyle {\int }_0^L{G}_{R^{\prime }}\left(x,t;\xi ,\tau \right)\frac{{W^{\prime }}_b{C}_b(\overline{T_a}-\overline{T})+\overline{W_b}{C}_b{T^{\prime }}_a+{Q^{\prime }}_m}{\rho C}\cdot \exp \left(\frac{\overline{W_b}{C}_b}{\rho C}\tau \right)d\xi }}$$E99

where, the mean temperature $\overline{T(x,t)}$ was given by Equation (86). Substituting Equation (100) into Equation (94), the temperature fluctuation $T^{\prime }(x,t)$ can be obtained. As illustration, two calculation examples using the above analytical solutions were presented in this section to investigate the temperature fluctuation phenomena in living tissues. In the following illustration calculations, the typical values for tissue properties were taken from [27]. For clarity, only effect of the pulsative blood perfusion ${W^{\prime }}_b$ alone will be considered while the pulsative arterial temperature and metabolic heat generation were not taken into account, namely, ${T^{\prime }}_a=0\begin{array}{cc},& {Q^{\prime }}_m\end{array}=0$, and a constant mean blood perfusion $\overline{W_b}$ was assumed. Here, the pulsative blood perfusion ${W^{\prime }}_b$ was far less than the mean value$\overline{W_b}$, and its simplest form can be expressed as a periodic quantity with constant frequency and oscillation amplitude such as $W_{bm}\cos \omega t$ (where $\omega$ is frequency, $W_{bm}$the amplitude, and$W_{bm}<<\overline{W_b}$). However, to be more general, a stochastic pulsative perfusion form as ${W^{\prime }}_b=0.1\overline{W_b}\left(0.5-ran\right)$ (where, $ran$is random number generated by Fortran function) was adopted for calculations.

Figure 4.

Temperature fluctuation due to pulsative blood perfusion ($\overline{W_b}=0.5kg/m^3$°C)

Fig. 4 depicted both the skin surface temperature fluctuation (the mean perfusion $\overline{W_b}=0.5kg/m^3$°C) and the blood perfusion perturbation ${W^{\prime }}_b$ causing this behavior. As a comparison, Fig. 4(a) and Fig. 4(b) gave temperature fluctuations spanning two different time scope. It can be seen that small perturbation on blood perfusion resulted in an evident and observable temperature fluctuation in the living tissues. This result accords with the commonly observed phenomenon that the measured tissue temperature appears as fluctuating even when the measured animal is under physiologically basal state. It can also be found that the frequency of temperature fluctuation appears much smaller than that of the blood perfusion fluctuation, which implies that intravital biological tissue tends to keep its temperature stable. This result indicates that the stochastic fluctuation of blood perfusion in intravital biological tissue may also contribute to the tissue temperature oscillations, and the internal relations between blood perfusion fluctuation and the temperature oscillation need further clarification.

In this section, the perturbation model for characterizing the temperature fluctuation in living tissues was illustrated and its exact analytical solution was obtained which has wide applicability. One of the most important results in this section is perhaps that small perturbation in blood perfusion result in evidently observable temperature fluctuation in the living tissues. And the larger blood perfusion, the more liable for the living tissues to keep its temperature stable. This model provides a new theoretical foundation for better understanding the thermal fluctuation behavior in living tissues.

4. Bioheat transfer problems with phase change

4.1. Analytical solutions to 3-D phase change problems during cryopreservation

Derivation such solutions was based on the moving heat source method, in which all the thermal properties were considered as constants, and phase transition was assumed to occur in a single temperature [28]. The density, specific heat and heat conductivity of solid phase were considered to be the same as those of the liquid phase, respectively. To simplify the problem, only computation in a regular geometry characterized by Cartesian coordinates was considered, as shown in Fig. 5. According to the geometrical symmetry, only 1/8 of the whole cubic tissue was chosen as the study object, whose center is set as the origin point, and l, s₁, s₂ represent the distances between the origin and the boundaries of the tissue along x, y, z directions, respectively.

Figure 5.

Schematic of 1/8 cuboidal tissue subject to cryopreservation [29]

The energy equations for different phase regions were then written. For the liquid phase:

$$k{\nabla }^2{T}_l(x,y,z,t)=\rho c\frac{\partial T_l(x,y,z,t)}{\partial t},\text{in region}{R}_l,t>0$$E100

For the solid phase:

$$k{\nabla }^2{T}_s(x,y,z,t)=\rho c\frac{\partial T_s(x,y,z,t)}{\partial t},\text{in region}{R}_s,t>0$$E101

where$$k$$, $$\rho$$, c denote thermal conductivity, density and specific heat of tissue, respectively;$$T_s$$, $$T_l$$are the temperature distributions for solid and liquid phase, respectively; t is time;$$R_l$$, $$R_s$$denote solid and liquid region; the subscript l, s represent the liquid and solid phase, respectively. To obtain an analytical solution, all the parameters were assumed as uniform and remain constant.

It should be pointed out that the physical properties for the biological tissues would change during the phase change process. Therefore it may cause certain errors when assuming both the frozen and unfreezing regions take the same physical parameters. According to existing measurements, the density changes little and thus can be used as a constant. However, the other parameters, especially the thermal conductivity and the specific heat, would change significantly. For such case, one can choose to adopt an equivalent physical property to represent the original parameter, i.e. the parameters could take into concern contributions from both frozen and unfreezing phase, such that $$k=\overline{k}=\left(k_l+k_s\right)/2$$ and$$c=\overline{c}=\left(c_l+c_s\right)/2$$. This will minimize the errors via a simple however intuitive way. The solid and liquid phases are separated by an obvious interface, which can be expressed as

$$S(x_0,y_0,z_0,t)=z_0-s(x_0,y_0,z_0)=0$$E102

where, s denotes the moving solid-liquid interface, $$(x_0,y_0,z_0)$$represents position of any point in this specific interface and each of them is time dependent.

In the solid-liquid interface, conservation of energy and continuum of temperature read as

$$k\frac{\partial T_s}{\partial n}-k\frac{\partial T_l}{\partial n}=\rho L{\frac{\partial s}{\partial t}|}_n\text{at}\text{S}(x,y,z,t)=0,\begin{array}{cc}& \end{array}t>0$$E103

$$T_s=T_l{\text{=T}}_m\text{at}\text{S}(x,y,z,t)=0,\begin{array}{cc}& \end{array}t>0$$E104

where n is the unit normal vector; L is the latent heat of tissue; $$T_m$$is the phase change temperature; $$v_n(t)={\frac{\partial s}{\partial t}|}_n$$represents the velocity of the interface. Based on the moving heat source method [28, 29], the above phase change problem can be equivalently combined to a heat conduction problem with an interior moving heat source term, i.e.

$${\nabla }^{2}T+\frac{\rho L}{k}\frac{\partial s(x,y,t)}{\partial t}\delta (z-z_0)=\frac{1}{\alpha }\frac{\partial T}{\partial t},\text{in}\text{region}{R}_s+R_l,t>0$$E105

where, $$\alpha =k/\rho c$$is the thermal diffusivity; $$z_0$$is z coordinate of the arbitrary point p in the interface; $$\delta (z-z_0)$$is the delta function. For brief, Equation (106) can be rewritten as

$${\nabla }^{2}T+\frac{q_r(x,y,z,t)}{k}=\frac{1}{\alpha }\frac{\partial T}{\partial t}$$E106

where, a generalized volumetric heat source has been expressed as

$$q_r(x,y,z,t)=\rho L\frac{\partial s(x,y,t)}{\partial t}\delta (z-z_0)$$E107

The typical cooling situations most encountered in a cryopreservation [8, 10] include the following cases: (a) convective cooling at all boundaries by liquid nitrogen; (b) fixed temperature cooling at all boundaries through contacting to copper plate with very low temperature; (c) fixed temperature cooling at upside and underside surface of tissues and convective cooling at side faces; (d) convective cooling at upside and underside surface of tissues and fixed temperature cooling at side faces. Usually, the boundary types as (c) and (d) were adopted to increase the cooling rate. In this section, the analytical solutions will be presented according to the above four cooling cases, respectively.

Case (a): Convective cooling at all boundaries

Clinically, one of the most commonly used cooling approaches is to immerse the processed tissue into liquid nitrogen and frequently shift it up and down so as to enhance the heat exchange between the tissue and the liquid. Such boundary conditions can be defined as

$$\begin{array}{l}\frac{\partial T}{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=0;\begin{array}{cc}& \end{array}-k\frac{\partial T}{\partial x}=h(T-T_f)\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=l;\\ \frac{\partial T}{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=0;\begin{array}{cc}& \end{array}-k\frac{\partial T}{\partial y}=h(T-T_f)\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=s_1;\\ \frac{\partial T}{\partial z}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=0;\begin{array}{cc}& \end{array}-k\frac{\partial T}{\partial z}=h(T-T_f)\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=s_2;\end{array}$$E108

where, $$T_f$$is temperature of the cooling liquid; $$h$$is the convective heat transfer coefficient. Here, the three planes, i.e.$$x=0$$，$$y=0$$，$$z=0$$ are defined as adiabatic boundaries, and the other three planes, i.e. x=l, y=s₁, z=s₂, are subjected to forced convective cooling conditions. Without losing any generality, the initial temperature is defined as

$$T(x,y,z,0)=T_0(x,y,z)$$E109

Using transformation$$\theta =T-T_f$$, Equations (107), (109) and (110) can be rewritten as

$${\nabla }^2\theta +\frac{q_r}{k}=\frac{1}{\alpha }\frac{\partial \theta }{\partial t}$$E110

$$\begin{array}{l}\frac{\partial \theta }{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=0;\begin{array}{cc}& \end{array}-k\frac{\partial \theta }{\partial x}=h\theta \begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=l;\\ \frac{\partial \theta }{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=0;\begin{array}{cc}& \end{array}-k\frac{\partial \theta }{\partial y}=h\theta \begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=s_1;\\ \frac{\partial \theta }{\partial z}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=0;\begin{array}{cc}& \end{array}-k\frac{\partial \theta }{\partial z}=h\theta \begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=s_2;\end{array}$$E111

$$\theta (x,y,z,0)=F(x,y,z)=T_0(x,y,z)-T_f$$E112

To solve for the Green’s function of the above equations, the following auxiliary problem needs to be considered for the same region:

$${\nabla }^{2}G=\delta (x-\xi )\delta (y-\eta )\delta (z-\lambda )\delta (t-\tau )+\frac{1}{\alpha }\frac{\partial G}{\partial t}$$E113

$$\begin{array}{l}\frac{\partial G}{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=0;\begin{array}{cc}& \end{array}-k\frac{\partial G}{\partial x}=hG\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=l;\\ \frac{\partial G}{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=0;\begin{array}{cc}& \end{array}-k\frac{\partial G}{\partial y}=hG\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=s_1;\\ \frac{\partial G}{\partial z}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=0;\begin{array}{cc}& \end{array}-k\frac{\partial G}{\partial z}=hG\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=s_2;\end{array}$$E114

$$G(x,y,z,0)=0$$E115

The final expression for the Green’s function of Equations (111) and (112) can be obtained as:

$$\begin{array}{l}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )=-{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \sum _{j=0}^{\infty }\frac{\cos ({\beta }_{j}x)\cos ({\gamma }_{n}y)\cos ({\mu }_{k}z)\cos ({\beta }_j\xi )\cos ({\gamma }_n\eta )\cos ({\mu }_k\lambda )}{N({\beta }_j)N({\gamma }_n)N({\mu }_k)}}}}\times \\ \begin{array}{ccc}\begin{array}{ccc}\begin{array}{ccc}& & \end{array}& & \end{array}& & \end{array}H(t-\tau )\exp [-\alpha ({\beta }_j^2+{\gamma }_n^2+{\mu }_k^2)(t-\tau )]\\ \\ j=0,1,2;n=0,1,2;k=0,1,2,\end{array}$$E116

where,$$N({\beta }_j)=\frac{l[{\beta }_j^2+{(h/k)}^2]+h/k}{2[{\beta }_j^2+{(h/k)}^2]}$$ ,$${\beta }_j$$ are positive roots of the equation$${\beta }_j\tan ({\beta }_{j}l)=h/k$$;$$N({\gamma }_n)=\frac{s_1[{\gamma }_n^2+{(h/k)}^2]+h/k}{2[{\gamma }_n^2+{(h/k)}^2]}$$, $${\gamma }_n$$are positive roots of equation$${\gamma }_n\tan ({\gamma }_n{s}_1)=h/k$$; and$$N({\mu }_k)=\frac{s_{\text{2}}[{\mu }_k^2+{(h/k)}^2]+h/k}{2[{\mu }_k^2+{(h/k)}^2]}$$, $${\mu }_k$$are positive roots of$${\mu }_k\tan ({\mu }_k{s}_2)=h/k$$.

Finally, according to the above results and expression for the heat source term in Equation (108), the analytical solution to the temperature field under totally convective cooling conditions can then be obtained as:

$$\begin{array}{l}\theta (x,y,z,t)=T(x,y,z,t)-T_f=-{\displaystyle {\int }_0^{t}d\tau }{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}\frac{q_r}{k}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )d\xi }}}d\eta d\lambda -\\ \begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& \end{array}{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}\frac{1}{\alpha }G(x,y,z,t;\xi ,\eta ,\lambda ,0)\cdot F(\xi ,\eta ,\lambda )d\xi }}}d\eta d\lambda \end{array}$$E117

From the first term containing time in the above analytical solution, it can be seen that the thermal diffusivity $$\alpha$$ appears in the exponential term, i.e.$$\exp [-\alpha ({\beta }_j^2+{\gamma }_n^2+{\mu }_k^2)(t-\tau )$$, which indicates that the time for the temperature to reach the thermal equilibrium state depends exponentially on the thermal diffusivity.

Case (b): Fixed temperature cooling at all boundaries

Clinically, direct cooling the tissues through contacting it to copper plate pre-cooled by liquid nitrogen has been proved to be more effective than cooling by convection [28]. Therefore, it is very essential to get the temperature field of the tissue under totally fixed temperature cooling boundary conditions. For this problem, the form of the control equations still remain the same, so did the solution procedures of the Green’s function method, since only the boundary conditions were slightly changed. Assuming that $$T_p$$ is the temperature of the cooling plate, using transformation$$\theta =T-T_p$$, and solving the auxiliary problem via the similar way as that used in the convective cooling case, one can obtain the Green’s function solution for the present case as:

$$\begin{array}{l}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )=-{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \sum _{j=0}^{\infty }\frac{8\cos ({\beta }_{j}x)\cos ({\gamma }_{n}y)\cos ({\mu }_{k}z)\cos ({\beta }_j\xi )\cos ({\gamma }_n\eta )\cos ({\mu }_k\lambda )}{s_1{s}_{2}l}}}}\times \\ \begin{array}{ccc}\begin{array}{ccc}\begin{array}{ccc}& & \end{array}& & \end{array}& & \end{array}H(t-\tau )\exp [-\alpha ({\beta }_j^2+{\gamma }_n^2+{\mu }_k^2)(t-\tau )]\\ j=0,1,2;n=0,1,2;k=0,1,2,\end{array}$$E118

Then the transient temperature field can be constructed as:

$$\begin{array}{l}\theta (x,y,z,t)=T(x,y,z,t)-T_p=-{\displaystyle {\int }_0^{t}d\tau }{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}\frac{q_r}{k}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )d\xi }}}d\eta d\lambda \\ \begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& -\end{array}{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}\frac{1}{\alpha }G(x,y,z,t;\xi ,\eta ,\lambda ,0)\cdot F(\xi ,\eta ,\lambda )d\xi }}}d\eta d\lambda \end{array}$$E119

where,

$$F(x,y,z)=\theta (x,y,z,0)=T_0(x,y,z)-T_p$$E120

.

In practice, demanded by certain specific cooling rate and mechanical factors, sometimes one has to apply different cooling strategies on each side of the tissue surfaces. Thus, it is essential to take into account the complex hybrid boundary conditions. For brief, we assume that the temperature of the cooling plate $$T_p$$ is equal to that of the cooling fluid$$T_f$$. Then the Green’s function solutions for the following two boundaries can be obtained accordingly.

Case (c): Fixed temperature cooling at upside and underside surface and convective cooling at side faces

$$\begin{array}{l}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )=-{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \sum _{j=0}^{\infty }\frac{2\cos ({\beta }_{j}x)\cos ({\gamma }_{n}y)\cos ({\mu }_{k}z)\cos ({\beta }_j\xi )\cos ({\gamma }_n\eta )\cos ({\mu }_k\lambda )}{N({\beta }_j)N({\gamma }_n)s_2}}}}\times \\ \begin{array}{ccc}\begin{array}{ccc}\begin{array}{ccc}& & \end{array}& & \end{array}& & \end{array}H(t-\tau )\exp [-\alpha ({\beta }_j^2+{\gamma }_n^2+{\mu }_k^2)(t-\tau )]\\ j=0,1,2;n=0,1,2;k=0,1,2,\end{array}$$E121

Case (d): Convective cooling at upside and underside and fixed temperature cooling at side faces

$$\begin{array}{l}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )=-{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \sum _{j=0}^{\infty }\frac{4\cos ({\beta }_{j}x)\cos ({\gamma }_{n}y)\cos ({\mu }_{k}z)\cos ({\beta }_j\xi )\cos ({\gamma }_n\eta )\cos ({\mu }_k\lambda )}{N({\mu }_k)s_{2}l}}}}\times \\ \begin{array}{ccc}\begin{array}{ccc}\begin{array}{ccc}& & \end{array}& & \end{array}& & \end{array}H(t-\tau )\exp [-\alpha ({\beta }_j^2+{\gamma }_n^2+{\mu }_k^2)(t-\tau )]\\ j=0,1,2;n=0,1,2;k=0,1,2,\end{array}$$E122

Considering that expressions for the transient temperature field for the above two cases still remain similar to that of Equation (118), they have not been rewritten here for brief.

It should be pointed that there still exist many difficulties to calculate the exact temperature field from the above analytical solutions. However, the solution forms can still be flexibly applied to analyze certain special problems. As indicated in [28], in the freezing or warming process there must exist a maximum cooling or warming rate at some places of the tissue, which is varying with the time. Theoretically, this transient position can be predicted by using Equations (118) and (120) or other equations for corresponding processes. For example, one can obtain $\frac{\partial }{\partial x_i}\left[\frac{\partial T\left(x,y,z,t\right)}{\partial t}\right]$ (where x_i may represent x, y, z direction) easily by utilizing any of the above Green’s function solutions. After making it equal to zero, one can solve for x_i which is just the position where the maximum cooling or rewarming rate occurs. Knowing such position is of importance for the operation of a successful cryopreservation procedure, since the maximum cooling or warming rate is the crucial factor to cause injury to biological materials. Here, computation of $\partial T\left(x,y,z,t\right)/\partial t$ can eliminate the time integral term in Equations (118) and (120), which would simplify the solution form. Overall, the present analytical method in virtue of its straightforward form is of great significance to evaluate the phase change problem in cryobiology.

4.2. Analytical solutions to 3-D phase change problems during cryosurgery

Cryosurgery is very different from cryopreservation, since living tissue has to be considered. Consequently, it must take into account the effects of blood perfusion and metabolic heat generation into bioheat equation. Here, the Pennes equation is applied to characterize the heat transfer process in the living tissue. To avoid the complex boundary conditions, the calculation tissue domain is chosen as a whole cuboid as shown in Fig. 6, where a cryoprobe with length l₁ was settled in the center. Then the location of the cryoprobe is at$$x_{\text{0}}=l/2$$,$$y_{\text{0}}=s_1/2$$ , and the range for z coordinate is$$\text{0}\text{.5(}{s}_2-l_1)\le \text{z}\le \text{0}\text{.5(}{s}_2+l_1)$$.

Figure 6.

Schematic of the living tissue domain subject to cryosurgery

The energy equations for the tissue before and after it was frozen are respectively as:

For the liquid phase

$$\rho c\frac{\partial T_l}{\partial t}=k{\nabla }^2{T}_l+{\omega }_b{\rho }_b{c}_b(T_a-T_l)+Q_m+Q_{\text{c}},\text{in region}{R}_l,t>0$$E123

For the solid phase

$$\rho c\frac{\partial T_s(x,y,z,t)}{\partial t}=k{\nabla }^2{T}_s(x,y,z,t)+Q_{\text{c}},\text{in region}{R}_s,t>0$$E124

where, $${\rho }_b$$, $$c_b$$are density, specific heat of blood; $${\omega }_b$$is blood perfusion; $$T_a$$is supplying arterial temperature; $$Q_m$$is volumetric metabolic heat; $$Q_{\text{c}}$$is heat sink, and the other parameters represent the same meanings as before.

The control equations in the solid-liquid interface are the same as before. The above phase change problem can then be equivalently transformed to a heat conduction problem, i.e.

$$\rho c\frac{\partial T}{\partial t}=k{\nabla }^{2}T+Q_r$$E125

where, $$Q_r$$is the moving heat source term, which consists of three parts in a cryosurgical process: the metabolic heat source term$$Q_m$$, the heat sink term $$Q_{\text{c}}$$and the phase change heat source term$$q_r$$, i.e.

$$Q_r=q_r+Q_m+Q_c=\rho L\frac{\partial s(x,y,t)}{\partial t}\delta (z-z_0)+q_m^{\text{'}}H(z-z_0)+Q_c$$E126

where, $$q_m^{\text{'}}=w_b{c}_b{\rho }_b(T_a-T)+q_m$$, reflecting contributions of the blood heat transfer and the metabolic heat generation in unfrozen region. Other parameters and functions have the same definitions with those of cryopreservation.

To simplify the problem, the cryoprobe inserted into the deep tissue is treated as a linear heat sink [29] and assumed to supply a constant cold amount$$Q_{\text{c}}$$, i.e.$$Q_c=B$$, which can also include many different discrete terms representing cooling effects from multiple cryoprobes and expressed as $$Q_{\text{c}}={\displaystyle \sum Q_i\delta \left(x-x_i,y-y_i,z-z_i\right)}$$ (which reflects the amount of cold at the location of$$\left(x_i,y_i,z_i\right)$$, where i is the sequence number of the cryoprobe; $$Q_i$$is the amount of cold released by the ith probe; $$\delta$$is the Delta function). The solution expressed below indicated that the Green’s function method can deal with the multi-probe freezing problem. This is however not easy to be dealt with even by numerical computation. In this side, the analytical solution embodies its particular theoretical significance.

Equation (125) can be rewritten as:

$$\frac{\partial T}{\partial t}=\alpha {\nabla }^{2}T-\frac{{\omega }_b{c}_b{\rho }_b}{\rho c}H(z-z_0)T+g(z,t)$$E127

where,$$g(z,t)=[({\omega }_b{c}_b{\rho }_b{T}_a+q_m)H(z-z_0)+q_r+Q_c]/\rho c$$. In a cryosurgical procedure, the cryoprobe is inserted into the target tissue, which will subject to a specific temperature drop due to cooling of heat sink effect of the probe. Since the skin surface is exposed to the ambient environment, the boundary conditions can thus be treated as:

$$\begin{array}{l}-k\frac{\partial T}{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=0;\begin{array}{cc}& \end{array}-k\frac{\partial T}{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=l\\ -k\frac{\partial T}{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=0;\begin{array}{cc}& \end{array}-k\frac{\partial T}{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=s_1\\ -k\frac{\partial T}{\partial z}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=0;\begin{array}{cc}& \end{array}-k\frac{\partial T}{\partial z}=h(T-T_f)\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=s_2\end{array}$$E128

As shown in Fig. 6, the origin (x=0, y=0, z=0) is defined differently from the cryopreservation. The adiabatic boundaries are applied on the side surfaces along the x, y axis. This is because the side surfaces are far from the heat sink and are not influenced by the heat sink in the deep tissue and the external convective conditions. Finally, the initial condition is defined as$$T(x,y,z,0)=T_0(x,y,z)$$, and every boundary is defined as adiabatic at this initial state. Although the control equations and boundary conditions in cryosurgery are very different from that in cryopreservation, the solution procedures still remain similar if certain transformations are introduced. In the following, only those steps different from the above will be addressed. To make solution of the problem feasible, we have adopted before the following specific transformation:

$$T(x,y,z,t)=T_0(x,y,z)+W(x,y,z,t)\exp [-\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}H(z-z_0)t]$$E129

Substituting it into Equation (127), one obtains

$$\frac{\partial W}{\partial t}=\alpha {\nabla }^{2}W+g_1(z,t)\exp [\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}H(z-z_0)\text{t]}$$E130

where, $$g_1(z,t)$$=$$\{[{\omega }_b{c}_b{\rho }_b(T_a-T_0)+q_m]H(z-z_0)+q_r+Q_c\}/\rho c$$. To simplify the equation, the volumetric moving heat source term in Equation (130) can be expressed as

$$q_r^{*}=g_1(z,t)\exp [\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}H(z-z_0)\text{t]}$$E131

The boundary conditions are rewritten as:

$$\begin{array}{l}-k\frac{\partial W}{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=0;\begin{array}{cc}& \end{array}-k\frac{\partial W}{\partial x}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}x=l\\ -k\frac{\partial W}{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=0;\begin{array}{cc}& \end{array}-k\frac{\partial W}{\partial y}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}y=s_1\\ -k\frac{\partial W}{\partial z}=0\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=0;\begin{array}{cc}& \end{array}-k\frac{\partial W}{\partial z}=h[\text{W}-f(t)]\begin{array}{cc}& \end{array}at\begin{array}{cc}& \end{array}z=s_2\end{array}$$E132

where, $$f(t)=(T_f-T_0)\exp [\frac{{\omega }_b{\rho }_b{c}_b}{\rho c}H(z-z_0)t]H(t)$$, and the initial condition is defined as$$W(x,y,z,t)=0$$. Applying the same theoretical strategies as used in solving cryopreservation into the above problem, one can obtain the Green’s function expression for a 3-D cryosurgical process as:

$$\begin{array}{l}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )=\frac{\text{1}}{s_1{s}_{2}l}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{n=0}^{\infty }{\displaystyle \sum _{j=0}^{\infty }{R}_1}}}{R}_2\frac{2[{\mu }_k^2+{(h/k)}^2]}{[{\mu }_k^2+{(h/k)}^2]+h/k}\exp [-\alpha ({\beta }_j^2+{\gamma }_n^2+{\mu }_k^2)(t-\tau )]\times \\ \begin{array}{ccc}& & \end{array}\begin{array}{c}\end{array}H(t-\tau )\cos ({\beta }_{j}x)\cos ({\gamma }_{n}y)\cos ({\mu }_{k}z)\cos ({\beta }_j\xi )\cos ({\gamma }_n\eta )\cos ({\mu }_k\lambda )\\ j=0,1,2;n=0,1,2;k=0,1,2,\end{array}$$E133

where,$$R_1=\{\begin{array}{l}1,\begin{array}{c}\end{array}j=0\\ 2,\begin{array}{c}\end{array}j=1,2,3,\cdot \cdot \cdot \end{array}$$;$$R_2=\{\begin{array}{l}1,\begin{array}{c}\end{array}n=0\\ 2,\begin{array}{c}\end{array}n=1,2,3,\cdot \cdot \cdot \end{array}$$;$${\beta }_j=\frac{j\pi }{l}$$;$${\gamma }_n=\frac{n\pi }{s_1}$$; $${\mu }_k$$is the positive root of$${\mu }_k\tan ({\mu }_k{s}_2)=h/k$$. Taking into account of the moving source term $$q_r^{*}$$ and the boundary conditions, i.e. Equation (132), one can obtain the transient temperature field of the tissue corresponding to Equation (130) as:

$$\begin{array}{l}W(x,y,z,t)=-{\displaystyle {\int }_0^{t}d\tau }{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}\frac{q_r^{*}}{k}G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )d\xi }}}d\eta d\lambda -\\ \begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& \end{array}{\displaystyle {\int }_0^{t}d\tau }{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}\frac{1}{k}{G(x,y,z,t;\xi ,\eta ,\lambda ,\tau )|}_{\lambda =s_2}\cdot h\cdot f(t)d\xi d\eta }}-\\ \begin{array}{cc}\begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}& \end{array}{\displaystyle {\int }_0^l{\displaystyle {\int }_0^{s_1}{\displaystyle {\int }_0^{s_2}\frac{1}{\alpha }G(x,y,z,t;\xi ,\eta ,\lambda ,0)F(\xi ,\eta ,\lambda )d\xi }}}d\eta d\lambda \end{array}$$E134

The above procedures illustrate the basic strategy to exactly solve the three dimensional phase change problem of biological tissues in vivo, which involves the blood perfusion and metabolic heat etc. However, the integral equation is so complex due to moving phase change front inherited in the integral term, that calculating the equation based on the above analytical expressions is still a challenge. This requests certain development of the applied mathematics. However, a simplified form for the present solution can be utilized to analyze some specific one dimensional heat transfer problems.

4.3. Analytical solutions to 3-D temperature distribution in tissues embedded with large blood vessels during cryosurgery

From the viewpoint of heat transfer, a large blood vessel (also termed a thermally significant vessel) denotes a vessel larger than 0.5 mm in diameter [30]. Anatomically, tumors are often situated close to or embedded with large blood vessels, since a tumor’s quick growth ultimately depends on nutrients supplied by its blood vessel network. During cryosurgery, the blood flow inside a large vessel represents a source which heats the nearby frozen tissues and, thereby, limits freezing lesions during cryosurgery. Under this condition, a part of the vital tumor cells may remain in the cryolesion and lead to recurrence of tumors after cryosurgical treatment. More specifically, tumor cell survival in the vicinity of large blood vessels is often correlated with tumor recurrence after treatment [30]. Consequently, it is difficult to implement an effective cryosurgery when a tumor is contiguous to a large blood vessel. To better understand the effect of blood flow to the temperature distribution of living tissues subject to freezing, a conceptual model for characterizing the heat transfer in 3-D cylindrical tissues embedded with a single blood vessel was illustrated in this section. And a closed form analytical solution to this model was provided to explore different factors’ thermal influences to the freezing mechanism of living tissues.

The geometry used for the analysis is depicted as Fig. 7, which is consisted of three distinct concentric cylinders: the most interior region representing a large blood vessel, the intermediate for unfrozen liquid-phase tissue and the outer the frozen tissue. In Fig. 7, symmetrical condition in θ direction can be applied. Then the 3-D bioheat transfer will degenerate to a 2-D problem. The cylinder is long enough so that its end effects to the heat transfer can be neglected. For simplicity in analytical solution, only steady state temperature field was assumed in both the vessel and the surrounding tissue. And the same constant thermal properties for different tissue area were considered. Therefore heat conduction in the regions of both liquid and frozen tissue can be described by only a single equation.

Figure 7.

Schematic of cylindrical tissues embedded with a single blood vessel [30]

To carry out the theoretical analysis, additional assumptions for tissues were made as follows: there is no heat flow across the boundaries at z=0 and z=L. In reality, if the vessel length is long enough, these adiabatic boundary conditions can be closely satisfied. And the cylindrical tissue surface (r=b) was assumed to be immersed into a cooling medium or subjected to a circular cryoprobe with constant freezing temperature T₀. This situation also reflects the cooling of an interior tissue inside the biological body, which is frozen on the surface. After introducing the dimensionless parameter$$T_t^{*}=(T_t-T_0)/(T_a-T_0)$$, where T_a stands for the body core temperature usually fixed at 37°C, the temperature distribution for the cylindrical tissue area can be described by the following equations [30]

$$\{\begin{array}{l}\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial T_t^{*}(r,z)}{\partial r})+\frac{\partial T_t^{*}(r,z)}{\partial z^2}=0\\ T_t^{*}=0\begin{array}{cccc}& & ,& \begin{array}{ccc}& at& r=b\end{array}\end{array}\\ T_t^{*}=f(z)\begin{array}{cccc}& ,& & \begin{array}{cc}at& r=a\end{array}\end{array}\\ \frac{\partial T_t^{*}(r,z)}{\partial z}=0\begin{array}{ccccc},& & at& z=0,L& \end{array}\end{array}$$E135

where, f(z) is temperature along the vessel wall, which is unknown here and needs to be calculated later using temperature continuity condition on the vessel wall.

Blood flow velocity profile in the vessel can be obtained as $$u=u_0(1-\frac{r^2}{a^2})$$ by solving the Navier-Stokes equation for an incompressible steady-state fluid, where a stands for the radius of blood vessel, and $$u_0$$ the average blood flow velocity over the vessel cross-section. Both temperature and heat flux on the vessel wall ($$r=a$$) satisfy continuity conditions. Defining the following dimensionless parameters:$$T_b^{*}=\frac{T_b^{}-T_0}{T_a-T_0}$$,$$T_c^{*}=\frac{T_c^{}-T_0}{T_a-T_0}$$ , the governing equations for the vessel read as

$$\{\begin{array}{l}\frac{K_b}{{\rho }_b{C}_b}\frac{1}{r}(\frac{\partial }{\partial r}(r\frac{\partial T_b^{*}}{\partial r}))=u_0(1-\frac{r^2}{a^2})\frac{d{T}_c^{*}}{dz}\\ \frac{\partial T_b^{*}}{\partial r}=0\begin{array}{ccccccc}& & ,& & & & at\begin{array}{c}\end{array}r=0\end{array}\\ T_b^{*}=f(z)\begin{array}{ccc}& ,& \begin{array}{ccc}& & \end{array}\begin{array}{ccc}& at\begin{array}{c}\end{array}r=a& \end{array}\end{array}\\ T_b^{*}=1\begin{array}{ccccccc}& & ,& & & & \end{array}at\begin{array}{c}\end{array}z=0\end{array}$$E136

where, $$T_c^{*}$$is the vessel bulk temperature within cylindrical symmetry, and defined as

$$T_c^{*}=\frac{4}{a^2}{\displaystyle {\int }_0^a{T}_b^{*}(1-\frac{r^2}{a^2}})rdr$$E137

In the above equations, subscript t, b designate tissue and blood respectively, * means the dimensionless parameters. The variable separation method was applied to solve the Equation (135), whose solution can be rewritten as product of an axial term Z(z) and a radial one R₀(r), i.e.$$T_t^{*}(r,z)=R_0(r)Z(z)$$, with Z(z) satisfying

$$\{\begin{array}{l}\frac{d^{2}Z(z)}{d{z}^2}+{\beta }_m^{2}Z(z)=0\begin{array}{ccc}\begin{array}{cc}& \end{array}& & (0<z<L)\end{array}\\ \frac{dZ(z)}{dz}=0\begin{array}{cccc}& & & \begin{array}{ccc}& & \begin{array}{cc}& (z=0,L)\end{array}\end{array}\end{array}\end{array}$$E138

and R₀(r) satisfying

$$\{\begin{array}{l}\frac{d^2{R}_0(r)}{d{r}^2}+\frac{1}{r}\frac{d{R}_0(r)}{dr}-{\beta }_m^2{R}_0(r)=0\begin{array}{ccc}& & (a<r<b)\end{array}\\ R_0=0\begin{array}{cccc}& & & \begin{array}{cccc}& & & \begin{array}{cccc}& & & \begin{array}{cc}\begin{array}{cc}& \end{array}& (r=b)\end{array}\end{array}\end{array}\end{array}\end{array}$$E139

The solution to Equation (138) is thus obtained as$$Z(z)=C_1\cos ({\beta }_{m}z)+C_2\sin ({\beta }_{m}z)$$. $${\beta }_m$$is the positive roots of$$\sin (\beta L)=0$$. Then

$${\beta }_m=\frac{m\pi }{L},m=0,1,2\cdot \cdot \cdot$$E140

Therefore

$$Z(z)=\cos ({\beta }_{m}z)$$E141

with normal as

$$N({\beta }_m)=\{\begin{array}{l}L/2\begin{array}{ccc}& ,& \end{array}m=1,2,3,\cdot \cdot \cdot \\ L\begin{array}{ccc}& ,& \end{array}m=0\end{array}$$E142

The solution to Equation (139) can be obtained as

$$R_0=\{\begin{array}{l}{C}_3[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}r)+K_0({\beta }_{m}r)],{\beta }_m\ne 0\\ C_4\ln (r/b),\begin{array}{ccccccc}& & & & & & \end{array}{\beta }_0=0\end{array}$$E143

Then the solution to Equation (135) can be expressed as

$$T_t^{*}(r,z)={\displaystyle \sum _{m=1}^{\infty }{A}_m}[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}r)+K_0({\beta }_{m}r)]Z({\beta }_m,z)+A_0\ln (r/b)Z({\beta }_0,z)$$E144

where, $$A_m$$,$$A{}_0$$ are coefficients. Substituting Equation (135) into Equation (144) leads to

$$f(z)={\displaystyle \sum _{m=1}^{\infty }{A}_m}[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}a)+K_0({\beta }_{m}a)]Z({\beta }_m,z)+A_0\ln (r/b)Z({\beta }_0,z),0<z\le L$$E145

where, $$A_m$$, $$A_0$$are obtained as

$$A_m=\frac{1}{[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}a)+K_0({\beta }_{m}a)]N({\beta }_m)}{\displaystyle {\int }_0^{L}Z({\beta }_m,z)f(z)dz}$$E146

$$A_0=\frac{1}{\ln (a/b)N({\beta }_0)}{\displaystyle {\int }_0^{L}Z({\beta }_m,z)f(z)dz}$$E147

Substituting Equations (146-147) into Equation (144) leads to

$$\begin{array}{l}{T}_t^{*}(r,z)=(T_a-T_0)(\begin{array}{l}\\ \\ \end{array}{\displaystyle \sum _{m=1}^{\infty }\frac{2[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}r)+K_0({\beta }_{m}r)]}{L[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}a)+K_0({\beta }_{m}a)]}\cos ({\beta }_{m}z){\displaystyle {\int }_0^L\cos ({\beta }_{m}z)f(z)dz}}\\ \begin{array}{cccc}\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& & & \end{array}& & & +\frac{\ln (r/b)}{L\ln (a/b)}{\displaystyle {\int }_0^{L}f(z)dz}\begin{array}{l}\\ \\ \end{array})+T_0\end{array}\end{array}$$E148

As to the exact temperature profile within blood vessel, if defining$$U=\frac{u_0{\rho }_b{C}_b}{4K_b}\frac{d{T}_c^{*}}{dz}$$, one has from Equation (136)

$$T_b^{*}=U[r^2-\frac{r^4}{4a^2}+C_5(z)\ln r+C_6(z)]$$E149

Applying Equation (136) to this equation, $$C_5(z)=0$$was obtained. Substituting Equation (149) into Equation (137) leads to

$$T_c^{*}=[\frac{7a^2{}_{}}{24}+C_6(z)]\frac{u_0{\rho }_b{C}_b}{4K_b}\frac{d{T}_c^{*}}{dz}$$E150

Further, it can be derived as

$$T_c^{*}=T_{c0}^{*}{\displaystyle {\int }_0^z\frac{1}{[\frac{7a^2{}_{}}{24}+C_6(z)]\frac{u_0{\rho }_b{C}_b}{4K_b}}}dz$$E151

where, $$T_{c0}^{*}$$is defined as$$T_{c0}^{*}=\frac{4}{a^2}{\displaystyle {\int }_0^a{T}_{b0}^{*}(1-\frac{r^2}{a^2})rdr}$$, $$T_{b0}^{*}$$is the vessel inlet temperature at z=0. From Equation (136), one obtains $$T_{b0}^{*}$$=1. Therefore

$$T_{c0}^{*}=1$$E152

Substituting Equation (152) into Equation (151) leads to

$$T_b^{*}(r,z)=\frac{[r^2-\frac{r^4}{4a^2}+C_6(z)]}{4[\frac{7a^2{}_{}}{24}+C_6(z)]}{e}^{{\displaystyle {\int }_0^z}\frac{1}{[\frac{7a^2{}_{}}{24}+C_6(z)]\frac{u_0{\rho }_b{C}_b}{4K_b}}dz}$$E153

At r=a, one has

$$f(z)=\frac{[\frac{3}{4}{a}^2+C_6(z)]}{4[\frac{7a^2{}_{}}{24}+C_6(z)]}{e}^{{\displaystyle {\int }_0^z}\frac{1}{[\frac{7a^2{}_{}}{24}+C_6(z)]\frac{u_0{\rho }_b{C}_b}{4K_b}}dz}$$E154

Using the continuity condition for heat flux on the blood vessel wall, one approximately has

$$k_t\frac{\partial T_t^{}}{\partial r}|{}_{r=a}=k_b\frac{\partial T_b^{}}{\partial r}|{}_{r=a}$$E155

The two terms in this equation can further be obtained from Equation (148) and Equation (153), respectively, which are

$$\begin{array}{l}\frac{\partial T_t^{}(r,z)}{\partial r}|{}_{r=a}=(T_a-T_0)\{\begin{array}{l}\\ \\ \\ \end{array}{\displaystyle \sum _{m=1}^{\infty }\frac{2{\beta }_m[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_1({\beta }_{m}a)+K_1({\beta }_{m}a)]}{L[-\frac{K_0({\beta }_{m}b)}{I_0({\beta }_{m}b)}{I}_0({\beta }_{m}a)+K_0({\beta }_{m}a)]}\cos ({\beta }_{m}z){\displaystyle {\int }_{\begin{array}{c}\end{array}0}^{\begin{array}{c}\end{array}L}\cos ({\beta }_{m}z)f(z)dz}}\\ \begin{array}{cccc}\begin{array}{cccc}\begin{array}{ccc}& & \end{array}& & & \end{array}& & & +\frac{1}{La\ln (a/R)}{\displaystyle {\int }_{{}^{}0}^{{}^{}L}f(z)dz}\begin{array}{l}\\ \\ \\ \end{array}\}\end{array}\end{array}$$E156

and

$$\frac{\partial T_b^{}}{\partial r}|{}_{r=a}=\frac{a(T_a-T_0)}{4[\frac{7a^2{}_{}}{24}+C_6(z)]}{e}^{{\displaystyle {\int }_0^z}\frac{1}{[\frac{7a^2{}_{}}{24}+C_6(z)]\frac{u_0{\rho }_b{C}_b}{4K_b}}dz}$$E157

Exact calculation on $$C_6(z)$$ from Equation (155) appears extremely difficult. However, if treating $$C_6$$ as a constant, and substituting Equation (156) and Equation (157) into Equation (155), then integrating it from z = 0 to L, one has

$${\overline{C}}_6=a^2[\begin{array}{l}\\ \end{array}\frac{K_b}{K_t}\ln (\frac{a}{b})-\frac{3}{4}\begin{array}{l}\\ \end{array}]$$E158

This constant $${\overline{C}}_6$$ does not exactly satisfy Equation (155) indeed. However it is a simple and intuitive approximation. According to former calculation [30], it can be found that heat fluxes calculated using temperature from both sides of the blood vessel wall were almost identical. Therefore, an approximate estimation on $$C_6$$ can also be made as follows. In Equation (155), if the axial position z is fixed, $$C_6$$can be obtained through numerical iteration. And the relative error $$[C_6(z)-{\overline{C}}_6]/{\overline{C}}_6$$was very small along the axial direction. For the case of aorta, the largest error is less than 1%, while for the case of terminal branches, it is 6%. Overall, the smaller vessel diameter, the larger error in the estimated value$${\overline{C}}_6$$. However, since the present discussion was mainly focused on the case of large blood vessel (for tissues with extremely small vessel, a collective model such as Pennes equation is often applicable), therefore treating $$C_6(z)$$ as a constant was an acceptable approximation. But for more complex situations where the above approximation cannot hold, exactly solving the Equation (158) for $$C_6(z)$$ is still very necessary in the future.

When analyzing the thermal effect of large blood vessel in cryosurgery, an important issue is when the blood vessel begins to freeze and how to control cryoprobe’s temperature to completely freeze the target tumor. Substituting Equation (158) into Equation (153), one can set up the relation for the blood vessel temperature $$T_b(a,L)$$ to reach the phase change point T_m, i.e.

$$T_0=\frac{\frac{\frac{3}{4}{a}^2+{\overline{C}}_6}{4(\frac{7}{24}{a}^2+{\overline{C}}_6)}{e}^{\frac{L}{4(\frac{7}{24}{a}^2+{\overline{C}}_6)(\frac{u_0{\rho }_b{C}_b}{4K_b})}}{T}_a-T_m}{\frac{\frac{3}{4}{a}^2+{\overline{C}}_6}{4(\frac{7}{24}{a}^2+{\overline{C}}_6)}{e}^{\frac{L}{4(\frac{7}{24}{a}^2+{\overline{C}}_6)(\frac{u_0{\rho }_b{C}_b}{4K_b})}}-1}$$E159

Up to now, the above analysis was based on a steady state assumption. And only a single blood vessel was considered for the sake of analytical solution, although it does provide certain important information for understanding the phase change heat transfer in living tissues with blood vessel. Clinically, knowledge on the transient temperature response is still very necessary for the successful operation of a cryosurgery. However, such non-steady state problem cannot be dealt with by the present method. This needs further efforts in the near future using numerical approach.

5. Conclusion

This chapter has presented an overview on several typical closed form analytical solutions to 3-D bioheat transfer problems with or without phase change as developed before in the authors’ laboratory. In these solutions, relatively complex boundary conditions and heating/cooling on skin surface or inside biological bodies were addressed. In addition, the theoretical strategies towards analytically solving the complex 3-D bioheat transfer problems were outlined by the mathematical transformation, the Green’s function method, and the moving heat source model etc.

The analytical solutions introduced in this chapter can be used to predicate the evolution of temperature distribution inside the target tissues during tumor hyperthermia, cryosurgery, cryopreservation, thermal diagnostics, thermal comfort analysis, brain hypothermia resuscitation, and burn injury evaluation. Through fitting the predicted with the experimentally measured temperatures at the skin surface, some thermal parameters of biological tissues such as blood perfusion, thermal conductivity, and heat capacity, can be estimated non-invasively. Moreover, based on the requirements for freezing/heating necrosis temperature of tissue, an approach to optimize the parameters of cryosurgical/hyperthermic treatment can be obtained using the presented analytical solutions. Therefore, the presented analytical solutions are very useful for a variety of thermal-oriented biomedical studies. However, it should be pointed out that although such analytical solutions have some versatility in dealing many bioheat transfer problems, numerical approaches are still needed for more complex situations. In fact, the relation between analytical and numerical solutions should be complementary. On one hand, numerical approach can deal with more complex problem than analytical way. On the other hand, the analytical results can serve as benchmark solutions for numerical analyses on complex situations. In summary, it is believed that even the applications with some simplified conditions do not affect the applicability of the present analytical solutions.

Nomenclature

Acknowledgement

Part of the researches as presented in this chapter has been supported by the National Natural Science Foundation of China under grants Grant Nos. 51076161 and 81071255, the Specialized Research Fund for the Doctoral Program of Higher Education, and Research Fund from Tsinghua University under Grant No. 523003001.