Dual-phase lag model of heat transfer between blood vessel and biological tissue

Ewa Majchrzak; Mikołaj Stryczyński; Ewa Majchrzak; Mikołaj Stryczyński

doi:10.3934/mbe.2021081

Mathematical Biosciences and Engineering

2021, Volume 18, Issue 2: 1573-1589. doi: 10.3934/mbe.2021081

Previous Article Next Article

Research article Special Issues

Dual-phase lag model of heat transfer between blood vessel and biological tissue

Ewa Majchrzak ^,,
Mikołaj Stryczyński

Department of Computational Mechanics and Engineering, Silesian University of Technology Gliwice, Poland

Received: 15 December 2020 Accepted: 18 January 2021 Published: 01 February 2021

A single blood vessel surrounded by the biological tissue with a tumor is considered. The influence of the heating technique (e.g. ultrasound, microwave, etc.) is described by setting a fixed temperature for the tumor which is higher than the blood and tissue temperature. The temperature distribution for the blood sub-domain is described by the energy equation written in the dual-phase lag convention, the temperature distribution in the biological tissue with a tumor is described also by the dual-phase lag equation. The boundary condition on the contact surface between blood vessel and biological tissue and the Neumann condition are also formulated using the extended Fourier law. So far in the literature, the temperature distribution in a blood vessel has been described by the classical energy equation. It is not clear whether the Fourier's law applies to highly heated tissues in which a significant thermal blood vessel is distinguished, therefore, taking into account the heterogeneous inner structure of the blood, the dual-phase lag equation is proposed for this sub-domain. The problem is solved by means of the implicit scheme of the finite difference method. The computations were performed for various values of delay times, which were taken from the available literature, and the influence of these values on the obtained temperature distributions was discussed.

Keywords:

bioheat transfer,
thermal interaction between blood vessel and tissue,
dual-phase lag model,
finite difference method

Citation: Ewa Majchrzak, Mikołaj Stryczyński. Dual-phase lag model of heat transfer between blood vessel and biological tissue[J]. Mathematical Biosciences and Engineering, 2021, 18(2): 1573-1589. doi: 10.3934/mbe.2021081

Related Papers:

[1]	Hao Jia, Bin Chen, Dong Li, Yuzhen Jin . Strategy of boundary discretization in numerical simulation of laser propagation in skin tissue with vascular lesions. Mathematical Biosciences and Engineering, 2021, 18(3): 2455-2472. doi: 10.3934/mbe.2021125
[2]	J. A. López Molina, M. J. Rivera, E. Berjano . Electrical-thermal analytical modeling of monopolar RF thermal ablation of biological tissues: determining the circumstances under which tissue temperature reaches a steady state. Mathematical Biosciences and Engineering, 2016, 13(2): 281-301. doi: 10.3934/mbe.2015003
[3]	Oualid Kafi, Nader El Khatib, Jorge Tiago, Adélia Sequeira . Numerical simulations of a 3D fluid-structure interaction model for blood flow in an atherosclerotic artery. Mathematical Biosciences and Engineering, 2017, 14(1): 179-193. doi: 10.3934/mbe.2017012
[4]	Caixia Zheng, Huican Li, Yingying Ge, Yanlin He, Yugen Yi, Meili Zhu, Hui Sun, Jun Kong . Retinal vessel segmentation based on multi-scale feature and style transfer. Mathematical Biosciences and Engineering, 2024, 21(1): 49-74. doi: 10.3934/mbe.2024003
[5]	Mette S. Olufsen, Ali Nadim . On deriving lumped models for blood flow and pressure in the systemic arteries. Mathematical Biosciences and Engineering, 2004, 1(1): 61-80. doi: 10.3934/mbe.2004.1.61
[6]	Alejandro Rincón, Fabiola Angulo, Fredy E. Hoyos . Analysis of a generalized Fujikawa’s growth model. Mathematical Biosciences and Engineering, 2020, 17(3): 2103-2137. doi: 10.3934/mbe.2020112
[7]	EYK Ng, Leonard Jun Cong Looi . Numerical analysis of biothermal-fluids and cardiac thermal pulse of abdominal aortic aneurysm. Mathematical Biosciences and Engineering, 2022, 19(10): 10213-10251. doi: 10.3934/mbe.2022479
[8]	Nattawan Chuchalerm, Wannika Sawangtong, Benchawan Wiwatanapataphee, Thanongchai Siriapisith . Study of Non-Newtonian blood flow - heat transfer characteristics in the human coronary system with an external magnetic field. Mathematical Biosciences and Engineering, 2022, 19(9): 9550-9570. doi: 10.3934/mbe.2022444
[9]	Hu Dong, Gang Liu, Xin Tong . Influence of temperature-dependent acoustic and thermal parameters and nonlinear harmonics on the prediction of thermal lesion under HIFU ablation. Mathematical Biosciences and Engineering, 2021, 18(2): 1340-1351. doi: 10.3934/mbe.2021070
[10]	Adélia Sequeira, Rafael F. Santos, Tomáš Bodnár . Blood coagulation dynamics: mathematical modeling and stability results. Mathematical Biosciences and Engineering, 2011, 8(2): 425-443. doi: 10.3934/mbe.2011.8.425

Abstract

1. Introduction

Modeling of the tissue heating process during thermal therapy is most often based on the well-known Pennes equation ^[1]. It is the classical Fourier's equation in which the source functions describing the blood perfusion and the metabolic heat generation are introduced. The form of this equation results from the assumption of instantaneous propagation of the thermal wave in the domain considered. In recent years, this model is often replaced by the Cattaneo-Vernotte equation ^[2,3] or the dual-phase lag equation ^[4]. In the Cattaneo-Vernotte equation the phase lag called the relaxation time τ_q appears. This parameter takes into account the effect of the finite value of the thermal wave velocity and it is the lag time of the heat flux in relation to the temperature gradient. In the dual-phase lag equation the additional parameter called the thermalization time τ_T is introduced, which takes into account the delay time of the temperature gradient in relation to the heat flux. The thermalization time takes into account the small-scale response in space while the relaxation time takes into account the small-scale response in time ^[4].

The thermally significant blood vessels can affect the temperature distribution in the heated tissue during thermal therapy, e.g. ^[5,6,7,8,9]. On the other hand, the influence of a strongly heated tumor region on the temperature distribution of blood flowing through the vessel is also interesting. The mathematical models of this type of phenomena used so far are based on the Fourier-Kirchhoff equation for the blood sub-domain and one of the equations mentioned above for the tissue sub-domain. In this paper the energy equation for the blood sub-domain is assumed in the convention of the dual-phase lag model and the temperature distribution in the tissue domain is also described by the dual-phase lag equation. The model is supplemented by the appropriate boundary and initial conditions. These conditions are formulated here on the basis of the extended Fourier law in which the relaxation and thermalization times appear ^[10,11,12]. To the best of our knowledge, the extended energy equation with two lag times is presented here for the first time. The problem formulated is solved using the implicit scheme of the finite difference method. In the final part of the paper the results of computations are shown and the conclusions are presented.

2. Methods

2.1. Dual-phase lag equation

Heat conduction for the macroscale problems is described by the Fourier law

${\bf{q}}\, (X, \, \, t) = - \, {\rm{ \mathsf{ λ} }}\, \nabla \, T(X, \, t)$

(1)

where λ is the thermal conductivity of the material, $\nabla T(X, t)$ is the temperature gradient, q (X, t) is the heat flux, X, t are the geometrical co-ordinates and time. Introducing the formula (1) to the Fourier-Kirchhoff equation ^[13]

$c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot {\kern 1pt} \nabla T(X, t)} \right] = - \nabla \, \cdot \, {\bf{q}}{\kern 1pt} {\kern 1pt} (X, \, \, t) + Q(X, t)$

(2)

one obtains equation in the form

$c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right] = \nabla \left[ {{\rm{ \mathsf{ λ} }}{\kern 1pt} \nabla \, T(X, \, \, t)} \right]\, + Q(X, t)$

(3)

where c is the specific heat, ρ is the mass density, w is the velocity vector and Q (X, t) is the capacity of internal heat sources.

The generalization of the Fourier is presented by Tzou ^[4], in particular

${\bf{q}}\, (X, \, \, t + {{\rm{ \mathsf{ τ} }} _q}\, ) = - \, {\rm{ \mathsf{ λ} }}\, \nabla \, T(X, \, \, t + {{\rm{ \mathsf{ τ} }} _T}\, )$

(4)

where τ_q and τ_T are the phase lags.

Using the Taylor series expansions, the following first-order approximation of Eq (4) is obtained

${\bf{q}}(X, \, \, t)\, + \, {{\rm{ \mathsf{ τ} }}_q}\, \frac{{\partial \, {\bf{q}}\, (X, \, \, t)}}{{\partial \, t}}\, = - \, {\rm{ \mathsf{ λ} }}\, \, \left[ {\nabla \, T(X, \, \, t)\, + \, {{\rm{ \mathsf{ τ} }}_T}\, \frac{{\partial \, \nabla \, T(X, \, \, t)}}{{\partial \, t}}} \right]$

(5)

From Eq (5) it follows that

$- {\kern 1pt} {\bf{q}}(X, \, \, t)\, = \, {{\rm{ \mathsf{ τ} }}_q}\, \frac{{\partial \, {\bf{q}}\, (X, \, \, t)}}{{\partial \, t}}\, + \, {\rm{ \mathsf{ λ} }}\, \, \left[ {\nabla \, T(X, \, \, t)\, + \, {{\rm{ \mathsf{ τ} }}_T}\, \frac{{\partial \, \nabla \, T(X, \, \, t)}}{{\partial \, t}}} \right]$

(6)

Introducing this formula to the Eq (2) one has

$\begin{gathered} c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right] = \\ {{\rm{ \mathsf{ τ} }}_q}\frac{\partial }{{\partial \, t}}\left[ {\nabla \cdot \, {\bf{q}}(X, \, \, t)} \right] + \nabla \left[ {{\rm{ \mathsf{ λ} }}\, \nabla \, T(X, \, \, t)} \right] + {{\rm{ \mathsf{ τ} }}_T}\, \nabla \left[ {{\rm{ \mathsf{ λ} }}\, \frac{{\partial \, \nabla \, T(X, \, \, t)}}{{\partial \, t}}} \right] + Q(X, t) \\ \end{gathered}$

(7)

Because (c.f. Eq (2))

$\nabla \, \cdot \, {\bf{q}}{\kern 1pt} {\kern 1pt} (X, \, \, t) = - c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right] + Q(X, t)$

(8)

thus

$\begin{gathered} c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right] = {{\rm{ \mathsf{ τ} }}_q}\frac{\partial }{{\partial \, t}}\left\{ { - c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right] + Q(X, t)} \right\} + \\ \nabla \left[ {{\rm{ \mathsf{ λ} }}\, \nabla \, T(X, \, \, t)} \right] + {{\rm{ \mathsf{ τ} }}_T}\, \nabla \left[ {{\rm{ \mathsf{ λ} }}\, \frac{{\partial \, \nabla \, T(X, \, \, t)}}{{\partial \, t}}} \right] + Q(X, t) \\ \end{gathered}$

(9)

$\begin{gathered} c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right] + {{\rm{ \mathsf{ τ} }}_q}\frac{{\partial \, }}{{\partial \, t}}\left[ {c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right]} \right] = \\ \nabla \left[ {{\rm{ \mathsf{ λ} }}\, \nabla \, T(X, \, \, t)} \right] + {{\rm{ \mathsf{ τ} }}_T}\frac{{\partial \, }}{{\partial \, t}}\, \left[ {\nabla \left[ {{\rm{ \mathsf{ λ} }}\nabla \, T(X, \, \, t)} \right]} \right] + Q(X, t) + {{\rm{ \mathsf{ τ} }}_q}\frac{{\partial \, Q(X, t)}}{{\partial \, t}} \\ \end{gathered}$

(10)

For the constant thermophysical parameters the Eq (10) takes a form

$\begin{gathered} c{\rm{ \mathsf{ ρ} }}\left[ {\frac{{\partial \, T(X, t)}}{{\partial \, t}} + {\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t) + {{\rm{ \mathsf{ τ} }}_q}\frac{{{\partial ^2}\, T(X, \, \, t)}}{{\partial \, {t^2}}} + {{\rm{ \mathsf{ τ} }}_q}\frac{{\partial \, \left[ {{\bf{w}}{\kern 1pt} \cdot \nabla {\kern 1pt} T(X, t)} \right]}}{{\partial \, t}}} \right] = \\ {\rm{ \mathsf{ λ} }}{\nabla ^2}T(X, \, \, t) + {\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_T}\, \frac{{\partial \, \left[ {{\nabla ^2}T(X, \, \, t)} \right]}}{{\partial \, t}}\, + Q(X, t) + {{\rm{ \mathsf{ τ} }}_q}\frac{{\partial \, Q(X, t)}}{{\partial \, t}} \\ \end{gathered}$

(11)

It should be noted that in the case of the dual-phase lag model the boundary conditions should be formulated in a different way than in the Fourier model ^[10,11,12]. For example, the Neumann boundary condition has the following form

${q_b}\, (X, t) + \, {{\rm{ \mathsf{ τ} }}_q}\, \frac{{\partial {q_b}(X, t)\, }}{{\partial \, t}}\, = - \, {\rm{ \mathsf{ λ} }}\, \, \left[ {{\bf{n}} \cdot \nabla \, T(X, t)\, + \, {{\rm{ \mathsf{ τ} }}_T}\, \frac{{\partial \left[ {{\bf{n}} \cdot \nabla \, T(X, t)} \right]\, }}{{\partial \, t}}} \right]$

(12)

where n is the normal outward vector and q_b (X, t) is the known boundary heat flux.

Under the assumption that for t = 0: $\nabla \, T(X, t) = 0$ and $\nabla \, \cdot \, {\bf{q}}{\kern 1pt} {\kern 1pt} (X, \, \, t) = 0$ (Eq (2)), the initial conditions are the following ^[14]

$t = 0:\quad T(X, t) = {T_p}, \quad {\left. {\frac{{\partial \, T(X, t)}}{{\partial \, t}}} \right|_{\, \, t\, = \, 0}} = \frac{{Q(X, 0)}}{{c{\rm{ \mathsf{ ρ} }}}}$

(13)

where T_p is the initial temperature.

2.2. Formulation of the problem

A single blood vessel and surrounding biological tissue with a tumor region, as shown in Figure 1, is considered (axisymmetric problem). As in the paper ^[8], it is assumed that the influence of the heating technique (e.g. ultrasound, microwave, etc.) is described by setting a fixed temperature for the tumor which is higher than the blood and tissue temperature.

Figure 1. Domain considered.

DownLoad: Full-Size Img PowerPoint

Blood temperature is described by the equation (c.f. Eq (11))

${c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}\left( {\frac{{\partial \, {T_b}}}{{\partial \, t}} + v\frac{{\partial {T_b}}}{{\partial {\kern 1pt} z}} + {{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}\frac{{{\partial ^2}\, {T_b}}}{{\partial \, {t^2}}} + {{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}{\kern 1pt} v\frac{{{\partial ^2}{T_b}}}{{\partial \, t\partial {\kern 1pt} z}}} \right) = {{\rm{ \mathsf{ λ} }}_b}{\nabla ^2}{T_b} + {{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}{\kern 1pt} {{\rm{ \mathsf{ λ} }}_b}\frac{{\partial \left( {{\nabla ^2}{T_b}} \right)\, }}{{\partial \, t}}\, \, + {Q_{met{\kern 1pt} b}}$

(14)

where

${\nabla ^2}{T_b}\, = \frac{1}{r}\frac{{\partial \, }}{{\partial \, r}}\left( {{\kern 1pt} r\frac{{\partial \, {T_b}}}{{\partial \, r}}} \right){\rm{ + }}\frac{{{\partial ^2}{T_b}\, }}{{\partial \, {z^2}}} = \frac{1}{r}\frac{{\partial \, {T_b}}}{{\partial \, r}} + \frac{{{\partial ^2}{T_b}\, }}{{\partial \, {r^2}}} + \frac{{{\partial ^2}{T_b}\, }}{{\partial \, {z^2}}}$

(15)

and v is the blood flow velocity in the axial direction, Q_{met b} is the metabolic heat source (constant value of Q_{met b} is assumed here), while T_b = T_b(r, z, t).

Temperature field in the tissue with a tumor is described by equation

$c{\rm{ \mathsf{ ρ} }}\left( {\frac{{\partial \, T}}{{\partial \, t}} + {{\rm{ \mathsf{ τ} }}_q}\frac{{{\partial ^2}\, T}}{{\partial \, {t^2}}}} \right) = {\rm{ \mathsf{ λ} }}{\nabla ^2}T + {\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_T}\, \frac{{\partial \left( {{\nabla ^2}T} \right)\, }}{{\partial \, t}}\, + Q + {{\rm{ \mathsf{ τ} }}_q}\frac{{\partial \, Q}}{{\partial \, t}}$

(16)

where

$Q = \, {w_b}{c_b}\left( {{T_a} - T} \right) + {Q_{met}}$

(17)

while w_b is the blood perfusion rate, T_a is the arterial blood temperature, Q_met is the metabolic heat source (constant value of Q_met is assumed here) and T = T(r, z, t).

On the contact surface between blood vessel and biological tissue the continuity condition of temperature and heat flux is assumed

$(r, z) \in {\Gamma _c}\, :\, \;\left\{ \begin{array}{l} {T_b} = T \\ {q_b} = q{\kern 1pt} \, \, \\ \end{array} \right.{\kern 1pt} \,$

(18)

On the outer surface of the tissue the constant temperature (body core temperature) is accepted. For z = 0, z = Z and r = 0 the no-flux conditions are assumed.

The initial conditions are also known (c.f. Eq (13))

$\begin{gathered} t = 0:\quad {T_b} = T = {T_p}, \\ {\left. {\frac{{\partial \, {T_b}}}{{\partial \, t}}} \right|_{\, \, t\, = \, 0}} = \frac{{{Q_{met{\kern 1pt} b}}}}{{{c_{b{\kern 1pt} }}{{\rm{ \mathsf{ ρ} }}_b}}}, \quad {\left. {\frac{{\partial \, T}}{{\partial \, t}}} \right|_{\, \, t\, = \, 0}} = \frac{{{w_b}{c_b}\left( {{T_a} - {T_p}} \right) + {Q_{met}}}}{{c{\rm{ \mathsf{ ρ} }}}} \\ \end{gathered}$

(19)

It should be noted that in the condition (18) the dependence (12) must be taken into account. After mathematical manipulations one has ^[12]

$(r, z) \in {\Gamma _c}\, :\, \;\left\{ \begin{array}{l} {T_b} = T \\ - \, {{\rm{ \mathsf{ λ} }}_b}{\kern 1pt} {\bf{n}} \cdot \nabla \, {T_b} - \, {{\rm{ \mathsf{ λ} }}_b}\left( {{{\rm{ \mathsf{ τ} }}_{Tb}} + {{\rm{ \mathsf{ τ} }}_q}} \right)\frac{{\partial \left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)\, }}{{\partial \, t}} - {{\rm{ \mathsf{ λ} }}_b}{{\rm{ \mathsf{ τ} }}_q}{{\rm{ \mathsf{ τ} }}_{Tb}}\frac{{{\partial ^2}\, \left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)}}{{\partial \, {t^2}}}\, \, = \\ - \, {\rm{ \mathsf{ λ} }}\, {\bf{n}} \cdot \nabla \, T - \, {\rm{ \mathsf{ λ} }}\, \left( {{{\rm{ \mathsf{ τ} }}_T} + {{\rm{ \mathsf{ τ} }}_{qb}}} \right)\frac{{\partial \left( {{\bf{n}} \cdot \nabla \, T} \right)\, }}{{\partial \, t}} - {\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_{qb}}{{\rm{ \mathsf{ τ} }}_T}\frac{{{\partial ^2}\, \left( {{\bf{n}} \cdot \nabla \, T} \right)}}{{\partial \, {t^2}}}\, \, \\ \end{array} \right.{\kern 1pt} \,$

(20)

2.3. Method of solution

The problem formulated is solved using implicit scheme of the finite difference method. Let $T_{i, j}^f = T({r_i}, {z_j}, f\Delta t)$ where ∆t is the time step, r_i = ih, z_j = jk (Figure 2) and f = 0, 1, …, F. Taking into account the initial conditions (19) one has: $T_{i, j}^{{\kern 1pt} 0} = {T_p},$ for blood subdomain: $T_{i, j}^{{\kern 1pt} 1} = {T_p} + {Q_{met{\kern 1pt} b}}\Delta t/({c_{b{\kern 1pt} }}{{\rm{ \mathsf{ ρ} }}_b})$ , for tissue sub-domain: $T_{i, j}^{{\kern 1pt} 1} = {T_p} + \left[ {{w_{b{\kern 1pt} }}{c_b}\left( {{T_a} - {T_p}} \right) + {Q_{met}}} \right]\Delta t/\left( {c{\rm{ \mathsf{ ρ} }} + {{\rm{ \mathsf{ τ} }}_q}{w_{b{\kern 1pt} }}{c_b}} \right).$

Figure 2. Discretization.

DownLoad: Full-Size Img PowerPoint

For transition t^f → t^f+1(f ≥ 1) the following approximate form of Eq (14) resulting from the introduction of adequate differential quotients is used

$\begin{gathered} {c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}\frac{{\, T_{i, j}^{f + 1} - T_{i, j}^f}}{{\Delta \, t}} + {c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}\frac{{T_{i, j}^{f + 1} - 2T_{i, j}^f + T_{i, j}^{f - 1}}}{{{{\left( {\Delta t} \right)}^2}}} = \\ {{\rm{ \mathsf{ λ} }}_b}\left( {{\nabla ^2}{T_b}} \right)_{{\kern 1pt} {\kern 1pt} i, j}^{{\kern 1pt} f + 1} + \frac{{{{\rm{ \mathsf{ λ} }}_b}{{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}}}{{\Delta \, t}}\left[ {\left( {{\nabla ^2}{T_b}} \right)_{{\kern 1pt} {\kern 1pt} i, j}^{{\kern 1pt} f + 1} - \left( {{\nabla ^2}{T_b}} \right)_{{\kern 1pt} {\kern 1pt} i, j}^{{\kern 1pt} f}} \right] + {Q_b} - \\ {c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}v\frac{{T_{i, j + 1}^{f + 1} - T_{i, j - 1}^{f + 1}}}{{2k}} - \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}{\kern 1pt} v}}{{\Delta \, t}}\left( {\frac{{T_{i, j + 1}^{f + 1} - T_{i, j - 1}^{f + 1}}}{{2k}} - \frac{{T_{i, j + 1}^f - T_{i, j - 1}^f}}{{2k}}} \right) \\ \end{gathered}$

(21)

$\begin{gathered} {c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}\frac{{\, T_{i, j}^{f + 1} - T_{i, j}^f}}{{\Delta \, t}}{\kern 1pt} + {c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}\frac{{T_{i, j}^{f + 1} - 2T_{i, j}^f + T_{i, j}^{f - 1}}}{{{{\left( {\Delta t} \right)}^2}}} = \\ \frac{{{{\rm{ \mathsf{ λ} }}_b}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}} \right)}}{{\Delta \, t}}\left( {{\nabla ^2}{T_b}} \right)_{{\kern 1pt} {\kern 1pt} i, j}^{{\kern 1pt} f + 1} - \frac{{{{\rm{ \mathsf{ λ} }}_b}{{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}}}{{\Delta \, t}}\left( {{\nabla ^2}{T_b}} \right)_{{\kern 1pt} {\kern 1pt} i, j}^{{\kern 1pt} f} + {Q_b} - \\ \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}v\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}} \right)\left( {T_{i, j + 1}^{f + 1} - T_{i, j - 1}^{f + 1}} \right)}}{{2k\Delta \, t}} + \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}{\kern 1pt} v\left( {T_{i, j + 1}^f - T_{i, j - 1}^f} \right)}}{{2k\Delta \, t}} \\ \end{gathered}$

(22)

where

$\left( {{\nabla ^2}{T_b}} \right)_{{\kern 1pt} {\kern 1pt} i, j}^{{\kern 1pt} p}\, = \frac{1}{{{r_{i, j}}}}\frac{{T_{i + 1, j}^p - T_{i - 1, j}^p}}{{2h}} + \frac{{T_{i - 1, j}^p - 2T_{i, j}^p + T_{i + 1, j}^p}}{{{h^2}}} + \frac{{T_{i, j - 1}^p - 2T_{i, j}^p + T_{i, j + 1}^p}}{{{k^2}}}$

(23)

while p = f or p = f+1.

After mathematical manipulations one has

$\begin{gathered} T_{i, j}^{f + 1} = \frac{{{{\rm{ \mathsf{ λ} }}_b}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}} \right)}}{{D\Delta \, t}}\left( {\frac{1}{{{r_{i, j}}}}\frac{{T_{i + 1, j}^{f + 1} - T_{i - 1, j}^{f + 1}}}{{2h}} + \frac{{T_{i - 1, j}^{f + 1} + T_{i + 1, j}^{f + 1}}}{{{h^2}}} + \frac{{T_{i, j - 1}^{f + 1} + T_{i, j + 1}^{f + 1}}}{{{k^2}}}} \right) - \\ \frac{{{{\rm{ \mathsf{ λ} }}_b}{{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}}}{{{D_{\kern 1pt} }\Delta \, t}}\left( {\frac{1}{{{r_{i, j}}}}\frac{{T_{i + 1, j}^f - T_{i - 1, j}^f}}{{2h}} + \frac{{T_{i - 1, j}^f - 2T_{i, j}^f + T_{i + 1, j}^f}}{{{h^2}}} + \frac{{T_{i, j - 1}^f - 2T_{i, j}^f + T_{i, j + 1}^f}}{{{k^2}}}} \right) + \frac{{{Q_b}}}{D} - \\ \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}v\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}} \right)\left( {T_{i, j + 1}^{f + 1} - T_{i, j - 1}^{f + 1}} \right)}}{{2k\Delta \, tD}} + \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}{\kern 1pt} v\left( {T_{i, j + 1}^f - T_{i, j - 1}^f} \right)}}{{2k\Delta \, tD}} + \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}\left( {\Delta \, t + 2{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}} \right)}}{{D{{\left( {\Delta t} \right)}^2}}}T_{i, j}^f - \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}{{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}}}{{D{{\left( {\Delta t} \right)}^2}}}T_{i, j}^{f - 1} \\ \end{gathered}$

(24)

where

$D = \frac{{{c_b}{\kern 1pt} {{\rm{ \mathsf{ ρ} }}_b}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{q{\kern 1pt} b}}} \right)}}{{{{\left( {\Delta t} \right)}^2}}} + \frac{{{\rm{2}}{{\rm{ \mathsf{ λ} }}_b}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}} \right)}}{{{h^2}\Delta \, t}} + \frac{{{\rm{2}}{{\rm{ \mathsf{ λ} }}_b}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_{T{\kern 1pt} b}}} \right)}}{{{k^2}\Delta \, t}}$

(25)

Approximate form of Eq (16) is the following

$\begin{gathered} c{\kern 1pt} \, {\rm{ \mathsf{ ρ} }}\frac{{\, T_{i, j}^{f + 1} - T_{i, j}^f}}{{\Delta \, t}}\, + c\, {\rm{ \mathsf{ ρ} }}{\kern 1pt} {{\rm{ \mathsf{ τ} }}_q}\frac{{T_{i, j}^{f + 1} - 2T_{i, j}^f + T_{i, j}^{f - 1}}}{{{{\left( {\Delta t} \right)}^2}}} = {\rm{ \mathsf{ λ} }}\left( {{\nabla ^2}T} \right)_{i, j}^{f + 1} + \\ \frac{{{\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta \, t}}\left[ {\left( {{\nabla ^2}T} \right)_{i, j}^{f + 1} - \left( {{\nabla ^2}T} \right)_{i, j}^f} \right] + {w_b}{c_b}\left( {{T_a} - T_{i, j}^{f + 1}} \right) + {Q_{met}} - {{\rm{ \mathsf{ τ} }}_q}{w_b}{c_b}\frac{{\, T_{i, j}^{f + 1} - T_{i, j}^f}}{{\Delta \, t}} \\ \end{gathered}$

(26)

$\begin{gathered} \left( {c{\rm{ \mathsf{ ρ} }} + {{\rm{ \mathsf{ τ} }}_q}{w_b}{c_b}} \right)\frac{{\, T_{i, j}^{f + 1} - T_{i, j}^f}}{{\Delta \, t}}\, + c{\rm{ \mathsf{ ρ} }}\, {{\rm{ \mathsf{ τ} }}_q}\frac{{T_{i, j}^{f + 1} - 2T_{i, j}^f + T_{i, j}^{f - 1}}}{{{{\left( {\Delta t} \right)}^2}}} = \frac{{{\rm{ \mathsf{ λ} }}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_T}} \right)}}{{\Delta \, t}}\left( {{\nabla ^2}T} \right)_{i, j}^{f + 1} - \\ \frac{{{\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta \, t}}\left( {{\nabla ^2}T} \right)_{i, j}^f + {w_b}{c_b}\left( {{T_a} - T_{i, j}^{f + 1}} \right) + {Q_{met}} \\ \end{gathered}$

(27)

and after mathematical manipulations one has

$\begin{gathered} T_{i, j}^{f + 1} = \frac{{{\rm{ \mathsf{ λ} }}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_T}} \right)}}{{E\Delta \, t}}\left( {\frac{1}{{{r_{i, j}}}}\frac{{T_{i + 1, j}^{f + 1} - T_{i - 1, j}^{f + 1}}}{{2h}} + \frac{{T_{i - 1, j}^{f + 1} + T_{i + 1, j}^{f + 1}}}{{{h^2}}} + \frac{{T_{i, j - 1}^{f + 1} + T_{i, j + 1}^{f + 1}}}{{{k^2}}}} \right) - \\ \frac{{{\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_T}}}{{E\Delta \, t}}\left( {\frac{1}{{{r_{i, j}}}}\frac{{T_{i + 1, j}^f - T_{i - 1, j}^f}}{{2h}} + \frac{{T_{i - 1, j}^f - 2T_{i, j}^f + T_{i + 1, j}^f}}{{{h^2}}} + \frac{{T_{i, j - 1}^f - 2T_{i, j}^f + T_{i, j + 1}^f}}{{{k^2}}}} \right) + \\ \frac{{\left( {c{\rm{ \mathsf{ ρ} }} + {{\rm{ \mathsf{ τ} }}_q}{w_b}{c_b}} \right)\Delta t + 2c{\rm{ \mathsf{ ρ} }}\, {{\rm{ \mathsf{ τ} }}_q}}}{{E{{\left( {\Delta t} \right)}^2}}}T_{i, j}^f - \frac{{c\, {\rm{ \mathsf{ ρ} }}{{\rm{ \mathsf{ τ} }}_q}}}{{E{{\left( {\Delta t} \right)}^2}}}T_{i, j}^{f - 1} + \frac{{{w_b}{c_b}{T_a} + {Q_{met}}}}{E} \\ \end{gathered}$

(28)

where

$E = \frac{{\left( {c{\rm{ \mathsf{ ρ} }} + {{\rm{ \mathsf{ τ} }}_q}{w_b}{c_b}} \right)\Delta t + c\, {\rm{ \mathsf{ ρ} }}{{\rm{ \mathsf{ τ} }}_q}}}{{{{\left( {\Delta t} \right)}^2}}} + \frac{{{\rm{2 \mathsf{ λ} }}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_T}} \right)}}{{{h^2}\Delta \, t}} + \frac{{{\rm{2 \mathsf{ λ} }}\left( {\Delta \, t + {{\rm{ \mathsf{ τ} }}_T}} \right)}}{{{k^2}\Delta \, t}} + {w_b}{c_b}$

(29)

The second part of boundary condition (20) is approximated with respect to time in the following way

$\begin{gathered} - \, {{\rm{ \mathsf{ λ} }}_b}{\kern 1pt} {\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)^{f + 1}} - \, {{\rm{ \mathsf{ λ} }}_b}\left( {{{\rm{ \mathsf{ τ} }}_{Tb}} + {{\rm{ \mathsf{ τ} }}_q}} \right)\frac{1}{{\Delta t}}\left[ {{{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)}^{f + 1}} - {{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)}^f}} \right] - \\ {{\rm{ \mathsf{ λ} }}_b}{{\rm{ \mathsf{ τ} }}_q}{{\rm{ \mathsf{ τ} }}_{Tb}}\frac{1}{{{{\left( {\Delta t} \right)}^2}}}\left[ {{{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)}^{f + 1}} - 2{{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)}^f} + {{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)}^{f - 1}}} \right]\, \, = \\ - \, {\rm{ \mathsf{ λ} }}\, {\left( {{\bf{n}} \cdot \nabla \, T} \right)^{f + 1}}\, - \, {\rm{ \mathsf{ λ} }}\, \left( {{{\rm{ \mathsf{ τ} }}_T} + {{\rm{ \mathsf{ τ} }}_{qb}}} \right)\frac{1}{{\Delta t}}\left[ {{{\left( {{\bf{n}} \cdot \nabla \, T} \right)}^{f + 1}} - {{\left( {{\bf{n}} \cdot \nabla \, T} \right)}^f}} \right] - \\ {\rm{ \mathsf{ λ} }}{{\rm{ \mathsf{ τ} }}_{qb}}{{\rm{ \mathsf{ τ} }}_T}\frac{1}{{{{\left( {\Delta t} \right)}^2}}}\left[ {{{\left( {{\bf{n}} \cdot \nabla \, T} \right)}^{f + 1}} - 2{{\left( {{\bf{n}} \cdot \nabla \, T} \right)}^f} + {{\left( {{\bf{n}} \cdot \nabla \, T} \right)}^{f - 1}}} \right]\, \, \\ {\kern 1pt} \, \\ \end{gathered}$

(30)

$\begin{gathered} {A_1}{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)^{f + 1}} - {A_3}{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)^f} + {A_5}{\left( {{\bf{n}} \cdot \nabla \, {T_b}} \right)^{f - 1}}\, = \\ {A_2}{\left( {{\bf{n}} \cdot \nabla \, T} \right)^{f + 1}} - {A_4}{\left( {{\bf{n}} \cdot \nabla \, T} \right)^f} + {A_6}{\left( {{\bf{n}} \cdot \nabla \, T} \right)^{f - 1}}\, \, \\ \end{gathered}$

(31)

where

$\begin{gathered} {A_1}{\rm{ = }}{{\rm{ \mathsf{ λ} }}_b}\frac{{{{\left( {\Delta t} \right)}^2} + \Delta t\left( {{{\rm{ \mathsf{ τ} }}_{Tb}} + {{\rm{ \mathsf{ τ} }}_q}} \right) + {{\rm{ \mathsf{ τ} }}_q}{{\rm{ \mathsf{ τ} }}_{Tb}}}}{{{{\left( {\Delta t} \right)}^2}}}, \quad A{}_2 = {\rm{ \mathsf{ λ} }}\frac{{{{\left( {\Delta t} \right)}^2} + \Delta t\left( {{{\rm{ \mathsf{ τ} }}_T} + {{\rm{ \mathsf{ τ} }}_{qb}}} \right) + {{\rm{ \mathsf{ τ} }}_{qb}}{{\rm{ \mathsf{ τ} }}_T}}}{{{{\left( {\Delta t} \right)}^2}}} \\ {A_3} = {{\rm{ \mathsf{ λ} }}_b}\frac{{\Delta t\left( {{{\rm{ \mathsf{ τ} }}_{Tb}} + {{\rm{ \mathsf{ τ} }}_q}} \right) + 2{{\rm{ \mathsf{ τ} }}_q}{{\rm{ \mathsf{ τ} }}_{Tb}}}}{{{{\left( {\Delta t} \right)}^2}}}, \quad {A_4} = {\rm{ \mathsf{ λ} }}\frac{{\Delta t\left( {{{\rm{ \mathsf{ τ} }}_T} + {{\rm{ \mathsf{ τ} }}_{qb}}} \right) + 2{{\rm{ \mathsf{ τ} }}_{qb}}{{\rm{ \mathsf{ τ} }}_T}}}{{{{\left( {\Delta t} \right)}^2}}}, \quad {A_5} = {{\rm{ \mathsf{ λ} }}_b}\frac{{{{\rm{ \mathsf{ τ} }}_q}{{\rm{ \mathsf{ τ} }}_{Tb}}}}{{{{\left( {\Delta t} \right)}^2}}}, \quad {A_6} = {\rm{ \mathsf{ λ} }}\frac{{{{\rm{ \mathsf{ τ} }}_{qb}}{{\rm{ \mathsf{ τ} }}_T}}}{{{{\left( {\Delta t} \right)}^2}}}\, \\ \end{gathered}$

(32)

Let (s, j) is the node located on the contact surface. Taking into account the first part of boundary condition (20), one obtains

$\begin{gathered} {A_1}\frac{{T_{s, j}^{f + 1} - T_{s - 1, j}^{f + 1}}}{h} - {A_3}\frac{{T_{s, j}^f - T_{s - 1, j}^f}}{h} + {A_5}\frac{{T_{s, j}^{f - 1} - T_{s - 1, j}^{f - 1}}}{h}\, = \\ {A_2}\frac{{T_{s + 1, j}^{f + 1} - T_{s, j}^{f + 1}}}{h} - {A_4}\frac{{T_{s + 1, j}^f - T_{s, j}^f}}{h} + {A_6}\frac{{T_{s + 1, j}^{f - 1} - T_{s, j}^{f - 1}}}{h}\, \, \\ \end{gathered}$

(33)

$\begin{gathered} T_{s, j}^{f + 1} = \frac{{{A_1}}}{{{A_1} + {A_2}}}T_{s - 1, j}^{f + 1} + \frac{{{A_2}}}{{{A_1} + {A_2}}}T_{s + 1, j}^{f + 1} + \frac{{{A_3}}}{{{A_1} + {A_2}}}\left( {T_{s, j}^f - T_{s - 1, j}^f} \right) - \frac{{{A_4}}}{{{A_1} + {A_2}}}\left( {T_{s + 1, j}^f - T_{s, j}^f} \right) \\ - \frac{{{A_5}}}{{{A_1} + {A_2}}}\left( {T_{s, j}^{f - 1} - T_{s - 1, j}^{f - 1}} \right) + \frac{{{A_6}}}{{{A_1} + {A_2}}}\left( {T_{s + 1, j}^{f - 1} - T_{s, j}^{f - 1}} \right)\, \, \\ \end{gathered}$

(34)

The no-flux boundary condition (c.f. Eq (12))

${\bf{n}} \cdot \nabla \, T(X, t)\, + \, {{\rm{ \mathsf{ τ} }}_T}\, \frac{{\partial \left[ {{\bf{n}} \cdot \nabla \, T(X, t)} \right]\, }}{{\partial \, t}} = 0$

(35)

is approximated with respect to time, namely

$\left( {{\bf{n}} \cdot \nabla \, T} \right)_{}^{f + 1} + \, \frac{{{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t}}\left[ {\left( {{\bf{n}} \cdot \nabla \, T} \right)_{}^{f + 1} - \left( {{\bf{n}} \cdot \nabla \, T} \right)_{}^f} \right] = 0$

(36)

For the nodes with the coordinate z = 0 one has

$\frac{{\Delta t + {{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t}}\frac{{T_{i, 1}^{f + 1} - T_{i, 0}^{f + 1}}}{k} - \, \frac{{{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t}}\frac{{T_{i, 1}^f - T_{i, 0}^f}}{k} = 0$

(37)

and then

$T_{i, 0}^{f + 1} = T_{i, 1}^{f + 1} - \frac{{{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t + {{\rm{ \mathsf{ τ} }}_T}}}\left( {T_{i, 1}^f - T_{i, 0}^f} \right)$

(38)

Similarly, for nodes with the coordinate z = Z one obtains

$\frac{{\Delta t + {{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t}}\frac{{T_{i, n}^{f + 1} - T_{i, n - 1}^{f + 1}}}{k} - \, \frac{{{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t}}\frac{{T_{i, n}^f - T_{i, n - 1}^f}}{k} = 0$

(39)

and then

$T_{i, n}^{f + 1} = T_{i, n - 1}^{f + 1} + \frac{{{{\rm{ \mathsf{ τ} }}_T}}}{{\Delta t + {{\rm{ \mathsf{ τ} }}_T}}}\left( {T_{i, n}^f - T_{i, n - 1}^f} \right)$

(40)

For each transition t^f → t^f+1(f ≥ 1) the system of Eqs (24), (28), (34), (38) and (40) is solved using the iterative method.

3. Results

3.1. Temperature distribution in the tissue with a tumor

The thermally significant blood vessel of dimensions R = 0.001 m, Z = 20R = 0.02 m is considered. The blood flow velocity is equal to v = 0.2 m/s ^[6,15]. The outer radius of the domain considered is R₁ = 5R. The heating target volume (tumor) is specified as 1.5R ≤ r ≤ 2.5R, 1.5Z/4 ≤ z ≤ 3Z/4 (Figure 1). The following values of thermophysical parameters are assumed: thermal conductivities λ_b = 0.488 W/(m K), λ = 0.5 W/(m K), specific heats c_b = 3770 J/(kgK), c = 4000 J/(kgK), mass densities ρ_b = 1060 kg/m³, ρ = 1000 kg/m³, metabolic heat sources Q_metb = 245 W/m³, Q_met = 245 W/m³, blood perfusion rate w_b = 0.5 kg/(m³ s), arterial blood temperature T_a = 37 ℃ ^[16].

The values of phase lags reported in literature vary greatly from one to another. For example, in the papers ^[9,17] the relaxation and thermalization times are of the order of 10 s. On the other hand, in the paper ^[18] it is found that the phase lag times for heat flux and temperature gradient for the living tissues are very close to each other and are within the range [0.464 s, 6.825 s]. Due to the lack of precise data, in this work the extreme values of time delays, namely τ_qb = τ_q = τ_Tb = τ_T = 0.464 s and τ_qb = τ_q = τ_Tb = τ_T = 6.825 s are considered and as suggested by Zhang ^[18], the same values of these delay times are accepted.

The initial temperature of tissue and blood is equal to T_p = 37 ℃, on the outer surface of the tissue the boundary temperature equals 37 ℃, while in the tumor region the constant temperature 45 ℃ is accepted.

The computations have been done under the assumption that m = 50, n = 200, m₁ = 15, m₂ = 25, n₁ = 50, n₂ = 150 (Figure 2) and Δt = 0.005 s.

In Figures 3 and 4 the temperature distribution after 5 second is shown. Figures 5 and 6 illustrate the temperature history at the points A (0.0014 m, 0.01 m), B (0.003 m, 0.01 m) obtained using the DPL model and the classical ones. As one can see, for small values of time delays (Figure 5) the temperature values in the domain stabilize quite quickly and the system reaches the steady state conditions, while for large values of time delays (Figure 6), this process takes much longer. The differences between the courses of thermal processes are great especially in the initial stage of heating. After reaching the steady state, the curves coincide, of course. As expected, the heating process is slower for the model with lag times and the longer delay times prolong the duration of the tissue heating process.

Figure 3. Temperature distribution after 5 seconds (τ_qb = τ_q = τ_Tb = τ_T = 0.464 s).

DownLoad: Full-Size Img PowerPoint

Figure 4. Temperature distribution after 5 seconds (τ_qb = τ_q = τ_Tb = τ_T = 6.825 s).

DownLoad: Full-Size Img PowerPoint

Figure 5. Temperature history at the points A (0.0014 m, 0.01 m), B (0.003 m, 0.01 m), solid lines-DPL model (τ_qb = τ_q = τ_Tb = τ_T = 0.464 s), dashed lines-classical model.

DownLoad: Full-Size Img PowerPoint

Figure 6. Temperature history at the points A(0.0014 m, 0.01 m), B(0.003 m, 0.01 m), solid lines-DPL model (τ_qb = τ_q = τ_Tb = τ_T = 6.825 s), dashed lines-classical model.

DownLoad: Full-Size Img PowerPoint

3.2. Temperature distribution in the radial direction

Figures 7–9 illustrate the temperature distribution in the radial direction for z = Z/2 = 0.01 m for DPLM and models based on the Fourier-Kirchhoff equations, respectively. The curves correspond to times of 1, 2, 3, 4 and 5 s. In the blood sub-domain (r ≤ 0.001 m), a slight increase in temperature is observed only near the blood vessel wall. The introduction of delay times slows down the changes of instantaneous and local temperatures in the sub-domains considered.

Figure 7. Temperature distribution in the radial direction (z = Z/2)-DPL model (τ_qb = τ_q = τ_Tb = τ_T = 0.464 s).

DownLoad: Full-Size Img PowerPoint

Figure 8. Temperature distribution in the radial direction (z = Z/2)-DPL model (τ_qb = τ_q = τ_Tb = τ_T = 6.825 s).

DownLoad: Full-Size Img PowerPoint

Figure 9. Temperature distribution in the radial direction (z = Z/2)-Fourier-Kirchhoff model.

DownLoad: Full-Size Img PowerPoint

3.3. Temperature of the blood vessel wall

As can be seen in the previous Figures, the temperature in the blood vessel does not change much. Figures 10 and 11 show the changes in the temperature of the blood vessel wall along the z axis for times 1, 2, 3, 4 and 5 s. The influence of the increased temperature in the tumor region on the changes in the wall temperature is clearly visible. Outside of this region, its temperature varies to a small extent. For low values of delay times (Figure 7), after about five seconds the system reaches steady state conditions, while for large values of delay times (Figure 8) it is just the beginning of the heating process.

Figure 10. Temperature of blood vessel wall (τ_qb = τ_q = τ_Tb = τ_T = 0.464 s) after 1, 2, 3, 4 and 5 s.

DownLoad: Full-Size Img PowerPoint

Figure 11. Temperature of blood vessel wall (τ_qb = τ_q = τ_Tb = τ_T = 6.825 s) after 1, 2, 3, 4 and 5 s.

DownLoad: Full-Size Img PowerPoint

It should be noted that the influence of the blood flow direction on the computations results is clearly visible. As it approaches the heated tumor, the blood warms up more and more, and it cools more slowly as it leaves the heated sub-domain.

3.4. Comparison of DPL model with the classical model for the blood sub-domain

As mentioned before, in this paper the energy equation for the blood sub-domain is assumed in the convention of the dual-phase lag model. Therefore, it should be checked whether the introduction of delay times to the equation describing the temperature field in the blood sub-domain significantly changes the computations results. For this purpose, calculations were made for τ_qb = τ_q = 3 s, τ_Tb = τ_T = 1 s ^[19] (the same relaxation time and thermalization time for both sub-domains) and τ_qb = τ_Tb = 0, τ_q = 3 s, τ_T = 1 s, respectively (mixed model).

The results of the comparison are presented in the form of heating curves at the points A and B (Figure 12) and additionally at the point C (0.0009 m, 0.01 m) located in the blood sub-domain and at the point D (0.001m, 0.01 m) located on the blood vessel wall (Figure 13).

Figure 12. Heating curves at the points A(0.0014 m, 0.01 m) and B(0.003 m, 0.01 m), solid lines-DPL model (τ_qb = τ_q = 3 s, τ_Tb = τ_T = 1 s), dashed lines - mixed model (τ_qb = τ_Tb = 0, τ_q = 3 s, τ_T = 1 s).

DownLoad: Full-Size Img PowerPoint

Figure 13. Heating curves at the points C(0.0009 m, 0.01 m) and D(0.001 m, 0.01 m), solid lines - DPL model (τ_qb = τ_q = 3 s, τ_Tb = τ_T = 1 s), dashed lines - mixed model (τ_qb = τ_Tb = 0, τ_q = 3 s, τ_T = 1 s).

DownLoad: Full-Size Img PowerPoint

As can be seen, the differences are slight and do not exceed 0.5 ℃ on the blood vessel wall (point D) and are practically invisible in the blood vessel sub-domain (point C) as well as in the tissue sub-domain (points A and B). Summing up, the introducing the delay times into the equation describing the temperature distribution in the blood vessel has a relatively small impact on the results (for the assumed values of relaxation and thermalization times, of course).

4. Discussion

The single blood vessel surrounded by the biological tissue with a tumor has been considered. A constant temperature of 45 ℃ has been assumed in the heated tumor region. The temperature fields have been described by dual-phase lag equation both for the blood region as well as for the tissue with a tumor. For comparison, the classical models have been also considered. Additionally, computations were made for the 'mixed' model, in which the temperature field in the tissue sub-domain is described by the DPL equation, while the temperature field in the blood vessel by the classical model. For such formulated mathematical description, an algorithm based on FDM in an implicit form was presented. The numerous numerical computations concerning the different variants of the basic model have been made. It turned out (which could of course be expected) that the visible differences between the solutions can be observed. The limit results of all variants for the steady state coincide, which confirms the correctness of the assumed numerical algorithm.

It should be noted that for all variants of computations and assumed delay times, the temperature in the blood vessel did not rise much (about 1 ℃), mostly near the blood vessel wall. Comparison of the calculation results obtained by the DPL model and mixed model in which the classical equation for the blood sub-domain has been assumed showed slight changes in temperature.

Obviously, the results are closely related to the adopted relaxation and thermalization times, and there is a need for further experimental studies to more accurately estimate these delay times.

Here, a constant temperature of the tumor region has been assumed. In future, the heating schemes should be taken into account, this means the external heating power density and heating duration ^[6,7].

Acknowledgments

The research is financed from financial resources from the statutory subsidy of the Faculty of Mechanical Engineering, Silesian University of Technology in 2021.

Conflict of interest

Authors declare no conflicts of interest in this paper.

References

[1]	H. H. Pennes, Analysis of tissue and arterial blood temperatures in the resting human forearm, J. Appl. Physiol., 1 (1948), 3–23. doi: 10.1152/jappl.1948.1.1.3
[2]	M. C. Cattaneo, A form of heat conduction equation which eliminates the paradox of instantaneous propagation, C.R. Acad. Sci. I - Math., 247 (1958), 431–433.
[3]	P. Vernotte, Les paradoxes de la theorie continue de l'equation de la chaleur, Comp. Rend., 246 (1948), 3154–3155.
[4]	D. Y. Tzou, A unified field approach for heat conduction from macro- to micro- scales, J. Heat Transfer, 117 (1995), 8–16. doi: 10.1115/1.2822329
[5]	J. Crezee, J. W. Lagendijk, Temperature uniformity during hyperthermia: the impact of large vessels, Phys. Med. Biol., 37 (1992), 1321–1337. doi: 10.1088/0031-9155/37/6/009
[6]	T. C. Shih, H. L. Liu, T. L. Horng, Cooling effect of thermally significant blood vessels in perfused tumor tissue during thermal therapy, Int. Commun. Heat Mass Transfer, 33 (2006), 135–141. doi: 10.1016/j.icheatmasstransfer.2005.08.003
[7]	T. C. Shih, T. L. Horng, H. W. Huang, K. C. Ju, T. C. Huang, P. Y. Chen, et al., Numerical analysis of coupled effects of pulsatile blood flow and thermal relaxation time during thermal therapy, Int. J. Heat Mass Transfer, 55 (2012), 3763–3773. doi: 10.1016/j.ijheatmasstransfer.2012.02.069
[8]	K. Khanafer, J. L. Bull, I. Pop, R. Berguer, Influence of pulsatile blood flow and heating scheme on the temperature distribution during hyperthermia treatment, Int. J. Heat Mass Transfer, 50 (2007), 4883–4890. doi: 10.1016/j.ijheatmasstransfer.2007.01.062
[9]	M. Jamshidi, J. Ghazanfarian, Blood flow effects in thermal treatment of three-dimensional non-Fourier multilayered skin structure, Heat Transfer Eng., 2020 (2020), 1–18.
[10]	J. R. Ho, C. P. Kuo, W. S. Jiaung, Study of heat transfer in multilayered structure within the framework of dual-phase-lag heat conduction model using lattice Boltzmann method, Int. J. Heat Mass Transfer, 46 (2003), 55–69. doi: 10.1016/S0017-9310(02)00260-0
[11]	E. Majchrzak, G. Kałuża, Analysis of thermal processes occurring in the heated multilayered metal films using the dual-phase lag model, Arch. Mech., 69 (2017), 275–287.
[12]	E. Majchrzak, B. Mochnacki, Dual-phase lag model of thermal processes in a multi-layered microdomain subjected to a strong laser pulse using the implicit scheme of FDM, Int. J. Therm. Sci., 13 (2018), 240–251.
[13]	G. Hauke, An Introduction to Fluid Mechanics and Transport Phenomena, Springer, New York, 2008.
[14]	M. Ciesielski, B. Mochnacki, E. Majchrzak, Integro-differential form of the first-order dual phase lag heat transfer equation and its numerical solution using the Control Volume Method, Arch. Mech., 72 (2020), 415–444.
[15]	J. C. Chato, Heat transfer to blood vessels, J. Biomech. Eng., 102 (1980), 110–118. doi: 10.1115/1.3138205
[16]	D. A. Torvi, J. D. Dale, A finite element model of skin subjected to a flash fire, J. Biomech. Eng., 116 (1994), 250–255. doi: 10.1115/1.2895727
[17]	W. Kaminski, Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure, J. Heat Transfer, 112 (1990), 555–560. doi: 10.1115/1.2910422
[18]	Y. Zhang, Generalized dual-phase lag bioheat equations based on nonequilibrium heat transfer in living biological tissues, Int. J. Heat Mass Transfer, 52 (2009), 4829–4834. doi: 10.1016/j.ijheatmasstransfer.2009.06.007
[19]	B. Mochnacki, E. Majchrzak, Numerical model of thermal interactions between cylindrical cryoprobe and biological tissue using the dual-phase lag equation, Int. J. Heat Mass Transfer, 108 (2017), 1–10. doi: 10.1016/j.ijheatmasstransfer.2016.11.103

This article has been cited by:

1.	Stanisław Kukla, Urszula Siedlecka, Mariusz Ciesielski, Fractional Order Dual-Phase-Lag Model of Heat Conduction in a Composite Spherical Medium, 2022, 15, 1996-1944, 7251, 10.3390/ma15207251
2.	Mohammed Sobhy, Ashraf M. Zenkour, Refined Lord–Shulman Theory for 1D Response of Skin Tissue under Ramp-Type Heat, 2022, 15, 1996-1944, 6292, 10.3390/ma15186292
3.	Mahmoud Ragab, Ahmed E. Abouelregal, Huda F. AlShaibi, Rasha A. Mansouri, Heat Transfer in Biological Spherical Tissues during Hyperthermia of Magnetoma, 2021, 10, 2079-7737, 1259, 10.3390/biology10121259
4.	Ashraf M. Zenkour, Tareq Saeed, Amal M. Aati, Refined Dual-Phase-Lag Theory for the 1D Behavior of Skin Tissue under Ramp-Type Heating, 2023, 16, 1996-1944, 2421, 10.3390/ma16062421
5.	Joana F. Soeiro, Filipa L. Sousa, Maria V. Monteiro, Vítor M. Gaspar, Nuno J. O. Silva, João F. Mano, Advances in screening hyperthermic nanomedicines in 3D tumor models, 2024, 9, 2055-6756, 334, 10.1039/D3NH00305A
6.	Hany H. Sherief, Mohamed F. Zaky, Mohamed F. Abbas, Samar A. Mahrous, Muhammad Mubashir Bhatti, Mathematical modeling of heat transfer in tissues with skin tumor during thermotherapy, 2024, 19, 1932-6203, e0298256, 10.1371/journal.pone.0298256
7.	Jyoti Pal, Panchatcharam Mariappan, S. Sundar, Application of Finite Pointset Method to Study Two‐Way Coupled Transient Bio‐Thermoelastic Effects in Skin Tissue, 2025, 4, 2702-4288, 10.1002/appl.70000
8.	Md. Alamgir Hossain, R. Nasrin, Eid S. Alatawi, Human Skin Burn Intensity Resulting From Various Incidents Utilizing Bioheat Transfer Model: A Comparative Analogy, 2025, 2688-4534, 10.1002/htj.23291

Reader Comments

Your name:*

Email:*
© 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

3.9

Metrics

Article views(3435) PDF downloads(249) Cited by(8)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(13)

Mathematical Biosciences and Engineering

Dual-phase lag model of heat transfer between blood vessel and biological tissue

Related Papers:

Abstract

1. Introduction

2. Methods

2.1. Dual-phase lag equation

2.2. Formulation of the problem

2.3. Method of solution

3. Results

3.1. Temperature distribution in the tissue with a tumor

3.2. Temperature distribution in the radial direction

3.3. Temperature of the blood vessel wall

3.4. Comparison of DPL model with the classical model for the blood sub-domain

4. Discussion

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

Dual-phase lag model of heat transfer between blood vessel and biological tissue

Related Papers:

Abstract

1. Introduction

2. Methods

2.1. Dual-phase lag equation

2.2. Formulation of the problem

2.3. Method of solution

3. Results

3.1. Temperature distribution in the tissue with a tumor

3.2. Temperature distribution in the radial direction

3.3. Temperature of the blood vessel wall

3.4. Comparison of DPL model with the classical model for the blood sub-domain

4. Discussion

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog