Introduction

The thermal transfer through a living skin tissue is extremely essential to a wide range of therapies and their applications1. Improvements in the microwave, laser, and other technologies also extended bioheat transfer work. Several researchers comment on biothermal transmission analyzes in biological living tissues, including Pennes2. He optimized the parabolic model for biological tissues for the first bioheat transfer. Also, the Pennes biological transmission method is used to model the nature and consistency of thermal behavior in biological tissues and living bodies. External similarities to such unpredictable results that display non-Fourier or hyperbolic conduction behavior. Cattaneo3 and Vernotte4 provided different modulations of Fourier's law of heat conduction as a linear extension form of the Fourier law to describe the hyperbolic equation type. They suggested a model of thermal waves measure the effect of microwave and thermal flux activity. Different forms of work have been undertaken to address various skin tissue conditions without damaging the balanced tissue surrounding them. Xu et al.5 solved analytical the Pennes bioheat transfer equation (PBTE), reviewed skin bioheat transfer, skin structure, thermal damage, skin biomechanics, and bio-thermomechanics. The non-Fourier thermomechanical activity of skin tissues under various surface thermal loading limits was analyzed in DPL, hyperbolic and parabolic models of biomass transport, Xu and al. studied and noticed substantial variations between Pennes, thermal wave and DPL anticipations models6. Also, Rossmann et al.7 reviewed the temperature dependence of the dielectric properties, thermal properties, and perfusion of biological tissues at hyper-thermic and ablation temperatures. Poor et al.8 focused on the temperature disturbance of skin tissue due to time-dependent surface heating. Tzou's model has been extended by the dual-phase-lag principle, with the delayed activity in a high rate of response taken into consideration. Although the process is lagging, the small-scale response is caught in time9,10. Tzou11 introduced a phase-lag for a temperature gradient. Askarizadeh et al.12 utilized the dual-phase-lag (DPL) model in treating the transient heat transfer problems in skin tissue. Dutta and Kundu established a two-dimensional thermal model of malignant tissues, focused on the bi-dimensional local thermal non-balance therapy model for bio-heat13.

Liu et al. first address the overall bioheat transmission paradigm in living tissue and analyze14. Liu et al. are tested during the hyperthermia procedure in the two-stage model for heat transfer problems in the biological tissue15,16. The Liu and Chen thermal conductor model (DPL) used the non-Fourier thermal activity for the diagnosis of hyperthermia to describe the thermal transportation occurring in biological tissue17. In terms of the properties of blood and tissue and the interphase charge of the heat transfer parameter and perfusion rate, Zhang articulated phase lag or relaxation periods18. He observed that the lap periods are very similar together for living tissues. Dutta and Kundu performed a thermal wave propagation study for a constant and variable flux of heat on a skin tissue surface in biological tissue, connected with the diagnosis of hyperthermia.19.

A finite element approach to the temperature distribution of deep skin tissues in the human limb was suggested by Agrawal et al.20. Kumar et al. also researched the thermal behavior, using a dual-phase-lag (DPL) model of time-fractional21. To address the effect of injecting nanoparticles in irregularly-shaped liver tumors on cancer cells, Shao et al. have developed an RFA model22. A dual-phase-lag (DPL) bioheat transfer model analysis approach was used by Poor et al. on the skin tissue as the finite domain of continuous, cosine, and pulse heat flow conditions on the surface of the skin8.

Liu and Xu developed a closed-form analytical approach for the Pennes bioheat equation solution for changes in skin tissue temperature induced by sinusoidal heat flow23. Shih et al. have overcome the consequence of a semi-infinite biological skin tissue temperature reaction study undergoing sine flux24. Regarding periodic constant and pulsed train heat flux boundary conditions, Ahmadika et al. shall have the analytical solution for the bioheat (parabolize) and hyperbolized bioheat transfer models of the Penne25. The statistical model was used by Horng et al. to explain the effects of pulsatile blood flow in the thermal disruption26. He found out that the influence on the thermal area on tumor and the slight variation only on a blood continuum range from standard to parabolic, velocity profiles on the blood flow. Shih et al. investigated the coupled effects of the pulsatile blood flow and thermal relaxation time in living tissues during thermal treatments27. Youssef and Alghamdi28 have solved a one-dimensional application of thermoelastic dual-phase-lag skin tissue under specific thermal loading. Kundu and Dewanjee29 introduced a new method for non-Fourier thermal response in a single layer skin tissue.

For one calming period in deformable structures, Youssef updated Chen and Gurtin's principle of thermal actions. Two specific temperatures are the thermodynamic temperature and the conductive temperature in this device30,31. The disparity in the meaning of these two temperature forms is commensurate with the importance of the material's thermal supply. Youssef with El-Bary and Alghamdi applied the two-temperature heat conduction in many applications32,33,34,35,36,37.

In the present paper, a mathematical model of bioheat skin tissue will be constructed in the context of a two-temperature dual-phase-lag heat conduction model. The surface of the skin tissue is subjected to a constant heat flux. The separation of variables as a finite domain will be applied to obtain the exact analytical solution. The effects of the dual-phase-lag time parameter, value of the heat flux, and two-temperature parameter will be studied and discussed.

Bioheat transfer models for biological tissues

The bioheat transmitting model was developed to measure the time-dependent temperature change as a heat reaction due to the usage of either heat source or thermal heating. The first model of biological tissues was found out by Pennes based on Fourier's law of heat conduction as follows28,36,37:

$$ K\nabla^{2} T = \rho C\frac{\partial T}{{\partial t}} - \rho_{b} w_{b} C_{b} \left( {T_{b} - T} \right) - \left( {Q_{met} + Q_{ext} } \right), $$
(1)

where \(\rho ,\,C\) and \(K\) are the density, specific heat, and thermal conductivity of the tissue, respectively.\(C_{b} ,\,w_{b} \,\), \(\rho_{b}\) and \(T_{b}\) denote the specific heat of the blood, blood perfusion rate, blood density, and blood temperature, respectively. T is the absolute temperature function. \(Q_{met}\) denotes the metabolic heat generated by the chemical reaction inside the tissue, and it is constant, and the external heat source is given by \(Q_{ext} = Q_{ext} \left( t \right)\). \(\nabla^{2}\) is the well-known Laplace operator.

The modified Fourier's law of heat conduction

Ventott and Cattaneo modified the classical Fourier thermal drive law by postulating the concept of finite thermal wave propagation speed and taking the following hyperbolic form of the thermal wave28,36,37:

$$ K\nabla^{2} T = \rho C\left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\frac{\partial T}{{\partial t}} + w_{b} C_{b} \rho_{b} \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\left( {T - T_{b} } \right) - \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\left( {Q_{met} + Q_{ext} } \right), $$
(2)

where \(\tau_{q} = \frac{\alpha }{{c_{0}^{2} }} > 0\) is a material property and is called the relaxation time parameter, and \(\alpha\) is the thermal diffusivity while \(c_{0}\) is the speed of the thermal wave inside the medium.

Dual-phase-lag model (DPL) of bioheat transfer

The dual-phase-lag (DPL) model based on the dual reaction between the gradient of the temperature \(\nabla T\) and the heat flux \(q\) which modified the well-known classical Fourier's law of heat conduction, thus, we have the following heat conduction equation28,36,37:

$$ K\left( {1 + \tau_{T} \frac{\partial }{\partial t}} \right)\nabla^{2} T = \rho C\left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\frac{\partial T}{{\partial t}} + w_{b} C_{b} \rho_{b} \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\left( {T - T_{b} } \right) - \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\left( {Q_{met} + Q_{ext} } \right) $$
(3)

where \(\tau_{T} \ge 0\) is the second parameter of the relaxation time, which is the phase-lag of the temperature gradient passing through the medium.

Two-temperature dual-phase-lag (TTDPL) bioheat transfer model on skin tissues

According to Youssef's model of heat conduction formulation, we have two equations of heat conduction on the biological tissues as following36,37 :

$$ K\left( {1 + \tau_{T} \frac{\partial }{\partial t}} \right)\nabla^{2} T_{C} = \rho C\left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\frac{{\partial T_{D} }}{\partial t} + w_{b} C_{b} \rho_{b} \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)T_{D} - \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\left( {Q_{met} + Q_{ext} } \right) $$
(4)

and

$$ T_{D} = T_{C} - \beta \nabla^{2} T_{C} $$
(5)

where \(\beta\) is a non-negative constant, which is called the two-temperature parameter36,37. \(T_{D}\) and \(T_{C}\) are the dynamical temperature and the conductive temperature, respectively. The value \(\beta = 0\) represents the one-temperature model, then, we obtain \(T_{C} = T_{D}\).

We consider the function of dynamical temperature increment takes the form

$$ \theta = T_{D} - T_{b} $$
(6)

and the function of conductive temperature increment takes the form

$$ \varphi = T_{C} - T_{b} $$
(7)

Hence, we have

$$ K\left( {1 + \tau_{T} \frac{\partial }{\partial t}} \right)\nabla^{2} \varphi = \rho C\left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\frac{\partial \theta }{{\partial t}} + w_{b} C_{b} \rho_{b} \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\theta - \left( {1 + \tau_{q} \frac{\partial }{\partial t}} \right)\left( {Q_{met} + Q_{ext} } \right) $$
(8)

and

$$ \theta = \varphi - \beta \,\nabla^{2} \varphi $$
(9)

Method

We consider the region \(0 \le x \le L\) is filled with skin tissue as in Fig. 1 and obeys the two-temperature dual-phase-lag (TTDPL) as in Eqs. (8) and (9).

Figure 1
figure 1

The skin tissue affected by constant heat flux on the surface.

The medium is considered to be quiescent initially and without any external heat source \(\left( {Q_{ext} = 0} \right)\) while \(Q_{met}\) is constant. The skin is subjected to a constant heat flux \(q_{0}\) on its bounding surface \(x = 0\), while the bounding surface \(x = L\) has no heat flux (see Fig. 1).

Hence, the heat conduction equations take the forms:

$$ \left[ {\beta \rho C\tau_{q} \frac{{\partial^{2} }}{{\partial t^{2} }} + \left( {K\tau_{T} + \beta \rho C + \beta \tau_{q} w_{b} C_{b} \rho_{b} } \right)\frac{\partial }{\partial t} + \left( {\beta w_{b} C_{b} \rho_{b} + K} \right)} \right]\frac{{\partial^{2} \varphi }}{{\partial x^{2} }} = \left( {\tau_{q} \rho C\frac{{\partial^{2} }}{{\partial t^{2} }} + \left( {\rho C + \tau_{q} w_{b} C_{b} \rho_{b} } \right)\frac{\partial }{\partial t} + w_{b} C_{b} \rho_{b} } \right)\varphi - Q_{met} $$
(10)

and

$$ \theta = \varphi - \beta \frac{{\partial^{2} \varphi }}{{\partial x^{2} }} $$
(11)

The initial conditions are:

$$ \left. {\varphi \left( {x,t} \right)} \right|_{t = 0} \left. { = \theta \left( {x,t} \right)} \right|_{t = 0} = \left. {\frac{{\partial \varphi \left( {x,t} \right)}}{\partial t}} \right|_{t = 0} = \left. {\frac{{\partial \theta \left( {x,t} \right)}}{\partial t}} \right|_{t = 0} = 0 $$
(12)

The boundary conditions are:

$$ \left. { - K\frac{{\partial \varphi \left( {x,t} \right)}}{\partial \,x}} \right|_{x = 0} = q_{0} \quad \left. {\frac{{\partial \varphi \left( {x,t} \right)}}{\partial \,x}} \right|_{x = L} = 0 $$
(13)

The boundary value problem (10)–(13) consists of non-homogeneous partial differential equations with non-homogeneous boundary conditions on the surface of the skin tissue. Thus, the equations must be formulated in terms of a steady part and a transient part as follows38,39:

$$ \varphi \left( {x,t} \right) = \phi_{1} \left( {x,t} \right) + \phi_{2} \left( x \right) $$
(14)

Hence, we obtain the transient part in the form

$$ \left[ {\beta \rho C\tau_{q} \frac{{\partial^{2} }}{{\partial t^{2} }} + \left( {K\tau_{T} + \beta \rho C\, + \beta \tau_{q} w_{b} C_{b} \rho_{b} } \right)\frac{\partial }{\partial t} + \left( {\beta w_{b} C_{b} \rho_{b} + K} \right)} \right]\frac{{\partial^{2} \varphi_{1} \left( {x,t} \right)}}{{\partial x^{2} }} - \left( {\tau_{q} \rho C\frac{{\partial^{2} }}{{\partial t^{2} }} + \left( {\rho C + \tau_{q} w_{b} C_{b} \rho_{b} } \right)\frac{\partial }{\partial t} + w_{b} C_{b} \rho_{b} } \right)\varphi_{1} \left( {x,t} \right) = 0 $$
(15)

and the steady-state part in the form

$$ \left( {\frac{{d^{2} }}{{dx^{2} }} - \lambda^{2} } \right)\phi_{2} \left( x \right) = - \gamma \psi $$
(16)

where \(\lambda^{2} = \frac{{w_{b} C_{b} \rho_{b} }}{{\left( {\beta w_{b} C_{b} \rho_{b} + K} \right)}} > 0,\quad \gamma = \frac{1}{{\left( {\beta w_{b} C_{b} \rho_{b} + K} \right)}} > 0\), and \(\psi = Q_{met}\).

The boundary conditions for the steady-state equation are

$$ - K\left. {\frac{{d\phi_{2} \left( x \right)}}{dx}} \right|_{x = 0} = q_{0} ,\quad - K\left. {\frac{{d\phi_{2} \left( x \right)}}{dx}} \right|_{x = L} = 0 $$
(17)

Then, regarding the boundary conditions on (13), the solution of the Eq. (10) is in the form (See “Appendix”):

$$ \phi_{2} \left( x \right) = \frac{{q_{0} \cosh \lambda \left( {L - x} \right)}}{\lambda K\sinh \lambda L} + \frac{\gamma }{{\lambda^{2} }}\psi $$
(18)

The initial and boundary conditions of the transient Eq. (12) are:

$$ \left. {\phi_{1} \left( {x,t} \right)} \right|_{t = 0} = - \phi_{2} \left( x \right),\quad \left. {\frac{{\partial \phi_{1} \left( {x,t} \right)}}{\partial t}} \right|_{t = 0} = 0 $$
(19)

and

$$ \left. {\frac{{\partial \phi_{1} \left( {x,t} \right)}}{\partial x}} \right|_{x = 0} = 0,\quad \left. {\frac{{\partial \phi_{1} \left( {x,t} \right)}}{\partial x}} \right|_{x = L} = 0 $$
(20)

To solve the Eq. (15) we expand the function \(\phi_{1} \left( {x,t} \right)\) in the following Fourier series expansion38,39

$$ \phi_{1} \left( {x,t} \right) = \sum\limits_{n = 0}^{\infty } {\vartheta_{n} \left( t \right)\cos \left( {\frac{n\pi }{L}x} \right)} $$
(21)

which satisfies the boundary conditions (20).

Substitute from Eqs. (21) into (15), we get

$$ \left[ {\frac{{\partial^{2} }}{{\partial t^{2} }} + A_{1n} \frac{\partial }{\partial t} + A_{2n} } \right]\vartheta_{n} \left( t \right) = 0,\quad n = 0,1,2,3 \ldots $$
(22)

where \(A_{1n} = \frac{{\omega_{n}^{2} \left( {\eta \tau_{T} + \beta + \beta \tau_{q} \varepsilon } \right) + \left( {1 + \tau_{q} \varepsilon } \right)}}{{\tau_{q} \left( {\beta \omega_{n}^{2} + 1} \right)}},\;\;A_{2n} = \frac{{\omega_{n}^{2} \left( {\beta \varepsilon + \eta } \right) + \varepsilon }}{{\tau_{q} \left( {\beta \omega_{n}^{2} + 1} \right)}},\;\;\varepsilon = \frac{{w_{b} C_{b} \rho_{b} }}{\rho C},\;\;\eta = \frac{K}{\rho C}\) and \(\omega_{n} = \frac{n\pi }{L}\).

The general solutions to (22) take the following forms38,39:

$$ \vartheta_{n} \left( t \right) = a_{n} f_{1} \left( {k_{1n} ,t} \right) + b_{n} f_{2} \left( {k_{2n} ,t} \right),\quad n = 0,1,2, \ldots $$
(23)

which give

$$ \vartheta_{1} \left( {x,t} \right) = \sum\limits_{n = 0}^{\infty } {\left[ {a_{n} f_{1} \left( {k_{1n} ,t} \right) + b_{n} f_{2} \left( {k_{2n} ,t} \right)} \right]} \cos \left( {\omega_{n} x} \right) $$
(24)

where \(k_{1n} ,k_{2n}\) is the solution of the following characteristic equation

$$ k^{2} + A_{1n} k + A_{2n} = 0 $$
(25)

The roots of the characteristic Eq. (25) are in the forms:

$$ k_{1n} = \frac{{ - A_{1n} + \sqrt {\Delta_{n} } }}{2},\quad k_{2n} = \frac{{ - A_{1n} - \sqrt {\Delta_{n} } }}{2} $$
(26)

where \(\,\Delta_{n} = A_{1n}^{2} - 4A_{2n}\).

To apply the initial conditions (19), we have to expand \(\phi_{2} \left( x \right)\) in Fourier series expansion to be in the form:

$$ \phi_{2} \left( x \right) = \frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }} + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{1}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) $$
(27)

Hence, the initial conditions (19) give the following system of algebraic equations:

$$ \left[ {a_{0} f_{1} \left( {k_{10} ,0} \right) + b_{0} f_{2} \left( {k_{20} ,0} \right)} \right]\cos \left( {\omega_{0} x} \right) + \sum\limits_{n = 1}^{\infty } {\left[ {a_{n} f_{1} \left( {k_{1n} ,0} \right) + b_{n} f_{2} \left( {k_{2n} ,0} \right)} \right]} \cos \left( {\omega_{n} x} \right) = - \phi_{2} \left( x \right) $$
(28)

and

$$ \left[ {a_{0} \frac{{\partial f_{1} \left( {k_{10} ,0} \right)}}{\partial t} + b_{0} \frac{{\partial f_{2} \left( {k_{20} ,0} \right)}}{\partial t}} \right]\cos \left( {\omega_{0} x} \right) + \sum\limits_{n = 1}^{\infty } {\left[ {a_{n} \frac{{\partial f_{1} \left( {k_{1n} ,0} \right)}}{\partial t} + b_{n} \frac{{\partial f_{2} \left( {k_{2n} ,0} \right)}}{\partial t}} \right]} \cos \left( {\omega_{n} x} \right) = 0 $$
(29)

For the case of \(\Delta_{n} > 0\)

$$ f_{1} \left( {k_{1n} ,t} \right) = {\text{e}}^{{k_{1n} \,t}} \;\; {\text{and}}\;\;f_{2} \left( {k_{2n} ,t} \right) = {\text{e}}^{{k_{2n} \,t}} $$
(30)

Thus, when \(t = 0\) we have

$$ f_{1} \left( {k_{1n} ,0} \right) = f_{2} \left( {k_{2n} ,0} \right) = 1,\quad \frac{{\partial f_{1} \left( {k_{1n} ,0} \right)}}{\partial t} = k_{1n} ,\quad \frac{{\partial f_{2} \left( {k_{2n} ,0} \right)}}{\partial t} = k_{2n} ,\quad n = 0,1,2,3 \ldots $$

Then, Eqs. (28) and (29) introduce the following system

$$ \left[ {a_{0} + b_{0} } \right]\cos \left( {\omega_{0} x} \right) + \sum\limits_{n = 1}^{\infty } {\left[ {a_{n} + b_{n} } \right]} \cos \left( {\omega_{n} x} \right) = - \frac{{q_{0} }}{{2LK\lambda^{2} }} - \frac{\gamma \psi }{{2\lambda^{2} }} - \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{1}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) $$
(31)

and

$$ \left[ {k_{10} a_{0} + k_{20} b_{0} } \right]\cos \left( {\omega_{0} x} \right) + \sum\limits_{n = 1}^{\infty } {\left[ {k_{1n} a_{n} + k_{2n} b_{n} } \right]} \cos \left( {\omega_{n} x} \right) = 0 $$
(32)

when \(n = 0\) we have

$$ \omega_{0} = 0,\;\cos \left( {\omega_{0} x} \right) = 1,\;A_{10} = \frac{{\varepsilon \tau_{q} + 1}}{{\tau_{q} }},\,\,\,A_{20} = \frac{\varepsilon }{{\tau_{q} }},\,\Delta_{0} = \left( {\frac{{\varepsilon \tau_{q} - 1}}{{\tau_{q} }}} \right)^{2} ,\,k_{10} = \frac{ - 1}{{\tau_{q} }},\;k_{20} = - \varepsilon . $$

Hence, the system in Eqs. (31) and (32) will be reduced to the following system

$$ a_{0} + b_{0} + \sum\limits_{n = 1}^{\infty } {\left[ {a_{n} + b_{n} } \right]} \cos \left( {\omega_{n} x} \right) = - \frac{{q_{0} }}{{2LK\lambda^{2} }} - \frac{\gamma \psi }{{2\lambda^{2} }} - \frac{{q_{0} }}{LK}\sum\limits_{n = 1}^{\infty } {\frac{1}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) $$
(33)

and

$$ - \frac{{a_{0} }}{{\tau_{q} }} - \varepsilon b_{0} + \sum\limits_{n = 1}^{\infty } {\left[ {k_{1n} a_{n} + k_{2n} b_{n} } \right]} \cos \left( {\omega_{n} x} \right) = 0 $$
(34)

The equations in (33) and (34) lead to the following two systems of algebraic equations

$$ a_{0} + b_{0} = - \left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right),\quad a_{0} + \tau_{q} \varepsilon b_{0} = 0 $$
(35)

and

$$ a_{n} + b_{n} = \frac{{ - q_{0} }}{{LK\left( {\lambda^{2} + \omega_{n}^{2} } \right)}},\quad k_{1n} a_{n} + k_{2n} b_{n} = 0,\quad n = 1,2,3, \ldots $$
(36)

By solving the above two systems, we get

$$ a_{0} = \left( {\frac{{\varepsilon \tau_{q} }}{{1 - \varepsilon \tau_{q} }}} \right)\left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right),\quad b_{0} = \left( {\frac{ - 1}{{1 - \varepsilon \tau_{q} }}} \right)\left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right) $$
(37)

and

$$ a_{n} = \left( {\frac{{k_{2n} }}{{k_{1n} - k_{2n} }}} \right)\frac{{q_{0} }}{{LK\left( {\lambda^{2} + \omega_{n}^{2} } \right)}},\quad b_{n} = \left( {\frac{{ - k_{1n} }}{{k_{1n} - k_{2n} }}} \right)\frac{{q_{0} }}{{LK\left( {\lambda^{2} + \omega_{n}^{2} } \right)}},\quad n = 1,2,3 \ldots $$
(38)

Now, the final solution of the heat conduction temperature increment, in this case, is in the form

$$ \varphi \left( {x,t} \right) = \frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }} + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{1}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) + \left( {\frac{1}{{1 - \varepsilon \tau_{q} }}} \right)\left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right)\left[ {\varepsilon \tau_{q} {\text{e}}^{{\frac{ - t}{{\tau_{q} }}}} - {\text{e}}^{ - \varepsilon t} } \right] + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\left( {\frac{{k_{2n} {\text{e}}^{{k_{1n} \,t}} - k_{1n} {\text{e}}^{{k_{2n} \,t}} }}{{\left( {k_{1n} - k_{2n} } \right)\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \right)} \cos \left( {\omega_{n} x} \right) $$
(39)

By using Eqs. (39) and (11), we get the dynamical temperature increment as follows:

$$ \theta \left( {x,t} \right) = \frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }} + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{{\left( {1 + \beta \omega_{n}^{2} } \right)}}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) + \left( {\frac{1}{{1 - \varepsilon \tau_{q} }}} \right)\left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right)\left[ {\varepsilon \tau_{q} {\text{e}}^{{\frac{ - t}{{\tau_{q} }}}} - {\text{e}}^{ - \varepsilon t} } \right] + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\left( {1 - \beta \omega_{n}^{2} } \right)\left( {\frac{{k_{2n} {\text{e}}^{{k_{1n} \,t}} + k_{1n} {\text{e}}^{{k_{2n} \,t}} }}{{\left( {k_{1n} - k_{2n} } \right)\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \right)} \cos \left( {\omega_{n} x} \right) $$
(40)

For the case of \(\Delta < 0\), we have

$$ f_{1} \left( {k_{1n} ,t} \right) = {\text{e}}^{{\left( {k_{1n} \,t} \right)}} \cos \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right),\quad f_{2} \left( {k_{2n} ,t} \right) = {\text{e}}^{{\left( {k_{2n} \,t} \right)}} \sin \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right), $$
(41)
$$ \frac{{\partial f_{1} \left( {k_{1n} ,t} \right)}}{\partial \,t} = {\text{e}}^{{\left( {k_{1n} \,t} \right)}} \left[ {k_{1n} \,\cos \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right) - \frac{{\sqrt { - \Delta } }}{2}\sin \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right)} \right] $$
(42)
$$ \frac{{\partial f_{2} \left( {k_{2n} ,t} \right)}}{\partial t} = {\text{e}}^{{\left( {k_{2n} \,t} \right)}} \left[ {k_{2n} \sin \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right) + \frac{{\sqrt { - \Delta_{n} } }}{2}\cos \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right)} \right] $$
(43)
$$ f_{1} \left( {k_{1n} ,0} \right) = 1,\quad f_{2} \left( {k_{2n} ,0} \right) = 0,\quad n = 0,1,2,3 \ldots $$
(44)
$$ \frac{{\partial f_{1} \left( {k_{1n} ,0} \right)}}{\partial t} = k_{1n} ,\frac{{\partial f_{2} \left( {k_{2n} ,0} \right)}}{\partial t} = \frac{{\sqrt { - \Delta_{n} } }}{2},\quad n = 0,1,2,3 \ldots $$
(45)
$$ A_{10} = \frac{1}{{\tau_{q} }},\;A_{20} = \frac{\varepsilon }{{\tau_{q} }},\;\Delta_{0} = \left( {\frac{{\varepsilon \tau_{q} - 1}}{{\tau_{q} }}} \right)^{2} ,\,k_{10} = \frac{ - 1}{{\tau_{q} }},\;k_{20} = - \varepsilon $$

Applying the initial conditions (19), we obtain

$$ a_{0} + \sum\limits_{n = 1}^{\infty } {a_{n} } \cos \left( {\omega_{n} x} \right) = - \frac{{q_{0} }}{{2LK\lambda^{2} }} - \frac{\gamma \psi }{{2\lambda^{2} }} - \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{1}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) $$
(46)

and

$$ \left( { - \frac{1}{{\tau_{q} }}} \right)a_{0} + \left( {\frac{{\tau_{q} - 1}}{{2\tau_{q} }}} \right)b_{0} + \sum\limits_{n = 1}^{\infty } {\left[ {a_{n} k_{1n} + b_{n} \frac{{\sqrt { - \Delta_{n} } }}{2}} \right]} \cos \left( {\omega_{n} x} \right) = 0 $$
(47)

Solving the above system, we obtain

$$ a_{0} = - \left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right),\quad b_{0} = - \left( {\frac{1}{{\tau_{q} - 1}}} \right)\left( {\frac{{q_{0} }}{{LK\lambda^{2} }} + \frac{\gamma \psi }{{\lambda^{2} }}} \right) $$
(48)
$$ a_{n} = - \frac{{q_{0} }}{{KL\left( {\lambda^{2} + \omega_{n}^{2} } \right)}},\quad b_{n} = \frac{{2k_{1n} q_{o} }}{{KL\left( {\lambda^{2} + \omega_{n}^{2} } \right)\sqrt { - \Delta_{n} } }} $$
(49)

Hence, the final solution of the heat conduction temperature increment in this case is:

$$ \begin{aligned} & \varphi \left( {x,t} \right) = \frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\psi }{2} + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{1}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) - \left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right){\text{e}}^{{\frac{ - t}{{\tau_{q} }}}} \cos \left( {\frac{{\sqrt { - \Delta_{0} } }}{2}\,t} \right) - \left( {\frac{1}{{\tau_{q} - 1}}} \right)\left( {\frac{{q_{0} }}{{LK\lambda^{2} }} + \frac{\gamma \psi }{{\lambda^{2} }}} \right){\text{e}}^{ - t} \sin \left( {\frac{{\sqrt { - \Delta_{0} } }}{2}\,t} \right) \\ & \quad \quad \quad \quad + \frac{{q_{0} }}{KL}\,\,\sum\limits_{n = 1}^{\infty } {\left[ {{\text{e}}^{{\left( {k_{1n} \,t} \right)}} \cos \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right) + \frac{{2k_{1n} }}{{\sqrt { - \Delta_{n} } }}{\text{e}}^{{\left( {k_{2n} \,t} \right)}} \sin \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right)} \right]} \frac{{\cos \left( {\omega_{n} x} \right)}}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}} \\ \end{aligned} $$
(50)

By using Eqs. (50) and (11), we get the dynamical temperature increment in this case as follows:

$$ \begin{aligned} & \theta \left( {x,t} \right) = \frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\psi }{2} + \frac{{q_{0} }}{KL}\sum\limits_{n = 1}^{\infty } {\frac{{\left( {1 + \omega_{n}^{2} \beta } \right)}}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}}} \cos \left( {\omega_{n} x} \right) - \left( {\frac{{q_{0} }}{{2LK\lambda^{2} }} + \frac{\gamma \psi }{{2\lambda^{2} }}} \right){\text{e}}^{{\frac{ - t}{{\tau_{q} }}}} \cos \left( {\frac{{\sqrt { - \Delta_{0} } }}{2}\,t} \right) - \left( {\frac{1}{{\tau_{q} - 1}}} \right)\left( {\frac{{q_{0} }}{{LK\lambda^{2} }} + \frac{\gamma \psi }{{\lambda^{2} }}} \right){\text{e}}^{ - t} \sin \left( {\frac{{\sqrt { - \Delta_{0} } }}{2}\,t} \right) \\ & \quad \quad \quad \quad + \frac{{q_{0} }}{KL}\,\,\sum\limits_{n = 1}^{\infty } {\left[ {{\text{e}}^{{\left( {k_{1n} \,t} \right)}} \cos \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}\,t} \right) + \frac{{2k_{1n} }}{{\sqrt { - \Delta_{n} } }}{\text{e}}^{{\left( {k_{2n} \,t} \right)}} \sin \left( {\frac{{\sqrt { - \Delta_{n} } }}{2}t} \right)} \right]} \frac{{\left( {1 + \omega_{n}^{2} \beta } \right)\cos \left( {\omega_{n} x} \right)}}{{\left( {\lambda^{2} + \omega_{n}^{2} } \right)}} \\ \end{aligned} $$
(51)

Results

In this study, the temperature distribution through skin tissue is investigated for two models (One temperature, two temperature) of bioheat transfer for constant heat flux condition on the skin surface. The values of the relevant thermal parameters used in the present calculations are in Table 1 as follows:

Table 1 The material properties of the skin tissue.

Discussions

Figures 2, 3, 4, 5 represent the conductive and dynamical temperature increment with respect to a wide range of skin distance x \((0 \le x \le 0.006\,\,{\text{m}})\) and at constant instance time \(t = 100({\text{s}})\).

Figure 2
figure 2

The temperature increment when \(t = 100\,({\text{s}}),\,\,\tau_{T} = 10\,({\text{s}}),\,\,\tau_{q} = 20\,({\text{s}})\) and \(q_{0} = 100\;({\text{W/m}}^{{2}} )\).

Figure 3
figure 3

The temperature increment when \(t = 100\,\left( {\text{s}} \right),\,\,\tau_{T} = 10\,\left( {\text{s}} \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and various values of the surface heat flux \(q_{0} \,\left( {{\text{W/m}}^{{2}} } \right)\) based on the two-temperature model.

Figure 4
figure 4

The temperature increment when \(t = 100\,\left( {\text{s}} \right),\,\,q_{0} = 100\,\,\left( {{\text{W/m}}^{{2}} } \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and various values \(\tau_{T} \left( {\text{s}} \right)\) based on the two-temperature model.

Figure 5
figure 5

The temperature increment when \(t = 100\,\left( {\text{s}} \right),\,\,q_{0} = 100\,\,\left( {{\text{W/m}}^{{2}} } \right),\,\,\tau_{T} = 10\,\left( {\text{s}} \right)\) and various values \(\tau_{q} \left( {\text{s}} \right)\) based on the two-temperature model.

Figure 2 shows that the two temperature parameter \(\beta\) has a significant effect on the temperature increment distribution. The red line represents the conductive temperature distribution, while the blue line represents the dynamical temperature distribution. The conductive temperature increment and dynamical temperature increment have the same behavior but with different values because of the effect of the two-temperature parameter. In the context of the one-temperature model \(\beta = 0\), the conductive and dynamical temperature increments are the same, and they appear in the solid blue curve. It is noted that an increase in the value of the two-temperature parameter leads to an increase in the value of conductive temperature increment and a decrease in the value of dynamical temperature increment. We noticed that thermal waves in the context of the one-temperature model have vanished before the thermal waves in the context of the two-temperature model.

Figure 3 shows that the conductive and dynamical temperature increment distributions when \(t = 100\;({\text{s)}},\,\,\tau_{T} = 10\;({\text{s}}),\,\,\tau_{q} = 20\;({\text{s}})\) and various values of the surface heat flux \(q_{0} \,({\text{W/m}}^{{2}} ) = \left( {100,200} \right)\) to stand on its effects. It is noted that the value of the heat flux on the surface of the skin tissue has a significant impact on the conductive and dynamical temperature increment distributions. An increase in the value of the heat flux leads to an increase in the value of the conductive and dynamical temperature increment distributions. A smaller value of the surface heat flux makes the temperature more closed to zero value at the end of the skin length.

Figure 4 shows that the conductive and dynamical temperature increment distributions when \(t = 100\;({\text{s}}),\,q_{0} = 100\;({\text{W/m}}^{{2}} ),\,\,\tau_{q} = 20\;({\text{s)}}\) and various values of the second relaxation time parameter \(\tau_{T} ({\text{s}}) = (0.0,10.0)\) to stand on its effects. It is noted that the value of the second relaxation time parameter has a significant impact on the conductive and dynamical temperature increment distributions. An increase in the value of the second relaxation time parameter leads to an increase in the value of the conductive and dynamical temperature increment distributions. The zero value of the first relaxation time parameter \(\tau_{T}\) (temperature lag-time) makes the thermal waves more closed to the zero value at the end of the skin length.

Figure 5 shows that the conductive and dynamical temperature increment distributions when \(t = 100\,{\text{(s}}),\,\,q_{0} = 100\,\,({\text{W/m}}^{{2}} ),\,\,\tau_{T} = 10\,({\text{s}})\) and various values of the first relaxation time parameter \(\tau_{q} \left( {\text{s}} \right) = \left( {10,20} \right)\) to stand on its effects. It is noted that the value of the first relaxation time \(\tau_{q}\) parameter has noticeable effects on the conductive and dynamical temperature increment distributions. An increase in the value of the second relaxation time parameter \(\tau_{q}\) leads to a decrease in the value of the conductive and dynamical temperature increment distributions. Increasing the value of the second relaxation time parameter \(\tau_{q}\)(temperature gradient lag-time) makes the thermal waves more closed to the zero value at the end of the skin length.

Figures 6, 7, 8, 9 represent the conductive and dynamical temperature increment with respect to a wide range of time x \(\left( {0 \le t \le 100\,\,{\text{s}}} \right)\) and at a constant distance \(x = 0.003\left( {\text{m}} \right)\).

Figure 6
figure 6

The temperature increment when \(x = 0.003\,\left( {\text{m}} \right),\,\,\tau_{T} = 10\,\left( {\text{s}} \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and \(q_{0} = 100\,\,\left( {W/m^{2} } \right)\).

Figure 7
figure 7

The temperature increment when \(x = 0.003\,\left( {\text{m}} \right),\,\,\tau_{T} = 10\,\left( {\text{s}} \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and various values of the surface heat flux \(q_{0} \,\left( {{\text{W/m}}^{{2}} } \right)\) based on the two-temperature model.

Figure 8
figure 8

The temperature increment when \(x = 0.003\,\left( {\text{m}} \right),\,\,q_{0} = 100\,\,\left( {{\text{W/m}}^{{2}} } \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and various values \(\tau_{T} \left( {\text{s}} \right)\) based on the two-temperature model.

Figure 9
figure 9

The temperature increment when \(x = 0.003\,\left( {\text{m}} \right),\,\,q_{0} = 100\,\,\left( {{\text{W/m}}^{{2}} } \right),\,\,\tau_{T} = 10\,\left( {\text{s}} \right)\) and various values \(\tau_{q} \left( {\text{s}} \right)\) based on the two-temperature model.

Figure 6 shows that the two temperature parameter \(\beta\) has a significant effect on the temperature increment distribution. The conductive temperature increment and dynamical temperature increment have the same behavior but with different values because of the impact of the two-temperature parameter. In the context of the one-temperature model \(\beta = 0\), the conductive and dynamical temperature increments are the same, and the solid blue curve represents them. It is noted that an increase in the value of the two-temperature parameter leads to an increase in the value of conductive temperature increment and a decrease in the value of dynamical temperature increment.

Figure 7 shows that the conductive and dynamical temperature increment distributions when \(x = 0.003\,\left( {\text{m}} \right),\,\,\tau_{T} = 10\,\left( {\text{s}} \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and various values of the surface heat flux \(q_{0} \,\left( {W/m^{2} } \right) = \left( {100,200} \right)\) to stand on its effects. It is noted that the value of the heat flux on the surface of the skin tissue has a significant impact on the conductive and dynamical temperature increment distributions. An increase in the value of the heat flux leads to an increase in the value of the conductive and dynamical temperature increment distributions.

Figure 8 shows that the conductive and dynamical temperature increment distributions when \(x = 0.003\,\left( {\text{m}} \right),\,\,q_{0} = 100\,\,\left( {{\text{W/m}}^{{2}} } \right),\,\,\tau_{q} = 20\,\left( {\text{s}} \right)\) and various values of the second relaxation time parameter \(\tau_{T} \left( {\text{s}} \right) = \left( {0.0,10.0} \right)\) to stand on its effects. It is noted that the value of the second relaxation time parameter has a significant impact on the conductive and dynamical temperature increment distributions. An increase in the value of the second relaxation time parameter leads to an increase in the value of the conductive and dynamical temperature increment distributions.

Figure 9 shows that the conductive and dynamical temperature increment distributions when \(x\left( {\text{m}} \right) = 0.003\,,\,\,q_{0} \left( {{\text{W/m}}^{{2}} } \right) = 100\,\,,\,\,\tau_{T} \left( {\text{s}} \right) = 10\,\) and various values of the second relaxation time parameter \(\tau_{q} \left( {\text{s}} \right) = \left( {0.0,10.0} \right)\) to stand on its effects. It is noted that the value of the first relaxation time parameter has a significant impact on the conductive and dynamical temperature increment distributions. An increase in the value of the first relaxation time parameter leads to a decrease in the value of the conductive and dynamical temperature increment distributions.

Figures 10 and 11 show the differences between the conductive temperature increment and dynamical temperature increment with respect to a wide range of time t \(\left( {0 \le t \le 100\left( {\text{s}} \right)\,} \right)\) and distance \(x\left( {0 \le x \le 0.006\left( {\text{m}} \right)\,} \right)\). Figure 8 represents the conductive temperature increment and dynamical temperature increment when \(q_{0} \left( {{\text{W/m}}^{{2}} } \right) = 100\) and Fig. 9 when \(q_{0} \left( {{\text{W/m}}^{{2}} } \right) = 200\). These two figures agree with the results in Figs. 3 and 6 and confirm that the value of the heat flux on the surface of the skin tissue has significant effects on the conductive and dynamical temperature increment distributions.

Figure 10
figure 10

The conductive and dynamical temperature increment when \(q_{0} = 100\left( {{\text{W/m}}^{{2}} } \right)\) based on the two-temperature model.

Figure 11
figure 11

The conductive and dynamical temperature increment when \(q_{0} = 200\left( {{\text{W/m}}^{{2}} } \right)\) based on the two-temperature model.

Figures 12 and 13 show the differences between the conductive temperature increment and dynamical temperature increment with respect to a wide range of time t \(\left( {0 \le t \le 100\left( {\text{s}} \right)\,} \right)\) and distance \(x\left( {0 \le x \le 0.006\left( {\text{m}} \right)\,} \right)\) for various values of the second relaxation time parameter \(\tau_{T} \left( {\text{s}} \right) = \left( {0.0,10.0} \right)\) This two figures agree with the results in the Figs. 4 and 7 and confirm that the value of the second relaxation time parameter has significant effects on the conductive and dynamical temperature increment distributions.

Figure 12
figure 12

The conductive and dynamical temperature increment when \(\tau_{T} \left( {\text{s}} \right) = 0.0\) based on the two-temperature model.

Figure 13
figure 13

The conductive and dynamical temperature increment when \(\tau_{T} \left( {\text{s}} \right) = 15\) based on the two-temperature model.

Figures 14 and 15 show the differences between the conductive temperature increment and dynamical temperature increment with respect to a wide range of time t \(\left( {0 \le t \le 100\left( {\text{s}} \right)\,} \right)\) and distance \(x\left( {0 \le x \le 0.006\left( {\text{m}} \right)\,} \right)\) for various values of the first relaxation time parameter \(\tau_{q} \left( {\text{s}} \right) = \left( {10,20} \right)\). These two figures agree with the results in Figs. 5 and 9 and confirm that the value of the first relaxation time parameter has significant effects on the conductive and dynamical temperature increment distributions.

Figure 14
figure 14

The conductive and dynamical temperature increment when \(\tau_{q} \left( {\text{s}} \right) = 10\).

Figure 15
figure 15

The conductive and dynamical temperature increment when \(\tau_{q} \left( {\text{s}} \right) = 20\).

Here, we should refer to the agreement of the results of this research with the results of other researches in order to confirm the validity of the current results and the validity of the proposed mathematical model as well as the method of solution.

The current results agree with the results of Youssef and Alghamdi28,37. Moreover, the current resluts agree with the result in figuer 6 of Kundu and Dewanjee29 for the case of boundary condition of cas (b). To verify our results with the results in Kundu and Dewanjee29, we represented Fig. 16 in which we use the same values of the parameters in Kundu and Dewanjee29 as follows:

$$ \rho_{b} W_{b} = 0.5\,{\text{kg/m}}^{{3}} {\text{s}},\;\tau_{T} = 0,\,\,\tau_{q} = 20\,{\text{s,}}\;{\text{K } = \text{ 0}}{\text{.2 W/m }}^{{\text{o}}} {\text{C,}}\;t = 160\,{\text{s,}}\,\,{\text{q}}_{0} = 500\,{\text{W/m}}^{{2}} $$
Figure 16
figure 16

The absolute conductive and dynamical temperatures of the current model against Kundu and Dewanjee29 when \(\rho_{b} W_{b} = 0.5\,\,{\text{kg/m}}^{{3}} {\text{s}},\,\,\tau_{T} = 0,\,\,\tau_{q} = 20\,{\text{s}}\), \({\text{K } = \text{ 0}}{\text{.2 W/m }}^{{\text{o}}} {\text{C}}\), \(\,t = 160\,{\text{s}}\), \(\,{\text{q}}_{0} = 500\,{\text{W/m}}^{{2}}\).

and calculating the absolut values of the dynamical and conductive temperature by adding \(T_{b} = 37\,^{{\text{o}}} {\text{C}}\) to the increments as follwos:

$$ T_{D} = \left( {\theta + T_{b} } \right),\,\,\,T_{C} = \left( {\varphi + T_{b} } \right) $$

Figure 16 show that the current result of the absolute dynamical and conductive temperature are close to the values of the absolute temperature in Kundu and Dewanjee29. Thus, we have a guarantee that the results of the current research are correct and verified.

Conclusion

  • The dual-phase lag time parameters \(\tau_{T} ,\,\,\tau_{q}\) have significant effects on the conductive and the dynamical temperature increment of the skin tissue.

  • The initial heat flux applying on the surface of the skin tissue has significant on the conductive and the dynamical temperature increment.

  • The two-temperature parameter has a significant effect on the conductive and the dynamical temperature increment.

  • The Two-temperature dual-phase-lag (TTDPL) bioheat transfer model is a successful model to simulate the thermal behavior of the skin tissue.

  • The results in the current work agree with the results of Kundu and Dewanjee29. Therefore, the two-temperature dual-phase-lag (TTDPL) bioheat transfer model be a succsessful model to study the bioheat transfer of a skin tissue and supplies us with exact solutions.