Bio-Heat Models Revisited: Concepts, Derivations, Nondimensalization and Fractionalization Approaches

The heat transfer in living tissues is an evergreen problem in mathematical modeling with great practical importance starting from the Pennes equation postulation. This study focuses on concept in model building, the correct scaling of the bio-heat equation (one-dimensional) by appropriate choice of time and length scales, and consequently order of magnitude analysis of effects, as well as fractionalization approaches. Fractionalization by different constitutive approaches, leading to application of different fractional differential operators, modeling the finite speed of the heat wave, is one of the principle problems discussed in the study. The correct choice of the damping (relaxation) function of the heat flux is of primarily importance in the formulation of the bio-heat equation with memory. Moreover, this affects the consequent scaling, order of magnitude analysis and solutions.

Model building is the art of selecting those aspects of a process that are relevant to the question being asked - J.H Holland.

(Holland, JH (1995) Hidden Order.
Addison-Wesley, New York, USA.)

1. Introduction

All biological bodies live in space-depended temperature fields (environments), that is even in single organism there are no uniform temperatures of organs and tissues [1]. The non-uniformity of the temperature fields causes energy transfer between the organs (and tissues) with the perfused blood as wells at the body interface trough the skin [1].

Heat transfer in living tissues is a complex process and the mathematical modeling needs simplifications allowing to apply its methods. Hereafter, we will consider, in general, human bio-heat problems where the core body temperature is maintained at a level of 36–37°C. Cases related to organisms which body temperature varies as the external temperature changes (the so-called poikilothermic organisms [2]) are out of the scope of this work. In general, the heat exchange between a living body and the environment, taking place through the skin, involves mechanisms of conduction, convection,radiation and perspiration (consisting of separation and evaporation of sweat) [3, 4]. In medical applications, for instance, heat transfer procedures such as cryosurgery, radio-frequency and laser thermal ablations, thermal resections, ultrasound, etc. [5–7] need not only experimental data but also deep understanding of the basic mechanisms of heat transport and adequate mathematical modeling.

The human body temperature is affected by many factors, such as: environmental conditions, temperatures of media surrounding tissues, metabolism in muscles and blood circulation. We will briefly comment two heat transfer processes involved in the basic models discussed further in this work:

(i) Heat transfer through the skin (section 1.1), and

(ii) Blood circulation and its contribution to the overall heat exchange in the body (section 1.2). Moreover, basic steps in formalization of the heat fluxes due to different causes such as conduction, convection and metabolic heat generation are commented prior to the analysis of the bio-heat models (section 1.3).

1.1. Heat Exchange Through the Skin

Skin is the largest single organ of the body [1] enabling protection from the surrounding environment and playing significant role in thermoregulation, sensory and host defense functions. The primary role of thermoregulation compromises heat generation, absorption, radiation, transmission and vaporization [1]. That is, the primary process occurring at the skin (and through the skin) is the thermal energy exchange with the environment. In humans the skin properties (structure and thickness) depend on the location and its function [1, 8]. The skin includes two layers epidermis and dermis [7, 9]. The former allows heat release by convection, radiation and perspiration (see Figure 1). In this context, the typical forearm skin temperature is about 31°C [7, 10] while the bone temperature is about 33°C. Hence, the increase in the temperature of a deep tissue requires higher temperatures to be imposed at the skin surface. The thermoregulation of the body, by help of the skin functions, is mainly by changing the blood flow rate in the micro-circulatory bed consisting of arterioles, arterial and vein capillaries, and venules [1]. In this context, an increase in the skin temperature to 42°C results in increased blood circulations (see below) in order to dissipate a portion of the thermal energy and to permit increase in the heat conduction mechanism through the deep tissue beneath [10, 11]. Last but not least, body fat affects strongly both the skin and deep tissue temperatures when the human body is at rest [11].

FIGURE 1

Figure 1. Schematic presentations of human body layers: see explanation in the text.

The common thermophysical parameters of the human skin are (collected from different sources) [1, 12–14]: density (ρ = 1.000 − 1.005[kg/m³]), heat conductivity (k = 0.5 − 0.628[W/mK]), heat capacity (c = 4.187 − 4.200[Jkg/K]). As it is shown by many authors (see data collected in [15]) the interior tissue temperature attains a constant value at a short distance from the surface, commonly this happens at a depth of 2–3 cm. Because of that in many studies the finite length of the tissue volume under consideration, commonly denoted as L, is within this range.

1.2. Blood Circulation and Heat Exchange at a Glance

Blood perfused tissues structures include layers of skin, fats, muscles and bones (see the schematic presentation vessels (see Figure 2A) (after [9]). Moreover, the blood circulation, trough the two principle sets of arteries and veins, is the principle mechanism in the body thermal regulation (see Figure 2B): Blood leaving the heart via the aorta (the largest blood vessel with a diameter of about 5000 μm is transported to the muscles through the arteries and the veins with diameters within the range of 300–1,000 μm: primary arteries with diameters of 100–300 μm are supplying blood to the secondary arteries 50–100 μm in diameter. Then, the blood is delivered to arterioles 20–40 μm connected to the smallest transport vessels, the capillaries [9]. The returning blood loop to the heart consists of veins.

FIGURE 2

Figure 2. Schematic presentations of human body circulatory system (A) and blood temperature distributions through the vessels of the circulatory system (B). See comments in the text.

The blood is leaving the aorta (with a temperature T_a) practically does not exchange energy with the tissue up to the point where it enters the primary arteries where thermal equilibration with the surrounding tissue begins; this process completes before the points where the blood enters the arterioles and the capillaries. Further, the blood flow attains the temperature of the perfused tissue T_t and its value does not vary too much (if it is not mixed with other blood streams) [9]. The chilled blood from the peripheral tissues and the hot blood from the internal organs are mixed in the vena cava and the right atrium, and ventricle. The next heat transfer steps of the blood circulation are the pulmonary circulation and remixing in left heart chamber: this allows the blood to attain the temperature at which it left the heart at the beginning of its flow through the circulation system [9, 16].

The thermoregulation due to blood circulation and perfusion, in accordance with prevailing theories, is to maintain the body core temperature at a constant value required for the normal physiological functions [17]. Alternatively, the blood circulations should maintain the body energy balance and therefore the tissue temperatures are results, not the causes for the thermoregulation process [18, 19].

The common blood parameters involved in the bio-heat models are Torvi and Dale [12], Yamada et al. [13], and Askarizadeh and Ahmadikia [14]: temperature ($T_b=37^{o}C$), metabolic heat generation ($q_m=1.19\times 10^{3}W/m^3$) perfusion rate (volumetric rate per unit tissue volume, i.e., ${\omega }_b=Q_{blod}/V_{tissue}=m^3{s}^{-1}/m^3$ (${\omega }_b=1.87\times 10^{-3}{s}^{-1}$).

Note 1: Sometimes the blood perfusion rate is repressed as mass blood flux per unit mass of the tissue, that is ${\omega }_b=G_{blood}/V_{tissue}=kg/m^{3}s$ [15]. The reader should take care what is the dimension used in a specific study, that mainly affects the nondimensalization procedures.

1.3. Bio-Heat Modeling Formalization

Analyzes of thermal energy balances in various biological tissues are multifaceted due to differences in tissue structures and functions. The relative importance of a certain heat transfer mechanism, the relevant time scale of the deposited energy and as well as differences in boundary and initial conditions strongly affect the energy balance [20]. In modeling, commonly simplifying assumptions are necessary thus making possible to model the basic properties of the thermal status of the organism (or its parts) and the effects of the boundary and initial conditions. The build-up of such models, in general, begins with application of the energy balance law over a control volume [20]

$\begin{array}{ll}\frac{d}{dt}{Q}_{gain}=\frac{d}{dt}{Q}_{storage}+\frac{d}{dt}{Q}_{loss}+\frac{d}{dt}W& (1)\end{array}$

That is, the rate of heat gained is balanced by the heat stored and heat loss (by conduction or heat exchange with advanced fluids), and the rate of work W performed by the tissue [20]. Precisely, the heat gain (2) due to heat generated by unit tissue segment q(r, t) and heat energy storage (3) can be expressed as integral over the control volume V, namely

$\begin{array}{ll}\frac{d}{dt}{Q}_{gain}={\displaystyle \underset{V}{\int }}q(r,t)dV& (2)\end{array}$

$\begin{array}{ll}\frac{d}{dt}{Q}_{storage}={\displaystyle \underset{V}{\int }}\rho c\left[\frac{\partial }{\partial t}T(r,t)\right]dV& (3)\end{array}$

Further, we will restrict the heat transfer mechanisms to only conduction and convection in order to make clearer the following analyzes of the models, their scaling and consequently generated dimensionless groups, and the fractionalization when the thermal inertial should be taken into account. In bio-heat transfer problems the heat conduction mechanism is commonly modeled by the Fourier law when no inertia is taken into account, that is

$\begin{array}{ll}{Q}_{conduction}=-\frac{k_t}{\Delta L_t}{A}_t(T_1-T_2),T_1>T_2& (4)\end{array}$

where k_t is the tissue heat conductivity with dimension [W/m · K], which is usually temperature-dependent; this transport coefficient strongly depends on the microscopic structures of the biological materials. In differential form the conduction component of the heat flux through the issue (per unit area) [W/m²] is

$\begin{array}{ll}{q}_{cond}=-k_t\frac{\partial }{\partial x}{T}_t(x,t)& (5)\end{array}$

while the convection (perfusion) fluid-body term is

$\begin{array}{ll}{q}_{conv}=h(T_w-T_{fluid})& (6)\end{array}$

The heat transfer coefficient (fluid-solid surface) h [W/m² · K] is dependent on the fluid velocity. In biological tissue the fluid-related heat transfer should be contributed to the blood perfusion and flow distribution due the heterogeneity of the living structures. Mathematically representing the blood perfusion component of the bio-heat thermal flux (through unit area of the tissue) q_b [W/m²] the structure of the relationship is the same as in the case of heat transfer at the solid-fluid boundary (as in the case of simple blood flow in a cylindrical vessel commented in section 2.2). However, now the flux is proportional to the difference between the blood temperature in the artery (denoted as T_art or T_a) and the venous level temperature T_ven, namely

$\begin{array}{ll}{q}_b={\omega }_b{\rho }_b{c}_b(T_{art}-T_{ven})& (7)\end{array}$

The choice of this temperature difference is a matter of arguments as we will see in the sequel when model constructions will be discussed. However, at a glance, these are the only practically measurable temperatures when thermal sensors have to be introduced in the living organism, despite the physical models assumptions considering another temperature differences. However, it could be considered that in capillary beds the very slow blood flows attain complete thermal equilibria with the tissues; consequently the perfusion associated heat flux (7) could be approximated as

$\begin{array}{ll}{q}_b={\omega }_b{\rho }_b{c}_b(T_{art}-T_t)& (8)\end{array}$

Hence, under such circumstances, the energy transfer by the blood stream, over the entire control volume, is (in original notations)

$\begin{array}{ll}\frac{d}{dt}{Q}_b={\displaystyle \underset{V}{\int }}{\omega }_b{\rho }_b{c}_b\left[T_{art}(r,t)-T_t(r,t)\right]dV& (9)\end{array}$

The result (9) is a not unconditional, despite the fact that it is frequently encountered in model formulations; in larger arteries it would be not valid due to strong blood mixing [20].

After considering the different heat transport mechanisms we may construct the balance of thermal energy (neglecting the work performed by the tissue) over an arbitrary volume element as

$\begin{array}{ll}\begin{array}{c}\hfill {\displaystyle \underset{V}{\int }}{\rho }_t{c}_t\frac{\partial }{\partial t}{T}_t(r,t)dV={\displaystyle \underset{V}{\int }}{k}_t\nabla T_t(r,t)dV+\\ \hfill +{\displaystyle \underset{V}{\int }}{\omega }_b{\rho }_b{c}_b\left[T_{art}(r,t)-T_t(r,t)\right]dV+{\displaystyle \underset{V}{\int }}{Q}_m(r,t)dV\end{array}& (10)\end{array}$

In one dimensional case and homogeneous distributed source of metabolic heat Q_m we get

$\begin{array}{l}{\rho }_t{c}_t\frac{\partial }{\partial t}{T}_t(x,t)=k_t\frac{{\partial }^2}{\partial x^2}{T}_t(x,t)+{\omega }_b{\rho }_b{c}_b\\ \left[T_{art}(x,t)-T_t(x,t)\right]+Q_m& (11)\end{array}$

This is the standard form of the heat transfer through living organisms (tissues) or bio-heat transfer equation. We will discuss it components and assumptions behind, which may vary under circumstances imposed by the modeled object, as well the contribution of the transport mechanism involved by scaling and dimensional analysis. The study addresses three principle issues in bio-heat modeling:

(i) Models building concept,

(ii) Model scaling and nondimensalization and

(iii) Model fractionalization by involving fractional operators.

1.4. Motivation of This Study

Working in the modern direction of modeling involving fractional operators the author was intrigued by development of so-called fractional bio-heat models and how they are related to the classical ones such as Pennes's, Wulff's, Klinger's, Chen-Holmes' models, etc. [20, 21]. The literature search revealed that the dominating approach is to fractionalize formally the Pennes equation. Moreover, the correct solution of models either based on local or non-local operators needs correct scaling and nondimensalization. This formulates a special interest about the time and length scales, velocity scale and finally the dimensionless groups appearing as pre-factors in the modeling differential equations; as weighting factors and allowing to see the contributions of different physical effects accounted by the model of interest.

The classical models mentioned above are based on the continuum approach allowing easily to apply the local (integer-order) calculus and as fundamental base the Fourier law. Therefore, in all these models, parabolic by nature, the heat flux speed is infinite, which is unphysical. The flux damping in models involving local derivatives are based on extension of the Fourier law by applying first order approximation by Taylor series expansion; the Maxwell-Cattaneo and the Dual-phase-lag approaches resulting in hyperbolic models are commonly applied.

Therefore, the main motivation of this study is to revisit the existing bio-heat models, precisely thus considering blood perfused tissue, and to attempt representing them from a common basis thus allowing to see what are the links between them and how the modern trends, such as the fractional models should be correctly constructed. Therefore, to be precise, this study considers the following problems related to 1-D bio-heat transfer models:

(1) Basic concepts of model build-up and the relevant physical and mathematical assumptions.

(2) Scaling and nondimensalization, thus detecting the controlling dimensionless groups.

(3) Fractionalization: (a) Critical analysis of the existing fractional models, (b) A systematic approach based on fractional Taylor series expansion and the fading memory concept.

Analytic and numerical solutions of existing bio-heat transfer models are out of scope of this work, even though we will analyze them to some extent in the context of the three main directions declared above. Comprehensive analysis of analytical solutions of existing bio-heat models (integer-order) with/without flux damping are available in Xu et al. [1], Datta [16], Valvano [17], Orr and Eberhart [20], Charny [21], and Fan and Wang [22].

1.5. Aim and Paper Organization

The main focus of the article was formulated in the preceding sections of the Introduction considering skin properties (section 1.1), blood circulation in human body (section 1.2) and the general approach in construction of bio-heat models (section 1.3). Now after these initial words the following part of this article is organized in three main directions, with a common point of view, as follows:

1) Integer order models without flux damping (classical models of Pennes, Wulff, Klinger, Chen-Holmes, porous media concepts) in section 3. The models developed prior to that of Pennes are briefly presented in sections 2.1 and 2.2.

2) Integer order models with flux relaxation based on the Maxwell-Cattaneo and Dual-Phase-Lag approaches are discussed in section 5.

3) Bio-heat models with memory, that actually means heat flux damping modeled by fractional-order operators are especially analyzed in section 7.

4) A systematic approach in bio-heat heat conduction models is formulated in section 7.2 via application of two principle approaches: (a) Fractionalization via fractional Taylor series (section 8.1) developed in section 8, and (b) Fractionalization by application of the fading memory concept with different relaxation kernels (section 9).

For completeness of the exposition in this article, we will briefly present at the very beginning (section 2) some mechanistic ideas (models) that appeared in parallel to the continuum models forming the core of this study.

2. Formal Mechanistic Bio-heat Models

The mechanistic models presented briefly here have been developed by applying formal analogy between the blood vessels-tissue configurations and the design of the existing heat transfer equipments. This can be easily explained by the motivation to apply the existing knowledge, but at the same time this approach has inherent limitations which are obvious when more complex vascular architectures should be modeled. Two simple cases are commented: (a) Surface heat flux to the tissue by a convective blood flow streaming the blood-tissue interface, (b) Heat transfer between a cylindrical blood vessel of final length and the surrounding tissue.

2.1. Surface Heat Flux Model

Following the comprehensive analysis of Charny [21] the surface heat flux q_s through tissue (per unit area A_t, that is q_s has a dimension W/m²) and with missing flood flow is

$\begin{array}{ll}{q}_s={\omega }_b{\rho }_b{c}_b{L}_t(T_b-T_s)+\frac{k_t}{L_t}(T_b-T_s)=K_{eff}(T_b-T_s)& (12)\end{array}$

where ω_b is the volumetric blood flow to the tissue per unit volume of it, that is $\left[{(m^3/s)}_{blood}/m_{tissue}^3\right]$ while the heat conductivity of the tissue with a thickness L_t is k_t. Here, the first term in (12) is the thermal energy transported by the blood flow through the tissue layer of depth L_t and volume V_t = A_tL_t, while the second term is the Fourier law. This construction allows the effective thermal conductivity K_eff to be determined by experimentally measured surface (at the skin) heat flux q_s and temperature T_s, and the body core temperature T_b at a depth L_t. This formulation allows the effect of blood perfusion to be evaluated as an additional contribution to the dimensionless (relative) effective heat conductivity $K_{eff}/K_o=1+{\omega }_b(c_b{L}_t^2/k_t)$ [21] using the tissue conductivity K₀ in absence of blood perfusion as a reference scale. The value of K_eff at large differences (T_b − T_s) does not vary thus corresponding to the maximum degree of vasoconstriction in the tissue [21]. In addition, the value of L_t was estimated by early physiological measurements at L_f = 2 cm [23] (see a comprehensive analysis in [21]).

The first term in (12) is related to macroscopic heat removal by a fluid flowing through the tissue (blood perfusion), while the second one is the heat transport by conduction mechanism. Now, in order to make more visible the contributions of each transport mechanism we transform the heat flux Equation (12) in two equivalent dimensionless forms, namely

$\begin{array}{l}{\bar{q}}_s=\frac{q_s}{\frac{k_t}{L_t}(T_b-T_s)}=1+{\omega }_b\frac{L_t^2}{a_b}\frac{k_b}{k_t}=1+{\omega }_b{\tau }_b\left(\frac{k_b}{k_t}\right),\\ {\tau }_b=\frac{L_t^2}{a_b},a_b=\frac{k_b}{{\rho }_b{c}_b}& (13)\end{array}$

$\begin{array}{l}{\bar{q}}_s=\frac{q_s}{\frac{k_t}{L_t}(T_b-T_s)}=1+{\omega }_b\left(\frac{{\rho }_b{c}_b}{{\rho }_t{c}_t}\right)\frac{L_t^2}{a_t}=1+{\omega }_b{\tau }_t\left(\frac{{\rho }_b{c}_b}{{\rho }_t{c}_t}\right),\\ {\tau }_t=\frac{L_t^2}{a_t},a_t=\frac{k_t}{{\rho }_t{c}_t}& (14)\end{array}$

The first form (13) defines the diffusion time scale of the blood as τ_b, while the second form (14) determines the heat diffusion time scale τ_t. Taking in to account that ω_s has dimension [s⁻¹] (see the definition above) the products ω_bτ_b and ω_bτ_t are dimensionless. Hence, this dimensionless presentation easily demonstrates the contribution of the blood perfusion to the total heat flux through the tissue.

2.2. Heat Exchanger Formalistic Models for Flow in Vessels

The models considered in this section are not the core of this article but prior to continue with more sophisticated continuum (section 3) and fractional models (section 7) we have to mentioned some mechanistic models, quite simplified, but related to well-known solutions of fluid-flow heat exchange in simple systems of tubes (modeling blood vessels).

The mechanistic tube flow models are analyzed by Chato [24, 25] using the facts that for individual blood vessels the capillaries are important for attaining the thermal equilibrium with the surrounding tissue: the smaller arteries and the veins have intensive heat transfer with the tissue, in contrast to the larger vessels affecting relatively little the heat exchange.

2.2.1. Heat Transfer From Blood Flowing in a Vessel

Following Chato [24, 25] the blood behaves as non-Newtonian yield stress medium (Casson fluid) with a constitutive equation

$\begin{array}{ll}\sqrt{{\tau }_s}=\sqrt{{\tau }_y}+\mu \frac{\partial V}{\partial y}& (15)\end{array}$

relating the shear stress τ_s, the yield stress τ_y and the yield-stress viscosity μ. In (15) V is the axial blood velocity, while y is the direction normal to the vessel wall.

For steady flow the relation of the mean Nusselt number Nu = h_bD/k_b based on the vessel diameter D is related to the Graetz number $Gr=Re·Pr·\left(\frac{D}{L}\right)$ and can be approximated as

$\begin{array}{ll}N{u}_m=4+0.155exp(1.58lo{g}_{10}Gr),Gr<10^3& (16)\end{array}$

This approach invokes two principle questions [25]:

(1) At what distance from the entrance the blood-tissue will attain the equilibrium? That is, what is the heat transfer entrance length? This question was further discussed in construction of the more complex continuum model of Chen and Holmes (see section 3.4).

(2) What is the blood temperature in axial direction at a given length along the vessel wall?

Considering a simple model of concentric cylindrical vessel-surrounding tissue of final thickness R₂ − R₁ (see Figure 3), then applying what is well-known from every textbook on heat transfer, we have

FIGURE 3

Figure 3. Mechanistic representations of single blood vessel-tissue heat transfer as a countercurrent energy exchange process: See the text about the studies in Chato [25].

Radial temperature distribution

$\begin{array}{ll}\begin{array}{l}\frac{T-T_0}{T_b-T_0}=\ln \left(\frac{R_2}{r}\right)\frac{1}{\frac{k_t}{(h_b{R}_1)}+\ln \left(\frac{R_2}{R_1}\right)},R_1\le r\le R_2\hfill \end{array}& (17)\end{array}$

Heat flux transferred to the tissue

$\begin{array}{ll}\begin{array}{l}{q}_1=2\pi k_t(T_b-T_0)\frac{1}{\frac{k_t}{(h_b{R}_1)}+\ln \left(\frac{R_2}{R_1}\right)}\hfill \end{array}& (18)\end{array}$

where h_b is the local heat transfer coefficient, while h_e is an effective heat transfer coefficient for the blood flow.

The axial temperature profile of the blood is simply modeled by one dimensional equation [25]

$\begin{array}{c}\frac{d{T}_b}{dx}+2\frac{h_e}{{\rho }_b{c}_{b}V{R}_1}(T_b-T_0)=0\Rightarrow \frac{d\theta }{dX}+{\Lambda }_x\theta =0,\\ \theta =\frac{T_b-T_0}{T_{b1}-T_0},X=\frac{x}{2R_1}=\frac{x}{D},{\Lambda }_x=4\frac{N{u}_e}{RePr}& (19)\end{array}$

leading to exponential decaying temperature distribution along the x axis with a rate depending on Λ_x.

With these two examples we close this direction in the blood-tissue heat exchange since they are out of the scope of this study (the continuum based models). Extensive analysis of systems based on mechanical analogies are available in Charny [21] and Chato [24, 25].

3. Energy Balance Based Continuum Models (Starting by Pennes)

Modeling of the blood flow effect on bioheat in living tissues due the complex architecture of the vascular system has been developed in two distinct directions [26]:

(1) An approach considering the vascular system and as a system of rigid tubes protruding the tissues. With this point of view only large vessels are considered while the smaller ones are neglected thus reducing the order of complexity (similar approach was commented briefly in section 2.2).

(2) The second approach eliminates the vascular flow effects and uses continuum models incorporating in average the blood flow effects on the temperature over a control volume of the tissue.

The continuum approach allows, depending on the problem that should be modeled by an additional source term in the modeling equation; or defining effective transport coefficients such as the heat conductivity and the blood-tissue heat transfer coefficient. The most popular models are the Pennes model (see section 3.1) where the blood perfusion is modeled by a fluid-surface heat transfer term with a heat transfer coefficient dependent on the blood flow rate, and the Weinbaum-Jiji models (see section 3.5) with heat conductivity dependent on the vascular structure. These models will be analyzed here together with other continuum models since they provide a sufficient basis to develop two modeling trends accounting relaxation effects in bio-heat transfer, namely:

(a) Heat-wave models (section 5) with simple heat flux lagging (section 5.1), and dual-phase lag models (section 5.2), both of them employing local derivatives.

(b) Non-local models modeling the heat flux and temperature gradient damping (lagging) by fractional differential operators (sections 7 and 8.1).

3.1. The Model of Pennes

Pennes [27] formulated the oldest and the most studied mathematical model based on temperature profiles in limb (forearms) assuming maximum of temperature at the axis of symmetry, thus reasonably formulating the model equation in cylindrical coordinates, namely

$\begin{array}{l}{\rho }_t{c}_t\frac{\partial T_t(r,t)}{\partial t}=\frac{k_t}{r}\frac{\partial }{\partial r}\left[r\frac{\partial T_t(r,t)}{\partial r}\right]+{\omega }_b{\rho }_b{c}_b\\ (T_{a0}-T_t(r,t))+Q_m& (20)\end{array}$

assuming a uniform distribution of the metabolic heat generation Q_m [W/m³] in the tissue layer, as well as uniform tissue heat conductivity k_t. With a convective boundary condition at the skin surface (r = R)

$\begin{array}{ll}{k_t\frac{\partial T_t}{\partial r}|}_{r=R}=h\left[T_t(R,t)-T_{\infty }\right]& (21)\end{array}$

The principle simplifying assumptions of the Pennes model are [4, 28]:

(a) Capillary bed equilibration: pre-arteriole heat transfer and pre-venule heat transfers are neglected. The only variable in the Pennes model is the tissue temperature.

(b) Isotropic blood flow in small capillaries, i.e., the blood flow direction is not considered.

(c) All blood vessels close to the capillary beds do not affect the heat exchange and the local vascular geometry is ignored.

(d) The blood reaches the capillary beds through the arterioles at deep body temperature (core temperature).

Actually, the Pennes model assumes that the arterial blood flow bathes the tissue at the same temperature everywhere in the limb cross section (see detailed comments in [21, 29]). This approach resembles the thermal problem of the longitudinal metal fin with constant cross section A_F (and cross section perimeter P_F), available in every textbook of heat transfer (see as good reference sources [30] and [31]), namely

$\begin{array}{ll}{\rho }_F{c}_F\frac{\partial T_F(x,t)}{\partial t}=k_t\frac{{\partial }^2{T}_F(x,t)}{\partial x^2}-h_F\frac{P_F}{A_F}\left[T_F(x,t)-T_{\infty }\right]& (22)\end{array}$

where h is the heat convective transfer coefficient and T_∞ equals the arterial blood temperature, as it is schematically illustrated in Figure 4 (Forearm-Fin analogy). This term in the Pennes' equation is a matter of arguments under different assumptions. In accordance with Wulff [32] the last term is not local in contrast to the other two (the time and space derivatives); this comment will appear in the criticism of the Wulff's continuum model [32] (see also [21, 29]) discussed further in this analysis. This standpoint would be true if the assumption of Pennes [27] about the thermal equilibrium between the tissue and venous blood is accepted. That is, if T_t (r, t) = T_b(ven) in this term is independent of the space coordinate x, then the last term in (3.6.1) is not local. However, if T_t (r, t) depends on the space coordinate r as T_F (x, t) in the fin Equation (22) then this term is local: it is a distributed heat sink such as the distributed metabolic term Q_m. Here, we can see the first contradiction in the Pennes equation formulation: in general the tissue temperature depends on both the space coordinate and time but the last cooling term it is assumed as space independent, even though actually it is; the source/sink term (termed also convective term) is homogeneously distributed along the x axis.

FIGURE 4

Figure 4. Forearm-Fin mechanical analogy of the logical construction used in the Pennes model. Present author conceived analogy and interpretations: see the explanations in the text.

The main deficiencies of the Pennes model can be briefly outlined as [28]:

(a) The thermal equilibration in the capillaries is a misconception since this process takes place in the pre-arteriole and post-venule vessels (with diameters within the range of 70–500 μ m).

(b) The blood flow spatial orientation is not taken into account, that is the countercurrent heat transfer between the artery-vein pairs is not considered, for example.

(c) The blood perfusion is overestimated since the arterial temperature changes continuously from the deep body temperature of the aorta to the secondary arteries and further to all branches of the system: (i) the pre-arteriole temperature is not equal to the body core temperature; (ii) The vein return temperature is not equal to the mean tissue temperature.

The nondimensalization of (22) yields

$\begin{array}{ll}\frac{\partial \theta }{\partial \bar{t}}=\frac{{\partial }^2\theta }{\partial x^2}-\frac{h}{k_t}{L}^2\frac{P_F}{A_F}\theta ,\Rightarrow \frac{\partial \theta }{\partial \bar{t}}=\frac{{\partial }^2\theta }{\partial x^2}-\bar{B}i\frac{P_F}{A_F}\theta & (23)\end{array}$

where the dimensionless variables are

$\begin{array}{ll}\theta =\frac{T-T_{\infty }}{T_t-T_{\infty }},\bar{t}=\frac{t}{t_0},t_0=\frac{L^2}{a_t},a_t=\frac{k_t}{{\rho }_t{c}_t}& (24)\end{array}$

The group $\bar{B}i=\frac{h}{k_t}{L}^2\frac{P_F}{A_F}$ is dimensionless and can be assumed as analog of the Biot number Bi. The product P_FL_F equals the total surface area of the fin; then $L_F^2(P_F/A_F)$ with a dimension of length defines the characteristic length scale of the fin l_F.

The solution of the Pennes model in radial coordinates [21] is

$\begin{array}{l}{T}_t=A_p{I}_0\left(r{a}_p\right)+\frac{Q_m}{{\omega }_b{\rho }_b{c}_b}+T_{a0},a_p=\sqrt{\frac{{\omega }_b{\rho }_b{c}_b}{k_t}},\\ A_p=\frac{T_{\infty }-\frac{Q_m}{{\omega }_b{\rho }_b{c}_b}-T_{a0}}{\frac{k_t{a}_p}{h}{I}_1\left(R{a}_p\right)+I_0\left(R{a}_p\right)}& (25)\end{array}$

where I₁ and I₀ are modified Bessel functions

First, the parameter a_p has a dimension [m⁻¹] because it can be represented as $a_p=\sqrt{{\omega }_b\frac{1}{a_b}\frac{k_b}{k_t}}$, where a_b is the thermal diffusivity of the blood. Taking into account that the dimension of ω_b is [s⁻¹] (volumetric blood flow rate per unit volume of the tissue) then the ratio (ω_b/a_b) has a dimension [1/m²]. Therefore, the products ra_p and Ra_p in (25) are dimensionless. Further, the ratio Q_m/(ω_bρ_bc_b) has a dimension of temperature per second [K/s] since the dimension of ω_b ((m³/s)/m³) is [s⁻¹], as it was commented above.

The scaling of the boundary condition in (21) yields

$\begin{array}{r}{\frac{\partial {\theta }_t}{\partial \bar{r}}|}_{r=R}=\frac{hR}{k_t}\theta (R,t)=Bi\left[\theta (R,t)\right],\theta =\frac{T_t(R,t)}{T_{art}(R,t)-T_{\infty }},\\ \bar{r}=\frac{r}{R}& (26)\end{array}$

where $Bi=\frac{hR}{k_t}$ is the thermal Biot number representing the ratio of the heat transported to the tissue by convection to the heat transported through the tissue by conduction.

It is well-known that for very small Biot numbers (Bi < 1) the temperature gradients in the body (tissue) are negligible and the model of Pennes reduces to the lumped heat capacitance model, namely

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial T_t}{\partial t}={\omega }_b{\rho }_b{c}_b(T_{a0}-T_t)+Q_m& (27)\end{array}$

or in a scaled version as

$\begin{array}{l}\frac{\partial T_t}{\partial t}=-{\omega }_b\frac{{\rho }_b{c}_b}{{\rho }_t{c}_t}(T_{a0}-T_t)+\frac{Q_m}{{\rho }_t{c}_t}\Rightarrow \frac{\partial {\theta }_t}{\partial \tau }={\omega }_b{t}_0\left(\frac{{\rho }_b{c}_b}{{\rho }_t{c}_t}\right)\\ +\frac{Q_m}{{\rho }_t{c}_t}{t}_0& (28)\end{array}$

For example, the steady-state version of the Pennes equation in rectangular coordinates (1-D) is

$\begin{array}{ll}\frac{{\partial }^2\theta }{\partial x^2}-\bar{B}i\frac{P_F}{A_F}\theta =0& (29)\end{array}$

and its solution is [21]

$\begin{array}{ll}\theta =\frac{sinh\left(\sqrt{\bar{B}i}x\right)}{sinh\left(\sqrt{\bar{B}i}\right)},0\le x\le 1& (30)\end{array}$

The Pennes equation was seriously criticized in consequent developed studies (see in the sequel), and we will stress the attention on this, but now we focus on nondimensalization and the physical meaning of the parameters in the solution (22). The principle shortcomings, in addition to the ones mentioned above, of the Pennes model can be outlined as:

(i) As mentioned above, the perfusion terms is assumed as that the tissue is bathed by the arterial flow, while the venous return from the distal part having an important effect on the tissue temperature is neglected.

(ii) The analysis of Shitzer and Kleiner [33] and Shitzer and Chato [34] reveals that the perfusion term may play a role of a sink in the warm core region of the tissue (precisely in the cylindrical limb model) and as source in the cooler peripheral area (see also analysis in [21, 29]. Further, with increase in the arterial bold flow, that is the increase in the perfusion rate, then the last term in the Pennes equation performs as a heat source in a large section of the living tissue [21]. Mathematically this means that the sign of the perfusion term should change (as source or as sink) depending on the location of the modeled tissue, a properties which is not encountered in the Pennes model.

(iii) The models do not consider the detailed structures of the vascular networks, especially the countercurrent arrangements as well as important anatomical characteristics of the blood circulatory system.

Despite these deficiencies the Pennes equation was extensively applied to model heat transfer in living tissues but mainly within body segments and with good very good agreements with experimental data. According to Charny [21] in most of cases the neglected heat transfer in the prearteriole and postvenule blood vessels of the circulatory systems can be used to explain discrepancies between the instrumental measurements and the data predicted by the Pennes model. We will meet again the Pennes model when problems related to its fractionalization will be discussed further in this article.

3.2. The Continuum Model of Wulff

The Wulff's approach [32] is a modeling of heat transfer through the tissue on a continuum basis and the main criticism is on the formulation of the perfusion terms. Actually, Wulff clearly states that Pennes' perfusion term is a global (integral term) that does not corresponds to the other terms (as derivatives) which are local derivatives. As mentioned earlier in this text, this is a matter of argument, but for better understanding what we mean, let us consider the basic idea of the Wulff's model. Following Wulff, since the blood is moving through the tissue it may remove or derive heat by convection in any directions but not only in the direction defined by of the local temperature gradient. Following this assumption, the convective energy transport through the tissue is Wulff [32] (see also [21, 29])

$\begin{array}{ll}{\rho }_b{\overline{h}}_b{v}_h=\frac{1}{4\pi }{\displaystyle \underset{\Omega }{\int }}{\rho }_b{h}_{b}vd\omega & (31)\end{array}$

In (31) ${\overline{h}}_b$ is the blood specific enthalpy formulated as

$\begin{array}{ll}{\overline{h}}_b={\displaystyle \underset{T_0}{\overset{T_b}{\int }}}{c}_p\left(T_b^{*}\right)d{T}_b^{*}+\frac{P}{{\rho }_b}+\Delta H_f(1-\epsilon )& (32)\end{array}$

Further, ω is the solid angle across the control surface of the blood vessel, while Ω is the entire solid angle, and v_h is the local mean apparent blood velocity. ΔH_f is the specific enthalpy of the metabolic reaction. With these assumptions the energy flux at any point of the tissue is [32, 35]

$\begin{array}{ll}q=-k_t\nabla T_t+{\rho }_b{h}_h{v}_h& (33)\end{array}$

Thus, the energy balance simply reads as

$\begin{array}{ll}\begin{array}{l}{\rho }_b{c}_b\frac{\partial T_t}{\partial t}=-\nabla ·q\hfill \end{array}& (34)\end{array}$

that in detailed form is [21, 29, 32]

$\begin{array}{ll}{\rho }_b{c}_p\frac{\partial T_t}{\partial t}=k_t\frac{{\partial }^2{T}_t}{\partial x^2}-{\rho }_b{v}_h{\displaystyle \underset{Metabolicheat}{{\displaystyle \underbrace{\left(c_p\frac{\partial T_b}{\partial x}-\Delta H_f\frac{\partial \epsilon }{\partial x}\right)}}}}& (35)\end{array}$

The last term in right-hand side of (35) is equivalent to the metabolic heat source Q_m. The main problem in application of the Wulff model is the determination of the local blood mass flux ρ_bv_h [35].

Further, Wulff assumed that the blood temperature gradient in (35) equals the gradient of the temperature gradient of the environment surrounding the tissue and this step permits to assume that T_b = T_t (the same assumption as in the Pennes equation !). Hence (35), can be written in a more convenient form, namely

$\begin{array}{ll}{\rho }_b{c}_p\frac{\partial T_t}{\partial t}=k_t\frac{{\partial }^2{T}_t}{\partial x^2}-{\rho }_b{v}_h{c}_p\frac{\partial T_t}{\partial x}+Q_m& (36)\end{array}$

which can be also re-written (with $a_t=\frac{k_t}{{\rho }_b{c}_p}$-thermal diffusivity of the tissue) as

$\begin{array}{ll}\frac{\partial T_t}{\partial t}=a_t\frac{{\partial }^2{T}_t}{\partial x^2}-v_h\frac{\partial T_t}{\partial x}+\frac{Q_m}{{\rho }_b{c}_b},a_t=\frac{k_t}{{\rho }_b{c}_b}& (37)\end{array}$

or in a more familiar form as 1–D convection-diffusion equation of heat transfer

$\begin{array}{ll}\frac{\partial T_t}{\partial t}+v_h\frac{\partial T_t}{\partial x}=a_t\frac{{\partial }^2{T}_t}{\partial x^2}+\frac{Q_m}{{\rho }_b{c}_b}& (38)\end{array}$

3.2.1. Nondimensalization of the Wulff's Model

The nondimensalization of (37) or (38) with $\theta =\frac{T_t(x)-T_t(0)}{T_t(L_t)-T_t(0)}$, $\bar{x}=x/L_f$, $\bar{t}=t/t_0$ and $t_0=L_t^2/a_t$ (as well as assuming k_t = k_b) yields

$\begin{array}{ll}\frac{\partial \theta }{\partial \bar{t}}+Pe\frac{\partial \theta }{\partial \bar{x}}=\frac{{\partial }^2\theta }{\partial {\bar{x}}^2}+Q_m^{*},Pe=\frac{v_h{L}_t}{a_t}=\frac{v_h{L}_t}{k_t}{\rho }_b{c}_b& (39)\end{array}$

The groups

$\begin{array}{ll}Pe=\frac{v_h{L}_t}{a_t}=\frac{v_h{L}_t}{k_t}{\rho }_b{c}_b,Q_m^{*}=\frac{Q_m}{{\rho }_b{c}_b}\frac{L_t^2}{a_t}=Q_m\frac{L_t^2}{k_t}& (40)\end{array}$

are dimensionless.

Moreover, the Peclet number Pe = advective transport/diffusion transport is a weighting coefficient showing the contribution of the perfusion to the total heat transfer in the living tissue. For low Peclet numbers (low perfusion rates) the heat conduction transport mechanism dominates.

3.3. The Continuum Model of Klinger

The model of Klinger [36, 37] is similar to that of Wulff [32] and focuses on correct description of experimental data from thermal clearance experiments. In the Klinger's model the tissue temperature is related to the rate at which the deep tissue temperature changed during point source heating [38, 39]. Taking into account the critique of Perl [38, 39] to the Pennes model that the unidirectional blood flow was postulated while the naturally existing nonunidirectional blood flow, Klinger [40] formulated the main problems emerging in the implantation of the Pennes equation, and what should be improved, that is

K1) The solution envisages a microscopic (local) temperature, that is the temperature at each point of the tissue.

K2) The measuring devices have finite dimensions and in most of the cases larger than the length scale of the tissue of interest and consequently, the measured temperature is an average value and cannot be accepted as a local one, i.e., it is not local, but do not misunderstand it as a result of a memory based differential operator since all derivatives used in the integer-order models are local operators !

K3) The correction of the missing nonunidirectional blood flow in Pennes' model is corrected by introducing the concept of convection multipoles related upon the in vivo vascular anatomy, thus accounting not only the magnitude of the blood flow but also its direction [21, 36, 37]. The Klinger model is based on the assumption of constant tissue properties (k_t, ρ_t and c_t) and incompressibility of the blood flow (that is ∇ · v_b = 0) and expressed as

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial T_t}{\partial t}+{\rho }_b{c}_b{\text{v}}_b·\nabla T_t=k_t{\nabla }^{2}T+Q& (41)\end{array}$

In a dimensionless form (41) becomes [21, 36, 37]

$\begin{array}{ll}\frac{\partial T_t}{\partial \tau }+Pe{\text{v}}^{*}·\bar{\nabla }{T}_t={\bar{\nabla }}^2{T}_t+\frac{Q{L}^2}{k_t}& (42)\end{array}$

Here v* denotes the blood velocity with respect to the characteristic velocity used in the formulation of the Peclet numberPe. In addition, ${\bar{\nabla }}^2$ is dimensionless Laplacian gradient operator based on a characteristic length scale L, while $\tau =L^2/a_t$ is the characteristic time (a_t is the tissue thermal diffusivity). As commented by Charny [21] this time characterization is equivalent to Fourier number of unity, that is when t = τ = L²/a. However, this is not true, as it possible to see from the preceding examples of nondimensalization of model equations. The main question arising here is: What is the definition of the characteristic length scale L? It should be determined arbitrary, that is depending on the region of the tissue considered and the experience of the modeler. It is noteworthy, to understand the main idea of this model, that Klinger introduced a dimensional velocity $\overline{v}=Pev*$ thus allowing the energy balance equation to presented as

$\begin{array}{ll}\psi T_t+\frac{Q{L}^2}{k_t}=0,\psi =\left[{\bar{\nabla }}^2-\left(Pe{\text{v}}^{*}·\bar{\nabla }\right)-\frac{\partial }{\partial \tau }\right]& (43)\end{array}$

Assuming the point heat source heating, the main idea in the studies of Klinger) is the replacement of the source term $Q{L}^2/k_t$ by a Dirac delta function (vanishing elsewhere but appearing at the position r₁); this immediately invoked the Green function based solutions (41) [21, 36, 37, 41]. In accordance with Klinger, in all vascular arrangements studied, the multiple moment is directly proportional to the total volumetric blood flow rate through the vessels with repeating tissue element. The effects of the multipole moments depends on the ratio (λ/L); here L is the macroscopic length scale, while λ is the length of the side of cube over which the average tissue temperature is determined, and the nondimensional volume is (λ/L)³. Since (λ/L) < 1 the results of Klinger's modeling experiments is that the countercurrent flow effect on the tissue heat transfer is a second order effect.

Hence, the concept of averaging the temperature around a certain position leads to [40]:

K4) The effect of the blood convection on the average temperature can be represented as superposition of terms, each of which mirrors a particular properties of local symmetry of the vascular system.

K5) The order of magnitude of the convection effect decreases with increasing in the symmetry.

K6) The Peclet number used by Klinger [36] is defined as a ratio of two characteristic times: $Pe=t_{conv}^{char}/t_{cond}^{char}$, where $t_{conv}^{char}={\mathrm{\text{v}}}_0/L$ and $t_{cond}^{char}=L^2/a_t$ (see the a comments after (refeq: Wulff-9) where Pe is defined as a ratio of the two heat fluxes (convection to conduction).

The model developed by Klinger is a serious contribution to bio-heat modeling since it address some important issues, among them [21, 36, 37, 40, 41]:

1) Since ε_b ≪ 1 [21] it follows that k_eff is practically independent of the blood flow, and consequently it is reasonable to accept k_eff ≈ k_{solid tissue}.

2) The Green function solution allows quantifying the vessel number density, blood perfusion rate and the vessel spatial construction. The principle outcome of this is that the temperature field are strongly influenced by the spatial geometry of the blood carrying vessels. More detailed analyzes are available elsewhere [21, 29].

3.4. The Model of Chen-Holmes

Chen and Holmes [42] have considered heat transfer applied in living tissue with hierarchical system of vessels with the continuum approach. The modeling technology is the same as in the models of Wulff and Klinger making energy balance of tissue-blood control volume. This physical hypothesis is based the existence of a large number of blood vessels in the tissue domain (volume) over which the energy balance is taken. As a rule, the characteristic length scale of this domain is much larger than the characteristic dimensions of the individual blood vessels, thus the continuum approach can be applied. Obviously, with this approach the blood flow through the small vessels protruding the tissue affects the definition of the effective heat transport coefficients such as tissue conductivity and the metabolic term. Hence, in terms of local mean temperature, the energy balance equation for both the tissue and vascular spaces are [21, 29, 42]:

In the solid tissue space

$\begin{array}{ll}d{V}_s\left[{\rho }_s{c}_s\frac{\partial T_s}{\partial t}\right]=\underset{conduction}{\underbrace{d{Q}_{ks}}}+\underset{\begin{array}{c}blood\\ compartment\end{array}}{\underbrace{d{Q}_{bs}}}+\underset{\begin{array}{c}metabolic\\ heating\end{array}}{\underbrace{d{Q}_m}}& (44)\end{array}$

In the vascular space

$\begin{array}{ll}d{V}_b\left[{\rho }_b{c}_b\frac{\partial T_b}{\partial t}\right]=\underset{conduction}{\underbrace{d{Q}_{kb}}}-\underset{\begin{array}{c}blood\\ compartment\end{array}}{\underbrace{d{Q}_{bs}}}+{\displaystyle \underset{S}{\int }({\rho }_b{c}_{b}T)\text{vds}}& (45)\end{array}$

where the last integral term is the contribution of the convective energy gain due to blood flow at velocity v across a surface area S. With help of (44) and (45) when dV (the bulk control volume) goes to zero the continuum energy balance of the tissue is

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial T_t}{\partial t}=\underset{conduction}{\underbrace{d{q}_k}}+\underset{metabolic}{\underbrace{d{q}_m}}+\underset{\begin{array}{c}blood\\ flow\end{array}}{\underbrace{d{q}_p}}& (46)\end{array}$

$\begin{array}{ll}d{q}_k=\frac{Q_{ks}+Q_{kb}}{dV}=\nabla ·k_k\nabla T_t,d{q}_p=\frac{1}{dV}{\displaystyle \underset{S}{\int }}({\rho }_b{c}_{b}T)\text{v}ds& (47)\end{array}$

Note 2: The temperature under the integral sign in (47) is not simply equal to T_s in contrast to the assumption used by Wulff [32]. Moreover, the tissue temperature T_t is defined as

$\begin{array}{ll}{T}_t=\frac{1}{\rho c}\left[\left(1-\frac{d{V}_b}{dV}\right){\rho }_s{c}_s{T}_s+\frac{d{V}_b}{dV}{\rho }_b{c}_b{T}_b\right]& (48)\end{array}$

and since the ratio dV_b/dV → 0 we get T_t → T_s.

Skipping details in the final formulation (see the analyzes in [21, 29] of the new equation) named Chen-Holmes bio-heat equation is

$\begin{array}{ll}\begin{array}{c}\hfill \rho c\frac{\partial T_t}{\partial t}={\displaystyle \underset{\begin{array}{c}overal\\ heat\\ conduction\end{array}}{{\displaystyle \underbrace{\nabla ·k_t\nabla T_t}}}}+{\displaystyle \underset{\begin{array}{c}effective\\ heatconduction\\ (withperfusion)\end{array}}{{\displaystyle \underbrace{\nabla ·k_p\nabla T_t}}}}+{\displaystyle \underset{\begin{array}{c}perfusion\\ heat\\ source\end{array}}{{\displaystyle \underbrace{{\omega }_j^{*}(T_a^{*}-T_t)}}}}-\\ \hfill -{\displaystyle \underset{\begin{array}{c}convective\\ heat\\ transfer\end{array}}{{\displaystyle \underbrace{{\rho }_b{c}_b{\text{v}}_p(\nabla ·\nabla T_t)}}}}+{\displaystyle \underset{\begin{array}{c}metabolic\\ heat\\ source\end{array}}{{\displaystyle \underbrace{\frac{d{q}_m}{dt}}}}}\end{array}& (49)\end{array}$

$\begin{array}{ll}\begin{array}{c}\hfill \rho c\frac{\partial T_t}{\partial t}={\displaystyle \underset{\begin{array}{c}overall(effective)\\ heat\\ conduction\end{array}}{{\displaystyle \underbrace{\nabla ·K_{eff}\nabla T_t}}}}+{\displaystyle \underset{\begin{array}{c}effective\\ heatconduction\\ (withperfusion)\end{array}}{{\displaystyle \underbrace{\nabla ·k_p\nabla T_t}}}}+{\displaystyle \underset{\begin{array}{c}perfusion\\ heat\\ source\end{array}}{{\displaystyle \underbrace{{\omega }_j^{*}(T_a^{*}-T_t)}}}}-\\ \hfill -{\displaystyle \underset{\begin{array}{c}convective\\ heat\\ transfer\end{array}}{{\displaystyle \underbrace{{\rho }_b{c}_b{\text{v}}_p(\nabla ·\nabla T_t)}}}}+{\displaystyle \underset{\begin{array}{c}metabolic\\ heat\\ source\end{array}}{{\displaystyle \underbrace{\frac{d{q}_m}{dt}}}}}\end{array}& (50)\end{array}$

Referring to the effective heat conduction of the tissue we address the formulated conductivity enhancement by the blood flow through the tissue. Hence, the effective heat conduction flux q_cond [the first term in the right-hand side of (50)] in the control volume of the perfused tissueδV is presented as

$\begin{array}{ll}\begin{array}{c}\hfill q_{cond}=\frac{q_{cond(tissue)}+q_{cond(blood)}}{\delta V}=\nabla ·(k_{eff}\nabla T_t),\\ \hfill k_{eff}={\epsilon }_b{k}_{blood}+(1-{\epsilon }_b)k_{solidtissue}\end{array}& (51)\end{array}$

Here, the effective heat conductivity is defined as in a two-phase (solid-liquid) medium where the blood volumetric fraction is

$\begin{array}{ll}{\epsilon }_b=\frac{\delta V_b}{\delta V}=\frac{controlvolumeintheblood}{controlvolumeintheperfusedtissue}& (52)\end{array}$

Therefore, since ε_b ≪ 1 [21] it follows that k_eff is practically independent of the blood flow, and consequently it is reasonable to accept k_eff ≈ k_{solid tissue}. Further, using the porosity of the tissue ε_b where blood flows, the effective density, heat capacity and local mean tissue temperature are defined as [29, 42]

$\begin{array}{lll}\rho & =& (1-{\epsilon }_b){\rho }_s+{\epsilon }_b{\rho }_b,c=\frac{1}{\rho }[(1-{\epsilon }_b)({\rho }_s{c}_s{T}_s)\\ & & +{\epsilon }_b({\rho }_b{c}_b{T}_b)]& (53)\end{array}$

$\begin{array}{ll}{T}_t=\frac{1}{\rho c}\left[(1-{\epsilon }_b)({\rho }_s{c}_s{T}_s)+{\epsilon }_b({\rho }_b{c}_b{T}_b)\right]& (54)\end{array}$

where the subscript b denotes blood.

3.4.1. Some Remarks on the Chen-Holmes Model

The Chen-Holmes final equations need some relevant sub-problems (of strong importance with respect to the physical background) to be commented briefly, that is:

Remark 1: The tissue temperature T_s can be replaced by the volume-weighted continuum temperature T_t as much as ε_b ≪ 1. Moreover, a perfusion term ${\omega }_j^{*}(T_a^{*}-T_t)$ similar to that in the Pennes equation also emerges in the Chen-Holmes equation. However, this is a formal similarity, since the analysis performed by Chen and Holmes about the thermal equilibration lengths in various blood transporting vessels revealed that the velocities are proportional to power of 1/3 of the vessel radii. Hence, for larger vessels (such as aorta and large arteries, as well as vein) the thermal equilibria require equilibration lengths to have order of meters [21, 42]. Otherwise, in vessel compromising microcirculations (arterioles, venules and capillaries) the thermal equilibration lengths are of order of microns. Therefore, the Chen-Holmes model and the relevant analysis clearly contradict the Pennes formulation that all tissue-blood heat transfer occurs in the capillary bed.

Remark 2: In the context of the preceding remark, Chen and Holmes emphasis on the fact that the formal perfusion heat source (54) differs in its physics from the perfusion term in the model of Pennes. There two principle differences albeit the formal similarities, that is Charny [21] and Chen and Holmes [42]: (1) The coefficient ${\omega }_j^{*}$ is the total perfusion bleed-off to the tissue only from the microvessels past the j^*th generation of branching. In contrast, the Pennes term ω_b considers the bleed-off from all generations of vasculature. (2) The perfusion term of Chen-Holmes (54) is proportional to the departure difference $(T_a^{*}-T_t)$ which is different from (T_a0 − T_t) (see Equation 3.6.1). The difference between $(T_a^{*}-T_t)$ and T_a0−Tt can attain an order of 10 − 50 percents using as a base for comparison the thermal equilibration lengths of the large vessels.

Remark 3: Chen and Holmes accepted the idea of Wulff [32] that the blood temperature equals the solid tissue temperature everywhere in the control volume. This assumption allows the perfusion through the control volume to be presented by a conventional convection term −ρ_bc_bv_p · ∇T_t, where ρ_bc_bv_p is the blood mass flux through the control volume of the tissue.

Remark 4: The effective (perfusion) conductivity [see the second term in the right-hand sides of (49) and (50)] simplifies the analysis assuming that the conductivity k_p is a vascular quantity, that is, it actually is related to the suggestion that the heat flux proportional to temperature gradient is normal to the differential surface area.

3.5. Weinbaum-Jiji Model

Weinbaum and Jiji [43–45] and Weinbaum et al. [46–49] developed a model starting from objections to the Pennes model addressing the missing directionality effects of the larger blood vessels on the blood-tissue-heat transfer. Additional critiques of Weinbaum and Jiji refereed to absence in the Pennes model characteristics describing the geometry of the vessel system, that is sizes, diameters, branching, etc. The series of studies began with the ideas presented in Weinbaum and Jiji [43].

The general step in the Weinbaum-Jiji model (WJM) is the modification of the Pennes equation with an effective thermal conductivity which is a function of the blood flow rate and the vessel architectures. This modification lays on the suggestion that small arteries and veins are parallel, thus the flow direction is countercurrent; hence, there are counterbalanced heating and cooling. The hypothesis is mainly applicable to intermediate skin layers (see comments in [21, 29]).

The governing equations for deep tissue layer are based on the vessel geometry and capillary bleed-off phenomena on a continuum basis. The continuity relationship related to the mass conservation in paired blood vessels is defined as $\frac{d}{ds}\left(n{a}^2ū\right)=-2na{g}_b$. Here, s is the location along the length of the countercurrent network which differs from the distance x due to the inclination angle of the vessel pair with respect to the normal axis to the skin surface [21]. Further, ū is the bulk mean electricity in the blood vessel, while g_b is the perfusion bleed-off per unit vessel surface area (the subscript b refers to blood).

Next, neglecting the axial heat conduction in arteries (with temperature T_a) and veins (with temperature T_v), the energy balances for the blood flows are

$\begin{array}{ll}({\rho }_b{c}_b)\frac{d}{ds}(n\pi a^2\overline{u}{T}_a)=-n{q}_a-({\rho }_b{c}_b)(2\pi an{g}_b)T_a& (55)\end{array}$

$\begin{array}{ll}({\rho }_v{c}_v)\frac{d}{ds}(n\pi a^2\overline{u}{T}_v)=-n{q}_a-({\rho }_v{c}_v)(2\pi an{g}_b)T_v& (56)\end{array}$

Here q_a is the flux of heat loss through the artery wall, while q_v is the heat gain by conduction (per unit length) into the vein through its wall; T_a and T_v are the blood bulk temperature in the arteries and the veins, respectively.

Further, the second term in (55) describes the heat loss from the artery due to perfusion, while in (56) the corresponding term is the gain from the atrial blood due to perfusion.

From (55) and (56), after substraction, it follows that

$\begin{array}{ll}({\rho }_b{c}_b)\left[\frac{d}{ds}(n\pi a^2\overline{u}{T}_a)-\frac{d}{ds}(n\pi a^2\overline{u})T_v\right]=\nabla ·k_t\nabla (T_t)+q_m& (57)\end{array}$

Applying the continuity equitation in the right-hand sides of these equations we get the heat fluxes, namely

$\begin{array}{ll}\begin{array}{c}{q}_a=-{(\rho c)}_b\left(n\pi a^2\overline{u}\right)\frac{d{T}_a}{ds},\\ q_v=-{(\rho c)}_b\left(n\pi a^2\overline{u}\right)\frac{d{T}_v}{ds},\\ q_a-q_v={(\rho c)}_b(\pi a^2\overline{u})\frac{d}{ds}\left[T_a-T_v\right]\end{array}& (58)\end{array}$

Taking into account that only T_a and T_v exist in the definitions of the fluxes (58) it is impossible to determine the tissue temperature T_t. To overcome this problem, it was suggested a simplification that [44]

$\begin{array}{ll}{T}_t\cong \frac{T_a+T_v}{2}& (59)\end{array}$

In addition, the authors suggested that the heat transfer through the tissue surrounding the vessels is dominated by heat conduction mechanism. As a consequence they assumed the simplification defined as

$\begin{array}{ll}{q}_a\cong q_v=\sigma k_t(T_a-T_v),\sigma =\frac{\pi }{cos{h}^{-1}\left(\frac{L_s}{a}\right)}& (60)\end{array}$

where σ is a shape factor related to the ratio of vessel spacing to the vessel diameter L_s/a.

After all these assumptions and simplifications the Weinbaum-Jiji equation [44] is

$\begin{array}{ll}\frac{n{\pi }^{2}a{k}_b}{4k_t}Pe\left[\frac{2gPe}{\sigma \overline{u}}\frac{d{T}_t}{ds}-\frac{d}{ds}\left(\frac{aPe}{\sigma }\frac{d{T}_t}{ds}\right)\right]=\nabla ·\left(k_t\nabla T_t\right)+q_m& (61)\end{array}$

with Peclet number defined as $Pe=\frac{2a{(\rho c)}_b\overline{u}}{k_b}$

The Weinbaum-Jiji concept leads to a parabolic models with infinite speed of the heat flux. We will remember about this when the heat flux relation should be taken into account (in models with phase-lag and fractional models). Now, we stress the attention on some physical facts that might help us to understand the physics behind. As it is mentioned in Weinbaum et al. [46] the thermal relaxation of the blood in paired arteries and veins occurs through two principle mechanisms:

(a) An axial departure decay due to radial diffusion in blood vessels and the surrounding tissue,

(b) Countercurrent heat exchange when an arteries and veins are closely juxtaposed.

It is worthy to note that an important characteristic of the Weinbaum-Jiji model is its applicability to heat transfer in peripheral tissues only, where the basic assumptions used in its construction are valid. That is, the assumptions that the mean blood temperature (T_a + T_v)/2 equals the tissue temperature. In this models construction, the superposition used by Baish et al. [50] is used to validate the underlying assumptions.

The Weinbaum-Jiji simplified model based on parallel blood transporting vessels was commented thoroughly by Wissler [51] and some questionable assumptions, crucial for its derivation, were outlined, among them:

(a) The tissue and blood temperatures are related: therefore is not necessary the tissue temperature to be calculated as an arithmetic mean of the artery and venous temperatures (see Equation 59).

(b) The temperature at which the blood enters a large vessel is important because the thermal equilibration needs a certain distance along the vessel wall (see the assumptions of the Chen-Holmes model commented above).

(c) Obviously the tissue temperature will be intermediate between the arterial an venous blood temperatures but Equation (59) cast doubts about its general applicability.

(d) The assumed close thermal coupling between blood and tissue (the basis of the Weinbaum-Jiji simplified model-see Figure 5) leads to Equation (61) where the arterial blood departure is missing.

(e) Wissler commented that the first version of the Weinbaum-Jiji model (published in [47] and not commented here) is more realistic that the simplified one (61).

FIGURE 5

Figure 5. Blood bleed-off trough thermally significant parallel blood vessels only. See explanation in text about the works of Weinbaum et al. [46] and Nakayama et al. [52] and the ideas developed by Wissler [51]. Note: The area is subdivided by area of solid A_s and area of capillaries A_c.

For more details see Wissler [51] and the commented works of Weinbaum and Jiji [43–45] and Weinbaum et al. [46–49], the presentation in Chapter 10 of Jiji [28] as well as the analyzes in Charny [21] and Khanafer and Varfai [29].

3.6. Characteristic Length Scales of Temperature Distributions in Tissue

The temperature distributions in living tissues is related to number of spacial scales, which mainly are taken into account in the Chen-Holmes model (see Remark 1 in section 3.4.1). Let us stress the attention on some of them [53].

3.6.1. Range where the heat transfer is controlled by conduction

The Pennes equation , for instance defines the length scale

$\begin{array}{ll}{L}_D~\sqrt{\frac{a_t}{{\omega }_b}}& (62)\end{array}$

where a_t is the tissue thermal diffusivity, defined with assumptions that ρ_b ≈ ρ_t and c_b ≈ c_t. This length scale is associated with the time scale ${\tau }_{conduction}=L_D^2/a_t\approx {\omega }_b$.

3.6.2. Length scale of the blood-tissue thermal equilibrium

Another, characteristic length scale L_v is considered as the mean distance in a single vessel where the arterial flow achieves thermal equilibrium the surrounding tissue when moving from the large arteries to the arterioles. This length scale, actually, equals the length scale L_II characterizing the distance along vessel axis at which the blood and tissue become in thermal equilibrium. In accordance with [24, 46] the definition is

$\begin{array}{ll}{L}_{II}~\left(\frac{R_v^2}{a_t}{v}_b\right)ln\left(\frac{d_{mean}}{R_v}\right)& (63)\end{array}$

where R_v is the radius of vessel, d_mean is the mean distance between vessels of the same length L; v_b is the blood velocity (averaged over the vessel cross-section). Following Lubashevsky and Gafiychuk [53] and Witmore [54] commonly L ~ d_mean an therefore taking into account that L_v ~ L_II we get

$\begin{array}{ll}{L}_v~\frac{a_t}{{\omega }_b}\left[ln{\left(\frac{d_{mean}}{R_v}\right)}^{-1}\right],{\omega }_b~\frac{R_v^2{v}_b}{L^3}& (64)\end{array}$

since the total blood flow through the vessel is ${\pi }^2{R}_v^2{v}_b$ and consequently $d_{mean}^2{L}_v~L_v^3$.

The value related to L_v directly affects the tissue thermal diffusivities [53]: typically $a_t\approx 2\times 10^{-7}{m}^2/s$, while ${\omega }_b\approx 6\times 10^{-2}{s}^{-1}$ and d_mean/R_v ≈ 40. This yields length scales L_D ≈ 0.6 cm. and L_v ≈ 0.3 cm, while the time scale τ_cond ≈ 3 min.

3.7. Some Comments on the Continuum Bio-Heat Models

The continuum models permit different features to be modeled but under correctly defined circumstances. Here we stress the attention on the comments of Crezee and Lagendijk [55] (even though they repeat some previous standpoints), among them:

(1) The basic problem is that the structure of vascularized tissues is too complex consisting of solid and fluid compartments with architecture that cannot be easily described and modeled. In the solid section the conduction is the dominant heat transfer mode while convection takes part in the fluid filled compartments. The vast number of small blood transporting vessels the modeler task to incorporate the role of all of them, down to singular capillaries, in the structures of mathematical models.

(2) Albeit the fact that it would be possible to require exact information on particular vascular geometry at micro-level, it becomes difficult to model it and the most important issue is that the results are unpractical despite the huge amounts of computing resources involved in the simulation process.

In the context of the previous comments, the adequate models should compromise between what we know and what is reasonable in the model build-up. Following this standpoint and looking at the modeling process pragmatically we may point out some comments on the existing continuum models, that is

(3) The first attempt of Pennes is criticized, as commented above, focusing on the basic assumption incorporated in the perfusion term: the blood enter the tissue with the arterial temperature and leave it with the averaged local temperature. From this, the main critical points are:

(a) the heat transfer related to the blood mass transport is neglected;

(b) the actual blood temperature entering a particular tissue compartment is not accounted for.

(c) the individual heat transfer processes of individual tissue such as cooling or heating are not incorporated in the model.

(d) the entire venous vessel architecture is neglected with the assumption that there is entrance lengths to attain thermal equilibria.

Despite these critical notes the Pennes model is widely popular due to its simple mathematical structure [56–61].

The models created after Pennes are based on effective transport coefficients, particularly on effective heat conductivity k_eff:

(4) The Wulff's model is the first incorporating such a transport coefficient by a substitution of the Pennes' perfusion term by the product ρ_bc_bU · gradT, thus introducing the mean fluid velocity U as a variable. Actually, this is an attempt to create a model of heat transfer in living tissues, similar to the ones known from the convection heat transfer with unidirectional flow; but, the blood transport in tissues is not unidirectional.

(5) The models of Chen-Holmes and Weinbaum-Jiji also use k_eff, precisely an enhanced thermal conductivity (tensor). This approach requires details about every component in the vascularized structure, but assuming isotropy, for the sake of simplicity in modeling, the conductivity tensor reduces to scalar k_eff.

After these comments we may write the bio-heat transfer model in a generalized in the form proposed by Crezee and Lagendijk [55]

$\begin{array}{ll}{c}_t{\rho }_t\frac{\partial T}{\partial T}=\left(\nabla k_{eff}\nabla T\right)-f_c{c}_b{\rho }_b(T-T_a)+q_m& (65)\end{array}$

where k_eff and < f_c ≤ 1 are phenomenological parameters.

Thus, repeating the Pennes construction, the generalized equation (65) is only an instructive heat balance model allowing on its basis particular bio-heat models to be developed: after successive procedures allowing to establish the model structures by direct averaging microscopic governing equations (see more comments in [53]).

At the end of this section we may conclude that creations of continuum based models of bio-heat transfer need serious preliminary analyzes of the particular tissue of interest and the associate blood transporting architectures of vessels.

4. Bio-heat Models Using the Porous Media Concept

All continuum models commented above, in their origins, start with criticisms of the flaws of the Pennes models and attempts to improve its deficiencies. As we saw each models is applicable to a particular type of living tissue and a generalization is hard to be attained. Alternatively, to the continuum models, a concept considering the matrix of tissue, arteries, veins and capillary vessels a porous medium with specific porosity variations, effective heat conductivities and heat dispersion by blood flow have been developed [21, 29, 52, 62–66]. Following this concept the thermal energy transfer between the blood flow and the tissue is modeled with the basic assumption of thermal non-equilibrium [29, 67–69] described by two energy equations for the blood phase and the solid phase [29], namely (in original notations)

Blood phase

$\begin{array}{l}\epsilon {\left(\rho c\right)}_b\left[\frac{\partial {\langle T\rangle }^b}{\partial t}+{\langle U\rangle }^b·\nabla {\langle U\rangle }^b\right]=\nabla ·\left[k_b^a·\nabla {\langle T\rangle }^b\right]+\\ +h_{bs}\left[{\langle T\rangle }^s-{\langle T\rangle }^b\right]& (66)\end{array}$

Solid porous phase

$\begin{array}{l}\left(1-\epsilon \right){\left(\rho c\right)}_s\frac{\partial {\langle T\rangle }^s}{\partial t}=\nabla ·\left[k_s^a·\nabla {\langle T\rangle }^s\right]-\\ -h_{bs}\left[{\langle T\rangle }^s-{\langle T\rangle }^b\right]+(1-\epsilon )q_m& (67)\end{array}$

Here 〈T〉^b, 〈T〉^s, ${\text{k}}_b^a$, ${\text{k}}_s^a$, h_bs are local volume-averaged arterial blood temperature, local volume-averaged solid-tissue temperature, blood effective thermal conductivity tensors (for the blood and solid phases, respectively), and the interstitial convective heat transfer coefficient, correspondingly. Further, ε and 〈U〉^b are the solid phase porosity and the blood velocity vector (all vectors are denoted as bold letters), accordingly. The locally averaged interstitial convective heat transfer coefficient h_bs depends on the blood velocity and solid phase geometry; the heat transfer between the phases is proportional to the difference $h_{bs}\left[{\langle T\rangle }^s-{\langle T\rangle }^b\right]$. 〈T〉^b and 〈T〉^s are locally averaged blood and solid tissue temperatures, while ${\text{k}}_b^a$ and ${\text{k}}_s^a$ the blood and solid tissue effective conductivity tensors (averaged), as mentioned above.

With the assumption of isotropic heat conduction the effective conductivities are

$\begin{array}{ll}{k}_b^a=\epsilon k_b+k_b^t,k_s^a=(1-\epsilon )k_s& (68)\end{array}$

where ${\text{k}}_b^t$ denotes the thermal dispersion conductivity as a natural consequence of the already existing models of heat transfer in porous media [67, 68].

For tissues with small-sized blood vessels when ε << 1 the flowing blood is thermally equilibrated , then using (66) and (67) the following simplified model can be developed (in original notations)

$\begin{array}{l}{\left(\rho c\right)}_bϵ+(1-\epsilon ){\left(\rho c\right)}_s\frac{\partial \langle T\rangle }{\partial t}+\underset{(bloodperfusion)}{\underbrace{\epsilon {(\rho c)}_b{\langle U\rangle }^b·\nabla \langle T\rangle }}=\\ =\nabla \left(k_s^a+k_s^a\right)·\nabla \langle T\rangle +q_m(1-\epsilon )& (69)\end{array}$

The blood perfusion term in (69) corresponds to the ideas developed by Wulff [32], Charny [21], and Klinger [36, 37, 40, 41] in contrast to the Pennes model where a uniform blood perfusion ωρ_bc_b(T_a,in − T_v,out) is assumed in the model construction.

4.1. Nakayama-Kuwahara Porous Medium Approach

The porous media concept in bio-heat transfer uses the volume averaging theory (VAT)[52, 65, 66], a technique well-established in the area of fluid-saturated porous media. Looking at the anatomy of the living tissues, three principle compartments can be identified: blood vessels, cells and interstitium (see Figure 6). For the sake of simplicity the volume averaging theory applicable to living tissues [52, 65, 66] considers only two regions: the vascular region and extravascular region (cells and interstitium). This approach considers the entire biological structure as a fluid-saturated porous medium through which the blood infiltrates. The extravascular region is regarded as a solid matrix (despite the fact that exatravascular liquid exists).

FIGURE 6

Figure 6. A schematic presentation of a control volume V in a tissue with blood supply through arteries and veins in parallel relevant to the studies of porous media approach.

The volume averaging needs the control volume V to be representative as characteristics but enough small in order to apply the tools of integer-order calculus. In such a case the volume averaging of certain variable ϕ is $\langle \varphi \rangle \equiv \frac{1}{V}{\displaystyle \underset{V}{\int }\varphi dV}$ and therefore for the region occupied by the blood we have ${\langle \varphi \rangle }^b\equiv \frac{1}{V_b}\underset{V_b}{\int }\varphi dV$. The relationship between the two averaged values is 〈ϕ〉 = ε〈ϕ〉^b, where ε = V_b/V defines the local porosity as a ratio of the volume of the fluid (blood) to the entire volume in the vascular space. In the dominating cases ε < 0.1 [65].

Further, considering the macroscopic governing equations of continuity, Navier-Stokes equation and the energy equation, the Nakayama-Kuwahara model consists of two energy equations based on VAT [52, 65, 66], namely (in original notations-take into account the summation notations)

Blood phase

$\begin{array}{l}{\epsilon }_b{\rho }_b{c}_b\frac{\partial {\langle T\rangle }^b}{\partial t}+{\rho }_b{c}_b\frac{\partial }{\partial x_j}\langle U_j\rangle \\ {\langle T\rangle }^b=\frac{\partial }{\partial x_j}\left(\epsilon k_b\frac{\partial {\langle T\rangle }^b}{\partial x_j}+\epsilon k_{dis\left(j,k\right)}\frac{\partial {\langle T\rangle }^b}{\partial x_k}\right)-\\ -a_b{h}_b\left({\langle T\rangle }^b-{\langle T\rangle }^s\right)-{\rho }_b{c}_b\omega \left({\langle T\rangle }^b-{\langle T\rangle }^s\right)& (70)\end{array}$

where k_{dis(j, k)} is the dispersion related heat conductivity [70]

The solid tissue phase

$\begin{array}{l}\left(1-\epsilon \right){\rho }_s{c}_s\frac{\partial {\langle T\rangle }^s}{\partial t}=\frac{\partial }{\partial x_j}\left[\left(1-\epsilon \right)k_s\frac{\partial {\langle T\rangle }^s}{\partial x_j}\right]\\ +a_b{h}_b\left({\langle T\rangle }^b-{\langle T\rangle }^s\right)+{\rho }_b{c}_b\omega \left({\langle T\rangle }^b-{\langle T\rangle }^s\right)+\left(1-\epsilon \right)q_m& (71)\end{array}$

Equations (70) and (71) are correct for all cases of thermal non-equilibrium (for more details about their derivations we refer [65]). Now, the interesting part of this section is to see how these equations are related to the models discussed earlier and derived from a continuum point of view.

4.2. Relations of the Porous Media Concept to the Existing Continuum Models

4.2.1. Pennes Model in Porous Media Terms

In terms of (70) and (71) the Pennes model can be expressed as [52, 65]

$\begin{array}{l}\left(1-\epsilon \right){\rho }_s{c}_s\frac{\partial {\langle T\rangle }^s}{\partial t}=\frac{\partial }{\partial x_j}\left[\left(1-\epsilon \right)k_s\frac{\partial {\langle T\rangle }^s}{\partial x_j}\right]\\ +{\rho }_b{c}_b{\omega }_{Pennes}\left(T_{a0}-{\langle T\rangle }^s\right)+\left(1-\epsilon \right)q_m& (72)\end{array}$

where ω_Pennes is the mean blood perfusion rate; T_a0 is the mean branchial artery temperature Comparing to the solid phase energy equation (71) it is obvious that

$\begin{array}{l}{\rho }_b{c}_b{\omega }_{Pennes}\left(T_{a0}-{\langle T\rangle }^s\right)=a_b{h}_b\left({\langle T\rangle }^b-{\langle T\rangle }^s\right)\\ +{\rho }_b{c}_b\omega \left({\langle T\rangle }^b-{\langle T\rangle }^s\right)& (73)\end{array}$

where a_b is the specific surface area contacting with the blood.

Assuming that $T_{a0}\simeq {\langle T\rangle }^b$ for small vessels [65] it follows that

$\begin{array}{ll}{\omega }_{Pennes}\simeq \omega +\frac{a_b{h}_b}{{\rho }_b{c}_b}& (74)\end{array}$

Hence, the Pennes perfusion rate is an effective value including the surface heat transfer [65]. The second term in (74) is a ratio of two components of the perfusion heat flux [see (73)]. Therefore the assumption for complete thermal equilibrium is valid when the Peclet number is too small (the diffusion transport dominates and the advection could be neglected).

4.2.2. Klinger and Wulff Models in Porous Media Terms

In terms of (70) and (71) the Klinger model is (in original notations)

$\begin{array}{l}\left(1-\epsilon \right){\rho }_s{c}_s\frac{\partial {\langle T\rangle }^s}{\partial t}=\frac{\partial }{\partial x_j}\left[\left(1-\epsilon \right)k_s\frac{\partial {\langle T\rangle }^s}{\partial x_j}\right]-\\ -{\rho }_b{c}_b\frac{\partial \langle U_j\rangle {\langle T\rangle }^s}{\partial x_j}+\left(1-\epsilon \right)q_m& (75)\end{array}$

Suggesting a thermal equilibrium, that is 〈T〉^b = 〈T〉^s [65] it follows from (70) and (71) that

$\begin{array}{l}\left[{\epsilon }_b{\rho }_b{c}_b+\left(1-\epsilon \right){\rho }_s{c}_s\right]\frac{\partial {\langle T\rangle }^s}{\partial t}+{\rho }_b{c}_b\frac{\partial }{\partial x_j}\langle U_j\rangle {\langle T\rangle }^s=\\ =\frac{\partial }{\partial x_j}\epsilon k_b+\left(1-\epsilon \right)k_s\frac{\partial {\langle T\rangle }^s}{\partial x_j}+\\ +\epsilon k_{dis\left(j,k\right)}\frac{\partial {\langle T\rangle }^s}{\partial x_k}+\left(1-\epsilon \right)q_m& (76)\end{array}$

Equation (76) reduces to the Klinger model (41) when the porosity (that is the ratio of vascular volume to the total volume) is sufficiently small. Further, when the blood flow is strong and the macroscopic diffusion transport can be neglected, then energy equation for the blood phase (70) becomes

$\begin{array}{ll}{\rho }_b{c}_b\frac{\partial }{\partial x_j}\langle U_j\rangle {\langle T\rangle }^b=-a_b{h}_b\left({\langle T\rangle }^b-{\langle T\rangle }^s\right)-{\rho }_b{c}_b\omega \left({\langle T\rangle }^b-{\langle T\rangle }^s\right)& (77)\end{array}$

It is easy to see that with the approximation (77) the energy equation for the solid phase (71) reduces to the Klinger equation (76). Hence, the condition the Peclet number to be sufficiently larger than unity (i.e., dominating convective transport with respect the negligible conduction heat transfer) leads to the Klinger model; the high Peclet number condition relates the Klinger model to the heat transport in large vessels.

5. Integer-Order Models With Relaxations (Hyperbolic Models)

The simplest model of heat conduction is based on the Fourier law assuming that the heat flux is related to the temperature gradient

$\begin{array}{ll}q=-k\frac{\partial T}{\partial x}& (78)\end{array}$

via the heat conductivity k as a transport coefficient. This a reliable model confirmed by experiments, but there are serious limitations for its applications. Precisely (78), is an adequate model of the microscopic phenomenon of heat diffusion in cases when the length scales are much greater than the mean free path and the time scales are much greater then the thermal relaxation time (see the excellent comments in [71]). However, one unphysical outcome, such as with the diffusion models based on the Fick's law, is the infinite speed of temperature disturbances, that is the heat flux has no damping. Physically, temperature distribution is due to motion of particle carriers such as electrons and quanta such as photons [71] and reasonably the realistic picture is the thermal disturbances to propagate with finite speed.

The models of Pennes, Wulff, Klinger, Chen-Holmes and Nakayama are based on the Fourier law leading to parabolic equation with infinite speed of the heat flux. However, in real systems the flux propagates with finite speed. In metals, for example, the relaxation times are of orders from 10⁻⁸ to 10⁻¹⁴s [72–74], while in biological tissues the relaxation times are of order of 10–30 s [72–75]. In biological tissues, for example, in muscles under local strong heating the experiments reveal anisotropy in heat transport that cannot be explained by the Fourier law only [76–78]. This leads to formulation of thermal wave bio-heat models [79, 80] based on two main approaches:

(1) Maxwell-Cattaneo approach with heat flux time lag, known also as Single-Phase Lag (SPL) approach [81, 82].

(2) Double-Phase-Lag (DPL) approach with relaxations in both the heat flux and temperature gradient propagation [81–84].

Both approaches lead to hyperbolic equations with finite speeds of the solutions discussed in the sequel of this section.

5.1. Single-Phase Lag Heat Conduction (Maxwell-Cattaneo Approach)

5.1.1. Basic Approach in Modeling of Heat Flux Relaxation by SPL Approach

Considering the need of heat flux relaxation (damping) the simplest approach is to assume first order approximation (via Taylor series expansion) of the flux with respect to the time [81–84]. Hence,the extension of the Fourier law (see as an example [85]) as a model with a single time-lag of the flux is

$\begin{array}{ll}q(x,t+{\text{λ}}_q)=q(x,t)+{\text{λ}}_q\frac{\partial q(x,t)}{\partial t}=-k\frac{\partial T}{\partial x},{\text{λ}}_q=\frac{a}{C^2},{\text{λ}}_q>0& (79)\end{array}$

According to this postulation the temperature gradient established at a point x at time t gives rise to a heat flux vector at x at a later time t + λ_q, that is there is time-shift between the heat flux and the temperature gradient. The heat flux relaxation time λ_q, called also a phase-lag parameter has been interpreted from different points of view [1]:

(i) from the framework of the second law of the irreversible thermodynamics [81, 86],

(ii) λ_q comes from the phase-lag shift between the heat flux vector and the temperature gradient [86] in high-rate responses;

(iii) λ_q is a constant representative for the interaction of different inner structural element of a material with the heat flux [73] and may take large values within the range (10⁻³ − 10³s) [74].

By using the formal temporal translation $\acute{t}=t+{\lambda }_q$ (79) may be written in the equivalent delayed form [81]

$\begin{array}{ll}q(x,t)=-k\frac{\partial T(x,\acute{t}-{\text{λ}}_q)}{\partial x},\acute{t}>{\text{λ}}_q& (80)\end{array}$

Actually, SPL is a result of a first order Taylor expansion and this model is consistent with the Second law of Thermodynamics [87, 88].

Equations (79) and (80) are the forward and backward forms of SPL model, respectively [81]. For λ_q < 0 the role of cause and effect are exchanged in corresponding forward and backward forms [81].

Obviously, there are no restrictions to expand the flux (the forward SPL problem) with a truncated Taylor series of various orders. Thus, with second-order of Taylor expansion we have [81]

$\begin{array}{ll}q+{\text{λ}}_q\frac{\partial q}{\partial t}+\frac{{\text{λ}}_q^2}{w}\frac{{\partial }^{2}q}{\partial t^2}+O({\text{λ}}_q^3)=-k\frac{\partial T}{\partial x},{\text{λ}}_q>0& (81)\end{array}$

The divergence of (79) and by help of the energy balance (see Equation 84) yields a hyperbolic heat equation

$\begin{array}{ll}{\text{λ}}_{qp}\frac{{\partial }^{2}T}{\partial t^2}+{}_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial x^2}& (82)\end{array}$

The solution of (82) is a weakly-damped wave equation [71] where the temperature disturbance speed is $v_T=\sqrt{k/{\lambda }_q{c}_p}$.

Particularly, the first order approximation yields (79) which is equivalent to the constitutive equation with fading memory (see section 9)

$\begin{array}{ll}q(t)=-{\displaystyle {\int }_0^{\infty }\stackrel{\sim }{R}}(t,s)\frac{\partial }{\partial x}T(t,s)ds,\stackrel{\sim }{R}(t,s)=\frac{k}{{\text{λ}}_q}{e}^{-\frac{t-s}{{\text{λ}}_q}}& (83)\end{array}$

The hyperbolic heat equation (82) can be regarded as a singular perturbation of the Fourier heat conduction model, recovered from (82) for λ_q → 0. This is a formal limit, as mentioned in [71] because any solution of (82) reduces to a solution of the Fourier heat conduction equation for λ_q → 0. This is a very sensitive situation since if λ_q cannot be considered as mathematically small the resulting solutions are unphysical (see detailed analysis in [71] and the comments in [89, 90]).

Similar interrelations between Taylor series expansions describing short-memory effects and integral constitutive equations with rapidly fading memory are widely encountered in other physical applications displaying memory effects such as viscoelasticity and electromagnetic wave propagation [81].

In the following parts of this section we will use the forward forms of SPL and DPL because only these versions are related to the existing time-fractional bio-heat models (section 6.2.1) as well to the systematic approach in equation fractionalizations (developed in section 7.2).

5.1.2. SPL (Heat-Wave) Model Without Blood Perfusion Relaxation

With (79) and applying the energy balance equation

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial T}{\partial t}=-\frac{\partial q}{\partial x}& (84)\end{array}$

we get a hyperbolic version of the Pennes model [22, 91] where both the blood perfusion term and the metabolic heat source are not affected by the thermal relaxation.

$\begin{array}{ll}{\text{λ}}_q\frac{{\partial }^{2}T}{\partial t^2}+\left(1+{\omega }_b\frac{{\rho }_b{c}_b}{{\rho }_t{c}_t}\right)\frac{\partial T}{\partial t}+{\omega }_b\frac{{\rho }_b{c}_b}{{\rho }_t{c}_t}\left(T_t-T_b\right)=a\frac{{\partial }^{2}T}{\partial x^2}+\frac{q_m}{{\rho }_t{c}_t}& (85)\end{array}$

Hereafter, in this section, we will continue with the two-dimensional model of Bravo et. al.[92] in a dimensionless form, namely

$\begin{array}{ll}D{e}_T\frac{\partial {\theta }^2}{\partial {\tau }^2}+\frac{\partial \theta }{\partial \tau }=\left[\frac{{\partial }^2\theta }{\partial {\xi }^2}+\frac{{\partial }^2\theta }{\partial {\eta }^2}\right]+\sigma \left({\theta }_a-\theta \right)=\left({\Phi }_m+{\Phi }_r\right)& (86)\end{array}$

where the dimensionless variables are defined as

$\begin{array}{ll}\theta =\frac{T}{T_c},\xi =\frac{x}{L},\eta =\frac{y}{L},\tau =\frac{at}{L^2},\Phi =\frac{L^2{q}_m}{k{T}_c}& (87)\end{array}$

In the model (86) the relaxation time ${\lambda }_q=a/C^2$ is defined (as in the linear viscoelasticity [93–95] and heat waves [96, 97]) as a ratio of the thermal diffusivity a to the thermal velocity of propagation in the medium C. The relaxation time λ_q allows definition of the thermal version of the Deborah number $D{e}_T={\lambda }_q(a/L^2)$ as a ratio of λ_q to the characteristic time of the heat diffusion $L^2/a_t$. Moreover, Φ is the dimensionless metabolic heat, despite the fact that the authors termed this terms as diffusive flux of energy.

The contribution of the blood perfusion (as in the model of Pennes) is weighted by the dimensionless factor σ = Lω_b/a_t. Precisely, in [92] σ is presented as $\sigma =\frac{L^2{\omega }_b}{a_t}$ which actually is a ratio of the blood perfusion rate ${\omega }_b\left[s^{-1}\right]$ to a velocity scale (volume of fluid flow per volume of tissue) defined through the thermal diffusivity as $u_T=a_t/L\left[s^{-1}\right]$ without physical meaning). Comparing (85), (86), and (87) it is easy to see how the dimensionless groups can be defined.

In the context of dimensional consistency, let us look at Equation (85). The fist term in left-hand side has dimension s⁻¹ (if we assume conditionally that T is dimensionless only for these comments), while the second and third terms are dimensionless since ω_b in the prefactors has dimension s⁻¹. The dimension of the first term in left-hand side is balanced by the dimension of first term in right-hand side (because a_t has a dimension m²s⁻¹).

With the assumption to represent the solution as a sum of steady state θ₀ and transient terms θ_t

$\begin{array}{ll}\theta (\xi ,\eta ,\tau )={\theta }_0(\xi ,\eta )+{\theta }_t(\xi ,\eta ,\tau )exp(-a\tau )& (88)\end{array}$

the modeling dimensionless model and associated boundary conditions are

$\begin{array}{l}D{e}_T\frac{\partial {\theta }_{\tau }}{\partial {\tau }^2}+\left[1-2\sigma D{e}_T\right]\frac{\partial {\theta }_{\tau }}{\partial \tau }+{\sigma }^{2}D{e}_T=\frac{{\partial }^2\theta }{\partial {\xi }^2}+\frac{\partial \theta }{\partial {\eta }^2}\\ +{\Phi }_{r}exp(\sigma t)& (89)\end{array}$

$\begin{array}{ll}\eta =1,-\frac{d{\theta }_t(\xi ,1,\eta )}{d\eta }=\left[{\Phi }_s-B{i}_0{\theta }_{air}-{\theta }_{0s}\right]exp\left(\sigma \tau \right),& (90)\end{array}$

$\begin{array}{l}\eta =0,{\theta }_t(\eta ,0,\tau )=0\xi =0,\frac{d{\theta }_{\tau }}{d\xi }=0,\\ \xi =1,\frac{d{\theta }_{\tau }}{d\xi }=0& (91)\end{array}$

The Biot number at steady-state Bi₀ is a measure of the thermal resistance of the tissue when the heat flux is imposed at its surface, while θ_air is the dimensionless temperature of the ambient air. In this context, the last term in (90) allows to be altered with small perturbation of the atmosphere, i.e., as in an air stream exhibiting a stochastic behavior. Hence, the model is oriented to modeling thermal reactions of living tissues when external fluid (air) with varying temperature is supplying or removing thermal energy from bodies. From the analyzes above, it becomes clear that this is a model oriented to bio-heat problems appearing close to the skin surface. Second, in this model construction, the blood perfusion is modeled as in the Pennes model without detailed analysis of the vascular architecture. Third, the definition of the Deborah number needs experimentally defined values of the thermal relaxation time of the tissue λ and the speed of thermal wave in the same medium.

The solutions performed in Bravo [92] with Biot number Bi₀ less and almost equal to unity, indicate that they are related to cases with negligible thermal gradients inside the tissue. Furthermore, the range of variation in the Deborah number was chosen from De_T = 0 (this means classical Fourier law without flux lagging) up to De_T = 0.0013, that should mean fast thermal relaxation of the tissue. However, it is hard to judge these results since neither λ nor C are presented numerically in this work.

At the end, of this subsection, we may say that, actually, the only new element is the application of the SPL concept (Maxwell-Cattaneo approach) for heat flux relaxation. The definition of the Deborah number as it is done is not useful since by definition it should be a ratio of relaxation time λ_q to the current time t. Precisely, with a general definition of the time scale t₀ we have for the left-hand side (after nondimensalization with respect to the temperature θ = T/T_c and the space coordinates)

$\begin{array}{ll}{\text{λ}}_q\frac{\partial \theta }{\partial t^2}+\frac{\partial \theta }{\partial t}\Rightarrow {\text{λ}}_q\frac{1}{t_0^2}\frac{\partial \theta }{\partial {\tau }^2}+\frac{1}{t_0}\frac{\partial \theta }{\partial t}\Rightarrow \frac{1}{t_0}\left[\left(\frac{{\text{λ}}_q}{t_0}\right)\frac{{\partial }^2\theta }{\partial {\tau }^2}+\frac{\partial \theta }{\partial \tau }\right]& (92)\end{array}$

Then, multiplying both sides of the dimensionless form of the equation by $t_0=\frac{L^2}{a}$ we get (86). Hence, the thermal Deborah number De_T defined in Bravo [92] is a ratio of the characteristic heat flux relaxation time to the characteristic thermal diffusion time. Now, making dimensionless both the nominator and denominator of De_T we get

$\begin{array}{ll}D{e}_T=\frac{{\lambda }_q}{t_0}=\frac{{\lambda }_q/t}{t_0/t}=\frac{De}{1/Fo}& (93)\end{array}$

Hence, the dimensionless prefactor De_T is a product of two dimensionless times (dimensionless numbers), that is De_T = DeFo.

We can see, that for small De_T, which is important when λ_q < < t₀, the first term in the model (86) can be neglected while the parabolic term ∂θ/∂τ determines the model behavior with dominating diffusion and perfusion terms in the right-hand side of (86). This can be simply tested by dividing both sides of (86) by De_T. We have to stress the attention on the fact that the definition of De_T, which can be considered as a fixed value of the Fourier number Fo for t = λ_q, is by analogy, taking similarity of the relaxation process in viscoelastic transient flow, as mentioned above. It is not related to material flow, as mentioned in [92] since the heat conduction equation as (86) does not encounter a flow behavior of the material.

In the context of the SPL approach we have to mention the convolution form of the thermal flux damping (relaxation) with exponential memory kernel [96, 97], (related to 83), namely

$\begin{array}{ll}{q}_{relax}=-\frac{k}{\tau }{\displaystyle {\int }_0^t}{e}^{t-\tau }\nabla (x,\tau )d\tau & (94)\end{array}$

which is erroneously represented as

$\begin{array}{ll}{q}_{relax}(x,t)=-\frac{k}{\tau }{e}^{-t/\tau }{\displaystyle {\int }_0^t}{e}^{s/\tau }\nabla (x,s)ds& (95)\end{array}$

following [98], thus destroying (violating) the memory effect (damping) in the flux relaxation model. The coefficient k in (94) should has a dimension [(W/m · K) · s]. It will be of a special interest further in this article when the fading memory formalism will be applied as a technique allowing correct fractionalization of model equations [see the SPL relationship (94) and the comments about (83)].

The Maxwell-Cattaneo approach was also used by Tang et al. [99] in solution of surface heating of biological tissues, that is excluding the perfusion effect from the problem. Well-performed analytical-numerical study of skin burn injury induced by radiation heating by 3-D model applying the Maxwell-Cattaneo approach was carried out by Dai et al. [100]; see also the results of Liu et al. [79].

5.1.3. Heat Wave Model With a Blood Perfusion Relaxation

The Maxwell-Cattaneo approach uses only the first order approximation of the heat flux (79) but as commented earlier there are no restrictions more terms of the Taylor series expansion to be involved [101] in the extension of the Fourier law, namely

$\begin{array}{ll}q(x,t+{\tau }_q)\approx ={\displaystyle \sum \frac{{\tau }_q^n}{n!}}\frac{d^n}{d{t}^n}q(x,t)=-k\nabla T(x,t)& (96)\end{array}$

Then, using the Pennes construction of the heat flux and its component (conduction, perfusion and metabolic source), and the external heating q_r, we may write [101]

$\begin{array}{l}{\rho }_t{c}_t\frac{\partial T_t(x,t+{\tau }_q)}{\partial t}=k{\nabla }^2{T}_t+{\omega }_b{c}_b(x,t+{\tau }_q)\\ \left[T_b-T_t(x,t+{\tau }_q)\right]+\\ +q_m(x,t+{\tau }_q)+q_r(x,t+{\tau }_q)& (97)\end{array}$

For the sake of simplicity of the analysis let us assume that q_m = 0 and q_r = 0. Then, following [101] let us consider only the temperature elevation θ = T(x, t) − T(x, 0) above the steady-state version of (97). In this case (subtracting (97) from the steady-state version) the resulting equation is

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial {\theta }_t(x,t+{\tau }_q)}{\partial t}=k{\nabla }^2{\theta }_t-{\omega }_b{c}_b(x,t+{\tau }_q)\theta (x,t+{\tau }_q)& (98)\end{array}$

Assuming the traveling approach in solution of (98) (not presented here) Liu et al. [101] concluded in their analysis that the higher blood perfusion rate the higher frequency in osculation of the tissue temperature, an effect that cannot be accounted by the classical models of Pennes, Klinger, Wulff, and Chen-Holmes.

This is only an attempt in modeling the integer-order bio-heat transfer by application of the Maxwell-Cattaneo approach with a step ahead considering a relaxation of the blood perfusion term. No more studies in this direction cannot be detected in the literature, thus we stop comments on this approach. However, we will come back to this idea, in a different form, when the bio-heat models using fractional-differential operators will be used to model the finite speed of heat flux.

5.2. Dual-Phase-Lag Heat Conduction Concept

As an extension of the Maxwell-Cattaneo approach, which cannot work correctly when in the medium there are micro-structural interaction effects, the so-called Dual-Phase-Lag (DPL) was proposed by Tzou [84] (see also the analysis and solutions in [83, 102]) as

$\begin{array}{ll}q+{\lambda }_q\frac{\partial q}{\partial t}=-k\frac{\partial T}{\partial x}-k{\lambda }_T\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial x}\right)& (99)\end{array}$

This construction of the heat flux-temperature gradient relationship can be easily understand if the Taylor series expansion is simultaneously applied to both sides of the Fourier law (see the beginning section 5), that is in terms of first order of approximations we get

$\begin{array}{ll}\underset{Fourierlaw}{\underbrace{q=-k\frac{\partial T}{\partial x}}}\underset{\begin{array}{l}Taylor\\ series\end{array}}{\underbrace{\Rightarrow }}\underset{DPLmodel}{\underbrace{q+{\lambda }_q\frac{\partial q}{\partial t}=-k\left[\frac{\partial T}{\partial x}+\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial x}\right)\right]}}& (100)\end{array}$

In the DPL model (100) the relaxation times λ_q and λ_T can be interpreted as the periods arising from thermal inertia of the materials and the micro-structural interactions [1]. Especially λ_T (λ_q was discussed earlier) can be considered as heat diffusion with sharp wave front (the hyperbolic equations exhibit finite speed solutions) induced by the appearance of λ_q [83]. From (100) it follows that at any point of the heat conductor the heat flux vector at t + λ_q corresponds to a temperature gradient at time t + λ_T [1, 103].

The basic condition [84] (see also in [81]) is that both λ_q and λ_T are positive and the delay (the time-shift between λ_q and λ_T) should be also positive, that is

$\begin{array}{ll}{\lambda }_d\approx {\lambda }_q-{\lambda }_T>0& (101)\end{array}$

This condition allows to avoid violation of the common physical causality between ∂T/∂x and q [81]. In this context, if a temporal translation $\acute{t}=t+{\lambda }_T$ or $\acute{t}=t+{\lambda }_q$ is applied, then (99) can be presented as the SPL model with a time delay λ_d = λ_q − λ_T [81]. The stability of the DPL heat conduction models and the thermodynamic restrictions are analyzed in Fabrizio and Franchi [81], Quintanilla and Racke [103], Fabrizio and Lazzari [104], and Fabrizio et al. [105], but we skip these problems since they are out of the scope of this work.

It is noteworthy that when λ_q = λ_T the result of lagging disappears and we get a diffusion model of heat conduction [84]. For example, the experiments of Tang et al. [106] with IR irradiated tissue revealed that λ_q = 20 s and λ_T = 14 s (data obtained in the first 50 s of the tissue irradiation) and these data confirm the wave-like behavior of the heat transport [84].

Hence, by first and second order (or higher) approximations via Taylor series it is possible to construct different phase-lag heat conduction models with several relaxation times thus accounting for heat transport through materials with complex structures. Hereafter in this section, we will focus the attention only on simple phase-lag models (obtained by first order approximations). For more complex models of such a type we refer Xu et al. [1].

5.2.1. Dual-Phase Lag Bio-Heat Model

The approach used in (100) to model relaxations by mixed time-space derivative is well-known from the viscoelastic models [93–95] and was used also in construction of models with fractional derivatives [94, 95] (see further in this article).

Now, with the Pennes model [14] presented as

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial T_t}{\partial t}=-k\frac{\partial q}{\partial x}+{\omega }_b{\rho }_b{c}_b(T_b-T_t)& (102)\end{array}$

where T_b is the blood temperature in the artery, that is T_b − T_a0 in the Pennes' model. Then, expressing the heat flux as the second version of (100) we get [14]

$\begin{array}{ll}\begin{array}{c}{\lambda }_q\frac{{\partial }^2{T}_t}{\partial t^2}+\left(1+{\lambda }_q\frac{{\omega }_b{\rho }_b{c}_b}{{\rho }_t{c}_t}\right)\frac{\partial T_t}{\partial t}+\frac{{\omega }_b{\rho }_b{c}_b}{{\rho }_t{c}_t}\left(T_b-T_t\right)=\\ =a\left(1+{\lambda }_T\frac{\partial }{\partial t}\right)\frac{\partial T_t}{\partial x^2}+\frac{q_m}{{\rho }_t{c}_t}\end{array}& (103)\end{array}$

with initial conditions T(x, 0) = T₀ and T(x, 0) = 0 in accordance with the version of the model solved.

In (103) the product $\left({\lambda }_q\frac{{\omega }_b{\rho }_b{c}_b}{{\rho }_t{c}_t}\right)$ in the pre-factor of the second term in the left-hand side is completely dimensionless since the product λ_qω_b is dimensionless.

Two simple examples of DPL model solutions are briefly presented next.

5.2.1.1. Constant flux as boundary condition to a tissue of final depth

With a boundary condition $-k\frac{\partial T}{\partial x}=q_0$ to a tissue with final depth L and introducing generalized time (t₀), length (L_x), temperature (T_scale) and heat flux (q_scale) scales (and order of magnitude analysis, as it was already done above) the following dimensionless variables can be defined as [14]

$\begin{array}{ll}\begin{array}{c}\xi =x\sqrt{\frac{W_b{c}_b}{k}},{\xi }_L=L\sqrt{\frac{W_b{c}_b}{k}},\Rightarrow L_x=\sqrt{\frac{k}{W_b{c}_b}}\\ \eta =t\left(\frac{W_b{c}_b}{{\rho }_t{c}_t}\right)\Rightarrow t_0=\frac{{\rho }_t{c}_t}{W_b{c}_b}\\ \theta \left(\xi ,\eta \right)=\frac{T\left(x,t\right)-T_b}{T_{scale}},T_{scale}=\frac{q_0}{\sqrt{k{W}_b{c}_b}}\\ {\Lambda }_q=\frac{{\lambda }_q}{t_0},{\Lambda }_T=\frac{{\lambda }_T}{t_0}\\ \psi =q_m\left(\frac{1}{q_0}\frac{k}{W_b{c}_b}\right)\Rightarrow q_{scale}=q_o\sqrt{\frac{W_b{c}_b}{k}}\end{array}& (104)\end{array}$

Hence, the dimensionless models is

$\begin{array}{l}A\frac{{\partial }^2\theta }{\partial {\eta }^2}+B\frac{\partial \theta }{\partial \eta }+\theta =\left(1+D\frac{\partial }{\partial t}\right)\frac{{\partial }^2\theta }{\partial {\xi }^2}+\psi \\ 0\le A={\Lambda }_q\le 1,1\le B=1+{\Lambda }_q\le 2,\\ 0\le D=\frac{{\Lambda }_T}{{\Lambda }_q}=\frac{{\lambda }_T}{{\lambda }_q}\le 1& (105)\end{array}$

with dimensionless initial and boundary conditions

$\begin{array}{l}\theta (\xi ,\eta =0,{\theta }_{\eta }(\xi ,\eta )=0\\ {\theta }_{\xi }(0,\eta )=-1,{\theta }_{\xi }({\xi }_L,\eta )=0& (106)\end{array}$

The Laplace transform solution of (105), (106) in the s space is [14]

$\begin{array}{ll}\begin{array}{c}\overline{\theta }(\xi ,s)=\frac{cosh\left[\sqrt{\beta }(\xi -{\xi }_L)\right]}{s\sqrt{\beta }sinh\left(\sqrt{\beta }{\xi }_L\right)}+\frac{\psi }{s(1+Ds)\beta }\\ \beta (s)=\frac{E{s}^2+Fs+G}{1+Hs}\end{array}& (107)\end{array}$

The time-domain solution needs inversion via the Bromwich integral but we will skip these huge expressions [14] since this is out of the scope of present work where model constructions and nondimensalizations are at issue.

Let us now see how the relaxation times and how their ratios control the contributions of different heat transfer effects, as terms, in the model equation. Dividing all terms in (105) by A we get

$\begin{array}{ll}\begin{array}{c}\frac{{\partial }^2\theta }{\partial {\eta }^2}+\frac{B}{A}\frac{\partial \theta }{\partial \eta }+\theta =\frac{1}{A}\left(1+D\frac{\partial }{\partial t}\right)\frac{{\partial }^2\theta }{\partial {\xi }^2}+\frac{1}{A}\psi \\ \frac{1}{A}=\frac{1}{{\Lambda }_q},\frac{B}{A}=1+\frac{1}{{\Lambda }_q},\frac{D}{A}=\frac{{\Lambda }_T}{{\Lambda }_q}\end{array}& (108)\end{array}$

The ratio $\frac{A}{A}=\frac{{\Lambda }_T}{{\Lambda }_q}=\frac{{\lambda }_t}{{\lambda }_q}$ has order of magnitude of unity. Thus, the diffusion term in the right-hand side of (108) has a strong influence on heat transport (by heat conduction). Further, B/A = 1 + 1/Λ_q and for very small values of λ_q we have B/A >> 1, that is, there is a negligible inertia in the heat flux propagation; then the second term in the left-hand side of (108) dominates and determines a parabolic behavior of the model (while the first term $\frac{{\partial }^2\theta }{\partial {\eta }^2}$ can be neglected. Otherwise, for high values of Λ_q the model has strong hyperbolic behavior. In this case, the contribution of metabolic terms is low but cannot be neglected since this means eliminating of the metabolic heat from the model, which is unphysical.

5.2.1.2. Periodic surface heat flux

The boundary condition q = q₀exp(iωt) and dimensionless variables

$\begin{array}{ll}\begin{array}{c}\xi =x\sqrt{\frac{\omega }{a}},{\xi }_L=L\sqrt{\frac{\omega }{a}}\Rightarrow L_x=\sqrt{\frac{a}{\omega }}\\ \eta =t,,{\tau }_1=\omega {\lambda }_q,{\tau }_2=\omega {\lambda }_T\Rightarrow t_0=\frac{1}{\omega }\\ \theta (\xi ,\eta )=k\frac{T(x,t)-T_b}{q_0}\sqrt{\frac{\omega }{a}}\Rightarrow T_{scale}=\frac{q_0}{k}\sqrt{\frac{a}{\omega }}\\ {\Lambda }_0=\frac{W_b{c}_b}{{\rho }_t{c}_t\omega },\phi =\frac{q_m}{q_0}\sqrt{\frac{k}{{\rho }_t{c}_t\omega }}\end{array}& (109)\end{array}$

allow to express the governing energy equation as

$\begin{array}{ll}E\frac{{\partial }^2\theta }{\partial {\eta }^2}+F\frac{\partial \theta }{\partial t}+G\theta =\left(1+H\frac{\partial }{\partial \eta }\right)\frac{{\partial }^2\theta }{\partial {\xi }^2}+\phi & (110)\end{array}$

with initial conditions

$\begin{array}{ll}\theta (\xi ,0)={\theta }_{\xi }(\xi ,0)=0& (111)\end{array}$

where E = τ₁, F = 1 + Λ₀τ₁, G = Λ₀ and H = τ₁.

In the Laplace domain the solution of (110), (111) is

$\begin{array}{ll}\begin{array}{c}\overline{\theta }(\xi ,s)=\frac{cosh\left[\sqrt{\gamma }(\xi -{\xi }_L)\right]}{(s-i)\sqrt{\gamma }sinh(\sqrt{\gamma }{\xi }_L)}+\frac{\phi }{s(1+Hs)\gamma },\\ \gamma (s)=\frac{E{s}^2+Fs+G}{1+Hs}\end{array}& (112)\end{array}$

The simulations in Askarizadeh and Ahmadikia [14] used value of the surface heat flux $q_0=5\times 10^{3}W/m^2$ (a value typical for the heat flux released by fire [107]).

5.2.2. Closing Briefs on the DPL Models

At the end of this section we have to refer to some additional studies applying the dual-phase-lag approach to modeling bio-heat problem. The dominating studies are on the Pennes model, among them: laser radiation heat supply to the tissue [108], temperature-dependent perfusion [109], skin heat transfer [14], pulsed IR irradiation [106], cryosurgery of lung cancer [110], laser heating [111, 112], interaction between tissue and a cryoprobe [113, 114], human head relaxation times [115]. The main idea of the DPL, that is the Taylor series expansions of the integer-order heat conduction equation, thus creating damping and heat wave behavior, is fruitfully used in fractionalization of bio-heat models (see section 8.4).

6. Fractional Bio-heat Equations: The State-of-the-Art in Model Formulations

This section is devoted to the existing fractional models of bio-heat transfer and the emphasis is on the technology of fractionalization, not on solution techniques. The correct model build-up, with correct application of both physical laws (and hypothesizes) and mathematical tools is of primary importance prior to attack the equations with simple or sophisticated solution techniques. With modern mathematical tools, either analytical or numerical any equation could be solved to greater extent. However, the principle question is: Does this equation is the adequate and correct model of what you need to represent in the mathematical space by it? We try to answer this question when analyzing the existing fractional versions of bio-heat equations.

6.1. Scaling and Nondimensalization of Fractional Operators

Before starting any analyzes of time-fractional bio-heat models we have to see how scaling and nondimensalizations have to carried out when time-fractional operators are used. The next section presents the common fractional operators used and the techniques of nondimensalization (based on the recent work [116]).

6.1.1. Time-Fractional Operators With Singular (Power-Law) Memory

We start with same definitions which will be useful in reading the following text even though they are available in many texts devoted to fractional modeling and textbooks [117]

6.1.1.1. Fractional integral

The Riemann-Liouville fractional integral of order α > 0 is a natural result of the Cauchy multiple integral the m-fold primitive of a function f(t) expressed as a single integral convolution for arbitrary positive number α > 0 [117]

$\begin{array}{ll}{{}_{0}I}^{\alpha }f(t)=\frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t{(t-z)}^{\alpha -1}f(z)dz,t>0,n\in R^{+}}& (113)\end{array}$

The law of exponents for fractional integrals is: ${{}_{0}D}^{-\alpha }{{}_{0}D}^{-\gamma }f(t)={{}_{0}D}^{-\alpha -\gamma }f(t)={{}_{0}D}^{-\gamma }{{}_{0}D}^{-\alpha }f(t)$

The Laplace transform of the fractional integral through the convolution theorem yields

$\begin{array}{ll}L{[}_0{D}_t^{\alpha }f(t)]=L\left[\frac{t^{\alpha }}{\Gamma (\alpha )}\right]L[f(t);s]=s^{-\alpha }F(s)& (114)\end{array}$

where LT[R(s) > 0], and F(s) = LT[f(t)] are Laplace transforms.

6.1.1.2. Riemann-Liouville time-fractional derivative

The Riemann-Liouville time-fractional derivative is defined as

$\begin{array}{l}\begin{array}{l}{{}_{0}D}^{\alpha }f(t)={{}_{0}D}^m{{}_{0}I}^{m-\alpha }f(t)=\frac{1}{\Gamma (m-\alpha )}\frac{d^m}{d{t}^m}{\displaystyle {\int }_0^t\frac{f(z)}{{(t-z)}^{\alpha +1-m}}dz,}\\ m-1<\alpha <m,m\in \text{N}& (115)\end{array}\end{array}$

6.1.1.3. Caputo time-fractional derivative

The Caputo derivative of a casual function f(t) is defined [117] as

$\begin{array}{l}{{}^{C}D}_t^{\alpha }f(t)={{}_{0}I}^{m-\alpha }\frac{d^m}{t^m}f(t)={{}_{0}D}_t^{-(m-\alpha )}{f}^{(m)}=\\ \frac{1}{\Gamma (m-\alpha )}{{\displaystyle \int }}_0^t\frac{f^{(m)}(z)}{{(t-z)}^{\alpha +1-m}}dz,m-1<\alpha <m& (116)\end{array}$

For m = 1 we have the common definition

$\begin{array}{ll}\begin{array}{l}{{}^{C}D}_t^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}{{\displaystyle \int }}_0^t\frac{1}{{(t-s)}^{\alpha }}\frac{df(s)}{ds}ds\end{array}& (117)\end{array}$

Caputo derivative of a constant is zero, i.e., ${{}^{C}D}_t^{\alpha }C=0$ [117]. Further, when f(0) = f′(0) = f″(0) = … = f⁽ⁿ⁾(0) = 0, then the Riemann-Liouville and the Caputo derivatives coincide.

Let us see the nondimensalization procedure to ${{}^{C}D}_t^{\alpha }f(t)$, which is uncommon in the existing literature [117] but is strongly related to the analysis of the existing fractional bio-heat models (see section 6.2). The Caputo derivative of the function U(t) is

$\begin{array}{ll}{{}^{C}D}_t^{\alpha }U(t)=\frac{1}{\Gamma (1-\alpha )}\underset{0}{\overset{t}{{\displaystyle \int }}}\frac{1}{{(t-s)}^{\alpha }}\frac{d}{ds}U(s)ds& (118)\end{array}$

Changing the variables as Y = U/U₀ and $\bar{t}=t/T_R$, with U₀ and T_R as characteristic scales of the process, so that 0 < Y < 1, $0<\bar{t}<1$, yields: U = U₀ ⇒ dU/dt = U₀dY/dt ⇒ dU/dt = U₀T_R,

Further, the nondimensalization of (t − α)^{− α} yields $\frac{{(t-\alpha )}^{-\alpha }}{T_R^{\alpha }}$. Then, we get

$\begin{array}{ll}{{}^{C}D}_t^{\alpha }\overline{U}(\overline{t})=\frac{U_0}{{T_R}^{\alpha }}\frac{1}{\Gamma (1-\alpha )}\underset{0}{\overset{\overline{t}}{{\displaystyle \int }}}\frac{1}{{(\overline{t}-\overline{s})}^{\alpha }}\frac{dY(\overline{s})}{d\overline{s}}d\overline{s}& (119)\end{array}$

If the nondimensalization is only with respect to the time, than we get

$\begin{array}{ll}{{}^{C}D}_t^{\alpha }U(\overline{\overline{t}})=\frac{1}{{T_R}^{\alpha }}\frac{1}{\Gamma (1-\alpha )}\underset{0}{\overset{\overline{t}}{{\displaystyle \int }}}\frac{1}{{(\overline{t}-\overline{s})}^{\alpha }}\frac{dU(\overline{s})}{d\overline{s}}ds& (120)\end{array}$

The result of nondimensalization can be easy attributed to the scale-invariant property of the power-law function used as memory kernel as commented in Hristov [116].

Alternatively, when the power law is expressed through the ratio (t/τ)^{− α} where τ are the process relaxation time; the process time scale is T_R. The introduction of the relaxation time leads to τ^αt^−α and since τ = const. this is directly related to the scale invariance of the power-law function. With the nondimensalization procedure, as it was done above, we get

$\begin{array}{ll}{{}^{C}D}_t^{\alpha }U(\overline{t})={\left(\frac{\tau }{T_R}\right)}^{\alpha }\frac{1}{\Gamma (1-\alpha )}{\displaystyle \underset{0}{\overset{\overline{t}}{\int }}\frac{1}{{(\overline{t}-\overline{s})}^{\alpha }}\frac{dU(\overline{s})}{d\overline{s}}d\overline{s}}& (121)\end{array}$

that is, an expression similar to (120),

Hence, we got a dimensionless scaling coefficient ${(\tau /T_R)}^{\alpha }$ which could be considered as a weighting factor, too.

6.1.1.4. Caputo-Fabrizio fractional operator

The Caputo and Fabrizio [118] time-fractional operator with exponential kernel is defined as

$\begin{array}{l}{{}^{CF}D}_t^{\alpha }f(t)=\frac{M(\alpha )}{1-\alpha }{\displaystyle {\int }_0^{t}exp\left[-\frac{\alpha (t-s)}{1-\alpha }\right]\frac{df(t)}{dt}ds\Rightarrow }\\ LT{[}_{CF}^C{D}_t^{\alpha }f(t)]=\frac{p{L}_T[f(t)-f(0)]}{p+\alpha (1-p)}& (122)\end{array}$

where the normalizing function M(α) should obey the conditions M(0) = M(1) = 1, and $LT\left[{{}_{CF}{}^{C}D}_T^{\alpha }f(T)\right]$ is the operator Laplace transform.

From the definition (122) it follows that if f(t) = C = const., then ${{}_{CF}D}_t^{\alpha }f(t)=0$, an expected results as in the classical Caputo derivative [117]. Analyzes of applications of the Caputo-Fabrizio operator are available in Hristov [96] and Hristov [119–121] and will skip reference quotations here.

The constitution of the Caputo-Fabrizio operator the stretched time is multiplied by a dimensional factor α/(1 − α) and therefore it should have a dimension of inverse time (t⁻¹). Actually α is a dimensionless parameter and it is related to the relaxation time [119] via nondimensalization of the exponential relaxation kernel represented as exp[(t − s)/τ] [119–121] precisely nondimensalization of the exponential function, yields

$\begin{array}{ll}\begin{array}{c}exp\left(\frac{t-s}{\tau }\right)=exp\left(\frac{t/t_0-s/t_0}{\tau /t_0}\right)=exp\left(\frac{\bar{t}-\bar{s}}{\bar{\tau }}\right)\Rightarrow \\ \Rightarrow \frac{1-\alpha }{\alpha }=\bar{\tau }=\frac{\tau }{t_0}\Rightarrow \alpha =\frac{1}{1+\tau /t_0}\end{array}& (123)\end{array}$

where t₀ is the characteristic time of the relaxation process.

The nondimensalization of the kernel preserves the meaning of the exponential memory function and avoids doubts about the definition of the fractional order α as Hristov [119] (see also the next comments).

Now, let us repeat the basic steps of nondimensalization of the convolution integral, as it was done in section 6.1.1.3. With Caputo-Fabrizio operator applied to a function U(t) and dimensionless variables $Y=U/U_0,\bar{t}=t/T,0<Y<1,0<\bar{t}<1$ we have defined as [116]

$\begin{array}{ll}{}^{CF}{{\displaystyle D}}_t^{\alpha }U(t)=\frac{1}{1-\alpha }{\displaystyle \underset{0}{\overset{t}{\int }}\exp \left[-\frac{\alpha }{1-\alpha }\left(t-s\right)\right]}\frac{d}{ds}U(s)ds& (124)\end{array}$

where T_R is the characteristic time-scale of the process.

After nondimensalization we get [116]

$\begin{array}{l}{{}^{CF}D}_t^{\alpha }U(\overline{t})\Rightarrow U_0\left\{\frac{1}{1-\alpha }{\displaystyle \underset{0}{\overset{\overline{t}}{\int }}\exp \left[-\frac{\alpha }{1-\alpha }(\overline{t}-\overline{s})\right]\frac{dY(\overline{s})}{d\overline{s}}d\overline{s}}\right\}\Rightarrow \hfill \\ {\Rightarrow }^{CF}{D}_t^{\alpha }U(\overline{t})=U_0{[}^{CF}{D}_t^{\alpha }Y(\overline{t})]\hfill & (125)\end{array}$

and the Laplace transform is

$\begin{array}{ll}LT\left[{{}_{CF}{}^{C}D}_t^{\alpha }f(t)\right]=U_{0}LT\left\{\left[{{}_{CF}{}^{C}D}_t^{\alpha }U(t)\right]\right\}& (126)\end{array}$

When the nondimensalization is only with respect to the time by $\bar{t}=t/T$ we get

$\begin{array}{l}{}^{CF}{{\displaystyle D}}_t^{\alpha }U(\overline{t})\Rightarrow \frac{1}{1-\alpha }{\displaystyle \underset{0}{\overset{\overline{t}}{\int }}\exp \left[-\frac{\alpha }{1-\alpha }\left(\overline{t}-\overline{s}\right)\right]}\frac{dU\left(\overline{s}\right)}{d\overline{s}}d\overline{s}\\ {\Rightarrow }^{CF}{D}_t^{\alpha }U(t)& (127)\end{array}$

Therefore, no additional prefactors emerge in front of the convolution integral unlike the case of the Caputo-derivative with a power-law kernel. This result can be attributed to the invariant properties of the exponential function used as memory kernel. The construction of Caputo-Fabrizio operator by definition, implicitly, contains a dimensionless kernel since it comes from the classical (exp − t/τ), which is dimensionless. Additional scaling inside the kernels such as (t/T) / (τ/T) does not change the ratio (exp − t/τ), but yields the product $(1/\overline{\tau })\left[(-t/T)\right]$. Now, $(1/\overline{\tau })\in [0,\infty )$ can be to be replaced by the ratio [α/(1 − α)] ∈ [0, ∞), where α ∈ [0, 1].

6.2. Existing Formally Fractionalized Versions: Analysis and Scaling

6.2.1. Ezzat's Models

Ezzat et al. [122, 123] considered the classical Maxwell-Cattaneo construction of the heat flux dependence on the temperature gradient q(x, t + τ) = −k∇T(x, t). Recall, that in the Maxwell-Cattaneo approach, a first order approximation $q\approx q+\tau \frac{\partial q}{\partial t}$ (local expansion in a Taylor series) is used to reach the hyperbolic heat conduction equation (section 5.1). Mimicking this approach, Ezzat et al. [122, 123] applied fractional Taylor series expansion (of order α) (see section 8.1 for details about fractional Taylor series) with modified Riemann-Liouville derivatives of Jumarrie [124] [see the definition (130)] and the heat flux-temperature gradient was expressed as

$\begin{array}{ll}q+\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }q}{\partial t^{\alpha }}=-k\nabla T,\frac{{\partial }^{\alpha }q}{\partial t^{\alpha }}=I^{\alpha -1}\left[\frac{\partial q}{\partial t}\right],0<\alpha <1& (128)\end{array}$

where the Jumarries' fractional Taylor series are defined as [124]

$\begin{array}{l}f(x+h)={\displaystyle \sum _{k=0}^{\infty }\frac{h^{\alpha k}}{\Gamma (1+\alpha k)}}{f}^{(\alpha k)}(x)={\displaystyle \sum _{k=0}^{\infty }\frac{h^{\alpha k}}{(\alpha k)!}}{f}^{(\alpha k)}(x),\\ 0<\alpha <1& (129)\end{array}$

through the modified Jummarie derivative (130)

$\begin{array}{ll}{{}_{J}f}^{\alpha }(x)=\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}{\displaystyle {\int }_0^t{(x-z)}^{-\alpha -1}[f(x)-f(0)]dz,0<\alpha <1}& (130)\end{array}$

The second version of (129) comes from the relationship Γ(1 + αk) = :(αk)!

Recall, that τ^α/α! should be equivalent to λ_q in the simple DPL model (79) and we may use the symbol ${\lambda }_q^{\alpha }$. The dimension of ${\lambda }_q^{\alpha }={\tau }^{\alpha }/\alpha !$ is s^α. The relationship (128) for α → 1 reduces to the integer-order Cattaneo construction (79) as mentioned above. Next, with the Fourier law construction (131) at hands, and taking the partial time-derivative of order α the result is (132)

$\begin{array}{l}\rho c\frac{\partial T\left(x,t\right)}{\partial t}=-\nabla q+Q\left(x,t\right),Q\left(x,t\right)=W_b{c}_b\\ \left[T_b-T\left(x,t\right)\right]+Q_m& (131)\end{array}$

$\begin{array}{ll}\rho c\frac{\partial }{\partial t}\left[\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}T\left(x,t\right)\right]=-\nabla \left(\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}q\right)+\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}Q\left(x,t\right),0<\alpha <1& (132)\end{array}$

Then, multiplying (132) by τ^α/α! both sides of (133) and applying (128) in the right-hand side, the final Ezzat's equation is

$\begin{array}{ll}\rho c\left(T+\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }T}{\partial t^{\alpha }}\right)=k{\nabla }^{2}T+\left(Q+\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }Q}{\partial t^{\alpha }}\right),0<\alpha <1& (133)\end{array}$

Actually, the Ezzat's equation is a result of a formal fractionalization of the Pennes' equation by fractional Taylor series (first order expansion) (about fractional Taylor series expansion see section 8.1) of the heat flux and accepting the concept of Wulff (see the preceding text) about the perfusion term (i.e., the reference temperature is the arterial blood temperature T_a). With the pre-factors τ^α/α! in both sides of (133) the equation is dimensionally consistent but why these pre-factors should exist and what they model are questions unanswered in the studies of Ezzat et al. [122, 123]. We will discuss this problem further in section 8.2.

The scaling by $\theta =\frac{T-T_a}{T_0}$ and the length, and time scales (L_x and t₀) yields

$\begin{array}{l}\left(\theta +\frac{1}{t_0^{\alpha }}\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }\theta }{\partial {\overline{t}}^{\alpha }}\right)=\frac{a}{L_x^2}{\nabla }^2\theta +\left(1+\frac{1}{t_0^{\alpha }}\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }}{\partial {\overline{t}}^{\alpha }}\right)\frac{Q}{\rho c{T}_0},\\ \overline{t}=\frac{t}{t_0},a=\frac{k}{\rho c}& (134)\end{array}$

Since the ratio $a/L_x^2$(a = k/ρc is the thermal diffusivity of the tissue) has a dimension of inverse time [1/s] it follows that with a time scale $t_0=L_x^2/a$ we get the dimensionless ${(\tau /t_0)}^{\alpha }$ relaxation time-delay (the term is in accordance with basic postulate of Cattaneo), or in other words a fractional Fourier number defined through τ, that is ${}_{\tau }{}^{\alpha }Fo=a\tau /L_x^2$. Moreover, ${}_{\tau }{}^{\alpha }Fo=(\frac{\tau }{t})(at/L_x^2)=DeFo$. Hence, the ratio ${(\tau /t_0)}^{\alpha }$ in the Ezzat model (see (134) is

$\begin{array}{ll}{\left(\frac{\tau }{t_0}\right)}^{\alpha }={\left(\frac{\tau }{t}\right)}^{\alpha }{\left(\frac{t}{t_0}\right)}^{\alpha }={\left(\frac{t}{t_0}\right)}^{\alpha }{\left(\frac{at}{L_x^2}\right)}^{\alpha }={\left(DeFo\right)}^{\alpha }& (135)\end{array}$

For α = 1 in (135) we get the integer-order relationship (93). In case of t₀ = τ we get that ${(\tau /t_0)}^{\alpha }=1$ and it follows from (135) that ${(\tau /t_0)}^{\alpha }={(Fo)}^{\alpha }$. This leads to simplified terms of the time-fractional derivatives in (134), that is

$\begin{array}{ll}\frac{1}{t_0^{\alpha }}\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }\theta }{\partial {\overline{t}}^{\alpha }}\Rightarrow \frac{1}{\alpha !}\frac{{\partial }^{\alpha }\theta }{\partial F{o}^{\alpha }}& (136)\end{array}$

This result answers the question previously raised about appearance of the pre-factors in the Ezzat model (134). The answer is simple: incomplete nondimensalization of the model equation. In the above formula we did not use the correct expression of ${\lambda }_q^{\alpha }$, omitting the dimensionless α!, only for the sake of simplicity. Actually, the left-hand-side of (128) can be presented as

$\begin{array}{ll}q+\frac{{\tau }^{\alpha }}{\alpha !}\frac{{\partial }^{\alpha }q}{\partial t^{\alpha }}=q+{\lambda }_q^{\alpha }\frac{{\partial }^{\alpha }q}{\partial t^{\alpha }}& (137)\end{array}$

, where the fractional relaxation time is ${\lambda }_q^{\alpha }={\tau }^{\alpha }/\alpha !$.

Now, the problem remaining is the definition of the length scale L_x. From the general requirement all terms to have order of magnitude of unity and dimensionless pre-factors we get that (Q/ρcT₀) is dimensionless. Hence, the length scale cannot be defined through the known parameters and it has to be taken from the experiments, depending on the precision of measurements. In this context, we have to remember that the magnitude of the time delay (relaxation time) of the tissue should be also determined from experiments parallel to adequate value of the fractional order α, thus recalling for accurate inverse problems taking into account their ill-posed natures.

6.2.2. Damor's Model

The model of Damor et al. [125] is a fractional version of the bio-heat equation by a simple replacement the time derivative with a fractional counterpart as Caputo derivative (order α ∈ (0, 1]) with singular (power-law) kernel and the spatial derivative by a Riesz-Feller fractional derivative of order β ∈ (0, 2], namely

$\begin{array}{l}{\rho }_t{c}_t\frac{{\partial }^{\alpha }{T}_t\left(x,t\right)}{\partial t^{\alpha }}=k_{\alpha ,\beta }\frac{{\partial }^{\beta }{T}_t\left(x,t\right)}{\partial x^{\beta }}+W_b{c}_b\left(T_a-T_t\right)+q_m,\\ 0<\alpha <1,0<x<L& (138)\end{array}$

where W_b and c_b are related to the blood; the subscript t is related to the tissue, as in the model of Pennes.

With initial and boundary conditions (139)

$\begin{array}{ll}T\left(x,0\right)=T_a,{T\left(x,t\right)|}_{x\to \pm \infty }=0& (139)\end{array}$

Otherwise, for α = 1 and a Riemann-Liouville space derivative of order γ ∈ (0, 2] [see (141)] the formal fractionalization results in [126]

$\begin{array}{ll}{\rho }_t{c}_t\frac{\partial T_t\left(x,t\right)}{\partial t}=k_{\gamma }\frac{{\partial }^{\gamma }{T}_t\left(x,t\right)}{\partial x^{\gamma }}+W_b{c}_b\left(T_a-T_t\right)+q_m& (140)\end{array}$

$\begin{array}{ll}\frac{{\partial }^{\gamma }{T}_t\left(x,t\right)}{\partial x^{\gamma }}=\frac{1}{\Gamma \left(2-\alpha \right)}\frac{d}{dx}{\displaystyle \underset{0}{\overset{x}{\int }}\frac{T_t\left(t,s\right)}{{\left(x-s\right)}^{\gamma }}}ds,0<\gamma \le <2& (141)\end{array}$

For β = 2 the model (138) reduces to [125]

$\begin{array}{l}{\rho }_t{c}_t\frac{{\partial }^{\alpha }{T}_t\left(x,t\right)}{\partial t^{\alpha }}=k\frac{{\partial }^2{T}_t\left(x,t\right)}{\partial x^2}+W_b{c}_b\left(T_a-T_t\right)+q_m,t>0,\\ 0<x<L,0<\alpha <1& (142)\end{array}$

With dimensionless variables [125]

$\begin{array}{l}\xi =x{\left(\frac{W_b{c}_b}{k_{\alpha }}\right)}^{1/2},\eta =t{\left(\frac{W_b{c}_b}{{\rho }_t{c}_t}\right)}^{1/\alpha },\theta =\frac{T-T_a}{T_0},\\ \varphi =\frac{Q_{met}}{T_0{W}_b{c}_b}& (143)\end{array}$

equation (142) takes a dimensionless form, namely

$\begin{array}{ll}\frac{{\partial }^{\alpha }\theta }{\partial t^{\alpha }}+\theta =\frac{{\partial }^2\theta }{\partial {\xi }^2}+\varphi & (144)\end{array}$

with initial and boundary conditions

$\begin{array}{ll}\theta \left(\xi ,0\right)=0,{\theta \left(\xi ,\eta \right)|}_{\xi \to \pm \infty }=0& (145)\end{array}$

The dimensionless variables (143) reveal that the length and time scales are

$\begin{array}{ll}{L}_x={\left(\frac{k_{\alpha }}{W_b{c}_b}\right)}^{1/2},t_0={\left(\frac{{\rho }_t{c}_t}{W_b{c}_b}\right)}^{1/\alpha },{\theta \left(\xi ,\eta \right)|}_{\xi \to \pm \infty }=0& (146)\end{array}$

The nondimensalization with the dimensionless variables (143) seems, to some extent, artificial (without explanations, performed in a typical mathematical manner without physical explanations). To be precise, if we use the generally defined length and time scales as (L_x and t₀), and with classical dimensionless temperature θ = (T − T_a)/T₀ or θ = T/T₀ (in this case there is only a shift in the temperature scale origin), and applying the classical techniques of nondimensalization of differential equations we get from (142)

$\begin{array}{l}\left[{\rho }_t{c}_t\frac{T_0}{{\left(t_0\right)}^{\alpha }}\right]\frac{{\partial }^{\alpha }\theta \left(\xi ,\eta \right)}{\partial \eta }=k_{\alpha }\left[\frac{T_0}{{\left(L_x\right)}^2}\right]\frac{{\partial }^2{T}_t\left(\xi ,\eta \right)}{\partial {\xi }^2}\\ +W_b{c}_b\theta T_0+q_m,{\theta \left(\xi ,\eta \right)|}_{\xi \to \pm \infty }=0& (147)\end{array}$

Now, dividing all terms by the prefactor of the time fractional derivative $\left[{\rho }_t{c}_t\frac{T_0}{{(t_0)}^{\alpha }}\right]$ we get

$\begin{array}{l}\frac{{\partial }^{\alpha }\theta \left(\xi ,\eta \right)}{\partial {\eta }^{\alpha }}=\left[\left(\frac{k_{\alpha }}{{\rho }_t{c}_t}\right)\frac{1}{{\left(L_x\right)}^2}{\left(t_0\right)}^{\alpha }\right]\frac{{\partial }^2{T}_t\left(\xi ,\eta \right)}{\partial {\xi }^2}\\ +\left[\frac{W_b{c}_b}{{\rho }_t{c}_t\frac{1}{{\left(t_0\right)}^{\alpha }}}\right]\theta +\frac{q_m}{{\rho }_t{c}_t{T}_0}{\left(t_0\right)}^{\alpha }& (148)\end{array}$

Since at this moment all terms in (148), including the variable θ, should be of order of magnitude of unity and applying the order of magnitude analysis requirement we may obtain the length and time scales of the process. Hence, setting the pre-factor of the second term in the right-hand side $(W_b{c}_b/{\rho }_t{c}_t){(t_0)}^{\alpha }=1$ we get the time scale ${(t_0)}^{\alpha }=({\rho }_t{c}_t/W_b{c}_b)\Rightarrow t_0={({\rho }_t{c}_t/W_b{c}_b)}^{1/\alpha }$. Further, from $(k_{\alpha }/{\rho }_t{c}_t)(1/{(L_x)}^2){(t_0)}^{\alpha }=1$ we have $t_0={\left[{(L_x)}^2/a_{\alpha }\right]}^{1/\alpha }$, where a_α = k_α/ρ_tc_t is the fractional thermal diffusivity of the tissue with a dimension [m²/s^α]. Consequently, the requirement all prefactors of the dimensionless equation to be equal to 1 yields $L_x={(k_{\alpha }/W_b{c}_b)}^{1/2}$.

Therefore, the characteristic time and length scales depend on preliminary known process parameters and there is nothing magic in formulation of the dimensionless variables (143) (do not believe that they are formulated ad hoc). Moreover, since $t_0={({\rho }_t{c}_t/W_b{c}_b)}^{1/\alpha }$, consequently η^α defines the fractional Fourier number. ${}^{\alpha }Fo=(t/t_0)=\left[t/{({\rho }_t{c}_t/W_b{c}_b)}^{1/\alpha }\right]$. Thus, the time-fractional derivative in (144) can be expressed as

$\begin{array}{ll}\frac{{\partial }^{\alpha }\theta \left(\xi ,\eta \right)}{\partial {\eta }^{\alpha }}=\frac{{\partial }^{\alpha }\theta \left(\xi ,\eta \right)}{\partial \left({}^{\alpha }Fo\right)}& (149)\end{array}$

However, we refer the dimension of the fractional heat conductivity k_α for correctness and physical relevance of the model. When we define the fraction thermal diffusivity a_α = k_α/ρ_tc_t it is obvious that it obeys the dimensional homogeneity of the time-fractional diffusion equation. Nevertheless, the problems is to the initial definition of k_α. In the Damor's model this was not clarified and on this basis Ferras et al. [127] criticized it (see the next section 6.2.3). We may say briefly, that the crucial problem is the correct constitutive equation of the heat flux in terms of fractional integral with adequate memory function and we will discuss this problem at large further in this article.

6.2.3. Ferras et al. Model

Ferras et al. [127] adapted the Pennes equation by using timer-fractional Caputo derivative in the form

$\begin{array}{l}\frac{{\partial }^{\alpha }T\left(x,t\right)}{\partial t^{\alpha }}=A\frac{\partial }{\partial x}\left(k\left(x\right)\frac{\partial T\left(x,t\right)}{\partial x}\right)-BT\left(x,t\right)+C,\\ 0<t<T^{*},0<x<L,0<\alpha <1& (150)\end{array}$

where

$\begin{array}{ll}A=\frac{1}{{\rho }_t{c}_t}\frac{1}{{\tau }^{\alpha -1}},B=\frac{W_b{c}_b}{{\rho }_t{c}_t}\frac{1}{{\tau }^{\alpha -1}},C=\frac{W_b{c}_b+q_m}{{\rho }_t{c}_t}\frac{1}{{\tau }^{\alpha -1}}& (151)\end{array}$

It is worth-mentioning that in this model a new parameter τ (this is not the past time in the fractional operator) was added to assure the dimensional consistency. Thus, this is the same step as in the approach of Ezzat commented above [see the operation transforming (132) into (133)] where the multiplication by τ^α/α! is applied.

It is worth noting that Ferras et al. [127] focused the attention on the fading memory approach following Gurtin [128] and Gurtin and Pipkin [128] where the non-locality of the heat flux is expressed by a convolution integral

$\begin{array}{ll}q=-k{\displaystyle \underset{0}{\overset{\infty }{\int }}K\left(t-\tau \right)}\frac{\partial }{\partial x}T\left(x,\tau \right)d\tau & (152)\end{array}$

With this step Ferras et al. [127] got a heat conduction equation with memory.

$\begin{array}{ll}\frac{\partial T\left(x,t\right)}{\partial t}=a{\displaystyle \underset{t}{\overset{0}{\int }}K\left(t-\tau \right)}\frac{{\partial }^2}{\partial x^2}T\left(x,\tau \right)d\tau ,a=\frac{k}{\rho c}& (153)\end{array}$

Note 3: However, we only mention at this point that the model (153) is not complete, as it should follow from the fading memory approach of Gurtin [128] and Gurtin and Pipkin [128], because the right-hand side of (153) is only the elastic part of the heat flux [97] (see also the analysis in [96]).

The next step of Ferras et al. [127] is the assumption that the memory kernel in (153) is a of power-law type, that resulted in

$\begin{array}{ll}q=-\frac{k_{pk}}{\Gamma \left(\alpha \right)}\frac{d}{dt}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\left(t-\tau \right)}^{\alpha -1}\frac{\partial }{\partial x}T\left(x,\tau \right)}d\tau ,0<\alpha <1& (154)\end{array}$

which, in fact, is that the heat flux is proportional to a Riemann-Liouville derivative of the temperature gradient. This is ad hoc formulation of a constitutive equation about heat flux, without any physical background, but in the sequel of this section we will see way this was done.

Then, with (154) the model (153) was transformed into

$\begin{array}{ll}\frac{{\partial }^{\alpha }T\left(x,t\right)}{\partial t^{\alpha }}=a_{\alpha }\frac{\partial }{\partial x}\left(\frac{\partial T\left(x,t\right)}{\partial x}\right)& (155)\end{array}$

The result (154) is, in fact, the time-fractional diffusion equation where the thermal diffusivity a_α should have a dimension [m²/s^α]. Hence, we may see why (154) was constituted; the explanation is simple: to get an equation that we know how to solve.

Ferras et al. [127] raised a principle question about the construction of their model (150) by introducing the new parameter in the coefficients [see (151)] which physical meaning was difficult to be explained. Moreover, Ferras et al. [127] claimed that a new relationship between the heat flux and temperature gradient is established (153) and (154), which actually is not true since the fading memory approach of Gurtin and Pipkin [128] is incompletely used . We will explain these critical standpoints further in this article but now let us turn on last result of Ferras et al. [127]

$\begin{array}{ll}{\rho }_t{c}_t\frac{{\partial }^{\alpha }T\left(x,t\right)}{\partial t^{\alpha }}=k_{\alpha }\frac{\partial }{\partial x}\left(\frac{\partial T\left(x,t\right)}{\partial x}\right),k_{\alpha }=\frac{k_{pk}}{{\tau }^{\alpha -1}}& (156)\end{array}$

defining in this way a “new fractional conductivity” $k_{\alpha }=k_{pk}/{\tau }^{\alpha -1}$.

Note 4: The approach of Ferras et al. [127] is very dubious because it uses only the pseudo-elastic part of the heat flux (see further the section devoted to fading memory concept) expressing it as a Riemann-Liouville derivative (154) that leads to a result which is a well-known, i.e., the time-fractional diffusion equation. Now, the principle question is: What would be the result is the Gurtin and Pipkin [128] approach is correctly applied ? This problem is discussed in the next section where a systematic approach in construction of bio-heat equation with memories involving time-fractional operators is drawn.

6.3. Briefs on the Existing Time-Fractional Bio-heat Models

Fractional models form a modern trend in modeling of transport processes with relaxations. In this context, the use of time-fractional derivatives in heat transfer model is a reasonable step since the integer-order models are related to infinite heat flux speed. However, formal replacement of integer-order derivatives by fractional counterparts (see the model of Damor), in general, is not a correct approach since it is not based on basic laws. Hence, striking the balance at the moment, the only correctly derived time-fractional bio-heat is that of Ezzat (136). The Ferras' model uses a completely different basis for fractionalization and we will comment it further in this text when the fading memory approach will be applied. Besides, the existing studies demonstrate only fractionalization techniques and model solutions but comparisons of the simulations outputs with experimental data are still missing.

7. Time-Fractional Bio-heat Models: A Systematic Approach

7.1. Short Preamble

After the analysis of the existing fractional versions of the bio-heat equations in the preceding section 6.2 we will try to draw some principle steps how fractional operators can be implemented in existing models. From the analysis of the continuum models it became clear that the main work how to account for the vascular geometries and the effects of either large or small blood transporting vessels was already done. At this moment we accept these models as granted and the focus is on the modeling technology accounting the fact that there is a heat flux damping when it propagates through the living tissue. We already saw how this could be done by the Maxwell-Cattaneo and the Dual-phase lag approaches. Now, we will use some ideas from them and will draw new ones toward correct formulation of heat transfer models applicable to living tissues. The principle model where we will test these approaches is the Pennes model due to its simplicity and taking into account that the already existing time-fractional bio-heat models are mainly based on it.

7.2. Fractionalization: Principle Directions

From the lessons learned we may suggest two principle directions in fractionalization of heat conduction models:

(A) By application of fractional Taylor series expansions concerning three examples:

(A1) Heat flux damping by application of the fractional Taylor series expansion (first order approximation).

(A2) Both heat flux and thermal gradient damping by application of the fractional Taylor series expansions (first order approximations).

(A3) Thermal gradient damping in a viscoelastic manner by applying time-fractional derivatives.

First of all, the fractional series expansions (first order approximations) will be discussed as an analog of the methods applied in integer-order models with relaxation (in section 8.2). All these methods use time-fractional derivatives with singular (power-law kernels).

(B) The fading memory concept in two versions:

(B1) Simple fading memory (SFM) and

(B2) Extended Fading Memory (EFM).

These two fading memory versions will be investigated with various relaxation kernels, some of them different from the common power-law. Now, we will investigate these ideas step-by-step.

8. Fractionalization by Series Expansions

8.1. Fractional Taylor Series: Necessary Information

With this short note we only explain what the fractional Taylor series expansion is meaning in the process of fractionalization of model equations, thus making the further discussions more readable.

Many authors suggested general forms of power series (we already commented an example with Jumarie's derivative), precisely Taylor series with fractional derivatives. Following Hardy [129] the Riemann general Taylor series is

$\begin{array}{ll}{}_{R}f(x+x_0)={\displaystyle \sum _{k=0}^{\infty }\frac{x_0^{k+r}}{\Gamma (k+r+1)}}\left(J_s^{k+r}f\right)(x)& (157)\end{array}$

where $J_s^{k+r}$ is the Riemann-Liouville Fractional integral of order k + r (see section 6.1.1.1).

Similarly, in 1931 Watanabe [130] developed the following expression [131]

$\begin{array}{l}{}_{W}f(x)={\displaystyle \sum _{m=-k}^{n-1}\frac{{\left(x-x_0\right)}^{\alpha +m}}{\Gamma (\alpha +m+1)}}\left({\widehat{D}}_s^{\alpha +m}f\right)(x)+R_{n,k}(x),k<\alpha ,\\ s\le x_0<x& (158)\end{array}$

with a residual term

$\begin{array}{l}{R}_{n,k}(x)=\left(J_s^{\alpha +n}{\widehat{D}}_s^{\alpha +n}f\right)(x)+\frac{1}{\Gamma (-\alpha -k)}{\displaystyle {\int }_0^{x_0}{(x-t)}^{-\alpha -k-1}}\\ \left({\widehat{D}}_s^{\alpha -k-1}f\right)(t)dt& (159)\end{array}$

where ${\hat{D}}_s^{\alpha +n}$ is Riemann-Liouville of order α + n.

A Taylor series expansion in terms of Riemann-Liouville fractional derivatives was proposed by Trujillo et al. [132]

$\begin{array}{l}Tf(x)={\displaystyle \sum _{m=0}^n\frac{\Gamma (\alpha ){(x-x_0)}^{m\alpha }}{\Gamma \left[(m+1)\alpha \right]}}\left({\widehat{D}}_s^{m\alpha }f\right)x_0^{+}+R_n(x,x_0),\\ 0<\alpha \le 1,s\le x_0<x& (160)\end{array}$

with a residual terms

$\begin{array}{ll}{R}_n\left(x,x_0\right)=\frac{{\left(x-x_0\right)}^{(n+1)\alpha }}{\Gamma \left[(n+1)\alpha +1\right]}\left({\widehat{D}}_s^{(n+1)\alpha }f\right)(\zeta ),x_0\le \zeta \le x& (161)\end{array}$

At the end, we refer to the work of Odibat and Shawagfeh [133] where a generalized Taylor series expansion in terms of Caputo derivative $D_{x_0}^{m\alpha }$ of order mα was proposed

$\begin{array}{l}{}_{OS}f(x)={\displaystyle \sum _{m=0}^n\left\{\frac{1}{\Gamma (m\alpha +1)}{D}_{x_0}^{m\alpha }\left[f(x_0\right]\right\}}){\left(x-x_0\right)}^{m\alpha }+R_n^{\alpha }(x),\\ 0<\alpha \le 1,x_0<x\le b\\ R_n^{\alpha }(x)=\frac{\left(D_{x_0}^{(m+1)}f\right)(\xi )}{\Gamma [(m+1)\alpha +1]}{\left(x-x_0\right)}^{\left(m+1\right)\alpha }\\ D_{x_0}^{m\alpha }=D_{x_0}^{\alpha }·D_{x_0}^{\alpha }·D_{x_0}^{\alpha }···\\ D_a^{k\alpha }f(x)\in C(a,b],k=0,1,\mathrm{...},n+1,0<\alpha \le 1& (162)\end{array}$

more informatively, for the purposes of the following sections, we have

$\begin{array}{l}{}_{OS}f(x)=f(x)+\left[\frac{\Delta x^{\alpha }}{\Gamma (1+\alpha )}\right]\frac{{\partial }^{\alpha }f}{\partial x^{\alpha }}+\left[\frac{\Delta x^{2\alpha }}{\Gamma (1+2\alpha )}\right]\\ \frac{{\partial }^{\alpha }f}{\partial x^{\alpha }}\left(\frac{{\partial }^{\alpha }f}{\partial x^{\alpha }}f(x)\right),\Delta x=x-x_0,0<\alpha \le 1& (163)\end{array}$

In terms of fractional relaxation times, we may write

$\begin{array}{l}{}_{OS}f(x)=f(x)+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }f(x)}{\partial x^{\alpha }}++{\lambda }^{2\alpha }\frac{{\partial }^{\alpha }f(x)}{\partial x^{\alpha }}\left(\frac{{\partial }^{\alpha }f}{\partial x^{\alpha }}f(x)\right),\\ 0<\alpha \le 1& (164)\end{array}$

where

$\begin{array}{ll}{\lambda }^{\alpha m}=\frac{{\left(\Delta x\right)}^{\alpha m}}{\Gamma (1+m\alpha )},0\le m<\infty & (165)\end{array}$

and λ^αm has a dimension ${(di{m}_x)}^{m\alpha }$, where dim_x is the dimension of x (time or length, depending on the physics of the modeled process).

Now, we will stop with definitions of fractional Taylor series since this is out of the scope of this work. We will concentrate the attention on the first-order approximations only, that is for m = 1. For more details and information we refer to El-Ajou et al. [131]. Hereafter in our discussion we may assume that the fractional Taylor series are based on Caputo derivative, if some other type is not especially specified. Recall, that in the above equations, x was used optionally and this is not the space coordinate used in the model. The relationships remain the same if we will use the time t instead x.

8.2. Heat Flux Damping by Fractional Taylor Series Expansion (Fractional SPL-Q)

This approach means application of the constitutive flux-gradient relationship in the form

$\begin{array}{ll}q\left(x,t\right)=q_0+{\lambda }_q\frac{{\partial }^{\alpha }q\left(x,t\right)}{\partial t^{\alpha }}& (166)\end{array}$

without specification what type of fractional derivative is used (from the previous section it is easy to see that using only 2 terms for k = 0 and k = 1 we get similar results). When the Caputo derivative is used, then ${\lambda }_q^{\alpha k}$ is defined by Equation (165).

Actually this a fractional version of the Maxwell-Cattaneo approach (see section 5) (which reduces to (79) for α = 1), an approach already used by Ezzat et al. [122, 123]. We already analyzed this model (see section 6.2.1) and therefore it will be more interesting to skip repetitions and explain the other ideas for fractionalizations. The case with α ≠ 1 is termed here Fractional SPL with respect to the heat flux, SPL-Q.

It is noteworthy to mention that this is the simplest way to achieve model fractionalization (excluding the formal replacement of the derivatives by fractional ones, which in many case results in unphysical models). However, there are no restrictions to apply high order approximation by of the heat flux as

$\begin{array}{ll}q\left(x,t\right)=q_0+{\displaystyle \sum _{k=1}^K{\lambda }_q^k\frac{{\partial }^{{\alpha }_k}q\left(x,t\right)}{\partial t^{{\alpha }_k}}}& (167)\end{array}$

Obviously, this will lead to more complex models with fractional derivatives of different orders. In this context, their applicability could be proved only by testing models to real data and answer sincerely to the question: What this flux damping approach will model and how it will improve the existing model to which it will be applied? Hence, we have an idea, and a guiding line how to perform it, but final formulations of the consequent models and solutions, and generally, answers to the emerging questions are not available yet. At this moment we have to take into account the comments at the end of section 5.1.1, precisely that λ_q should be mathematically small since if we like the model to be physically sound.

8.3. Thermal Gradient Damping by Fractional Taylor Series Expansions (Fractional SPL-T)

Alternatively to the approach used in section 8.2, we apply the fractional Taylor expansion to the temperature gradient only

$\begin{array}{ll}\frac{\partial T(x,t)}{\partial x}=\frac{\partial T(x,t)}{\partial x}+{\lambda }_T\frac{{\partial }^{\beta }}{\partial t^{\beta }}\left(\frac{\partial T(x,t)}{\partial x}\right),0<\beta \le 1& (168)\end{array}$

where

$\begin{array}{ll}{\lambda }_T=\frac{{\tau }^{\beta }}{\Gamma (1+\beta )}& (169)\end{array}$

with a dimension s^β. Then, the fractionalized heat flux is

$\begin{array}{ll}q(x,t,\beta )=k\left(1+{\lambda }_T\frac{{\partial }^{\beta }}{\partial t^{\beta }}\right)\left(\frac{\partial T(x,t)}{\partial x}\right)& (170)\end{array}$

or

$\begin{array}{ll}q(x,t,\beta )=k\frac{\partial T(x,t)}{\partial x}+k_T^{\beta }\frac{{\partial }^{\beta }}{\partial t^{\beta }}\left(\frac{\partial T(x,t)}{\partial x}\right)& (171)\end{array}$

The transport coefficient $k_T^{\beta }=k{\lambda }_T$ with a dimension (W/mK)s^β is a damping heat conduction (some times termed as elastic heat conduction) which has no effect when the memory kernel of ∂^β/t^β goes to zero. For β = 1 in (170) and (171), we get the right-hand side of the integer-order DPL model. By analogy, we term this version of the fractional SPL, FSPL-T (fractional SPL with respect to the temperature gradient).

The approach demonstrated here can be applied by application of Taylor series expansions to either the time t, as it as done here, which corresponds to the assumption of a homogenous heat conductor, or with respect to the space coordinate x if the heat conduction is assumed as non-homogeneous in space (see such examples in [134, 135]).

8.4. Fractional Dual-Phase Lag Approach

This approach is a fractional version of the Dual-Phase Lag concept [83, 84, 102] and by analogy we will term it Fractional Dual-Phase Lag (FDPL) concept. The time-fractional version of the constitutive equation relating the heat flux and the temperature gradient is

$\begin{array}{ll}q\left(x,t\right)=q_0+{\lambda }_q\frac{{\partial }^{\alpha }q\left(x,t\right)}{\partial t^{\alpha }}=k\left(1+{\lambda }_T\frac{{\partial }^{\beta }}{\partial t^{\beta }}\right)\left(\frac{\partial T(x,t)}{\partial x}\right)& (172)\end{array}$

Obviously, λ_q ≠ λ_T. In order to preserve the causality of the heat flux with respect to the temperature, we need λ_q > λ_T and λ_q − λ_T > 0. Consequently, the causality requirement mentioned above needs α < β, that is the flux should have stronger memory than the temperature gradient. There is no mandatory condition requiring α = β but the condition λ_q − λ_T > 0 remains. For α = β = 1 we get the integer order DPL model (section 5.2), and obviously we should have λ_q > λ_T. Further, the approach expressed by (166) could be applied to FDPL, but we will skip this point since a serious physical background to do that is missing, albeit the formal application of high order approximations by fractional Taylor series is a valid mathematical approach.

8.5. Thermal Gradient Damping in a Viscoelastic Manner

The idea for this damping of the heat flux (by high order time derivatives) comes from the viscoelastic models of non-Newtonian general second grade fluids where the relationship between the shear stress τ_shear(t) and strain ε_strain(t) is presented as [93–95, 136]

$\begin{array}{l}{\tau }_{shear}\left(t\right)=\mu {\epsilon }_{strain}\left(t\right)+{\lambda }_e{D}_t^{\gamma }\left[{\epsilon }_{strain}\left(t\right)\right],\\ \tau \left(t\right)=\mu \epsilon \left(t\right)+E\frac{d\epsilon \left(t\right)}{dt}& (173)\end{array}$

Here μ is a fluid dynamic viscosity, while λ_e is the first normal stress modulus (analog to the heat conduction relaxation time in the expansions by Taylor series commented above [94, 136]

$\begin{array}{ll}\rho \frac{\partial u}{\partial t}=\mu \frac{{\partial }^{2}u}{\partial y^2}+{\alpha }_1{D}_t^{\gamma }\frac{{\partial }^{2}u}{\partial y^2}& (174)\end{array}$

Alternatively, the second version of (173) leads to the governing integer-order equation of the unidirectional flow of a second grade fluid [93, 94, 136, 137], namely

$\begin{array}{ll}\frac{\partial u}{\partial t}=\nu \frac{{\partial }^{2}u}{\partial y^2}+{\beta }_1\frac{{\partial }^{3}u}{\partial t\partial y^2}& (175)\end{array}$

where ν is fluid kinematic viscosity with dimension m²/s (analog of the thermal diffusivity a).

The last term in (175) is the same as in the integer order DPL model [see (99) or (100)].

To be precise, and for the sake of clarity of the following expression, we mention that from the basic relationship between the velocity u(x, t) and shear stress τ_fric (the flux of momentum in hydrodynamics

$\begin{array}{ll}\frac{\partial u}{\partial t}=-\nu \frac{\partial {\tau }_{fric}}{\partial x}& (176)\end{array}$

and the integer-order relationship

$\begin{array}{ll}{\tau }_{fric}=-\left[\nu +{\beta }_1\frac{\partial }{\partial t}\right]\frac{\partial u}{\partial x}& (177)\end{array}$

leading to (175) we actually have the SPL-T analog in the fluid flow.

The fractional version corresponding to (174) is

$\begin{array}{ll}{\tau }_{fric}=-\left[\nu +{\beta }_1\frac{{\partial }^{\gamma }}{\partial t^{\gamma }}\right]\frac{\partial u}{\partial x}& (178)\end{array}$

Therefore, using analogy between the fluid velocity field and the temperature field distribution in the body it is possible to formulate a constitutive relationship between the heat flux and temperature gradient as

$\begin{array}{l}q\left(x,t\right)=-k\frac{\partial T\left(x,t\right)}{\partial x}-k{\lambda }_{eT}\frac{{\partial }^{\beta }}{\partial t}\left(\frac{\partial T\left(x,t\right)}{\partial x}\right)\\ =-k\left(1+{\lambda }_{eT}\frac{{\partial }^{\beta }}{\partial t^{\beta }}\right)\frac{\partial T\left(x,t\right)}{\partial x}& (179)\end{array}$

Equation (179) is only the right-hand side of the Fractional Dual-Phase Lag model formulated above. For β = 1 we get the integer order DPL model (actually, only its half part related to the damping of the temperature gradient).

Then, with the energy balance relationship $\rho c_p\frac{\partial }{\partial t}T(x,t)=-\frac{\partial q(x,t)}{\partial x}$ (without convection and volumetric heat sources: for the sake of simplicity in presentation of this fractionalization technique) we get

Fractional model

$\begin{array}{ll}\frac{\partial T(x,t)}{\partial t}=k\frac{{\partial }^{2}T(x,t)}{\partial x^2}+{\lambda }_{eT}\frac{{\partial }^{\beta }}{\partial t^{\beta }}\left(k\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right),0<\beta <1& (180)\end{array}$

where λ_eT has a dimension s^β and therefore the product kλ_eT could be considered as a fractional (in time) heat conductivity k^β with a dimension (W/mK)s^β.

Integer-order model for β = 1

$\begin{array}{ll}\frac{\partial T(x,t)}{\partial t}=k\frac{{\partial }^{2}T(x,t)}{\partial x^2}+{\lambda }_{eT}\frac{\partial }{\partial t}\left(k\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right)& (181)\end{array}$

Further, applying this approach to the Pennes model we get

$\begin{array}{ll}& {\rho }_t{c}_t\frac{\partial T_t(x,t)}{\partial t}=k\frac{{\partial }^2{T}_t(x,t)}{\partial x^2}+{\lambda }_{eT}\frac{{\partial }^{\beta }}{\partial t^{\beta }}\left(k\frac{{\partial }^2{T}_t(x,t)}{\partial x^2}\right)\\ +& {\omega }_b{\rho }_b{c}_b(T_{a0}-T_t)+q_m& (182)\end{array}$

Actually, the fractionalization applied here is equivalent to the one-side approach (fractionalization only of the temperature gradient) demonstrated in section 8.3.

Note 5: We like to stress the attention on the last two versions of fractionalization where the damping effect is applied to the temperature gradient only (sections 8.3 and 8.5). At a glance, it would be a violation of the causalty between the heat flux and the temperature gradient. However, we have to mention that the Taylor series expansion and fractionalization of the viscoelastic manner are more philosophically related to the fading memory concept developed in the next section 9 and there is no violation of causality mentioned above.

8.6. Comments on Fractionalization by Taylor Series Expansions

Fractionalization of integer-order models by application of fractional Taylor series expansions (first-order approximation) is a feasible technique resulting in simple time-fractional equations. The techniques demonstrated in this section dealt only with models assuming homogenous tissues as heat conductors. Logically, the fractional Taylor series technique could be extended toward space-fractional models but this needs serious physical motivation and information about the non-homogeneity of the tissue of interest. Regarding the existing time fractional bio-heat models, only the Ezzat' equation (136) [through (128)] follows this approach correctly, despite the fact that the obsolete Jummarie fractional derivative is used. However, there are no obstacles the same approach to be applied with Taylor series based on the Caputo derivative (this will change only the pre-factors of the time-fractional terms). It might be expected that the approaches demonstrated in this section would be successfully applied for creation of new bio-heat models.

9. The Fading Memory Concept to Heat Conduction

9.1. Basic Principle of Fading Memory in Heat Conduction

The fading memory concept relating the flux j(x, t) to its gradient, for simple materials [128, 138–142], is modeled by the following integro-differential equation

$\begin{array}{ll}j(x,t)=-D_0\nabla C(x,t)-D^{\prime }{\displaystyle \underset{-\infty }{\overset{t}{\int }}}R(t-\tau )\nabla C(x,\tau )d\tau & (183)\end{array}$

as a manifestation of the Boltzmann linear superposition functional expressing the flux history [143] through the function of influence (memory kernel) R(t, τ). In (183) D₀ and D′ are transport coefficients (diffusivities).

The appropriate history value problem is related to the hereditary integral in (183) [144]

$\begin{array}{ll}d(t)={\displaystyle \underset{-\infty }{\overset{0}{\int }}}R(t-\tau )\nabla C(x,\tau )d\tau & (184)\end{array}$

allowing to give a function C(x, t) on (−∞ < t ≤ 0)

From (183) and (184) it follows that

$\begin{array}{lll}\nabla ·j(x,t)& =& -D_0\Delta C(x,t)-D^{\prime }{\displaystyle \underset{-\infty }{\overset{t}{\int }}}R(t-\tau )\Delta C(x,\tau )d\tau \\ & +& \nabla ·d(t)& (185)\end{array}$

Since C(x, t) is a causal function (vanishing for t < 0) and considered only for 0 < t < ∞ it follows that d(t) = ∇ · d(t) = 0, thus (185) can be rewritten as

$\begin{array}{ll}\nabla ·j(x,t)=-D_0\Delta C(x,t)-D^{\prime }{\displaystyle \underset{0}{\overset{t}{\int }}}R(t-\tau )\Delta C(x,\tau )d\tau & (186)\end{array}$

Following Gurtin [128], Gurtin and Pipkin [128] the heat flux can be expressed in a shorter version by the hereditary integral only

$\begin{array}{ll}\nabla ·j(x,t)=-D^{\prime }{\displaystyle \underset{0}{\overset{t}{\int }}}R(t-\tau )\Delta C(x,\tau )d\tau & (187)\end{array}$

which comes from (186) for D₀ = 0. Hence we have two forms of the fading memory concept:

(B1) Simple fading memory concept expressed by (187).

(B2) Extended fading memory concept expressed by (186).

Let us now see how these two versions of the fading memory work in the heat conduction problem.

9.2. Simple Fading Memory Concept and Effect of the Kernel

As it was several times commented above the Fourier law trough the energy balance results in the simple parabolic heat conduction equation by the relationship

$\begin{array}{ll}\rho c_p\frac{\partial T}{\partial t}=-\frac{\partial q}{\partial x},q(x,t)=-k\frac{\partial T(x,t)}{\partial x}\Rightarrow \rho c_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial x^2}& (188)\end{array}$

with infinite speed of the heat flux.

A damping function related to a finite speed of heat diffusion in rigid conductors was conceived by Cattaneo [145] by a generalization of the Fourier law through a linear superposition of the heat flux and its time history (with a time-delay s by a hereditary (memory) integral).

$\begin{array}{ll}q(x,t)=-{\displaystyle {\int }_0^t}R(x,t-s)\frac{\partial }{\partial x}T(x,t-s)ds& (189)\end{array}$

Now, the principle question is: What function R(x, t) should serve as adequate memory kernel? It is noteworthy to mention that the type of the relaxation function R(x, t) strongly depends on the reaction (response) of the material where the heat transport takes place. It was demonstrated in Hristov [116, 120], in the field of the linear viscoelasticity of materials, the relaxation response defines the type of the memory function in the hereditary integral. Hence, there is a freedom in choosing R(x, t), but this cannot be done voluntary since the physics should be taken into account. In the sequel we will see how the final form of the heat conduction equation depends on choice of R(x, t).

9.2.1. Exponential Kernel

For homogeneous rigid heat conductors the memory function R(x, t) is space-independent and can be represented, for instance, by the Jeffrey kernel [96, 97] as

$\begin{array}{ll}R(t,s)=exp(-(t-s)/\tau )& (190)\end{array}$

where the relaxation time τ is finite, i.e., τ = const.

Then, the energy balance yields the Cattaneo equation

$\begin{array}{ll}\frac{\partial T(x,t)}{\partial t}=-\frac{a_2}{\tau }{\displaystyle \underset{0}{\overset{t}{\int }}}exp\left[-\left(\frac{t-s}{\tau }\right)\right]\frac{\partial T(x,s)}{\partial x}ds,a_2=\frac{k_2}{\rho c_p}& (191)\end{array}$

For t → ∞ the right-hand side is zero and no physical meaning comes from this. Here, the thermal diffusivity a₂ is elastic thermal diffusivity relevant only to the fading memory term (see below the comments about (194). For τ → 0 the limit of the Cattaneo equation is the Fourier law.

The first order approximation of the heat flux, in τ yields [146] is

$\begin{array}{ll}q(x,t+\tau )=-k_1\frac{\partial T(x,t)}{\partial x},q(x,t+\tau )\approx q(x,t)+\tau \frac{\partial q(x,t)}{\partial x}& (192)\end{array}$

a technique already discussed in section 5.1

This leads to a first order differential equation [146]

$\begin{array}{ll}\frac{1}{\tau }q(x,t)+\frac{\partial q(x,t)}{\partial t}=-\frac{k_1}{\tau }\frac{\partial T(x,t)}{\partial x}& (193)\end{array}$

The integration of (193) with help of (192) results in the Cattaneo equation, which is the simplest giving rise to finite speed of flux propagation. Further, continuing the integration and applying the energy balance we will get the classical telegraphic hyperbolic equation. However, now we look at results from a different viewpoint trying to interpret the constitutive equations involving hereditary integral from the position of fractional calculus.

9.2.2. Relaxation Function of Joseph and Preziosi

Considering the heat-wave nature of the conduction energy transfer Joseph and Preziosi [97] have considered a modified relaxation function replacing the exponential kernel in (189) by

$\begin{array}{ll}{R}_{JP}=k_1\delta (s)+(k_2/\tau )exp(-s/\tau )& (194)\end{array}$

, where δ(s) is Dirac delta function, while k₁ and k₂ are the effective thermal conductivity and the elastic conductivity, respectively. This constitutive equation defines the relaxation function in a way that matches exactly the philosophy of fading memory concept: the first term of R_JP models the instantaneous (memory-less) reaction of the heat conductor, that is this corresponds to the Fourier law, while the second term corresponds to the flux damping effect and finite speed of propagation.

In this case the Fourier law leads to a flux defined as [97, 146]

$\begin{array}{ll}q(x,t)=-k_1\frac{\partial T(x,t)}{\partial x}-\frac{k_2}{\tau }{\displaystyle \underset{0}{\overset{t}{\int }}}{e}^{-\left(\frac{t-s}{\tau }\right)}\frac{\partial T(x,s)}{\partial x}ds& (195)\end{array}$

Consequently, the conservation equation of the internal energy in the heat conductor [128] results in the Jeffrey-type integro-differential equation [146]

$\begin{array}{l}\hfill \frac{\partial T(x,t)}{\partial t}=a_1\frac{{\partial }^{2}T(x,t)}{\partial x^2}+\frac{a_2}{\tau }{\displaystyle \underset{0}{\overset{t}{\int }}}{e}^{-\left(\frac{t-s}{\tau }\right)}\frac{{\partial }^{2}T(x,s)}{\partial x^2}ds,\hfill \\ \hfill a_1=k_1/\rho c_p,a_2=k_2/\rho c_p\hfill & (196)\end{array}$

The second term in (196) goes to zero for large times, since the memory function fades and only the first term remains; despite this the time for existence of the second term is related to the finite heat flux speed (i.e., the flux is damped in time). Equation (196) corresponds to the extended fading memory concept with a simple exponential kernel, as it will be demonstrated in the sequel (see section 9.3.1).

Equation (196) can be expressed as [96] (see also the analysis in [119])

$\begin{array}{l}\frac{\partial T(x,t)}{\partial t}=a_1\frac{{\partial }^{2}T(x,t)}{\partial x^2}+a_2{(1-\alpha )}_{CF}{D}_t^{\alpha }\frac{{\partial }^{2}T(x,t)}{\partial x^2},\\ t>0,0<\alpha \le 1& (197)\end{array}$

with a damping term controlled by the Caputo-Fabrizio operator (derivative) [118]

$\begin{array}{ll}{{}_{CF}D}_t^{\alpha }f\left(t\right)=\frac{M\left(\alpha \right)}{1-\alpha }{\displaystyle \underset{0}{\overset{t}{\int }}\exp \left[-\frac{\alpha \left(t-s\right)}{1-\alpha }\right]\frac{df\left(t\right)}{dt}ds}& (198)\end{array}$

where M(α) is a normalization function such that M(0) = M(1) = 1.

9.2.3. Simple Power-Law Kernel

The exponential memory function (190) is a bonded kernel since for t → 0⁺ it approach unity. However, the power-law memory kernel t^α−1/Γ(α) (with 0 < α <1) is unbounded at t → 0⁺ and time-scale invariant (the exponential function is not time scale invariant).

Now, constructing the shorter equation of the fading memory concept (187) with a power-law function we get

$\begin{array}{ll}q(x,t)=-K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\frac{\partial }{\partial x}T(x,\tau )d\tau & (199)\end{array}$

We can see that the right-hand side in (199) is the Riemann-Liouville fractional integral [117]. Hence, the heat flux gradient dq/dx is

$\begin{array}{ll}\frac{\partial q(x,t)}{\partial x}=-K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\frac{{\partial }^{2}T(x,\tau )}{\partial x^2}d\tau & (200)\end{array}$

Now, applying the energy balance equation we get

$\begin{array}{ll}\rho c_p\frac{\partial T(x,t)}{\partial t}={\displaystyle \underset{dampedheatconduction(elasticterm)}{{\displaystyle \underbrace{K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\frac{{\partial }^{2}T(x,\tau )}{\partial x^2}d\tau }}}}& (201)\end{array}$

At this moment we may perform a small thought experiment. Let us suggest that the heat conduction is modeled by the Continuous Time Random Walk (CTRW) concept. Then, the right-hand side of (199) corresponds to this model and we get

$\begin{array}{ll}\frac{d}{dx}\left[q_{elastic}(x,t)\right]=-K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\frac{{\partial }^{2}T(x,\tau )}{\partial x^2}d\tau & (202)\end{array}$

Differentiation of (202) with respect to the time yields

$\begin{array}{l}\frac{\partial }{\partial t}\left[\frac{\partial }{\partial x}(q_{elastic}(x,t))\right]=-K_{\alpha }\frac{1}{\Gamma (1-\alpha )}\\ \frac{d}{dt}{\displaystyle {\int }_0^t}\frac{1}{{(t-\tau )}^{\alpha }}\frac{{\partial }^{2}T(x,\tau )}{\partial x^2}d\tau & (203)\end{array}$

The right-hand side of (203) is the core of the constitutive equation about the heat flux relaxation in the fractional model of Ferras et al. [127] (see Equation 154). Further, the suggestion of Ferras et al. [127] simply means that applying the heat balance we have to write

$\begin{array}{ll}\frac{\partial }{\partial t}\left[\rho c_p\frac{\partial T(x,t)}{\partial t}\right]=-\frac{\partial }{\partial t}\left(\frac{\partial }{\partial x}{q}_{elastic}\right)=-\frac{\partial }{\partial x}\left(\frac{\partial }{\partial t}{q}_{elastic}\right)& (204)\end{array}$

where $\frac{\partial q}{\partial t}$ and $\frac{\partial q}{\partial x}$ should be taken from the elastic part of Equation 201).

Therefore, the final result of the Ferras' constitutive equation should be

$\begin{array}{ll}\rho c_p\frac{{\partial }^{2}T(x,t)}{\partial t^2}=K_{\alpha }\left({{}^{RL}D}_t^{\alpha }\left[\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right]\right)& (205)\end{array}$

but not equation (154).

In this context, we can see, as it was mentioned earlier in section 6.2.3, that the quest to find something that is well-known and easily solvable may lead to use of wrong premises (the Ferras's et al. constitutive relationship heat flux-thermal gradient) violating basic thermodynamic principles such as the fading memory concept or naively applying it. We are aware about the commonly mentioned and published solutions of the time-fractional heat-wave equation [117] (the final result of Ferras et al. [127]) but applying it we need serious physical background. In the specific case of bio-heat models such physical reasons are still missing (or not published yet).

9.2.4. Extended Power-Law Kernel

Here, we suggest an extended power-law (EPL) memory kernel, by analogy of the extended exponential relaxation function (189) of Joseph and Preziosi [97], namely

$\begin{array}{ll}{R}_{EPL}(t,s)=k_1\delta (s)+k_2\frac{t^{\alpha -1}}{\Gamma (\alpha )}& (206)\end{array}$

With

$\begin{array}{ll}q=-k{\displaystyle {\int }_0^t}{R}_{EPL}(t,s)\frac{\partial T(x,s)}{\partial x}ds& (207)\end{array}$

we get

$\begin{array}{ll}q=-k_1\frac{\partial T(x,t)}{\partial x}-k_2\frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t}\frac{1}{{(t-s)}^{1-\alpha }}\frac{\partial T(x,s)}{\partial x}ds& (208)\end{array}$

That is, the second term in (208) is the Riemann-Liouville fractional integral. The result (208) coincides with the one obtained by the extended fading memory concept and simple power-law kernel (see the next section 9.3.2).

Then, the heat conduction equation becomes

$\begin{array}{l}\frac{\partial T(x,t)}{\partial t}=a_1\frac{{\partial }^{2}T(x,t)}{\partial x^2}-a_2\frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t}\frac{1}{{(t-s)}^{1-\alpha }}\frac{{\partial }^{2}T(x,s)}{\partial x^2}ds,\\ a_1=\frac{k_1}{\rho c_p},a_2=\frac{k_2}{\rho c_p}& (209)\end{array}$

For large times, when the memory kernel goes to zero (209), reduces to the Fourier equation. The dimension of a₁ is m²/s, while k₁ has dimension [W/m · K], that is the same as in the Fourier law. The dimension of k₂ (and a₂) needs a special attention and this will be discussed in section 9.4.

9.3. Extended Fading Memory Concept and Effect of the Kernel

9.3.1. Exponential Kernel

With the general construction (186) and exponential kernel (190) we get

$\begin{array}{l}\frac{\partial T(x,t)}{\partial t}=a_1\frac{{\partial }^{2}T(x,t)}{\partial x^2}-\frac{a_2}{\tau }{\displaystyle \underset{0}{\overset{t}{\int }}}{e}^{-(\frac{t-s}{\tau })}\frac{{\partial }^{2}T(x,s)}{\partial x^2}ds,\\ a_1=\frac{k_1}{\rho c_p},a_2=\frac{k_2}{\rho c_p}& (210)\end{array}$

The result (210) matches (196) obtained with the simple fading memory concept and the extended exponential kernel (194).

9.3.2. Power-Law Kernel

The exponential memory function (190) is a bonded kernel since for t → 0⁺ it approach unity. However, the power-law memory kernel t^α−1/Γ(α) (with 0 < α < 1) is unbounded at t → 0⁺ and time-scale invariant (the exponential function is not time scale invariant).

Constructing the basic equation of the fading memory concept (183) with a power-law function we get

$\begin{array}{ll}q(x,t)=-K_1\nabla C(x,t)-K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\nabla T(x,\tau )d\tau & (211)\end{array}$

The second term in (211) is the Riemann-Liouville fractional integral [117] and the heat flux gradient dq/dx is

$\begin{array}{ll}\frac{dq(x,t)}{dx}=-K_1\frac{{\partial }^{2}T(x,t)}{\partial x^2}-K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\frac{{\partial }^{2}T(x,\tau )}{\partial x^2}d\tau & (212)\end{array}$

Then, applying the energy balance equation we get

$\begin{array}{ll}\begin{array}{c}\rho c_p\frac{\partial T(x,t)}{\partial t}={\displaystyle \underset{instantaneousconduction}{{\displaystyle \underbrace{K_1\frac{{\partial }^{2}T(x,t)}{\partial x^2}}}}}+\\ +{\displaystyle \underset{dampedheatconduction(elasticterm)}{{\displaystyle \underbrace{K_{\alpha }\frac{1}{\Gamma (\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{1}{{(t-\tau )}^{1-\alpha }}\frac{{\partial }^{2}T(x,\tau )}{\partial x^2}d\tau }}}}\end{array}& (213)\end{array}$

For large times, when the memory kernel vanishes, equation (213) reduces to the Fourier model.

9.4. On the Dimensionality in Constitutive Equations With Fading Memory Operators

The dimensions of the coefficients with temporal relaxations in both the fractional Taylor series and fading memory operators were mentioned several times, but now we like to stress the attention especially on this point. When the flux damping is modeled by either SPL or DPL in integer-order or fractional manner we have

$\begin{array}{ll}q=q_0+{\lambda }_q\frac{{\partial }^{\alpha }q(x,t)}{\partial t^{\alpha }},0<\alpha \le 1& (214)\end{array}$

the fading memory is incorporated in the construction of the fractional derivative. In this case the dimensionality of λ_q is [s^α], thus the second term in (214) has a dimension of heat flux, that is (W/m²).

Similarly, when the temperature gradient is damped by SPL or DPL, we have

$\begin{array}{ll}k\frac{\partial T}{\partial x}=k\frac{\partial T(x,t)}{\partial x}+k^{\alpha }\frac{{\partial }^{\alpha }T(x,t)}{\partial t^{\alpha }},0<\alpha \le 1& (215)\end{array}$

where k^α has dimension of (W/m · K)s^α, that is $k^{\alpha }=k·{\lambda }_T$.

Alternatively, we may write

$\begin{array}{ll}K=k+k^{\alpha }=k\left(1+\frac{k^{\alpha }}{k}\right)=k(1+{\lambda }_T)& (216)\end{array}$

The term (1 + λ_T) is valid only during the short time relaxation defined by λ_T and has a temporal dimension defined by λ_T. Therefore, for times coinciding with the relaxation (damping) period (0 < t ≤ λ_T) the dimension of K is (W/m · K)s^α in contrast to the long-time case (t > λ_T) when the dimension is (W/m · K). This is strongly related to the fact that the two terms in the right-side hand of (214) (the same in (215) are transient but:

(1) The local first term (with coefficient to k) is transient due the change in the boundary condition at t = 0, while

(2) The second one (with the memory) (with a lower terminal of the integral at t = 0 due to the causalty condition) is also transient but it is non-local and is related to the time history of the flux.

If the medium relaxation time is zero, that is an instantaneous reaction of the disturbed body is observed, this means that, for instance, with λ_T = λ_q we get K = k (see Equation 216). In this context, we may say something related to the heat conduction problem and the dimensionality of the coefficients involved.

In the classical heat conduction, since the time of Fourier, the dimension of the heat conductivity is time-independent and when a transient heat conductivity appears in constitutive equation with memory, this may cause ambiguities and cast doubts. The simple explanations are:

(a) The traditional time-independent conductivity occurs in problems when the transient processes are larger that the time duration of the action at the boundary (it is well-known that the Fourier heat conduction equation works only for large time, and fails for short times, as well as is thermodynamical inconsistent). That is, there is no need to know what is the relationship between the time of the action at the boundary, the total transient time of the process and the specific relaxation time (the accommodation time) of the materials, i.e., now knowledge about the process history is needed (memoryless process). Recall also the comments just after Equations (78) and (79).

(b) However, when the time duration of the action at the boundary is comparable with the material thermal relaxation time, then the damping of the heat flux should be taken into account and this is well described by a hereditary integral (a technique know since the time of Boltzmann [143]).

(c) Better understanding of the thoughts above, could be attained if we consider the heat transfer by conduction when the heat flux is related to the temperature gradient by the thermal diffusivity a_T = k/ρc_p. Then, we may write the well-known time-fractional heat-wave equation

$\begin{array}{ll}\frac{{\partial }^{\alpha }T}{\partial t^{\alpha }}=a_T^{\alpha }\frac{{\partial }^{2}T}{\partial x^2}& (217)\end{array}$

where the dimensionality of $a_T^{\alpha }$ is m²/s^α, and the time fractional derivative is either from Riemann-Liouville or Caputo type.

From (217) it follows that of the heat balance equation is ∂T/∂t = −∂q/∂x then the heat flux should be related to the thermal gradient by the relation $q=-a_T^{\alpha }(\partial T/\partial x)$.

Hence, there is no contradiction with the classical heat condition, but only development of the existing knowledge. This is a simple example involving a time-fractional derivative, as a more familiar operator to the people working in fractional modeling, but the same philosophy remains when fractional integral of Riemann-Liouville type appears in application of the fading memory principle. Obviously, the physics of the heat flux relaxation should be correctly approximated by the type of the memory kernel and the definition of the temporal dimension of coefficient in front of the fractional integral.

This short section only marks the problem and comments issues relevant to the discussed topic since for many readers the appearance of time-fractional-dependent heat conductivity (or thermal diffusivity), for instance, will be quite unusual and heretic, or at least erroneous. For more studies and examples involving fractional (non-local problems) and fractal (local problems) dimensions we refer to Westerlund [147], Muslih and Baleanu [148], Ebaid et al. [149], Gomez-Aguilar et al. [150, 151], Dokoumetzidis and Macheras [152], and Liu [153].

10. Fading Memory Concept to Bio-heat Equations: Examples

Now, at the end of this long analysis we present two examples how the fading memory concept works in the formulation of bio-heat equations with memories.

10.1. Pennes' Equation With a Fading Memory

10.1.1. Extended Fading Memory Concept

Let us consider equation (22) in section 3.1. Assuming that the non-local perfusion term and the metabolic heat are not relaxing, and applying the fading memory principle we get

$\begin{array}{l}{\rho }_F{c}_F\frac{\partial T_F(x,t)}{\partial t}=k_t\frac{{\partial }^2{T}_F(x,t)}{\partial x^2}+K_{\alpha }{\displaystyle {\int }_0^t}R(t,\tau )\frac{{\partial }^2{T}_F(x,\tau )}{\partial x^2}d\tau \\ -h_F\frac{P_F}{A_F}\left[T_F(x,t)-T_{\infty }\right]& (218)\end{array}$

Now, applying different memory kernels R(t, τ) we may obtain equations in different forms with different fractional operators, but recall the choice of R(t, t) should be physically motivated. The application of the fading memory principle to the Pennes equation results in a nice-looking model, but it still retains its main deficiency: the already commented problems with the perfusion term.

10.1.2. Simple Fading Memory Concept

In this simple case there is no instantaneous heat flux component, that is first therm in (218) is zero. Then, we have

$\begin{array}{l}\frac{\partial T_F(x,t)}{\partial t}=\frac{K_{\alpha }}{{\rho }_F{c}_F}{\displaystyle {\int }_0^t}R(t,\tau )\frac{{\partial }^2{T}_F(x,\tau )}{\partial x^2}d\tau -\frac{h_F}{{\rho }_F{c}_F}\frac{P_F}{A_F}\\ \left[T_F(x,t)-T_{\infty }\right],a_T^{\alpha }=\frac{K_{\alpha }}{{\rho }_F{c}_F}& (219)\end{array}$

where $a_T^{\alpha }$ is the thermal diffusivity with a dimension m²/s^α. Hence, only the damping term remains. The relaxation function R(t, τ) strongly depends on the physical concept assumed in the heat conduction modeling, precisely what should be the heat conduction transfer mechanism. With the power-law kernel, for instance, as in (213), accepting the random walk concept, we get Riemann-Liouville integral and consequently the time-fractional heat-wave equation with a sink. It is easy to check different versions of (219) by change in the memory kernel, but this should be seriously motivated in regards what type of tissue will be considered and how the relaxation function is related to its real thermal behavior.

10.2. Wulff's Model With a Fading Memory

10.2.1. Extended Fading Memory Concept

Let us start with Wulff formulation of the heat flux (33) [21, 29, 32]. Applying the fading memory principle we have

$\begin{array}{l}q={\displaystyle \underset{instantaneousheattransfer}{{\displaystyle \underbrace{\left[-k_t\nabla T_t+{\rho }_b{h}_h{v}_h\right]}}}}\\ -{\displaystyle \underset{elasticheattransport}{{\displaystyle \underbrace{K_{\alpha }{\displaystyle {\int }_0^t}R(t,\tau )\nabla T_t(x,\tau )d\tau +{\displaystyle {\int }_0^t}R(t,\tau ){\rho }_b{h}_h{v}_{h}d\tau }}}}& (220)\end{array}$

Thus, the energy balance ${\rho }_b{c}_b\frac{\partial T_t}{\partial t}=-\nabla ·q$ reads

$\begin{array}{l}{\rho }_b{c}_b\frac{\partial T_t}{\partial t}=\left[k_t\frac{{\partial }^2{T}_t}{\partial x^2}-{\rho }_b{h}_h\frac{\partial v_h}{\partial x}\right]+K_{\alpha }{\displaystyle {\int }_0^t}R(t,\tau )\frac{{\partial }^2{T}_t(x,\tau )}{\partial x^2}d\tau \\ -K_v{\displaystyle {\int }_0^t}R(t,\tau ){\rho }_b{h}_h\frac{\partial v_h}{\partial x}d\tau +Q_m& (221)\end{array}$

which is the Wulff equation with a fading memory. The coefficients K_α and K_v depend on the type of the memory function used in the hereditary integrals. Applying different memory kernels (with strong physical motivation) equation (221) can be expressed in terms of different fractional operators. However, this draw new research lines which are out of the scope of this work.

11. Some Comments on the Fractionalization Approaches

Now, after the performed analyzes in sections 7 and 7.2 (and the related results in sections 8 and 9) we may draw some conclusions on feasibility of the fractionalization approaches and emerging problems, namely

(1) The simple replacement of time-derivatives in existing integer-order models by fractional counterparts, irrespective of the type of the memory kernels used is blind and erroneous approach without physical background.

(2) Fractionalization by fractional Taylor series expansion of both the heat flux and the temperature gradient is logical, mathematically correct and based on analogies already used in the Maxwell-Cattaneo approach resulting in the hyperbolic heat-wave model and in the Dual-Phase Lag (DPL) model.

(3) Fractionalization by additional mixed time-space derivative is useful approach based on analogies of the modeling equations of the non-Newtonian fluid flow (second order fluid) and the fact that lagging can be obtained by Taylor series (in time) expansion of the temperature gradient. The result is a simple equation with only one damping term in contrast the results of the preceding (point 1 and point 2) final models (see the model of Ezzat, for example).

(4) The fading memory approach is not restrictive in the initial step of fractionalization, as in the cases when fractional Taylor series are applied; it is based on rigorous thermodynamic principles and gives a freedom to apply different memory kernels depending on the experimentally exhibited responses (relaxation function) of the materials (tissues) where the heat conduction should be modeled. A good example is the Wulff models with a fading memory derived here since it demonstrates that both the heat diffusion (conduction) and the perfusion blood flux can exhibit relaxations (i.e., these fluxed have finite speeds of propagation). There are a lot attractive and challenging problems in this area that should be solved.

(5) The fractionalization problem discussed are oriented only to damping heat fluxes in the time domain, that is the tissues where the heat transfer should take place are assumed as homogenous. To some extent, the spacial non-homogeneity of the tissues are modeled by effective transport coefficients (such as in the models of Wulff, Klinger and Weinbaum-Jiji) thus facilitating the process of fractionalization. Space-fractional models are open problems [154] waiting correct build-ups with correct physical analyzes of the both the tissue structures and vascular architectures.

12. Final Comments

Finally, after the long route through the beautiful forest of bio-heat models we may say that the first goal of the analysis done is successfully attained. Precisely, a direct parallel between the existing integer-order models and the corresponding fractional versions was drawn. To some extent, the existing fractional bio-heat models are not systematically derived, excluding the attempt of Ezzat which falls into one of the systematic way how fractionalization could be done. Actually, the good understanding of the contributions of different terms of the discussed models needs nondimensalization since this allows to see the contributions of different processes involved in the heat transfer in tissues. This a standard procedure in mathematical modeling but sometimes the reported studies define the dimensionless variables ad hoc without deep physical analysis and consequent order of magnitude analysis. Last but not least, collecting dispersed studies in both bio-heat model formulations, their analyzes, solution and fractionalizations (if this is well motivated) is a difficult task. We hope the work done here and its results will serve as a good source of information and new ideas, and we may expect that the systematic approach in the formulation of the time-fractional models will draw more correct approaches, from mathematical and physical viewpoints, of course. The correct model formulation is of primary importance in mathematical modeling and this goal has a priority over calculation techniques.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Xu F, Lu TJ, Seffen KA, Ng EYK. Mathematical modelling of skin bioheat transfer Appl Mech Rev. (2009) 62:050801. doi: 10.1115/1.3124646

CrossRef Full Text | Google Scholar

2. Hildebrand M, Goslow GE Jr. Analysis of Vertebrate Structure. New York, NY: Wiley (2001).

Google Scholar

3. Gokul KC, Gurung DB, Adhikary PR. Effect of blood perfusion and metabolism in temperature distribution in human eye. Adv Appl Math Biosci. (2013) 4:13–23.

Google Scholar

4. Nowakowska O, Bulinski Z. Mathematical modelling of heat transport in section of human arm. Comp Assist Meth Eng Sci. (2015) 22:347–63.

5. Hayes KW. Heat and cold in management of rheumatoid arthritis. Arthrit Care Res. (1994) 6:103–54.

PubMed Abstract | Google Scholar

6. Michlovitz SL. Thermal Agents in Rehabilitation. Philadelphia, PA: FA Davis (1996).

7. Petrofsky JS, Lohman E III, Suh J, Garcia J, Anders A, Sutterfield C, et al. Determination of the conductive heat exchange of the skin in relation to environmental temperature. J Appl Res. (2006) 6:157–69.

Google Scholar

8. Parsons KC. Human Thermal Environment. London: Taylor and Francis (1993).

PubMed Abstract

9. Zolfaghari A, Maerefat M. Bioheat transfer. In: Marco Aurelio dos Santos Bernardes, editor. Developments in Heat Transfer Rijeka: IntechOpen (2011). doi: 10.5772/822

CrossRef Full Text | Google Scholar

10. Petrofsky JS, Lind AR. The insulative power of body fat on deep muscle temperature and isometric endurance. J Appl Physiol. (1975) 39:639–42. doi: 10.1152/jappl.1975.39.4.639

CrossRef Full Text | Google Scholar

11. Petrofsky JS, Lind AR. The relationship of body fat content to deep muscle temperature and isometric endurance in man. Clin Sci Mol Med. (1975) 48:405–12. doi: 10.1042/cs0480405

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Torvi DA, Dale JD. A finite element model of skin subjected to a flash fire ASME J Bimech Eng. (1994) 116:250–5. doi: 10.1115/1.2895727

CrossRef Full Text | Google Scholar

13. Yamada Y, Tien T, Ohta M. Theoretical analysis of temperature variation of biological tissue irradiated by light. ASME/JSME Ther Eng. Conf. (1995) 4:575–81.

14. Askarizadeh H, Ahmadikia H. Analytical analysis of the dual-phase-lag model of bioheat transfer equation during transient heating of skin tissue. Heat Mass Transf. (2014) 48:317–27. doi: 10.1007/s00231-014-1373-6

CrossRef Full Text | Google Scholar

15. Kengne E, Sayde M, Lakhssasi A. Effect of convective term on temperature distribution in biological tissue. Eur Phys J Plus. (2013) 128:98. doi: 10.1140/epjp/i2013-13098-8

CrossRef Full Text | Google Scholar

16. Datta A. Biological and Bioenviromental Heat and Mass Transfer. New York, NY: Marcel Dekker (2002).

Google Scholar

17. Valvano JW. Bioheat transfer. In: Webster JG, editor. Encyclopedia of Medical Devices and Instrumentation, Vol. 1, 2nd ed. Hoboken, NJ: Wiley (2005). p. 188–97.

18. Chappuis P, Pittet P, Jequier E. Heat storage regulation in exercise during thermal transients. J Appl Physiol. (1976) 40:384–92. doi: 10.1152/jappl.1976.40.3.384

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Webb P. The physiology of heat regulation. Am J Physiol. (1995) 268:R838–50. doi: 10.1152/ajpregu.1995.268.4.R838

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Orr C-S, Eberhart RC. Overview of bioheat transfer. In: Welch AJ, van Gemert MJC, editors. Optical-Thermal Response of Laser-Irradiated Tissue. New York, NY: Plenum Press (1995). p. 367–84.

Google Scholar

21. Charny CK. Mathematical model of bioheat transfer. Adv Heat Transf. (1992) 22:19–155. doi: 10.1016/S0065-2717(08)70344-7

CrossRef Full Text | Google Scholar

22. Fan J, Wang L. Analytical theory of bioheat transport. J Appl Phys. (2011) 109:104702. doi: 10.1063/1.3580330

CrossRef Full Text | Google Scholar

23. Baish JW, Foster KR, Ayyaswamy PS. Perfused phantom models of microwave irradiated issue. ASME J Biomech Eng. (1986) 108:239–45. doi: 10.1115/1.3138609

PubMed Abstract | CrossRef Full Text

24. Chato JC. Heat transfer to blood vessels. Trans ASME J Biomemech Eng. (1980) 102:110–8.

PubMed Abstract | Google Scholar

25. Chato JC. Fundamentals of bioheat transfer. In: Gautherie M, editor. Thermal Dosimetry and Treatment Planning. Heidelberg: Springer (1990). p. 1–56.

Google Scholar

26. Zhu L, Schappeler T, Cordero-Tumangday C, Rosengard AJ. Thermal interaction between blood and tissue: development of a theoretical approach in predicting body departure during blood cooling/rewarming. In: Minkowycz WJ, SparrowEM, editors. Advances in Numerical Heat Transfer, Vol. 3. Boca Raton, FL: CRC Press, Taylor and Francis Group (2009). p. 197–219.

Google Scholar

27. Pennes HH. Analysis of tissue and arterial blood temperatures in the resting forearm J Appl Phys. (1948) 1:93–122. doi: 10.1152/jappl.1948.1.2.93

CrossRef Full Text | Google Scholar

28. Jiji LM. Heat Conduction, 3rd ed. Berlin: Springer (2009).

29. Khanafer K, Varfai K. Synthesis of mathematical models representing bioheat transport. In: Minkowycz W, Sparrow E, Abraham J, editors. Advances in Numerical Heat Transfer, Vol. 3. CRC Press (2009). p. 1–28.

30. Kreith F, Black WZ. Basic Heat Transfer. New York, NY: Harpers and Row Publ. (1980).

Google Scholar

31. Aziz A, Kraus AD. Transient heat transfer in extended surfaces. Appl Mech Rev. (1995) 48:317–50. doi: 10.1115/1.3005105

CrossRef Full Text | Google Scholar

32. Wulff W. The energy conservation equation for living tissues. IEEE Trans Biomed Eng. (1974) 21:494–5. doi: 10.1109/TBME.1974.324342

CrossRef Full Text | Google Scholar

33. Shitzer A, Kleiner M. Thermal behavior of biological tissues- A General analysis Bull Math Biol. (1976) 38:369–86. doi: 10.1016/S0092-8240(76)80018-3

CrossRef Full Text | Google Scholar

34. Shitzer A, Chato JC. Analytic solutions to the problem of transient heat in living tissue. ASME J Biomech Eng. (1978) 100:202–10. doi: 10.1115/1.3426211

CrossRef Full Text | Google Scholar

35. Mincowycz WJ, Sparrow EM, Abraham JP. Advances in Numerical Heat Transfer. Boca Raton, FL: CRC Press (2009).

Google Scholar

36. Klinger HG. Heat transfer in perfused biological tissue. I., General theory, Bull Math Biol. (1974) 36:403–15. doi: 10.1016/S0092-8240(74)80038-8

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Klinger HG. Heat transfer in perfused biological tissue. I. The macroscopic temperature distribution. Bull Math Biol. (1976) 38:183–99. doi: 10.1007/BF02461434

CrossRef Full Text | Google Scholar

38. Perl W. Heat and master distribution in body tissue and the determination of tissue bllodflow by local clearance methods. J Theor Biol. (1962) 2:201–3. doi: 10.1016/0022-5193(62)90025-5

CrossRef Full Text | Google Scholar

39. Perl W, Hirsh RL. Local blood flow in kidney tissue by heat clearance measurements. J Theor Biol. (1966) 19:251–80. doi: 10.1016/0022-5193(66)90126-3

CrossRef Full Text

40. Klinger HG. The description of heat transfer in biological tissue. Ann N Y Acad Sci. (1980) 335:133–6. doi: 10.1111/j.1749-6632.1980.tb50741.x

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Klinger HG. A simple solution of the non-stationary heat transport problem in capillarized biological tissue. Bull Math Biol. (1984) 46:879–89. doi: 10.1016/S0092-8240(84)80008-7

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Chen MM, Holmes KR. Microvascular contributions in tissue heat transfer. Ann N Y Acad Sci. (1980) 335:137–50. doi: 10.1111/j.1749-6632.1980.tb50742.x

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Weinbaum S, Jiji LM. A two phase theory for the influence of circulation on the heat transfer in surface tissue. In: Wells MK, editor. Advances in Bioengineering. New York, NY: ASME (1979). p. 179–82.

44. Weinbaum S, Jiji LM. A new simplified equation for the effect of blood flow on local average tissue temperature. ASME J Biomech Eng. (1985) 107:131–9. doi: 10.1115/1.3138533

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Weinbaum S, Jiji LM. The matching if thermal fields surrounding countercurrent microvessels and the closure approximation in the Weinbaum-Jiji equation ASME J Biomech Eng. (1989) 111:271–5. doi: 10.1115/1.3168378

CrossRef Full Text | Google Scholar

46. Weinbaum S, Jiji LM, Lemmons DE. Theory and experiment for the effect of vascular microstructure on the surface tissue heat transfer: part 1: anatomical foundation and model conceptualization. ASME J Biomech Eng. (1984) 106:321–30. doi: 10.1115/1.3138501

CrossRef Full Text

47. Weinbaum S, Jiji LM, Lemmons DE. Theory and experiment for the effect of vascular microstructure on the surface tissue heat transfer: part 2: model formulation and solution ASME J Biomech Eng. (1984) 106:331–41. doi: 10.1115/1.3138502

CrossRef Full Text | Google Scholar

48. Song WJ, Weinbaum S, Jiji LM. A theoretical model for peripheral tissue heat transfer in living tissue using the bioheat equation of Weinbaum and Jiji. ASME J Biomech Eng. (1987) 109:72–8. doi: 10.1115/1.3138646

CrossRef Full Text

49. Song WJ, Weinbaum S, Jiji LM, Lemmons DE. A combined macro and microvascular model for whole limb heat transfer. ASME J Biomech Eng. (1987) 110:259–68. doi: 10.1115/1.3108440

PubMed Abstract | CrossRef Full Text

50. Baish JW, Ayyaswamy PS, Foster KR. Small-scale temperature fluctuations in perfused tissues due to local hyperthermia. ASME J Biomech Eng. (1986) 108:246–50. doi: 10.1115/1.3138610

CrossRef Full Text

51. Wissler EH. Comments on the new bioheat equation proposed by Weinbaum and Jiji. ASME J Biomech Eng. (1987) 109:226–33. doi: 10.1115/1.3138673

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Nakayama A, Sano Y, Yoshikawa K. A rigorous derivation of the bioheat equation for local tissue heat transfer based on a volume averaging theory. Heat Mass Transf. (2010) 46:739–46. doi: 10.1007/s00231-010-0619-1

CrossRef Full Text | Google Scholar

53. Lubashevsky IA, Gafiychuk VV. Mathematical Description of Heat Transfer in Living Tissue. Lviv, VNTL Publishers (1999).

Google Scholar

54. Witmore RL. Rheology of Circulation London: Pergamon Press (1968).

55. Crezee J, Lagendgijk JUW. Experimental verification of bioheat transfer theories: measurement of temperature profiles around large artificial vessels in perfused tissue. Phys Med Biol. (1990) 35:905–23. doi: 10.1088/0031-9155/35/7/007

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Talaee MR, Kabiri A. Exact analytical solution of bioheat equation subjected to intensive moving heat source. J Mech Med Biol. (2017) 17:1750081. doi: 10.1142/S0219519417500816

CrossRef Full Text | Google Scholar

57. Kengne E. Akhssasssi, bioheat transfer for one-dimensional spherical biological tissues. Math Biosci. (2015) 269:1–9. doi: 10.1016/j.mbs.2015.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Liu KC. Thermal propagation analysis for living tissue with surface heating. Int J Therm Sci. (2008) 47:507–13. doi: 10.1016/j.ijthermalsci.2007.04.005

CrossRef Full Text | Google Scholar

59. Shen W, Zhang J, Yang F. Modeling and numerical simulations of bioheat transfer and biomechanics in soft tissue. Math Comput Model. (2005) 41:1251–65. doi: 10.1016/j.mcm.2004.09.006

CrossRef Full Text | Google Scholar

60. Lin M, Xu F, Lu TJ, Bai BF. A review of heat transfer in human tooth-Experimental characterization and mathematical modeling. Dental Mat. (2010) 26:501–13. doi: 10.1016/j.dental.2010.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Deng ZS, Liu J. Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies. Trans ASME J Biomech Eng. (2002) 124:638–49. doi: 10.1115/1.1516810

CrossRef Full Text | Google Scholar

62. Bejan A, Dincer I, Lorente S, Miguel A, Reis A. Porous and Complex Flow Structures in Modern Technologies. New York, NY: Springer-Verlag (1984).

Google Scholar

63. Khaled ARA, Vafai K. The role of porous medium in modeling flow and heat transfer in biological tissues. Int J Heat Mass Trans. (2003) 46:4989–5003. doi: 10.1016/S0017-9310(03)00301-6

CrossRef Full Text | Google Scholar

64. Khanafer K, Vafai K. The role of porous media in biomedical engineering as related to magnetic resonance imaging in drug delivery. Heat Mass Trans. (2006) 42:939–53. doi: 10.1007/s00231-006-0142-6

CrossRef Full Text | Google Scholar

65. Nakayama A, Kuwahara F. A general bioheat transfer model based on the theory of porous media. Int J Heat Mass Trans. (2008) 51:3190–9. doi: 10.1016/j.ijheatmasstransfer.2007.05.030

CrossRef Full Text | Google Scholar

66. Nakayama A, Kuwahara F, Liu W A general set of bioheat transfer equations based on the volume averaging theory. In Vafai K, editor. Porous Media, Applications in Biological Systems and Biotechnology. Boca Raton, FL: CRC Press (2010). p. 1–44. doi: 10.1201/9781420065428-2

CrossRef Full Text | Google Scholar

67. Amiri A, Vafai K. Analysis of dispersion effects and nonthermal equilibrium, non-Darcian, variable porosity incompressible-flow through porous media. Int J Heat Mass Trans. (1994) 37:939–54. doi: 10.1016/0017-9310(94)90219-4

CrossRef Full Text | Google Scholar

68. Amiri A, Vafai K. Transient analysis of incompressible flow through a packed bed Int J Heat Mass Trans. (1998) 41:4529–79. doi: 10.1016/S0017-9310(98)00120-3

CrossRef Full Text | Google Scholar

69. Lee DY, Vafai K. Analytical characterization and conceptual assessment of solid and fluid temperature differentials in porous media. Int J Heat Mass Trans. (1999) 42:423–35. doi: 10.1016/S0017-9310(98)00185-9

CrossRef Full Text | Google Scholar

70. Nakayama A, Kuwahara F, Hayashi T. Numerical modelling for three-dimensional heat and fluid flow through a bank of cylinders in yaw. J Fluid Mech. (2004) 498:139–59. doi: 10.1017/S0022112003006712

CrossRef Full Text | Google Scholar

71. Fabrizio M, Giorgi C, Morro A. Modeling heat conduction via fractional derivatives Heat Mass Transf. (2017) 53:2785. doi: 10.1007/s00231-017-1985-8

CrossRef Full Text | Google Scholar

72. Liu KC, Cheng PJ. Finite propagation of heat transfer in multilayer tissue. J Thermopphys Hear Transf. (2008) 22:775–82. doi: 10.2514/1.37267

CrossRef Full Text | Google Scholar

73. Kaminski W. Hyperbolic heat conduction equation for materials with a nonhomogenous inner structure. J Heat Transf. (1990) 112:555–60. doi: 10.1115/1.2910422

CrossRef Full Text | Google Scholar

74. Luikov AV. Analytical Heat Diffusion Theory. New York, NY: Academic Press (1968).

Google Scholar

75. Mitra K, Kumar S, Vedavarz A, Moalemi MK. Experimental evidence of hyperbolic heat conduction in processed meat. J Heat Transfer. (1995) 117:568–73. doi: 10.1115/1.2822615

CrossRef Full Text

76. Davidov EV, Lubashevsky IA, Milyaev IA, Musin RF. Nondiffusive heat transfer in muscle tissue. Preliminary results. arXiv: cond-mat/0102006. I(1). (2001).

Google Scholar

77. Banerjee A, Ogale AA, das C, Mitra K, Subramanian C. Temperature distribution in different material due to short pulse laser. Heat Transf Eng. (2005) 26:41–9. doi: 10.1080/01457630591003754

CrossRef Full Text | Google Scholar

78. Xu F, Lin M, Lu TJ. Modelling skin thermal pain sensation: role of non-Fourier thermal behaviour in transduction process of nociceptor Comp Biol Med. (2010) 40:478–86. doi: 10.1016/j.compbiomed.2010.03.002

CrossRef Full Text | Google Scholar

79. Liu J, Chen X, Xu LX. New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating. IEEE Trans Biomed Eng. (1999) 46:420–8. doi: 10.1109/10.752939

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Yang WH. Thermal (Heat) shock biothermomechanical viewpoint. J Biomed Eng. (1993) 115:617–21. doi: 10.1115/1.2895549

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Fabrizio M, Franchi F. Delayed thermal models: stability and thermodynamics. J Ther Stress. (2014) 37:160–73. doi: 10.1080/01495739.2013.839619

CrossRef Full Text | Google Scholar

82. Cheng L, Xu M, Wang L. From Boltzmann transport equation to single-phase-lagging heat equation. Int J Heat Mass Transf. (2008) 51:6018–23. doi: 10.1016/j.ijheatmasstransfer.2008.04.004

CrossRef Full Text | Google Scholar

83. Tzou DY. Experimental support for the lagging behaviour in heat propagation. J Thermophys Heat Transf. (1995) 9:686–93. doi: 10.2514/3.725

CrossRef Full Text | Google Scholar

84. Tzou DY. Macro- to Microscale Heat Transfer: The Lagging Behaviour. Washington, DC: Taylor and Francis (1997).

Google Scholar

85. Antaki PJ. Analysis of hyperbolic heat conduction in a semi-infinite slab with surface convection. Int J Heat Mass Transf. (1997) 40:3247–50. doi: 10.1016/S0017-9310(96)00351-1

CrossRef Full Text | Google Scholar

86. Tzou DY. An engineering assessment to the relaxation time in thermal wave propagation. Int J Heat Mass Transf. (1993) 36:1845–51. doi: 10.1016/S0017-9310(05)80171-1

CrossRef Full Text | Google Scholar

87. Fabrizio M, Morro A. Viscoelastic relaxation functions compatible with thermodynamics. J Elasticity. (1999) 19:63–75. doi: 10.1007/BF00041695

CrossRef Full Text | Google Scholar

88. Fabrizio M, Morro A. Dissipativity and irreversibility of electromagnetic systems. Math Mod Meth Appl Sci. (2000) 10:217–46. doi: 10.1142/S0218202500000148

CrossRef Full Text | Google Scholar

89. Marin E. On thermal waves. Eur J Phys. (2013) 34:L85. doi: 10.1088/0143-0807/34/5/L83

CrossRef Full Text | Google Scholar

90. Jaisaardsuetrong J, Straughan B. Thermal waves in a rigid heat conductor Phys Lett A. (2007) 366:433–6. doi: 10.1016/j.physleta.2007.02.058

CrossRef Full Text | Google Scholar

91. Kundu B, Dewanjee D. A new method for non-Fourier thermal response in single layer skin tissue. Case Stud. Ther Eng. (2015) 5:79–88. doi: 10.1016/j.csite.2015.02.001

CrossRef Full Text | Google Scholar

92. Bravo ME, De Jesus Sanchez P, Vargas Aguilar RO, Chavez AE. Heat transfer in biological tissues. In: Klapp J, Chavarria GR, Ovando AM, Villa AL, editors. Selected Topics of Computational and Experimental Fluid Mechanics. Cham: Spinger International Publishing (2015). p. 313–20. doi: 10.1007/978-3-319-11487-3_21

CrossRef Full Text | Google Scholar

93. Chen CI, Chen Ch K, Yang YT. Unsteady unidirectional flow of second grade fluid between the parallel plates with different given volume flow rate conditions. Appl Math Comput. (2003) 137:437–50. doi: 10.1016/S0096-3003(02)00149-2

CrossRef Full Text | Google Scholar

94. Hristov J. Integral-balance solution to the stokes first problem of a viscoelastic generalized second grade fluid. Therm Sci. (2012) 16:395–410. doi: 10.2298/TSCI110401077H

CrossRef Full Text | Google Scholar

95. Hristov J. A transient flow of a non-newtonian fluid modelled by a mixed time-space derivative: an improved integral-balance approach. In: Taş K, Baleanu D, Machado J, editors. Mathematical Methods in Engineering. Nonlinear Systems and Complexity, Vol. 24. Cham: Springer (2019). p. 153–74. doi: 10.1007/978-3-319-90972-1_11

CrossRef Full Text | Google Scholar

96. Hristov J. Transient heat diffusion with a non-singular fading memory: from the cattaneo constitutive equation with Jeffrey's kernel to the Caputo-Fabrizio time-fractional derivative. Ther Sci. (2016) 20:765–70. doi: 10.2298/TSCI160112019H

CrossRef Full Text | Google Scholar

97. Joseph DD, Preciozi L. Heat waves. Rev Mod Phys. (1989) 61:41–73. doi: 10.1103/RevModPhys.61.41

CrossRef Full Text | Google Scholar

98. Ozisik MN, Tzou DY. On the wave theory in heat conduction. J Heat Transf. (1994) 116:526–35. doi: 10.1115/1.2910903

CrossRef Full Text | Google Scholar

99. Tang MM, Trujillo M, Lopez Molina JA, Rivera MJ, Berjan EJ. Modeling the heating of biological tissue on the hyperbolic transfer equation. Math Comput Model. (2009) 50:665–72. doi: 10.1016/j.mcm.2008.12.023

CrossRef Full Text | Google Scholar

100. Dai W, Wang H, Jordan PM, Mickens RE, Bejan A. A mathematical model for skin burn injury induced by radiation heating. Int J Heat Mass Transf. (2008) 51:5497–510. doi: 10.1016/j.ijheatmasstransfer.2008.01.006

CrossRef Full Text | Google Scholar

101. Liu J, Zhang X, Wang CC, Lu WQ, Ren Z. Generalized time delay bioheat equation and preliminary analysis on its wave nature. Chinese Sci. Bull. (1997) 42:289–92. doi: 10.1007/BF02882462

CrossRef Full Text | Google Scholar

102. Tzou DY, Zhang YS. An analytical study on the fast-transient process in small scales. Int J Eng Sci. (1995) 33:1449–63. doi: 10.1016/0020-7225(94)00130-C

CrossRef Full Text | Google Scholar

103. Quintanilla R, Racke R. A note on stability in dual-phase-lag heat conduction. Int J Heat Mass Transf. (2006) 49:1209–13. doi: 10.1016/j.ijheatmasstransfer.2005.10.016

CrossRef Full Text | Google Scholar

104. Fabrizio M, Lazzari B. Stability and second law of thermodynamics in dual-phase-lag heat conduction. Int J Heat Mass Transf. (2014) 74:484–9. doi: 10.1016/j.ijheatmasstransfer.2014.02.027

CrossRef Full Text | Google Scholar

105. Fabrizio M, Lazzari B, Tibullo V. Stability and thermodynamic restrictions for dual-phase-lag thermal model. J Non Equilib. Thermodyn. (2016) 42:243–52. doi: 10.1515/jnet-2016-0039

CrossRef Full Text | Google Scholar

106. Tang DW, Araki N, Yamagishi N. Transient temperature responses in biological materials under pulsed IR irradiation. Heat Mass Transf. (2007) 43:579–85. doi: 10.1007/s00231-006-0125-7

CrossRef Full Text | Google Scholar

107. Shih TC, Yuan P, Lin WL, Kou HS. Analytical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux conditions on the skin surface. Med Eng Phys. (2007) 29:946–53. doi: 10.1016/j.medengphy.2006.10.008

CrossRef Full Text | Google Scholar

108. Hooshmand P, Moradi A, Kherzry B. Bioheat transfer analysis of biological tissues by laser radiation. Int J Therm Sci. (2015) 90:214–23. doi: 10.1016/j.ijthermalsci.2014.12.004

CrossRef Full Text | Google Scholar

109. Kengne E, Lakhssasi A. Analytico-numerical study of bio-heat transfer problems with temperature-dependent perfusion. Eur Phys J Plus. (2015) 130:89. doi: 10.1140/epjp/i2015-15089-1

CrossRef Full Text | Google Scholar

110. Kumar A, Kumar S, Katiyar VK, Telles S. Dula phase lag bio-heat transfer during cryosurgery of lung cancer: comparison of three heat transfer models. J Ther Biol. (2017) 69:228–37. doi: 10.1016/j.jtherbio.2017.08.005

CrossRef Full Text | Google Scholar

111. Zhou J, Zhang Y, Chen JK. An axisymmetric dual-phase lag bioheat model for laser heating of living tissues. Int J Term Sci. (2009) 48:1477–85. doi: 10.1016/j.ijthermalsci.2008.12.012

CrossRef Full Text | Google Scholar

112. Zhou J, Zhang Y, Chen JK. Dual-phase lag effects on thermal damage to biological tissues caused by laser irradiation. Comput Biol Med. (2009) 39:286–93. doi: 10.1016/j.compbiomed.2009.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

113. Majchrzak E. Numerical solution of dual phase lag model of bioheat transfer using the general boundary element method. Comput Model Eng Sci. (2010) 69:43–60.

Google Scholar

114. Mochnacki B, Majchrzak E. Numerical model of thermal interaction between cylindrical cryoprobe and biological tissue using the dual-phase lag equation. Int J Heat Mass Transf. (2017) 108:1–10. doi: 10.1016/j.ijheatmasstransfer.2016.11.103

CrossRef Full Text | Google Scholar

115. Tullius TK, Bayazitoglu Y. Analysis of relaxation times on the human head using thermal wave model. Int J Heat Mass Transf. (2013) 67:1007–13. doi: 10.1016/j.ijheatmasstransfer.2013.08.097

CrossRef Full Text | Google Scholar

116. Hristov J. Linear viscoelastic responses and constitutive equations in terms of fractional operators with non-singular kernels: pragmatic approach, memory kernel correspondence requirement and analyses. Eur Phys J Plus. (2019) 134:283. doi: 10.1140/epjp/i2019-12697-7

CrossRef Full Text | Google Scholar

117. Podlubny I. Fractional Differential Equations. New York, NY: Academic Press (1999).

Google Scholar

118. Caputo M, Fabrizio M. A new definition of fractional derivative without singular kernel. Progr Fract Differ Appl. (2015) 1:73–85.

119. Hristov J. Derivatives with non-singular kernels: from the Caputo-Fabrizio definition and beyond: appraising analysis with emphasis on diffusion models. In: Bhalekar S, editor. Frontiers in Fractional Calculus. Sharjah: Bentham Science Publishers (2017). p. 269–342. doi: 10.2174/9781681085999118010013

CrossRef Full Text | Google Scholar

120. Hristov J. Response functions in linear viscoelastic constitutive equations and related fractional operators. Math Model Nat Phenom. (2019) 14:305. doi: 10.1051/mmnp/2018067

CrossRef Full Text | Google Scholar

121. Hristov J. Linear viscoelastic responses: the Prony decomposition naturally leads into the Caputo-Fabrizio fractional operator. Front Phys. (2018) 6:135. doi: 10.3389/fphy.2018.00135

CrossRef Full Text | Google Scholar

122. Ezzat MA, Al-Sowayan NS, Al-Muhiameed ZIA, Ezzat SM. Fractional modelling of Pennes' bioheat transfer equation. Heat Mass Transf. (2014) 50:907–14. doi: 10.1007/s00231-014-1300-x

CrossRef Full Text | Google Scholar

123. Ezzat MA, El-Barry, Al-Sowayan NS. Tissue responses to fractional transfer heating with sinusoidal heat flux condition on skin surface. Anim Sci J. (2016) 10:1304–11. doi: 10.1111/asj.12568

CrossRef Full Text | Google Scholar

124. Jumarrie GD. Derivation and solutions of some fractional Black-Sholes equations in coarse-grained space and time, applications of Merton's optimal portfolio. Comput Math Appl. (2010) 59:1142–64. doi: 10.1016/j.camwa.2009.05.015

CrossRef Full Text | Google Scholar

125. Damor RS, Kumar S, Shukla AK. Solution of fractional bioheat equation in terms of Fox's H-function. SpingerPlus. (2016) 5:111. doi: 10.1186/s40064-016-1743-2

PubMed Abstract | CrossRef Full Text | Google Scholar

126. Damor RS, Kumar S, Shukla AK. Temperature distribution in living tissue with fractional bioheat model in thermal therapy. In: Patel H, Deheri G, Patel H, Mehta S, editors. Proceedings of the International Conference on Advances in Tribology and Engineering Systems. Lecture Notes in Mechanical Engineering. Springer India (2014). p. 493–8.

Google Scholar

127. Ferras LL, Ford NJ, Morgado ML, Nóbrega JM, Regbelo MS. Fractional Pennes' bioheat eqaution: theoretical and numerical studies. Frac Cal Appl Anal. (2015) 18:10180–1106. doi: 10.1515/fca-2015-0062

CrossRef Full Text | Google Scholar

128. Gurtin ME. A general theory of heat conduction with finite wave speed. Arch Rat Mech Anal. (1968) 31:113–26. doi: 10.1007/BF00281373

CrossRef Full Text | Google Scholar

129. Hardy G. Riemann's form of Taylor series. J London Math Soc. (1945) 20:48–57. doi: 10.1112/jlms/s1-20.1.48

CrossRef Full Text | Google Scholar

130. Watanabe Y. Notes on generalized derivatives of Riemann-Liouville and its application to Leibniz's formula. Tohoku Math J. (1931) 24:8–41.

131. El-Ajou A, Arqub OA, Al Zhour Z, Momani S. New results on fractional power series: theories and applications. Entropy. (2013) 15:5305–23. doi: 10.3390/e15125305

CrossRef Full Text | Google Scholar

132. Trujillo J, Rivero M, Bonilla B. On a Riemann-Liouville generalized Taylor's formula J Math Anal Appl. (1999) 231:255–65. doi: 10.1006/jmaa.1998.6224

CrossRef Full Text | Google Scholar

133. Odibat Z, Shawagfeh N. Generalized Taylor's formula. Appl Math Comput. (2007) 186:286–93. doi: 10.1016/j.amc.2006.07.102

CrossRef Full Text | Google Scholar

134. Mehdinejadiani B, Jafari H, Baleanu D. Derivation of fractional Boussinesq equation for modelling unconfined groundwater. Eur Phys J. (2013) 222:1805–12. doi: 10.1140/epjst/e2013-01965-1

CrossRef Full Text | Google Scholar

135. Mehdinejadiani B, Naseri AA, Jafari H, Chanbarzadeh A, Baleanu D. A mathematics model for simulation of a water table profile between two parallel subsurface drains using fractional derivatives. Comp Math Appl. (2013) 66:785–94. doi: 10.1016/j.camwa.2013.01.002

CrossRef Full Text | Google Scholar

136. Kang J, Xu M. Exact solutions for unsteady unidirectional flows of a generalized second -order fluid through a rectangular conduit. Acta Mech Sin. (2009) 25:181–6. doi: 10.1007/s10409-008-0209-3

CrossRef Full Text | Google Scholar

137. Bandelli R, Rajagopal KR. Start-up flows of second grade fluids in domains with one finite dimension. J Non-Linerar Mech. (1995) 30:817–39. doi: 10.1016/0020-7462(95)00035-6

CrossRef Full Text | Google Scholar

138. Coleman B, Noll W. Foundations of linear Viscoelasticity. Rev Mod Phys. (1961) 33:239–49. doi: 10.1103/RevModPhys.33.239

CrossRef Full Text | Google Scholar

139. Coleman B, Gurtin ME. Equipresence and constitutive equations for rigid heat conductors. Z Angew Math Phys. (1967) 18:188–208. doi: 10.1007/BF01596912

CrossRef Full Text

140. Nachlinger RR, Weeler L. A uniqueness theorem for rigid heat conductors with memory Q Appl Math. (1973) 31:267–73. doi: 10.1090/qam/436794

CrossRef Full Text

141. Nunciato JW. On heat conduction in materials with memory. Q Appl Math. (1971) 29:187–273. doi: 10.1090/qam/295683

CrossRef Full Text

142. Nunciato JW. On uniqueness in the linear theory of heat conduction with finite wave speds. SIAM J Appl Math. (1973) 25:1–4. doi: 10.1137/0125001

CrossRef Full Text | Google Scholar

143. Boltzmann L. Zur Theorie der Elastischen Nachwirkung. Sitzungsber Akad Wiss. Wien Mathem Naturwiss. (1874) 70:275–300.

144. Miller RK. An integrodifferential equation for rigid heat conductors with memory. J Math Anal Appl. (1978) 66:313–32. doi: 10.1016/0022-247X(78)90234-2

CrossRef Full Text | Google Scholar

145. Cattaneo C. On the conduction of heat (In Italian). Atti Sem Mat Fis Univ Modena. (1948) 3:83–101.

146. Araujo A, Ferreira JA, Oliveira P. The effect of memory terms in diffusion phenomena. J Comp Math. (2000) 24:191–102.

Google Scholar

147. Westerlund S. Dead matter has memory? Phys Scrip. (1991) 43:174–9. doi: 10.1088/0031-8949/43/2/011

CrossRef Full Text | Google Scholar

148. Muslih SI, Baleanu D. Fractional multipoles in fractional space. Nonlin Anal Real World Appl. (2007) 8:198–203. doi: 10.1016/j.nonrwa.2005.07.001

CrossRef Full Text | Google Scholar

149. Ebaid A, Masaedeh B, El-Zahar E. A new fractional model for the falling body problem. Chin Phys Lett. (2012) 34:029201. doi: 10.1088/0256-307X/34/2/020201

CrossRef Full Text | Google Scholar

150. Gomez-Aguilar JF, Rosales-Garcia JJ, Bernal-Alvarado JJ. Fractional mechanical oscillators. Rev Mex Fis. (2012) 58:348–52.

Google Scholar

151. Gomez-Aguilar JF. Space-time fractional diffusion equation using a derivative with non-singular and regular kernel. Physica A. (2017) 465:562–72. doi: 10.1016/j.physa.2016.08.072

CrossRef Full Text | Google Scholar

152. Dokoumetzidis A, Macheras P. Fractional kinetics in drug absorption and disposition processes. J Pharmacokinet Pharmacodyn. (2009) 3:165–78. doi: 10.1007/s10928-009-9116-x

CrossRef Full Text | Google Scholar

153. Liu J. Preliminary survey on the mechanism of the wave-like behaviors of heat transfer in living tissues. Forsch. Ingenieurwesen. (2000) 66:1–10. doi: 10.1007/s100100000031

CrossRef Full Text | Google Scholar

154. Damor RS, Kumar S, Shukla AK. Parametric study of fractional bioheat equation in skin tissue with sinusoidal heat flux. Frac Diff Cal. (2015) 5:43–53. doi: 10.7153/fdc-05-04

CrossRef Full Text | Google Scholar