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

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.


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.

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. where, ρ , c are the density and the specific heat of tissue, respectively; b ρ and b c denote the density and the specific heat of blood, respectively; X contains the Cartesian coordinates x , y and z ; Ω denotes the analyzed spatial domain; ( ) k X is the space dependent thermal conductivity; and 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. a T is the blood temperature in the arteries supplying the tissue and is often treated as a constant at 37°C; ( , ) T t X is the tissue temperature; m ( , ) Q t X is the metabolic heat generation; and r ( , ) Q t X 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.

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][22][23][24][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 s1 = s2 = 0.08m and height L = 0.03m was depicted as the shadowed region in Fig. 1, where s1 and s2 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.  [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: The generalized boundary conditions (BCs) often encountered in a practical clinical situation can be written as: where, ( ) The body core temperature was regarded as a constant ( c T ) on considering that the biological body tends to keep its core temperature to be stable, i.e.
The BCs at y and z directions can be expressed as 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 where, ( ) 0 , , T x y z 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: where, 0 h 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 f T the surrounding air temperature.
The solution to Equation (11) is: Through using the following transformation [13] Equation (2) was transformed to the following form: , , ; exp where, k c α ρ = is the thermal diffusivity of tissue. The corresponding boundary and initial conditions are: where, Using Green function method, ( ) , , ; W x y z t can be solved from the combined Equations (14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24). The Green's functions 1 G and 2 G to the second and third BCs can finally be obtained: The Eigen-values p β are the positive roots of the following equation For the third BC, the solution is 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, 1 ( ) P t is the point-heating power, δ is the Dirac function, 0 0 0 ( , , ) x y z 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 pointheating 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.

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 ∞ 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]: where Q is the applied microwave power, Ct is a proportional constant, ε is the microwave attenuation constant in tissue, and z0 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 T ∞ . Therefore, one has: where, R0 and R2 are the urethra and prostate radius, respectively; Tf 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-43) as: The Green's function for the above equation sets can be obtained as: , 0,1,2,3...
Here Jm and Nm are the Bessel functions of the first kind and the second kind, respectively.
is a heavy-side unit step function which has the following properties: Finally, the temperature field was constructed as:  This analytical solution has been applied to perform parametric studies on the bioheat transfer problems involved in prostate hyperthermia [26].

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.

T T T′
where, symbol " -" represents the mean value, and " ' " the fluctuation one. Here, temporal averaging is adopted and defined as: Compared with the mean value, the fluctuation value is generally a small quantity. Then one has the following statistical relation: Substituting Equations (64-67) into Equation (63) leads to: Using Equation (69), Equation (71) was simplified as: Subtracting Equation (72) from Equation (70) leads to: Equations (72) and (73) where ( ) O δ stands for the value far less than (1) O , since the pulsative physical quantities investigated in this study is a small value.
Omitting those terms less than (1) O in Equation (72) and those less than ( ) O δ in Equation (73), the Equations (72) and (73) can be respectively simplified as: For the interpretation of temperature fluctuation in living tissues, it is reasonable to apply the 1-D degenerated forms of Equations (75) and (76). The boundary condition at the skin surface can be chosen as convective case which is often encountered in reality, i.e. ( ), 0 At the body core, a symmetrical or adiabatic condition can be used, namely Then using the relations in Equations (64-67), the boundary and initial conditions of Equations (75) and (76) can be respectively obtained as: where, 0 x = is defined as the skin surface while x L = the body core; f h is the apparent heat convection coefficient between the skin and the environment which is the contribution from the natural convection and radiation; f T the environment temperature; and 0 ( ) T x the initial temperature distribution. The following transformation was introduced [27] ( , ) Then Equations (75) and (79-81) were respectively rewritten as: where, K C α ρ = is the diffusivity of tissue, and 0 , 0 ( ) 1, 0 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 (86-89), the Green's function G can finally be obtained as: Equations (76) and (82-84) 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 R′ in Equation (94) can be obtained as where, the mean temperature ( , ) T x t was given by Equation (85) ran is random number generated by Fortran function) was adopted for calculations. k g m =°C) and the blood perfusion perturbation b W′ 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.

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, s1, s2 represent the distances between the origin and the boundaries of the tissue along x, y, z directions, respectively. 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 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. 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 where, f T is temperature of the cooling liquid; h is the convective heat transfer coefficient.
Here, the three planes, i.e. 0 x = ， 0 y = ， 0 z = are defined as adiabatic boundaries, and the other three planes, i.e. x=l, y=s1, z=s2, are subjected to forced convective cooling conditions. Without losing any generality, the initial temperature is defined as To solve for the Green's function of the above equations, the following auxiliary problem needs to be considered for the same region:  Finally, according to the above results and expression for the heat source term in Equation (107), the analytical solution to the temperature field under totally convective cooling conditions can then be obtained as: From the first term containing time in the above analytical solution, it can be seen that the thermal diffusivity α appears in the exponential term, i.e.
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 p T is the temperature of the cooling plate, using transformation 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 p T is equal to that of the cooling fluid f T . Then the Green's function solutions for the following two boundaries can be obtained accordingly.

Case
Considering that expressions for the transient temperature field for the above two cases still remain similar to that of Equation (117), 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 (117)

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 l1 was settled in the center. Then the location of the cryoprobe is at 0 / 2 x l = , 0 1 / 2 y s = , and the range for z coordinate is 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 c Q , i.e. c Q B = , which can also include many different discrete terms representing cooling effects from multiple cryoprobes and expressed as 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 (124) can be rewritten as: . 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: 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 0 ( , , ,0) ( , , ) T x y z T 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 where,  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.

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.
where, * c T is the vessel bulk temperature within cylindrical symmetry, and defined as 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 (134), whose solution can be rewritten as product of an axial term Z(z) and a radial one R0(r), i.e. * 0 ( , ) and R0(r) satisfying The solution to Equation (137)  The solution to Equation (138) can be obtained as Then the solution to Equation (134) can be expressed as where, m A , 0 A are coefficients. Substituting Equation (134) into Equation (143) leads to where, m A , 0 A are obtained as As to the exact temperature profile within blood vessel, if defining This constant 6 C does not exactly satisfy Equation (154) 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 6 C can also be made as follows. In Equation (154), if the axial position z is fixed, 6 C can be obtained through numerical iteration. And the relative error 6 6 6 [ ( ) ]/ C z C C − 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 6 C .
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 6 ( ) C z as a constant was an acceptable approximation. But for more complex situations where the above approximation cannot hold, exactly solving the Equation (157) for 6 ( ) C 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 (157)  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.

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,