advertisement

Journal of Thermal Biology 84 (2019) 274–284 Contents lists available at ScienceDirect Journal of Thermal Biology journal homepage: www.elsevier.com/locate/jtherbio Numerical simulation of fractional non-Fourier heat conduction in skin tissue T P. Goudarzi∗, A. Azimi Department of Mechanical Engineering, Shahid Chamran University of Ahvaz, Ahvaz, Iran ARTICLE INFO ABSTRACT Keywords: Fractional non-Fourier-non-Fourier conduction-multilayer skin-dual-phase lag In this paper, a fractional non-Fourier heat conduction model is employed to simulate the heat diffusion through the skin tissue, as a biological system, upon immediate contact with a heat source. In order to study skin models and different boundary aspects, two problems: the three-layer skin tissue in contact with a hot water source and a single-layer skin tissue exposed suddenly to a heat source generated by a laser are investigated. In both cases, the super-diffusion fractional non-Fourier model is used to simulate the heat transfer diffused through the skin tissue. In the first case, the governing equation is solved using an implicit method, and in the second problem, its governing equation is solved using a finite volume method. In the fractional non-Fourier model, the effect of the model's essential parameters ( and ) on the prediction of temperature distribution in skin tissue is studied as well as the effect of other parameters such as the blood rate is studied. In addition, grid study has been investigated and the most efficient and appropriate gird is obtained. The results are validated against the DPL (Dual-Phase Lag) model's results. The fractional single-phase-lag model's results indicate that this model is highly precise and encompasses all the results of the dual-phase-lag model. The results also show the high precision of the model, taking into account both the microstructure interactions and the lags. 1. Introduction The Pennes heat transfer model is commonly used for the simulation of thermal behavior in bodies and biological tissues due to its simplicity and reliability. The Pennes bio-heat equation describes thermal behavior based on the Fourier classic law which shows the heat signals diffusion's infinite velocity. In fact, the living tissues are highly inhomogeneous and require a relaxation time in order to reserve sufficient energy and transfer it to the nearest cell. As a result, the SinglePhase Lag (SPL) biological heat transfer was proposed for the first time to solve the inconsistency in Pennes model, and to study the physical mechanisms and thermal wave behavior in living tissues (Jing Liu et al., 1999). The characteristics of wave diffusion in biological heat transfer have attracted the researchers’ interest. Liu (2008) and Özen et al. (2008) studied the biological heat transfer of thermal wave diffusion respectively in multilayer and inhomogeneous tissues. Shih et al. (2005) investigated the effect of thermal wave characteristics on temperature distribution during thermal treatment. However, to take into account the microstructure effects in the rapid transient heat transfer process, a phase lag was introduced for the absent temperature gradient in SPL (Cattaneo, 1958; Vernotte, 1958) which was named the Dual Phase Lag ∗ (DPL) model (Tzou, 1997). Antaki (2005) utilized DPL model to describe heat conduction in processed meat. Antaki (2005) estimated the time phase lag for heat flux to be 14–16 s, and the temperature gradient phase lag to be −0.043–0.056 s for processed meat. Xu et al. (2008) addressed the use of DPL model in skin tissue thermal behavior. Liu et al. (2012) employed DPL model to analyze temperature distribution in the three-layer skin. Hooshmand et al. (2015) develop an analytical solution for a generalized DPL model based on the nonequilibrium heat transfer in biological tissues during laser irradiation. Lin and Li (2016) studied heat transfer of skin subjected to the pulse laser heating and fluid cooling and they use an analytical method for solving the general problem of heat conduction (DPL method). On the other hand, in the past three decades, the efficiency of the fractional calculus, where the derivative and integral are rational, is proved in the simulation of abnormal behaviors observed in physical phenomena. However, until recent decades, no significant use of fractional calculus has been found in scientific and engineering fields, and it has only been studied in mathematics. Today, due to its inherent capability, the fractional calculus is highly promising in the simulation of anomalous behaviors in highly complicated processes in different scientific fields. In fact, the Fractional Single-Phase Lag (FSPL) is, in comparison to SPL and DPL, a Corresponding author. E-mail addresses: pirangoodarzi@yahoo.com (P. Goudarzi), a.azimi@scu.ac.ir (A. Azimi). https://doi.org/10.1016/j.jtherbio.2019.05.021 Received 23 February 2019; Received in revised form 4 May 2019; Accepted 20 May 2019 Available online 23 May 2019 0306-4565/ © 2019 Elsevier Ltd. All rights reserved. Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi newfound field in describing non-Fourier heat conduction, and is capable of describing the rapid transient effects of heat transfer (time effects) and the effects of microstructural interactions (space effects), simultaneously. Some major applications of fractional calculus are simulations of anomalous diffusion, viscoelasticity, the control theory, and bio-engineering. Ghazizadeh et al. (2012) employed FSPL for the first time to model non-Fourier conduction in processed meat. Povstenko (2011) generalizes theory of thermal stresses based on the fractional Cattaneo-type heat conduction equation. Povstenko and Povstenko (2013) studied fractional heat conduction in a composite medium consisting of a spherical inclusion. The heat conduction is described by the time-fractional heat conduction equation with the Caputo derivative of fractional order and used Laplace transform to solve equations. Zhang et al. (2014) model one-dimensional heat conduction problem where two cold waves collide in a layer and compare the solutions of the parabolic equation, the hyperbolic heat equation, the fractional heat equation, and the generalized Cattaneo equations. They investigated effect of fractional parameter on results. Povstenko (2015) published a book as fractional thermoelasticity consists of each generalization of heat conduction equation results in formulation of the corresponding generalized theory of thermal stresses. This book devoted to fractional thermoelectricity. Ezzat et al. (2016) used a fractional model of bioheat equation for describing quantitatively the thermal responses of skin tissue under sinusoidal heat flux conditions on skin surface. They use analytical method to solve problem. Yu et al. (2016) formulate a fractional thermal wave model for a bi-layered spherical tissue. They estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model by means of inverse analysis with Levenberg–Marquardt method. The present study is innovative compared to previous studies, specifically on biological tissue, in that it utilizes and analyzes FSPL results for the first time on skin tissue under the influence of transient heat flux due to the laser. This investigation is indicative of the efficiency and capability of such model in laser treatment and medicine. In this paper, two skin tissue samples as biological systems have been studied, and the results reveal the model's efficiency and high precision for investigating heat behavior of skin tissue. Thus, two problems of single-layer and three-layer skin models have been investigated under various conditions in order to expand the use and capability of FSPL. In both cases, the superdiffusion FSPL is used so as to simulate heat transfer diffused in skin tissue. In FSPL model, the effect of essential parameters ( and ) on prediction of heat and temperature history in skin tissue has been studied along with the effect of other parameters such as blood perfusion on temperature distribution. The grid study has also been done and the most efficient and appropriate grid has been obtained. The results have been validated against the results of DPL model. The FSPL results indicate that the model is highly precise and encompasses all the results of the Fourier model, the hyperbolic model, and the Dual Phase Lag model. The results also show the precision of the model, taking into account both the microstructure interactions and the times lags. c c b cb (Tb T ) + qmet + qext (2) T =k t 2T + wb b cb (Tb T ) + qmet + qext (3) Where T skin temperature q is heat flux, c is specific heat capacity of the tissue, wb is Perfusion rate of blood, b is Density of blood, cb is specific heat capacity of the blood, Tb is temperature of blood, qmet is metabolic heating source and qext is external heat generation. 2.2. SPL Non-Fourier conduction model Since the classic Fourier model is a thoroughly diffusive nature model and does not take into account the limited speed of propagating heat stimulations in transient conditions, Cattaneo and Vernotte (Cattaneo, 1958; Vernotte, 1958) independently introduced the following model based on relaxation between temperature gradient and heat flux to avoid such physical (Cattaneo, 1958; Vernotte, 1958), q (r , t + q) = (4) k T (r , t ) where q represent heat flux phase lags. By combining Eq. (4) and Eq. (2) and eliminating the heat flux q, the bio-heat transfer equation of SPL model will be obtained: k 2T = 1+ q t c T t wb b cb (Tb T) qmet qext (5) 2.3. DPL Non-Fourier conduction model Since the microstructure interactions are not considered in SPL model, Tzou presented DPL model to take into account such effects by adding the temperature phase lag term with SPL model. After that, DPL equation is stated as follows (Tzou, 1997): q (r , t + q) = k T (r , t + (6) T) where q and T represent heat flux and, temperature gradient phase lags, respectively. When q < T , the heat flux is generated before the temperature gradient, and if T < q , the temperature gradient is generated before the heat flux. If it is assumed T = 0 , Eq. (6) turns into SPL model, and if it is assumed q = T , Eq. (6) turns into the Fourier heat conduction equation. By combining Eq. (6) and Eq. (2) and eliminating the heat flux q from both equations, the general form of DPL heat conduction model will be obtained as follows (Zhou et al., 2009): 1+ T t k 2T = 1+ q t c T t wb b cb (Tb T) qmet qext (7) As it can be seen from the above equation, using DPL model increases the number of terms and creates compound derivatives on the right-hand side of the DPL heat conduction equation. These added terms and derivatives complicate the numerical and analytical solutions of the DPL equation with respect to those of SPL model. Tzou and Dai (2009) extended the concept of the time lag for simulating heat transfer in the media involving multiple heat carriers. In his study, he showed that as the number of heat carriers increases, the order of the time lag parameters also increases in the final equation. Increasing the order, in turn, complicates the physical response of the obtained equation due to creation many and higher order compound derivatives. The presence of such higher order and compound derivatives results in a limitation of mathematical stability in the equation and it makes analytical solution and numerical solution difficult. In the past two decades, extensive efforts have been made in order Bio-heat conduction models have been addressed in this section. 2.1. Pennes Fourier conduction model The diffusion term in Pennes bio-heat transfer equation is based on the Fourier classic law (Pennes, 1948), implying the infinite speed of the diffusion of heat signals, introduced by Joseph Fourier in his studies on heat conduction. k T (r , t ) div (q (r, t )) + wb And eliminating the heat flux q, the general bio-heat transfer equation will result: 2. Bio-heat conduction models q (r , t ) = T = t (1) By combining the Fourier law Eq. (1) in energy balance equation as: 275 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi Fig. 1. Schematic view of physical model. to study anomalous diffusion phenomenon. Here, fractional calculus has successfully proved high capability at offering a more precise description of mid-behaviors observed in most physical phenomena (Djorev et al., 2003). Non-Fourier heat transfer phenomenon, one of the anomalous diffusion categories, is also an active practical field for fractional calculus. c f t = r f (t r t) , r r=0 =( 1)r r Dt f (t ) = f t = 1 (n < n, n dn ) dt n q (r , t ) + (8) )n 1f ( )d ,n 1 (9) c f (t ) = f t = 1< )n (t (n ) <n ,n 1f (n) ( )d n (12) k T (r , t ) (13) + c b q 1+ T t1+ cb (Tb =k 2 T (r , qgen t) + t + qgen T ) + qmet + qext (14) 4. Governing equations, discretization and grid generation a (10) N 4.1. Test case 1 The skin single-layer tissue of thickness L = 1 mm and initial temperature T0 = 37 0C are shown in Fig. 1. At t = 0+, a laser beam with a constant intensity (Qin ) is shone onto the left skin surface boundary for 5 s. After 5 s, the laser stops. Fig. 1 depicts one-dimensional computational grid for the tissue. The skin single-layer model is used in this problem (Zhou et al., 2009). Considered the behavior of the problem under analysis, the fractional equations super-diffusion version is used, and the governing equation of the problem for FSPL is as follows: 3.1. Fractional Cattaneo model q + t Compte and Metzler (1997) generalized the Cattaneo model to the time fractional Cattaneo models and studied the properties of the corresponding fractional Cattaneo equations and introduced a fractional time derivative operator in the equation as follows: q (r , t ) = t + qgen t 1 Caputo Definition has the ability to consider the initial condition and boundary conditions. In this paper, initial conditions and boundary conditions are very important. Therefore, we use a fractional derivative based on Caputo definition. For more details on basic fractional theory and fractional mathematics, see reference Podlubny (1999). Among these, two heat transfer equations can be referred to as fractional models introduced in heat transfer. q (r , t ) + T t qgen = wb Caputo. c D a t t T ) + qmet + qext q (r , t ) = t q a N qgen t) + where 0 < < 1 is the fractional order of the derivative. is the t fractional order derivative based on Caputo definition. By combining Eq. (13) and Eq. (2), and then eliminating the heat flux q, FSPL equation results: t (t 2 T (r , In SPL and DPL models, the heat conduction equation was obtained by generalizing the first order of Taylor series in energy equation and eliminating heat flux between the two equations. Recently, Ghazizadeh (2010) generalized SPL model to FSPL model. They stated FSPL model by applying fractional Taylor series formula to SPL as follows: Riemann—Liouville: a cb (Tb =k 3.2. Fractional Non-Fourier conduction model [t / t ] 1 t b 1 + T (r , t ) t1+ c Through changing fractional derivative order, , in Eq. (12) from 0 to 1, heat transfer mechanism shifts from Fourier classic model to heat wave model. For 0 < < 1, Eq. (12) models mid-behaviors of heat transfer between pure diffusion and heat wave. In this section, we present basic equations of fractional calculus. Assume f (t ) continues at [a, t]. There are Several definitions of a fractional derivative such as Grunwald-Letnikov's definition, RiemannLiouville's definition, and Caputo's definition. These three definitions are defined as (Podlubny, 1999). Grunwald—Letnikov: Dt f (t ) = + qgen = wb 3. Basic definition of fractional calculus a T (r, t ) t k 2 T (r , t) 1+ q q t1 + =D 2q x2 + Dwb b cb T x (15) where D = is the thermal diffusion coefficient. The problem's boundary conditions are (Zhou et al., 2009): k c q = Qin (1 q=0 (11) The Operator in the above equation is formulated according to t Riemann-Liouville definition. By combining the above relation (11) with Eq. (2) and eliminating heat flux q, the generalized Cattaneo equation (CV) is obtained as: Rd) for x = 0 , for x = L , 0<t<t , 0 < t < tf , t = 5s t f = 40 s (16) where Rd represents diffuse reflectance of the laser shone on the surface, Qin laser intensity and t duration of laser shone on the surface. Values of these parameters are presented in Table 1. The initial conditions are as follows (Zhou et al., 2009): 276 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi ice mixture for 30 s. In order to model skin physics, a three-layer skin model of thickness L = 6 mm is used which consists of Epidermis, Dermis, and Fat, as schematically depicted in Fig. 2 (Kuo-Chi et al., 2011). As in the previous problem, the super-diffusion fractional equation based on caputo definition is employed considering the behavior of the problem under analysis, and for FSPL fractional model, the governing equation is written: Table 1 Values of the laser's parameters (Zhou et al., 2009). Parameter Value Diffusion reflection Laser intensity ( 0.05 20 W ) cm2 Laser duration (s) 5 T t + + q 1+ T t1+ qext + qmet i ci + q + wb b cb q i ci T t (qext + qmet ) i ci t 2T = Di + x2 q i ci + wb b cb (Tb i ci (wb b cb Tb) t T) (20) k where Di = ic , i = 1,2,3 is the thermal diffusion coefficient for three i i skin layers in Table 3and is the fractional order. The boundary conditions are as follows (Kuo-Chi et al., 2011):where u (t ) is unit step function that defines: T (0, t ) = 100 100u (t T (L , t ) = Tb t > 0 q (x ) = 0, u (t (17) q + wb x b c b (Tb T ) + qmet + qext t = TPt + t c qEt+ t qWt+ 2 x t + wb b cb (Tb (18) TPt ) + qmet 15) = T (x , 0) = Tb, (21) 1 0 t > 15 t 15 (22) T (x , 0) = t 2 T (x , t2 0) =0 (23) The boundary conditions have been obtained at the interfaces of two layers (ED, DF) assuming that temperature and heat flux are continuous. The governing equations are discretized using an implicit method. For fractional equation discretization, some fractional terms require being discretized. Based on Caputo fractional derivative, they are discretized in appendix B. After that, the resulted algebraic system of equations is solved using TDMA method, and temperature distribution through the skin is obtained. where qext is zero. In appendix A, used finite volume method for discretization governing equation. Eventually, equation (A-7) is written for all grid elements. After that, by solving the resulted equations system using TDMA method, the heat flux is obtained, which can after that be employed to yield temperature distribution by the following equation: TPt+ t>0 Moreover, the initial conditions are: After calculating heat flux by FSPL model, energy balance is employed in order to obtain temperature distribution (Zhou et al., 2009), T c = t 45 Tb = 37 0C Fig. 2. Three-layer skin coordinate and geometry. q = 0 for 0 < x < L when t = 0 t 15) (19) 4.3. Grid generation 4.2. Test case 2 For both problems, the grid study has been done both in temporal and in spatial dimensions. The most efficient grid obtained for the first problem is 120 spatial grids and 90-time step in Fig. 3, and for the second problem, it equals 900 spatial grids and 500-time step in Fig. 4. At t = 0 , the skin surface is exposed to an immediate contact with a heat source such as hot water at a constant 100 0C for 15 s. Then, the heat source is removed, and the skin surface is cooled using 0 0C water- Fig. 3. Grid study for case 1. 277 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi Fig. 4. Grid study for case 2. 5. Numerical results and discussion in the skin. However, as can be seen, in smaller values of that heat transfer mechanism tends towards pure diffusion, fluctuations decrease, which means less energy accumulation. Fig. 6 shows the effect of phase lag on temperature history at the boundary affected by laser at a constant fractional order, = 0.98. Here, the effect of phase lag is the same as the effect of fractional order in Fig. 5, i.e. its increase results in temperature increase. Figs. 5 and 6 conclude that and both take non-Fourier effects into account and are dependent on each other. Here, in order to indicate the accuracy of the FSPL model, its results are compared with DPL model results available in reference (Zhou et al., 2009). Fig. 7 shows the comparison of FSPL model's results with DPL model's results with high accuracy. DPL model's results were obtained for q = 16 and T = 0.05. In FSPL model, the necessary parameters, = 0.9985 and = 16, have been calculated through trial and error and they well predicted DPL model's results. According to the above results, these two parameters ( and ) in FSPL model predict both the environment transient heat flux effects ( q ) and microstructure interactions effects ( T ) in DPL model. In addition to contributing factors analyzed on heat transfer in skin tissue, other factors such as blood rate, laser intensity, and diffuse reflectance, have In the fractional non-Fourier model, the two essential parameters ( and ) are unknown. Through trial and error, and analyzing their effect on the predicted thermal behavior, values of these parameters have been obtained, and they predict DPL's results with high precision. In addition, the results have been validated against DPL results indicating a high level of accuracy. 5.1. Test case 1 In this case, the single-layer skin model is employed, and its physiological properties are presented in Table 2 (Zhou et al., 2009). Governing equations discretized by FVM method that given equation (A-6), by solving equation (A-6) using TDMA method, the heat flux is obtained and finally temperature distribution obtained with equation (19). Fig. 5 shows the effect of the fractional order parameter, , on temperature history on the boundary affected by laser, at a constant time lag, = 16 s . In addition, this figure shows that increase in results in increasing temperature, further approaching wave mode and widening fluctuations. It shows that decrease in to zero, results show Fourier model. The smaller the fractional order is, the more temperature distribution approaches that of diffusion or Fourier case. In addition, due to present phase lag and high transient flux, there are fluctuations in temperature history. These temperature fluctuations occur due to energy accumulation at the position imposed high heat flux. After a relaxation time, this energy is released and transferred to another position Table 2 The skin single-layer tissue's physiological properties (Zhou et al., 2009), test case 1. Value Properties 1000 Density ( 4187 Specific Heat Capacity ( 0.628 Thermal Conductivity ( kg ) m3 J ) kg K W m 0K ) 1060 Blood Density 3860 Blood Specific Heat Capacity ( 0.00187 kg ( 3) m Blood Perfusion ( m3 m3 tissue s 1190 Metabolic Heat ( W3 ) 37 Blood Temperature (0C ) J ) kg K ) m Fig. 5. The effect of fractional order ( ) on temperature history on the boundary affected by laser ( x = 0 ) in skin tissue, = 16, test case 1. 278 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi Fig. 8. The effect of blood perfusion on temperature history on the boundary affected by laser (x = 0 ) in skin tissue in FSPL mode with = 0.9985and = 16, test case 1. Fig. 6. The effect of on temperature history on the boundary affected by laser ( x = 0 ) in skin tissue, = 0.98, test case 1. an influence on heat transfer in skin tissue. Fig. 8 indicates the effect of blood perfusion parameter on heat transfer in skin tissue. As a result of higher blood rate, skin transfers more heat through to blood due to a convection process, and thus skin temperature drops. Also, by looking at Fig. 8, an increase in blood perfusion predicts lower temperature distribution. Fig. 9 depicts the effect of laser intensity shone on skin surface on temperature history. With an increase in laser intensity, the affected surface experiences temperature increases and the temperature significantly drops when the heat source is removed. Higher or lower laser intensity is employed depending on the thermal treatment type. Fig. 10 indicates the effect of skin surface's diffuse reflectance parameter on temperature history. According to the diagram, it can be stated that as the skin tissue diffuse reflectance is higher, it will diffuse a significant amount of the laser heat in the environment and absorb a little amount; therefore, increasing the diffuse reflectance leads to the lower temperature distribution in the skin tissue. Fig. 7. FSPL model temperature history with = 0.9985 and = 16on the boundary affected by laser (x = 0 ) and its comparison with that of DPL model q = 16 and T = 0.05 (Zhou et al., 2009), test case 1. Fig. 9. The effect of laser intensity on the temperature history on the boundary affected by (x = 0 ) laser in skin tissue in FSPL mode with = 0.9985and = 16, test case 1. 5.2. Test case 2 Now, the results of test case 2 will be assessed. In this case, skin tissue is assumed to consist of three layers, Epidermis, Dermis and Fat, and its surface is exposed to a hot water source. The physiological properties of each tissue layer are different as seen in Table 3 (Kuo-Chi et al., 2011). Governing equations discretized by FDM method that given equation (B-5), by solving equation (B-5) using TDMA method, temperature distribution obtained. 279 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi Fig. 12. The effect of fractional parameter ( ) on temperature history at ED interface in skin tissue, = 0.98, test case 2. Fig. 10. The effect of surface's diffuse reflectance on temperature history on the boundary affected by laser ( x = 0 ) in skin tissue in FSPL mode with = 0.9985 and = 16, test case 1. Fig. 11 shows the effect of the fractional order parameter, , in a constant phase lag, = 10 s , on temperature history at the skin first and second layer interface. The increase in fractional order parameter predicts a higher temperature. Also, Fig. 12 illustrates the effect of phase lag parameter, , in a constant fractional order, = 0.98, on temperature history. As concluded in the previous problem, an increase in phase lag value predicts a higher temperature and these two essential parameters are highly correlated. As can be seen, there is a significant temperature increase when skin surface is contacted with the hot water source for 15 s, and when the heat source is removed, a huge temperature drop occurs. If attention is paid to the moment of adding and removing hot water, no sudden increase or drop in temperature will occur at that moment. In fact, this increase or drop in temperature happens with a lag time which FSPL model has predicted with a significant accuracy taking into con- Table 3 Physiological properties of skin multilayer tissue (Kuo-Chi et al., 2011), test case 2. Fat Dermis Epidermis Blood Properties 971 1116 1190 1060 Density ( 2700 3300 3600 3770 Specific Heat Capacity ( 0.185 0.445 0.235 – Thermal Conductivity ( 368.3 368.1 368.1 – Metabolic heat 0.0044 0.0015 0.0001 – Thickness (m) kg ) m3 J kg 0 K W m 0K ) ) ( W3 ) m Fig. 13. The effect of fractional order ( ) on temperature history at DF interface in skin tissue, = 10 , test case 2. Fig. 11. The effect of fractional order ( ) on temperature history at ED interface in skin tissue, = 10 , test case 2. 280 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi Fig. 14. The effect of fractional order ( ) on FSPL model temperature distribution at the 15th second of the solution for = 10 and various , test case 2. Fig. 16. FSPL model temperature history with = 0.98and face and its comparison with that of DPL model q = 10 and et al., 2011), test case 2. sideration all of its effects. Fig. 13 illustrates the effect of fractional order in FSPL model on the tissue's history distribution at DF interface. As evident in the diagram, varying from 0 to approximately 1 would change the temperature distribution from that of diffusion mode to that of nearing wave mode. On the other hand, comparing the results in Figs. 11 and 13 clearly reveals that at the second interface (DF), which is further away from the skin surface, a lower temperature is predicted with a less sudden increase in a longer period. Figs. 14 and 15 respectively show the effects of on temperature distribution across skin tissue at 15 and 45 s after imposing the water heat source on the surface's skin tissue. As in the previous figures, the ‘s effects are clearly visible. The results indicate the fractional model's capability at giving a non-Fourier description under certain circumstances. It can prove the environment's diffused nature, its wavy nature and diffusion-wave mid modes by varying the values of and T = 10 at ED inter= 0.005 (Kuo-Chi parameters in the FSPL model. Also, it is observed in Figs. 14 and 15 that at the first 15 s, when the skin is in contact with hot water, the temperature approaches equilibrium farther away from the skin surface. For validation, the results of the FSPL model are compared to the results of DPL conducted in reference (Kuo-Chi et al., 2011). In Fig. 16, the temperature history results at the ED layers’ boundary, and in Fig. 17, the results of temperature distribution in skin tissue at the 15th second of the solutions are compared with those of DPL model. DPL model results are for q = 10 and T = 0.005, and the FSPL model results are for = 10 and = 0.98 through trial and error. Fig. 17. FSPL model temperature distribution with = 0.98and = 10 at the 15th second of the solution and its comparison with that of DPL model q = 10 and T = 0.005 (Kuo-Chi et al., 2011), test case 2. Fig. 15. The effect of fractional order ( ) on FSPL model temperature distribution at the 45th second of the solution for = 10 and various , test case 2. 281 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi The FSPL model can make a highly accurate prediction of DPL model results. regard. The FSPL model's results indicate that this model is highly precise and encompasses all the results of the dual-phase-lag model. The results also show the high precision of the model, taking into account both the microstructure interactions and the lags. It can be concluded that fractional non-Fourier model is highly efficient in describing heat in such phenomena as laser therapy in medicine. Application of this model can open a new field in medical advances and thermal aspects of the biological systems. For instance, with an empirical experiment conducted on a tissue, this FSPL model can help to obtain the required parameters ( and ) unique to that biological system. Prediction of temperature distribution and its effects can improve and enhance the precision of laser therapy or similar therapies in that biological system. 6. Conclusion Because of the sensitivity of biological tissues and the advances of science to cure diseases, it is required to consider the thermal aspects of these biological systems. In this paper, after assessing different conditions and obtaining the results, the capability of the FSPL model at describing various physical conditions in skin tissue has been concluded, and this model can be considered a more comprehensive model for describing heat transfer in skin tissue as a biological system to include all Fourier and non-Fourier models so far evaluated with this Symbols and Abbreviations c Specific heat capacity of the tissue [ cb Specific heat capacity of blood [ D Thermal diffusion coefficient [ Density of tissue k J ] kg K W J.m3 Heat flux Position vector (radius vector) Dt fractional order derivative time derivative order Gamma function kg ] m3 Thermal conductivity of the tissue q r ] kg [ 3] m Density of blood [ b J ] kg K W [ 0 ] m K [ W2 ] m wj Weighted arithmetic mean f (t ) Continues function wr Weight function u (t ) Unit step function qext External heat generation [ W3 ] qgen Heat generation in skin qmet Metabolic heating source [ W3 ] div Divergence operator m3 DF Dermic-Fat interface ED FSPL Epidermis-Dermic interface Fractional Single-Phase Lag SPL Single Phase Lag m [ W3 ] m Gradient operator m wb Perfusion rate of blood [ L T Skin thickness [mm ] Tb Temperature of blood [0C ] m3 tissue s ] Temperature of tissue [0C ] T0 DPL Initial temperature [0C ] Time [s ] Time duration of start to finish [s ] t tf t CV FVM Time duration of laser shone on the surface [s ] Time lag Temperature gradient time lag Heat flux time lag T q FDM TDMA SPL Dual Phase Lag Cattaneo Vernet Finite Volume Method Finite Difference Method Tri-Diagonal Matrix Algorithm Single Phase Lag Diffusion reflection Rd Qin Laser intensity [ W ] cm2 Appendix A Finite volume method used for discretization governing equation. First, the heat flux equation (13) is integrated: e t+ t w t q t q t wj = q 1 (2 (j 2 e t+ t q dt dx = t 1+ t D w t+ t t q + D wb x = t D t 2q + D wb x2 b cb T dt dx x (A-1) e b cb T (A-2) w the term is replaced by the following relation (Ghazizadeh et al., 2010): n w j=1 j = = q q x q+ where 1+ q + t ) . 1 2 (j (qin . 1)2 j+1 2qin j + qin j 1 ) 1 t ) (A-3) where is Gamma function and w is Weighted arithmetic mean (Ghazizadeh et al., 2010). Insert equation (A-2) in equation (A-1): 282 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi x qPt+ t = t D qPt + qEt+ t t+ t j=1 q t+ t 2qPt+ t + qW w j (qPt+ + D wb x t j+1 TEt b cb qPt+ t j ) t w j=1 j q (qPt j+1 qPt j ) t TW (A-4) 2 By separating j = 1 from series: x [qPt+ + t qPt + t+ t j=2 q = t D (qPt+ q wj (qPt+ qEt+ t t qPt ) t j+1 t+ t 2qPt+ t + qW qPt+ t j + D wb x (qPt q ) t j=2 q TEt b cb qPt t ) wj (qPt j+1 qPt j ) t TW (A-5) 2 By simplifying the equation, the following equation is obtained: aP qPt+ t = aE qEt+ t + aW qWt+ t +b x q (A-6) where t x aE = D , t x aW = D aP = aE + aW + x + x b=[ x+2 x x ] q t+ t j =2 q q qPt wj (qPt+ qPt t j+1 t qPt+ + D wb t j TEt b cb t j=2 ) t TW t 2 wj (qPt j +1 qPt j ) (A-7) Appendix B Finite difference method used for discretization governing equation. the governing equation is written: T t + + 1+ T t1+ q qext + qmet 1+ q t1 + where T i ci t 1 (1 = ) 1 (1 1 (1 where = w1j + q t t 1 (1 ) (t n=2 d Ti j ) 2Ti j 1 + Ti j 2 t2 a tj 2Ti j 1 + Ti j 2 t2 tj 1 = 1 (1 ) Ti j n j=1 2Ti j 1 + Ti j 2 t2 Ti j n j =1 n j=1 1 t1 + 2Ti j 1 + Ti j 2 t2 (Ti j n j=1 2Ti j (Tin 1 j +1 × (tn j (tn <n dz z ) dz z )1 (tn × + Ti j 2 ) 2Tin 1<1+ n Ti j n j=1 . 1 t2 ) t1 t2 1 . 2T ) 1 (1 1 . n w1 + j=1 j 1+ ) 1 . 1 (j1 (j n w j=1 j = wj = ) 1 (1 = = 1 (t 1 t j )1 (n + Tin + (tn j)1 j 1 ) j1 t j 1)1 + (n 1))1 (j (j 1)1 (B-2) is Gamma function and w is Weighted arithmetic mean (Ghazizadeh et al., 2010). Similar to equation (A-3), = 1+ (B-1) = 1 . ) 1 1 (1 = 1+ q t1 + . ) i ci a = = T) (wb b cb Tb) t q + t wb b cb (Tb i ci + x2 the term is obtained from equation (8) by the following relation (Ghazizadeh et al., 2010): 1+ T t1 + = 2T = Di t (qext + qmet ) q + i ci wb b cb q i ci + 1 (2 (j 2 ) . 1 2 (j (qin 2qin j + qin and q t following relation: j 1 ) 1 . t1 + 1)1 ) (qin j + 1 . j+1 1+ q t1 + 2qin j + qin j 1 ) 1 t 1)2 ) (B-3) and in this problem, we assume that: 283 Journal of Thermal Biology 84 (2019) 274–284 P. Goudarzi and A. Azimi (qext + qmet ) t (wb = 0, b cb Tb) t = 0, qr = 0 (B-4) By simplifying the resulted equation, we have: Di t n T x 2 i+1 Di t Tn x2 i 1 1+ q wb b cb q i ci + ( + 1+2 t i ci (wb Di t x2 = (1 + 2 t n j=2 w1j + t n j=2 b c b Tb + q wb b cb i ci t+ 1+ (Tin j + 1 wj (Tin j+1 t+ w c 2 b cb b q i i 2Tin j + Tin j ) wb b cb q i ci t+ 1+ q t 1+ ) 1+ ) t Tin Tin 1 Tin j 1) q 1+ tTin 2 + qmet + qext ) (B-5) Özen, Ş., Helhel, S., Çerezci, O., 2008. Heat analysis of biological tissue exposed to microwave by using thermal wave model of bio-heat transfer (TWMBT). Burns 34, 45–49. https://doi.org/10.1016/J.BURNS.2007.01.009. Pennes, H.H., 1948. Analysis of tissue and arterial blood temperatures in the resting human forearm. J. Appl. Physiol. 1, 93–122. https://doi.org/10.1152/jappl.1948.1. 2.93. Podlubny, I., 1999. Fractional Differential Equations : an Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Academic Press. Povstenko, Y.Z., 2011. Fractional cattaneo-type equations and generalized thermoelasticity. J. Therm. Stress. 34, 97–114. https://doi.org/10.1080/01495739.2010. 511931. Povstenko, Y., 2015. Fractional Thermoelasticity, Solid Mechanics and its Applications. Springer International Publishing, Cham. https://doi.org/10.1007/978-3-31915335-3. Povstenko, Y., Povstenko, Yuriy, 2013. Fractional heat conduction in an infinite medium with a spherical inclusion. Entropy 15, 4122–4133. https://doi.org/10.3390/ e15104122. Shih, T.-C., Kou, H.-S., Liauh, C.-T., Lin, W.-L., 2005. The impact of thermal wave characteristics on thermal dose distribution during thermal therapy: a numerical study. Med. Phys. 32, 3029–3036. https://doi.org/10.1118/1.2008507. Tzou, D.Y., 1997. In: Macro- to Microscale Heat Transfer : the Lagging Behavior, Published in 1997. Taylor & Francis. Taylor & Francis, Washington DC. Tzou, D.Y., Dai, W., 2009. Thermal lagging in multi-carrier systems. Int. J. Heat Mass Transf. 52, 1206–1213. https://doi.org/10.1016/J.IJHEATMASSTRANSFER.2008. 08.029. Vernotte, P., 1958. Les paradoxes de la theorie continue de l’equation de la chaleur. Compt. Rendus Chem. 246, 3154–3155. Xu, F., Seffen, K.A., Lu, T.J., 2008. Non-Fourier analysis of skin biothermomechanics. Int. J. Heat Mass Transf. 51, 2237–2259. https://doi.org/10.1016/J. IJHEATMASSTRANSFER.2007.10.024. Yu, B., Jiang, X., Wang, C., 2016. Numerical algorithms to estimate relaxation parameters and Caputo fractional derivative for a fractional thermal wave model in spherical composite medium. Appl. Math. Comput. 274, 106–118. https://doi.org/10.1016/J. AMC.2015.10.081. Zhang, W., Cai, X., Holm, S., 2014. Time-fractional heat equations and negative absolute temperatures. Comput. Math. Appl. 67, 164–171. https://doi.org/10.1016/J. CAMWA.2013.11.007. Zhou, J., Chen, J.K., Zhang, Y., 2009. Dual-phase lag effects on thermal damage to biological tissues caused by laser irradiation. Comput. Biol. Med. 39, 286–293. https:// doi.org/10.1016/J.COMPBIOMED.2009.01.002. References Antaki, P.J., 2005. New interpretation of non-fourier heat conduction in processed meat. J. Heat Transf. 127, 189. https://doi.org/10.1115/1.1844540. Cattaneo, C., 1958. A form of heat-conduction equations which eliminates the paradox of instantaneous propagation. Compt. Rendus Chem. 248, 431–433. Compte, A., Metzler, R., 1997. The generalized Cattaneo equation for the description of anomalous transport processes. J. Phys. A Math. Gen. 30, 7277–7289. https://doi. org/10.1088/0305-4470/30/21/006. Djorev, P., Borisova, E.G., Avramov, L.A., 2003. Interaction of the IR laser radiation with human skin: Monte Carlo simulation. In: In: Atanasov, P.A., Serafetinides, A.A., Kolev, I.N. (Eds.), Proceedings of the SPIE, vol. 5226. pp. 403–407. (2003). p. 403. https://doi.org/10.1117/12.519582. Ezzat, M.A., El-bary, A.A., Al-sowayan, N.S., 2016. Tissue responses to fractional transient heating with sinusoidal heat flux condition on skin surface. Anim. Sci. J. 87, 1304–1311. https://doi.org/10.1111/asj.12568. Ghazizadeh, H., 2010. Modeling non-fourier thermal transport based on the theory of fractional calculus. Tarbiat Modares. Ghazizadeh, H.R., Maerefat, M., Azimi, A., 2010. Explicit and implicit finite difference schemes for fractional Cattaneo equation. J. Comput. Phys. 229, 7042–7057. https:// doi.org/10.1016/J.JCP.2010.05.039. Ghazizadeh, H.R., Azimi, A., Maerefat, M., 2012. An inverse problem to estimate relaxation parameter and order of fractionality in fractional single-phase-lag heat equation. Int. J. Heat Mass Transf. 55, 2095–2101. https://doi.org/10.1016/J. IJHEATMASSTRANSFER.2011.12.012. Hooshmand, P., Moradi, A., Khezry, B., 2015. Bioheat transfer analysis of biological tissues induced by laser irradiation. Int. J. Therm. Sci. 90, 214–223. https://doi.org/10. 1016/J.IJTHERMALSCI.2014.12.004. Kuo-Chi, L., Po-Jen, C., Yan-Nan, W., 2011. Analysis of non-Fourier thermal behavior for multi-layer skin model. Therm. Sci. 15, 61–67. https://doi.org/10.2298/ TSCI11S1061L. Lin, S.-M., Li, C.-Y., 2016. Analytical solutions of non-Fourier bio-heat conductions for skin subjected to pulsed laser heating. Int. J. Therm. Sci. 110, 146–158. https://doi. org/10.1016/J.IJTHERMALSCI.2016.06.034. Liu, K.-C., 2008. Thermal propagation analysis for living tissue with surface heating. Int. J. Therm. Sci. 47, 507–513. https://doi.org/10.1016/J.IJTHERMALSCI.2007.04.005. Liu, Jing, Chen, Xu, Xu, L.X., 1999. New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating. IEEE Trans. Biomed. Eng. 46, 420–428. https:// doi.org/10.1109/10.752939. Liu, K.-C., Wang, Y.-N., Chen, Y.-S., 2012. Investigation on the bio-heat transfer with the dual-phase-lag effect. Int. J. Therm. Sci. 58, 29–35. https://doi.org/10.1016/J. IJTHERMALSCI.2012.02.026. 284