InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Mechanical Engineering » "Heat Transfer Phenomena and Applications", book edited by Salim N. Kazi, ISBN 978-953-51-0815-3, Published: October 24, 2012 under CC BY 3.0 license. © The Author(s).

Chapter 8

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

By Zhong-Shan Deng and Jing Liu
DOI: 10.5772/52963

Article top


Calculation geometry for 3-D case [13]
Figure 1. Calculation geometry for 3-D case [13]
Temperature distribution at cross-section z=0.04m after 1200s’ heating
Figure 2. Temperature distribution at cross-section z=0.04m after 1200s’ heating
D configuration of the prostate under microwave heating [26]
Figure 3. D configuration of the prostate under microwave heating [26]
Temperature fluctuation due to pulsative blood perfusion (Wb¯=0.5kg/m3°C)
Figure 4. Temperature fluctuation due to pulsative blood perfusion (Wb¯=0.5kg/m3°C)
Schematic of 1/8 cuboidal tissue subject to cryopreservation [29]
Figure 5. Schematic of 1/8 cuboidal tissue subject to cryopreservation [29]
Schematic of the living tissue domain subject to cryosurgery
Figure 6. Schematic of the living tissue domain subject to cryosurgery
Schematic of cylindrical tissues embedded with a single blood vessel [30]
Figure 7. Schematic of cylindrical tissues embedded with a single blood vessel [30]

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

Zhong-Shan Deng1 and Jing Liu1, 2

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, ρ, c are the density and the specific heat of tissue, respectively; ρband cb denote the density and the specific heat of blood, respectively; Xcontains the Cartesian coordinatesx, yandz; Ω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. Tais the blood temperature in the arteries supplying the tissue and is often treated as a constant at 37°C; T(X,t)is the tissue temperature; Qm(X,t)is the metabolic heat generation; and Qr(X,t) the distributed volumetric heat source due to externally applied spatial heating.

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

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

3. Bioheat transfer problems without phase change

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

Derivation of the solutions was based on the Green’s function method, since the Green’s function obtained for the differential equation is independent of the source term. Therefore it can be flexibly used to calculate the temperature distribution for various spatial or temporal source profiles. Furthermore, the Green’s function method is capable of dealing with the transient or space-dependent boundary conditions. Up to now, quite a few studies have applied the Green’s function method to solve the bioheat transfer problems [21-25]. However, in most of the existing analytical studies, the available solutions to the bioheat transfer problem are for the cases with one dimensional geometry, steady state, infinite domain, constant heating, or heat conduction equations not considering blood perfusion, which may not be practical for some real bio-thermal situations. In this section, the generalized analytical solutions, which have incorporated relatively complex situations such as the 3-D tissue domain, the transient or space-dependent boundary conditions, and volumetric heating, were especially addressed. Such solutions are expected to be very useful in a variety of bio-thermal practices. The 3-D computational domain with widths 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.


Figure 1.

Calculation geometry for 3-D case [13]

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


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




where, f1(y,z;t)is the time-dependent surface heat flux, f2(y,z;t)is the time-dependent temperature of the cooling medium, and hf 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 (Tc) 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, T0(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, h0is 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 Tf the surrounding air temperature.

The solution to Equation (11) is:


where,A=ωbρbcb/k. Through using the following transformation [13]


Equation (2) was transformed to the following form:


where, α=k/ρcis the thermal diffusivity of tissue. The corresponding boundary and initial conditions are:




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




The Eigen-values βp 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, P1(t)is the point-heating power, δis the Dirac function, (x0,y0,z0)is the position of the point-heating source.

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


Figure 2.

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

Air temperature (Tf)°C25
Artery blood temperature (Ta)°C37
Blood perfusion of tissue (ωb)ml/s/ml0.0005
Body core temperature (Tc)°C37
Density of tissue (ρ)Kg/m31000
Density of blood (ρb)Kg/m31000
Heat convection coefficient (h0)W/m2· °C10
Heat convection coefficient (hf)W/m2· °C100
Metabolic heat generation of tissue (Qm)W/m333800
Specific heat of tissue (c)J/Kg· °C4200
Specific heat of blood (cb)J/Kg· °C4200
Temperature of cooling medium (f2)°C15
Thermal conductivity of tissue (k)W/m· °C0.5

Table 1.

Typical thermophysical properties of soft biological tissues [13].

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

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

The transurethral microwave catheter (T3 catheter) was used as the heating apparatus in this section. Geometric presentation of the prostate with the inserted T3 catheter was shown in Fig. 3. It was modeled as a cylinder of 3.4cm in diameter and 3cm in length with constant temperature T 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:


Figure 3.

D configuration of the prostate under microwave heating [26]

To obtain an analytical solution, all these parameters were assumed to be uniform throughout the prostate and remained constant except for ωb which varied with respect to the heat power. The cooling effect from the chilled water running inside the catheter was modeled by an overall convection coefficient h. The external boundary at the capsule was prescribed as the body core temperatureT. 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-44) as:




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




βlcan be found by


μncan be calculated from


μn*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, Q*(r,θ,z)and F(r,θ,z) were given in Equations (51-52), respectively.

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

3.3. Analytical solutions to bioheat transfer with temperature fluctuation

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


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, A(x,y,z;t)denotes the transient physical quantity at the vicinity of point (x,y,z), 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:


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




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: CbWbTa¯and CbWbT¯ represent the mean transferred energy due to perfusion perturbations and temperature fluctuations in tissue and arterial blood, respectively. It is from Equation (74) that the temperature fluctuation (T) and the pulsative blood perfusion, arterial blood temperature and metabolic heat generation (W, Ta, andQm) were correlated. Clearly, since the mean tissue temperature T¯ is a space dependent value, Tis 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 O(δ) stands for the value far less thanO(1), since the pulsative physical quantities investigated in this study is a small value.

Omitting those terms less than O(1) in Equation (73) and those less than O(δ) in Equation (74), the Equations (73) and (74) 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 (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, x=0is defined as the skin surface while x=L the body core; hfis the apparent heat convection coefficient between the skin and the environment which is the contribution from the natural convection and radiation; Tfthe environment temperature; and T0(x) the initial temperature distribution. The following transformation was introduced [27]


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


where, α=K/ρCis the diffusivity of tissue, and H(t)={0,t<01,t>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 (87-90), the Green’s function G can finally be obtained as:


where, the Eigen-values βn 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 temperatureT(x,t)¯. The above solution is applicable to any transient environmental temperature Tf(t) and space-dependent metabolic heat generationQm(x,t)¯. 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, GR¯(x,t;ξ,τ)was given in Equation (91).

Consequently, the fluctuation variable R in Equation (95) can be obtained as


where, the mean temperature T(x,t)¯ was given by Equation (86). Substituting Equation (100) into Equation (94), the temperature fluctuation T(x,t) can be obtained. As illustration, two calculation examples using the above analytical solutions were presented in this section to investigate the temperature fluctuation phenomena in living tissues. In the following illustration calculations, the typical values for tissue properties were taken from [27]. For clarity, only effect of the pulsative blood perfusion Wb alone will be considered while the pulsative arterial temperature and metabolic heat generation were not taken into account, namely, Ta=0,Qm=0, and a constant mean blood perfusion Wb¯ was assumed. Here, the pulsative blood perfusion Wb was far less than the mean valueWb¯, and its simplest form can be expressed as a periodic quantity with constant frequency and oscillation amplitude such as Wbmcosωt (where ω is frequency, Wbmthe amplitude, andWbm<<Wb¯). However, to be more general, a stochastic pulsative perfusion form as Wb=0.1Wb¯(0.5ran) (where, ranis random number generated by Fortran function) was adopted for calculations.


Figure 4.

Temperature fluctuation due to pulsative blood perfusion (Wb¯=0.5kg/m3°C)

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

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

4. Bioheat transfer problems with phase change

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

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


Figure 5.

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

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

k2Tl(x,y,z,t)=ρcTl(x,y,z,t)t,in regionRl,t>0

For the solid phase:

k2Ts(x,y,z,t)=ρcTs(x,y,z,t)t,in regionRs,t>0

wherek, ρ, c denote thermal conductivity, density and specific heat of tissue, respectively;Ts, Tlare the temperature distributions for solid and liquid phase, respectively; t is time;Rl, Rsdenote solid and liquid region; the subscript l, s represent the liquid and solid phase, respectively. To obtain an analytical solution, all the parameters were assumed as uniform and remain constant.

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


where, s denotes the moving solid-liquid interface, (x0,y0,z0)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

kTsnkTln=ρLst|n at S(x,y,z,t)=0,t>0

where n is the unit normal vector; L is the latent heat of tissue; Tmis the phase change temperature; vn(t)=st|nrepresents 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, α=k/ρcis the thermal diffusivity; z0is z coordinate of the arbitrary point p in the interface; δ(zz0)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, Tfis temperature of the cooling liquid; his the convective heat transfer coefficient. Here, the three planes, i.e.x=0y=0z=0 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θ=TTf, 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:

G(x,y,z,t;ξ,η,λ,τ)=k=0n=0j=0cos(βjx)cos(γny)cos(μkz)cos(βjξ)cos(γnη)cos(μkλ)N(βj)N(γn)N(μk)×H(tτ)exp[α(βj2+γn2+μk2)(tτ)] j=0,1,2; n=0,1,2; k=0,1,2,

where,N(βj)=l[βj2+(h/k)2]+h/k2[βj2+(h/k)2] ,βj are positive roots of the equationβjtan(βjl)=h/k;N(γn)=s1[γn2+(h/k)2]+h/k2[γn2+(h/k)2], γnare positive roots of equationγntan(γns1)=h/k; andN(μk)=s2[μk2+(h/k)2]+h/k2[μk2+(h/k)2], μkare positive roots ofμktan(μks2)=h/k.

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


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.exp[α(βj2+γn2+μk2)(tτ), 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 Tp is the temperature of the cooling plate, using transformationθ=TTp, 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:

G(x,y,z,t;ξ,η,λ,τ)=k=0n=0j=08cos(βjx)cos(γny)cos(μkz)cos(βjξ)cos(γnη)cos(μkλ)s1s2l×H(tτ)exp[α(βj2+γn2+μk2)(tτ)] j=0,1,2; n=0,1,2; k=0,1,2,

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 Tp is equal to that of the cooling fluidTf. 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

G(x,y,z,t;ξ,η,λ,τ)=k=0n=0j=02cos(βjx)cos(γny)cos(μkz)cos(βjξ)cos(γnη)cos(μkλ)N(βj)N(γn)s2×H(tτ)exp[α(βj2+γn2+μk2)(tτ)] j=0,1,2; n=0,1,2; k=0,1,2,

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

G(x,y,z,t;ξ,η,λ,τ)=k=0n=0j=04cos(βjx)cos(γny)cos(μkz)cos(βjξ)cos(γnη)cos(μkλ)N(μk)s2l×H(tτ)exp[α(βj2+γn2+μk2)(tτ)] j=0,1,2; n=0,1,2; k=0,1,2,

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 xi[T(x,y,z,t)t] (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 T(x,y,z,t)/t can eliminate the time integral term in Equations (118) and (120), which would simplify the solution form. Overall, the present analytical method in virtue of its straightforward form is of great significance to evaluate the phase change problem in cryobiology.

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

Cryosurgery is very different from cryopreservation, since living tissue has to be considered. Consequently, it must take into account the effects of blood perfusion and metabolic heat generation into bioheat equation. Here, the Pennes equation is applied to characterize the heat transfer process in the living tissue. To avoid the complex boundary conditions, the calculation tissue domain is chosen as a whole cuboid as shown in Fig. 6, where a cryoprobe with length l1 was settled in the center. Then the location of the cryoprobe is atx0=l/2,y0=s1/2 , and the range for z coordinate is0.5(s2l1)z0.5(s2+l1).


Figure 6.

Schematic of the living tissue domain subject to cryosurgery

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

For the liquid phase

ρcTlt=k2Tl+ωbρbcb(TaTl)+Qm+Qc,in regionRl,t>0

For the solid phase

ρcTs(x,y,z,t)t=k2Ts(x,y,z,t)+Qc,in regionRs,t>0

where, ρb, cbare density, specific heat of blood; ωbis blood perfusion; Tais supplying arterial temperature; Qmis volumetric metabolic heat; Qcis 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, Qris the moving heat source term, which consists of three parts in a cryosurgical process: the metabolic heat source termQm, the heat sink term Qcand the phase change heat source termqr, i.e.


where, qm'=wbcbρb(TaT)+qm, reflecting contributions of the blood heat transfer and the metabolic heat generation in unfrozen region. Other parameters and functions have the same definitions with those of cryopreservation.

To simplify the problem, the cryoprobe inserted into the deep tissue is treated as a linear heat sink [29] and assumed to supply a constant cold amountQc, i.e.Qc=B, which can also include many different discrete terms representing cooling effects from multiple cryoprobes and expressed as Qc=Qiδ(xxi,yyi,zzi) (which reflects the amount of cold at the location of(xi,yi,zi), where i is the sequence number of the cryoprobe; Qiis 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,g(z,t)=[(ωbcbρbTa+qm)H(zz0)+qr+Qc]/ρc. In a cryosurgical procedure, the cryoprobe is inserted into the target tissue, which will subject to a specific temperature drop due to cooling of heat sink effect of the probe. Since the skin surface is exposed to the ambient environment, the boundary conditions can thus be treated as:


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 asT(x,y,z,0)=T0(x,y,z), and every boundary is defined as adiabatic at this initial state. Although the control equations and boundary conditions in cryosurgery are very different from that in cryopreservation, the solution procedures still remain similar if certain transformations are introduced. In the following, only those steps different from the above will be addressed. To make solution of the problem feasible, we have adopted before the following specific transformation:


Substituting it into Equation (127), one obtains


where, g1(z,t)={[ωbcbρb(TaT0)+qm]H(zz0)+qr+Qc}/ρc. To simplify the equation, the volumetric moving heat source term in Equation (130) can be expressed as


The boundary conditions are rewritten as:


where, f(t)=(TfT0)exp[ωbρbcbρcH(zz0)t]H(t), and the initial condition is defined asW(x,y,z,t)=0. Applying the same theoretical strategies as used in solving cryopreservation into the above problem, one can obtain the Green’s function expression for a 3-D cryosurgical process as:

G(x,y,z,t;ξ,η,λ,τ)=1s1s2lk=0n=0j=0R1R22[μk2+(h/k)2][μk2+(h/k)2]+h/kexp[α(βj2+γn2+μk2)(tτ)]×H(tτ)cos(βjx)cos(γny)cos(μkz)cos(βjξ)cos(γnη)cos(μkλ) j=0,1,2; n=0,1,2; k=0,1,2,

where,R1={1,j=02,j=1,2,3,;R2={1,n=02,n=1,2,3,;βj=jπl;γn=nπs1; μkis the positive root ofμktan(μks2)=h/k. Taking into account of the moving source term qr* 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 [30]. Anatomically, tumors are often situated close to or embedded with large blood vessels, since a tumor’s quick growth ultimately depends on nutrients supplied by its blood vessel network. During cryosurgery, the blood flow inside a large vessel represents a source which heats the nearby frozen tissues and, thereby, limits freezing lesions during cryosurgery. Under this condition, a part of the vital tumor cells may remain in the cryolesion and lead to recurrence of tumors after cryosurgical treatment. More specifically, tumor cell survival in the vicinity of large blood vessels is often correlated with tumor recurrence after treatment [30]. Consequently, it is difficult to implement an effective cryosurgery when a tumor is contiguous to a large blood vessel. To better understand the effect of blood flow to the temperature distribution of living tissues subject to freezing, a conceptual model for characterizing the heat transfer in 3-D cylindrical tissues embedded with a single blood vessel was illustrated in this section. And a closed form analytical solution to this model was provided to explore different factors’ thermal influences to the freezing mechanism of living tissues.

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


Figure 7.

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

To carry out the theoretical analysis, additional assumptions for tissues were made as follows: there is no heat flow across the boundaries at z=0 and z=L. In reality, if the vessel length is long enough, these adiabatic boundary conditions can be closely satisfied. And the cylindrical tissue surface (r=b) was assumed to be immersed into a cooling medium or subjected to a circular cryoprobe with constant freezing temperature 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 parameterTt*=(TtT0)/(TaT0), 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 [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 u=u0(1r2a2) by solving the Navier-Stokes equation for an incompressible steady-state fluid, where a stands for the radius of blood vessel, and u0 the average blood flow velocity over the vessel cross-section. Both temperature and heat flux on the vessel wall (r=a) satisfy continuity conditions. Defining the following dimensionless parameters:Tb*=TbT0TaT0,Tc*=TcT0TaT0 , the governing equations for the vessel read as


where, Tc*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.Tt*(r,z)=R0(r)Z(z), with Z(z) satisfying


and R0(r) satisfying


The solution to Equation (138) is thus obtained asZ(z)=C1cos(βmz)+C2sin(βmz). βmis the positive roots ofsin(βL)=0. Then




with normal as


The solution to Equation (139) can be obtained as


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


where, Am,A0 are coefficients. Substituting Equation (135) into Equation (144) leads to


where, Am, A0are obtained as


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


As to the exact temperature profile within blood vessel, if definingU=u0ρbCb4KbdTc*dz, one has from Equation (136)


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


Further, it can be derived as


where, Tc0*is defined asTc0*=4a20aTb0*(1r2a2)rdr, Tb0*is the vessel inlet temperature at z=0. From Equation (136), one obtains Tb0*=1. Therefore

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




Exact calculation on C6(z) from Equation (155) appears extremely difficult. However, if treating C6 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 C¯6 does not exactly satisfy Equation (155) indeed. However it is a simple and intuitive approximation. According to former calculation [30], it can be found that heat fluxes calculated using temperature from both sides of the blood vessel wall were almost identical. Therefore, an approximate estimation on C6 can also be made as follows. In Equation (155), if the axial position z is fixed, C6can be obtained through numerical iteration. And the relative error [C6(z)C¯6]/C¯6was 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 valueC¯6. However, since the present discussion was mainly focused on the case of large blood vessel (for tissues with extremely small vessel, a collective model such as Pennes equation is often applicable), therefore treating C6(z) as a constant was an acceptable approximation. But for more complex situations where the above approximation cannot hold, exactly solving the Equation (158) for C6(z) is still very necessary in the future.

When analyzing the thermal effect of large blood vessel in cryosurgery, an important issue is when the blood vessel begins to freeze and how to control cryoprobe’s temperature to completely freeze the target tumor. Substituting Equation (158) into Equation (153), one can set up the relation for the blood vessel temperature Tb(a,L) 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.

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.




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.


1 - J. Liu, C. C. Wang, 1997Bioheat Transfer (in Chinese). Beijing: Science Press. 435 p.
2 - J. Liu, Z. S. Deng, 2008Physics of Tumor Hyperthermia (in Chinese). Beijing: Science Press. 381 p.
3 - J. Liu, X. Chen, L. X. Xu, 1999New Thermal Wave Aspects on Burn Evaluation of Skin Subjected to Instantaneous Heating. IEEE Transactions on Biomedical Engineering 46420428
4 - Y. G. Lv, J. Liu, J. Zhang, 2006Theoretical Evaluation on Burn Injury of Human Respiratory Tract due to Inhalation of Hot Gas at the Early Stage of Fires. Burns 32436446
5 - J. Liu, 2007Cooling Strategies and Transport Theories for Brain Hypothermia Resuscitation. Frontiers of Energy and Power Engineering in China l: 3257
6 - Z. S. Deng, J. Liu, 2004Computational Study on Temperature Mapping Over Skin Surface and Its Implementation in Disease Diagnostics. Computers in Biology and Medicine 34: 2004.
7 - Fanger PO1970Thermal Comfort- Analysis and Applications in Environmental Engineering. New York: McGraw-Hill. 244 p.
8 - J. Liu, 2007Principles of Cryogenic Biomedical Engineering (in Chinese). Beijing: Science Press. 338 p.
9 - Diller KR1992Modeling of Bioheat Transfer Processes at High and Low Temperatures. Advances in Heat Transfer 22157357
10 - Hua TC, Ren HS1994Cryogenic Biomedical Technology (in Chinese). Beijing: Science Press. 412 p.
11 - Roemer RB1999Engineering Aspects of Hyperthermia Therapy. Annual Review of Biomedical Engineering. 1347376
12 - AA Baust. J. Gage, 1998Mechanism of Tissue Injury in Cryosurgery. Cryobiology. 37171186
13 - Z. S. Deng, J. Liu, 2002Analytical Study on Bioheat Transfer Problems with Spatial or Transient Heating on Skin Surface or Inside Biological Bodies. ASME Journal of Biomechanical Engineering. 124638649
14 - Pennes HH1948Analysis of Tissue and Arterial Blood Temperatures in the Resting Human Forearm. Journal of Applied Physiology. 193122
15 - J. Liu, 2006Bioheat Transfer Model. In: Akay M, editor. Wiley Encyclopedia of Biomedical Engineering. John Wiley & Sons. 429438
16 - W. Wulff, 1974The Energy Conservation Equation for Living Tissues. IEEE Transactions on Biomedical Engineering. 21494497
17 - Chen MM, Holmes KR1980Microvascular Contributions in Tissue Heat Transfer. Annals of the New York Academy of Sciences. 335137150
18 - S. Weinbaum, L. M. Jiji, 1985A New Simplified Bioheat Equation for the Effect of Blood Flow on Local Average Tissue Temperature. ASME Journal of Biomechanical Engineering. 107131139
19 - J. Liu, Z. S. Deng, 2009Numerical Methods for Solving Bioheat Transfer Equations in Complex Situations. In: Minkowycz WJ, Sparrow EM, Abraham JP, editors. Advances in Numerical Heat Transfer (3Taylor & Francis. 75120
20 - Z. S. Deng, J. Liu, 2002Monte Carlo Method to Solve Multi-Dimensional Bioheat Transfer Problem. Numerical Heat Transfer, Part B: Fundamentals. 42543567
21 - R. Vyas, M. L. Rustgi, 1992Green’s Function Solution to the Tissue Bioheat Equation. Medical Physics. 1913191324
22 - B. Gao, S. Langer, P. M. Corry, 1995Application of the Time-Dependent Green’s Function and Fourier Transforms to the Solution of the Bioheat Equation. International Journal of Hyperthermia. 11267285
23 - Durkee JW, Antich PP, Lee CE1990Exact-Solutions to the Multiregion Time-Dependent Bioheat Equation- Solution Development. Physics in Medicine and Biology. 35847867
24 - Durkee JW, Antich PP1991Exact-Solution to the Multiregion Time-Dependent Bioheat Equation with Transient Heat-Sources and Boundary-Conditions. Physics in Medicine and Biology. 36345368
25 - Kou HS, Shih TC, Lin WL2003Effect of the Directional Blood Flow on Thermal Dose Distribution during Thermal Therapy: An Application of a Green’s Function Based on the Porous Model. Physics in Medicine and Biology. 4815771589
26 - J. Liu, L. Zhu, L. X. Xu, 2000Studies on the Three-Dimensional Temperature Transients in the Canine Prostate during Transurethral Microwave Thermal Therapy. ASME Journal of Biomechanical Engineering. 122372379
27 - Z. S. Deng, J. Liu, 2001Blood Perfusion Based Model for Characterizing the Temperature Fluctuation in Living Tissues. Physica A: Statistical Mechanics and Its Applications. 300521530
28 - J. Liu, Y. X. Zhou, 2002Analytical Study on the Freezing and Thawing Processes of Biological Skin with Finite Thickness. Heat and Mass Transfer. 38319326
29 - F. F. Li, J. Liu, K. Yue, 2009Exact Analytical Solution to Three-Dimensional Phase Change Heat Transfer Problems in Biological Tissues Subject to Freezing. Applied Mathematics and Mechanics. 306372
30 - Y. T. Zhang, J. Liu, Y. X. Zhou, 2002Pilot Study on Cryogenic Heat Transfer in Biological Tissues Embedded with Large Blood Vessels. Forschung im Ingenieurwesen- Engineering Research. 67188197