Typical thermophysical properties of soft biological tissues .
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 . 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 , burn injury evaluation [3, 4], brain hypothermia resuscitation , disease thermal diagnostics , thermal comfort analysis , cryosurgery planning [8, 9], and cryopreservation programming . 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 . 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 . 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 .
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 , 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 . 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, , c are the density and the specific heat of tissue, respectively; and denote the density and the specific heat of blood, respectively; contains the Cartesian coordinates, and; denotes the analyzed spatial domain; is the space dependent thermal conductivity; and 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. is the blood temperature in the arteries supplying the tissue and is often treated as a constant at 37°C; is the tissue temperature; is the metabolic heat generation; and 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 . 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 , Chen-Holmes model addressing both the flow and perfusion properties of blood , and the Weinbaum-Jiji three-layer model to characterize the heat transfer in the peripheral tissues . 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 . 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 . 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 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.
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, is the time-dependent surface heat flux, is the time-dependent temperature of the cooling medium, and 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 () 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 and 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 and directions. However, it should be mentioned that more generalized boundary conditions in and directions can also be dealt with by the present approach, but they were not listed here for brevity.
The initial temperature is
where, 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, 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 the surrounding air temperature.
The solution to Equation (11) is:
where,. Through using the following transformation 
Equation (2) was transformed to the following form:
where, is the thermal diffusivity of tissue. The corresponding boundary and initial conditions are:
Using Green function method, can be solved from the combined Equations (14-24). The Green’s functions and to the second and third BCs can finally be obtained:
The Eigen-values are the positive roots of the following equation
Then, the solution of Equation (14) can be easily obtained. For the second BC at the skin surface, one has
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, is the point-heating power, is the Dirac function, 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.
|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/m3||1000|
|Density of blood (ρb)||Kg/m3||1000|
|Heat convection coefficient (h0)||W/m2· °C||10|
|Heat convection coefficient (hf)||W/m2· °C||100|
|Metabolic heat generation of tissue (Qm)||W/m3||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|
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 . 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 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 :
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. 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
One can rewrite the above equations (Equations (39-44) as:
The Green’s function for the above equation sets can be obtained as:
can be found by
can be calculated from
are the positive roots of the following equation:
Here Jm and Nm are the Bessel functions of the first kind and the second kind, respectively.
H(t-τ) is a heavy-side unit step function which has the following properties:
Finally, the temperature field was constructed as:
where, and 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 .
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 . 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 . 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, denotes the transient physical quantity at the vicinity of point (), and the temporal averaging period.
Compared with the mean value, the fluctuation value is generally a small quantity. Then one has the following statistical relation:
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: and 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 () and the pulsative blood perfusion, arterial blood temperature and metabolic heat generation (, , and) were correlated. Clearly, since the mean tissue temperature is a space dependent value, 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:
where stands for the value far less than, since the pulsative physical quantities investigated in this study is a small value.
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
where, is defined as the skin surface while the body core; is the apparent heat convection coefficient between the skin and the environment which is the contribution from the natural convection and radiation; the environment temperature; and the initial temperature distribution. The following transformation was introduced 
Then Equations (76) and (80-82) were respectively rewritten as:
where, is the diffusivity of tissue, and the Heaviside function.
If the Green’s function for the above equation system is obtained, its transient solution can thus be constructed . Through introducing an auxiliary problem corresponding to Equations (87-90), the Green’s function G can finally be obtained as:
where, the Eigen-values are the positive roots of the following equation
Then, the solution to Equation (87) can be obtained as
Substituting Equation (93) into Equation (86) leads to the mean temperature. The above solution is applicable to any transient environmental temperature and space-dependent metabolic heat generation. For the fluctuation temperature, the solving process is as follows. Using the similar transformation as Equation (86)
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, was given in Equation (91).
Consequently, the fluctuation variable in Equation (95) can be obtained as
where, the mean temperature was given by Equation (86). Substituting Equation (100) into Equation (94), the temperature fluctuation 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 . For clarity, only effect of the pulsative blood perfusion alone will be considered while the pulsative arterial temperature and metabolic heat generation were not taken into account, namely, , and a constant mean blood perfusion was assumed. Here, the pulsative blood perfusion was far less than the mean value, and its simplest form can be expressed as a periodic quantity with constant frequency and oscillation amplitude such as (where is frequency, the amplitude, and). However, to be more general, a stochastic pulsative perfusion form as (where, is random number generated by Fortran function) was adopted for calculations.
Fig. 4 depicted both the skin surface temperature fluctuation (the mean perfusion °C) and the blood perfusion perturbation 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 . 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.
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;, are the temperature distributions for solid and liquid phase, respectively; t is time;, 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 and. 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
where, s denotes the moving solid-liquid interface, 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
where n is the unit normal vector; L is the latent heat of tissue; is the phase change temperature; 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.
where, is the thermal diffusivity; is z coordinate of the arbitrary point p in the interface; is the delta function. For brief, Equation (106) can be rewritten as
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
where, is temperature of the cooling liquid; is the convective heat transfer coefficient. Here, the three planes, i.e.，，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
Using transformation, Equations (107), (109) and (110) can be rewritten as
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:
where,,are positive roots of the equation;, are positive roots of equation; and, are positive roots of.
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:
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 . 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 is the temperature of the cooling plate, using transformation, 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:
Then the transient temperature field can be constructed as:
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 is equal to that of the cooling fluid. 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
Case (d): Convective cooling at upside and underside and fixed temperature cooling at side faces
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 , 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 (where xi 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 xi 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 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 l1 was settled in the center. Then the location of the cryoprobe is at,, and the range for 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, , are density, specific heat of blood; is blood perfusion; is supplying arterial temperature; is volumetric metabolic heat; 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.
where, is the moving heat source term, which consists of three parts in a cryosurgical process: the metabolic heat source term, the heat sink term and the phase change heat source term, i.e.
where, , 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  and assumed to supply a constant cold amount, i.e., which can also include many different discrete terms representing cooling effects from multiple cryoprobes and expressed as (which reflects the amount of cold at the location of, where i is the sequence number of the cryoprobe; is the amount of cold released by the ith probe; 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:
where,. 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, 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:
Substituting it into Equation (127), one obtains
where, =. To simplify the equation, the volumetric moving heat source term in Equation (130) can be expressed as
The boundary conditions are rewritten as:
where, , and the initial condition is defined as. 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:
where,;;;; is the positive root of. Taking into account of the moving source term and the boundary conditions, i.e. Equation (132), one can obtain the transient temperature field of the tissue corresponding to Equation (130) as:
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 . 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 . 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 T0. 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, where Ta 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 
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 by solving the Navier-Stokes equation for an incompressible steady-state fluid, where a stands for the radius of blood vessel, and the average blood flow velocity over the vessel cross-section. Both temperature and heat flux on the vessel wall () satisfy continuity conditions. Defining the following dimensionless parameters:,, the governing equations for the vessel read as
where, 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 (135), whose solution can be rewritten as product of an axial term Z(z) and a radial one R0(r), i.e., with Z(z) satisfying
and R0(r) satisfying
The solution to Equation (138) is thus obtained as. is the positive roots of. Then
with normal as
The solution to Equation (139) can be obtained as
Then the solution to Equation (135) can be expressed as
where, , are obtained as
As to the exact temperature profile within blood vessel, if defining, one has from Equation (136)
Further, it can be derived as
where, is defined as, is the vessel inlet temperature at z=0. From Equation (136), one obtains =1. Therefore
At r=a, one has
Using the continuity condition for heat flux on the blood vessel wall, one approximately has
Exact calculation on from Equation (155) appears extremely difficult. However, if treating as a constant, and substituting Equation (156) and Equation (157) into Equation (155), then integrating it from z = 0 to L, one has
This constant does not exactly satisfy Equation (155) indeed. However it is a simple and intuitive approximation. According to former calculation , 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 can also be made as follows. In Equation (155), if the axial position z is fixed, can be obtained through numerical iteration. And the relative error 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. 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 as a constant was an acceptable approximation. But for more complex situations where the above approximation cannot hold, exactly solving the Equation (158) for 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 to reach the phase change point Tm, i.e.
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.
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.
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.