## 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:

where,

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*_{1} = *s*_{2} = 0.08m and height *L* = 0.03m was depicted as the shadowed region in Fig. 1, where *s*_{1} and *s*_{2} 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*.

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:

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

or

where,

The body core temperature was regarded as a constant (

The BCs at *y* and *z* directions can be expressed as

The reason for adopting the adiabatic conditions in the two ends of the

The initial temperature is

where,

where,

The solution to Equation (11) is:

(12) |

where,

Equation (2) was transformed to the following form:

where,

where,

Using Green function method,

(25) |

(26) |

where,

The Eigen-values

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

(34) |

For the third BC, the solution is

(35) |

Then the tissue temperature field can be constructed as:

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:

where,

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.

### 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

where Q is the applied microwave power, *C*_{t} is a proportional constant, *ε* is the microwave attenuation constant in tissue, and *z*_{0} 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:

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

where, *R*_{0} and *R*_{2} 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:

The initial temperature is

Using transformation:

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

where,

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

(52) |

where

Here *J*_{m} and *N*_{m} are the Bessel functions of the first kind and the second kind, respectively.

(59) |

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

Finally, the temperature field was constructed as:

(62) |

where,

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

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.

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

where,

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

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

Further,

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

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

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:

where

Omitting those terms less than

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.

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

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

where,

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

where,

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:

where, the Eigen-values

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

(92) |

Substituting Equation (93) into Equation (86) leads to the mean temperature

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

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

where,

Consequently, the fluctuation variable

where, the mean temperature

Fig. 4 depicted both the skin surface temperature fluctuation (the mean perfusion

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*_{1}, *s*_{2} represent the distances between the origin and the boundaries of the tissue along *x*, *y*, *z* directions, respectively.

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

For the solid phase:

where*c* denote thermal conductivity, density and specific heat of tissue, respectively;*t* is time;*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

where, *s* denotes the moving solid-liquid interface,

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

where *n* is the unit normal vector; *L* is the latent heat of tissue;

where, *z* coordinate of the arbitrary point *p* in the interface;

where, a generalized volumetric heat source has been expressed as

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

(108) |

where, *x=l, y=s*_{1}*, z=s*_{2}*,* are subjected to forced convective cooling conditions. Without losing any generality, the initial temperature is defined as

Using transformation

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

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

(116) |

where,

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:

(117) |

From the first term containing time in the above analytical solution, it can be seen that 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

(118) |

Then the transient temperature field can be constructed as:

(119) |

where,

.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

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

(121) |

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

(122) |

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 *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

### 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*_{1} was settled in the center. Then the location of the cryoprobe is at*z* coordinate is

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

For the liquid phase

For the solid phase

where,

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.

where,

where,

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*i* is the sequence number of the cryoprobe; *i*th probe;

Equation (125) can be rewritten as:

where,

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

Substituting it into Equation (127), one obtains

where,

The boundary conditions are rewritten as:

where,

(133) |

where,

(134) |

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.

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*_{0}. 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*_{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]

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

where,

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*_{0}(*r*), i.e.*Z*(*z*) satisfying

and *R*_{0}(*r*) satisfying

The solution to Equation (138) is thus obtained as

Therefore

with normal as

The solution to Equation (139) can be obtained as

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

where,

where,

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

(148) |

As to the exact temperature profile within blood vessel, if defining

Applying Equation (136) to this equation,

Further, it can be derived as

where, *z*=0. From Equation (136), one obtains

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

At *r*=*a*, one has

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

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

(156) |

and

Exact calculation on *z* = 0 to *L*, one has

This constant *z* is fixed,

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*_{m}, i.e.

(159) |

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.