Thermo-viscoelastic analysis of biological tissue during hyperthermia treatment

Journal Pre-proof

Thermo-viscoelastic analysis of biological tissue during hyperthermia treatment Xiaoya Li , Qing-Hua Qin , Xiaogeng Tian PII: DOI: Reference:

S0307-904X(19)30667-5 https://doi.org/10.1016/j.apm.2019.11.007 APM 13132

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

15 April 2019 23 October 2019 4 November 2019

Please cite this article as: Xiaoya Li , Qing-Hua Qin , Xiaogeng Tian , Thermo-viscoelastic analysis of biological tissue during hyperthermia treatment, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.11.007

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights 

The dual phase lag thermo-viscoelastic model is developed in this work.



Transient thermo-mechanical responses of the tumor-normal tissues are investigated.



The effects of relaxation time and viscoelastic property are discussed.

1

Thermo-viscoelastic analysis of biological tissue during hyperthermia treatment

Xiaoya Li a

a, b

b

a

, Qing-Hua Qin , Xiaogeng Tian *

State Key Laboratory for Strength and Vibration of Mechanical Structures, Shaanxi Engineering Research Center of Nondestructive Testing and Structural Integrity Evaluation, School of Aerospace, Xi’an Jiaotong University, Xi’an, 710049, PR China;

b

Research School of Electrical, Energy and Materials Engineering, College of Engineering and Computer Science, Australian National University, Acton, ACT, 2601, Australia; * Corresponding author. Tel.: +86 029 82665420. E-mail address: [email protected](X.-G. Tian)

Abstract: Although viscoelastic properties of biological tissue has been reported in many articles, no effort has been made to investigate the coupled thermal and mechanical behavior of biological tissue based on the viscoelastic theory. This provides us a motivation to study the transient thermoelastic coupling response in the context

of

generalized

thermo-viscoelastic

model.

The

dual

phase

lag

thermo-viscoelastic model is established to capture the micro-scale responses of biological tissue. The governing equations are solved by Laplace transformation. The effects of relaxation times and viscoelastic property on the responses of the tumor and normal tissues are discussed and illustrated graphically. According to the numerical results, we can obtain (1) the viscoelastic parameter has a significant effect on the distributions of displacement and stress; (2) the lagging thermo-viscoelastic responses depend not only on the ratio of  t /  q , but also on the absolute values of  t and  q . 2

Keywords: Transient responses; Dual phase lag model; Viscoelastic theory; Relaxation time; Laplace transformation

1. Introduction Cancer is a major leading death cause in modern human life. It is characterized by the abnormal cells uncontrollably growing and spreading to the surrounding normal tissues. A lot of thermal therapies, such as hyperthermia, thermal ablation and microwave radiation have been widely applied in the cancer treatment [1]. One of the challenges in the thermal treatments is selectively destroying the tumor cells without causing damage to the surrounding healthy tissues. For this, understanding the temperature and stress distributions in biological tissue is very important for achieving a good quality in clinical treatment. Pennes’ bioheat transfer equation is the most commonly used to describe the temperature profile for its simplicity [2]. It is known that Pennes’ bioheat transfer equation was established based on the classical Fourier’s law, which predicts thermal signal propagating at an infinitely speed. However, some literatures [3-5] have shown that thermal behavior in non-homogenous media requires a longer relaxation time to accumulate enough thermal energy for transferring to the nearest element. Mitra et al. [6] carried out experiments on processed meat with different boundary conditions and observed the wave-like behavior in the heat conduction process. As a result, many researchers have been attracted to depict the wave-like behavior in biological tissue by using different bioheat transfer models. Liu et al. [7] used thermal wave model to describe the temperature distribution in living tissue. A dual phase lag (DPL) model was proposed by Tzou [8] to capture micro-scale response in both time and space. Liu et al. [9] investigated the temperature behavior of multi-layer skin tissue in the context of DPL bioheat transfer equation. Ahmadikia [10] used the Pennes model, thermal wave model and DPL bioheat transfer equations to investigate the effects of thermal relaxation times on the first- and second-degree burn time, respectively. 3

It is noted that even a small change in heat-induced stress can suppress immune response, alter the production of hormones and protein denaturation [11]. However, most studies mainly focus on the heat conduction [12-15], while the heat induced deformation is not considered. Shen et al. [16] studied the static thermo-mechanical responses of skin tissue at high temperature. Xu et al. [17-20] investigated the heat transfer, thermal damage and heat-induced stress of human skin. Nevertheless, it is assumed that the mechanical behavior has no effect on the temperature response in these studies. Iordana and Alexandru [21] have proposed that the elastic parameters of tumor tissue have an important effect on the magnetic nanoparticles distribution and temperature filed. The classical coupling theory of thermoelasticity was first established by Biot [22] and predicted an infinite speed propagation of thermal signal, which contradicts physical facts [23]. To eliminate such paradox, a number of generalized thermoelastic theories involving a finite speed of heat conduction have been drawn much attention in recent decades [24-29]. In the above work, the material parameters of biological tissue are assumed to be linear-elastic properties, which may not predict the real thermo-mechanical responses in biological tissue. Liu and Bilston [30] studied the linear viscoelastic properties of bovine liver via three different experiments. Ocal et al. [31] obtained the viscoelastic properties of tissue based on the general Maxwell model for viscoelastic materials. Majeed [32] analysed the heat transfer process in ferromagnetic viscoelastic fluid. Attar et al. [33] investigated the thermo-viscoelastic behaviour of tumorous and healthy bovine liver tissue without considering the inertia effect. In summary, all the previous studies were mainly focused on the temperature behaviour [12-15], even if there are themoelastic coupling response study, they only limit to the elastic theory [34, 35] or the thermoelastic equations based on the Pennes’ bioheat transfer equation [33, 36]. No report was found in the investigation of the transient thermoelastic coupling behaviour in biological tissue based on the generalized thermo-viscoelastic theory. Besides, to capture the micro-scale responses 4

in both time and space, the dual phase lag thermo-viscoelastic model is developed. The influences of relaxation times and viscoelastic property on the temperature, displacement, stress and tumor-normal tissue interface are obtained and presented graphically.

2. Problem formulation and governing equations In present work, a problem for superficial cancer via short-pulsed laser is studied [37, 38]. The coordinate system is chosen that the x-axis is perpendicular to the tumor surface, and the y-and z-axes is parallel to tumor surface, which is assumed to be infinitely long in y-and z-axes. In the therapy, the laser beam is perpendicular to the tumor surface and the size of spot is much larger than the thickness of the considered biological tissue for the time period of interest. In this process, the interface of tumor and normal tissues is assumed to be continuous. So the problem can be considered as one-dimensional case and all the functions in tumor and normal tissues will only depend on the location x and time t , as shown in Fig. 1. Additionally, the metabolic heat generation can be neglected compared with the high-intensity laser on the tumor surface. The sub-symbols i  1 and 2 indicate the variables in the tumor and normal tissues, respectively.

Fig. 1 Schematic geometry of tumor and normal tissues subjected to laser

The displacements have the following forms:

uxi  ui  x, t  ,u yi  uzi  0

(1)

Considering the geometrical relation, we can obtain: 5

ui , yyi   zzi  xyi   xzi   yzi  0 x u ei   xxi  i x

 xxi 



where εi   xx ,  yy ,  zz



T i

(2) (3)

is the strain, ei is the cubical dilatation.

It is assumed that the tumor and normal biological tissues are the linear isotropic viscoelastic media with constant material parameters in this study. No phase change and no chemical reactions occur within the biological tissue. The dual-phase-lag thermo-viscoelastic equations in biological tissue can be expressed as: (a) The constitutive equations [39, 40]:

 

t ei      xxi   d  0 Rvi  t     ei  atii  d 3  

(4)

t e 1 t  Ri (t   ) i d   Rvi  t     ei  atii  d  0 3 0  

(5)

t

 xxi   Ri (t   ) 0

 yyi   zzi  

i Si  ati  Rvi  t    t

0

where σi   xx ,  yy ,  zz 

T i

ei  c d  i i i  T0

(6)

is the stress; i and ci are density and specific heat;

Si is entropy density; ati is the linear thermal expansion coefficient; T0 is the reference temperature; Ti is the biological tissue temperature;  i is the temperature increment, i  Ti  T0 , i / T0  1 . And the relaxation functions Ri  t  and

Rvi  t  can be expressed as [41]: t Ri  t   2i 1  A1  e  t t  1dt  0  

(7)

t Rvi  t   Ki 1  A2  e  t t  1dt  0  

(8)

where Ki   i  2i / 3 , i ,i are Lame’s constants;   0    1 ,     0  ,

A1  A1  0  and A2  A2  0  are non-dimensional empirical constants; Ri  t   0,

6

Rvi  t   0,3Rvi  t   Ri  t   0 . Eqs. (4)-(6) will be reduced to linear thermoelastic equations when A1  0, A2  0 in Eqs. (7) and (8). (b) The energy conservation equation [2]: qxi S   iT0 i  bi wbi cbi (Tb  Ti ) x t

(9)

where qxi is heat flux; Tb is the blood temperature; bi ,wbi ,cbi are the blood mass density, blood perfusion rate and specific heat, respectively. It is assumed Tb  T0 in present study. The item bi wbi cbi (Tb  Ti ) describes the heat conduction between the blood and tissue. (c) The dual-phase-lag heat conduction equation [8]: qxi   qi

qxi     ki 1   ti  i t x  t 

(10)

where ki is thermal conductivity;  qi and  ti are the relaxation times in the heat flux and temperature gradient, respectively. When  qi   ti  0 , Eq. (10) will be reduced to the classical Fourier’s heat conduction equation; when  ti   qi the thermal wave tends to a diffuse wave; and the temperature distribution appears to be a wave-like behavior when  ti   qi . For the sake of brevity, qxi is replaced by qi in the flowing equations. (d) The motion equation without body force:  xxi  2ui  i 2 x t

(11)

Substituting Eqs. (4) and (5) into Eq. (11), we can obtain: t  2ui 2 t  3ui  i 2   Ri (t   ) 2 d   Rvi  t    0 0 t 3 x  

  2ui i   2  ati  d  x   x

(12)

Substituting Eqs. (6) and (10) into Eq. (9), the temperature governing equation can be expressed as: 7

 2ui 2     t ki 2 1   ti  i  T0 ati 1   qi   Rvi  t    d  x  t  t  t  0 x     i ci 1   qi  i  bi wbi cbi 1   qi  i t  t  t  

(13)

The Dirichlet temperature boundary condition is used in this study, which means the tumor is subjected to sudden heating and its surface

 x  0

is free of stress. The

heat and elastic waves cannot arrive at the right boundary of normal tissue during the considered time. So the boundary conditions are given as:

1  0, t   0 H  t  ; xx1  0, t   0

(14)

2  L, t   0;u2  L, t   0

(15)

The continuous interface conditions are used in the tumor-normal tissue interface ( h is the interface position):

1  h, t   2 (h, t )

(16)

u1  h, t   u2 (h, t )

(17)

 xx1  h, t    xx 2 (h, t )

(18)

q1  h, t   q2 (h, t )

(19)

To completely describe this problem, initial conditions should be given. The tumor and normal tissues are assumed to be in a quiescent initial state, thus the initial conditions have the following form: i ( x, 0)  0 t u ui  x, 0   0; i  x, 0   0 t

i ( x, 0)  0;

(20) (21)

Thus far, the problem considered in present work is proposed and all the theoretical models and conditions, such as the governing equations, initial and boundary conditions, as well as tumor-normal tissue interfacial conditions, are given systematically. For simplicity, the following non-dimensional variables are introduced: 8

 x , u , h , L   1l  x, u , h, L  ;t , *

* i

qi* 

*

*

*

i

l

qi ,Ri*  t  

k1 Tw  T0 

* i1

, i*2  

σ  v  t , i1, i 2  ;σ*i  i i*  i  l K1 Tw  T0

2 Ri  t  R t  ;Rvi*  t   vi ;v  3K1 K1

K1

1

where Tw is the temperature of thermal therapy on the tumor surface; l and v are the characteristic length and velocity, respectively.

The asterisk of the

non-dimensional variables is dropped in the following for the sake of brevity. Eqs. (4), (5), (12) and (13) can be rewritten in dimensionless forms as: t

 xxi   Ri (t   ) 0

 yyi   zzi  

t  2ui  d   Rvi  t    0 x 

 ui   x  Tw  T0  atii  d  

t  2ui 1 t   u  R ( t   ) d   Rvi  t     i  Tw  T0  atii  d i  0 2 0 x   x 

(22)

(23)

t t i  2ui  3ui      2ui  R ( t   ) d   R t    Tw  T0  ati i  d (24)   i vi  2 2 2   0 0 1 t x    x x 

ati vlT0    2ui  c vl 2    t 1      1   K R t   d  i i    ti qi  i   0 1 vi 2  x  t  ki Tw  T0  t  t  x ki

 w c l2       1   qi  i  bi bi bi 1   qi  i t  t  ki t  

where Ri  t  

(25)

t 4i K 1  A1 g  t   ,Rvi  t   i 1  A2 g  t  ,g  t    e  t t  1dt. 0 3K1 K1

The boundary and interfacial conditions (Eqs. (14)-(16)) can be expressed as:

1  0, t   H  t  ; xx1  0, t   0

(26)

2  L, t   0;u2  L, t   0

(27)

1  h, t   2 (h, t );u1  h, t   u2 (h, t )

(28)

 xx1  h, t    xx 2 (h, t );q1(h, t )  q2 (h, t )

(29)

The initial conditions can be rewritten as:

i ( x, 0)  0;

i ( x, 0)  0 t

(30) 9

ui  x, 0   0;

ui  x, 0   0 t

(31)

3. Solution of the governing equations The

Laplace

transformation

is

often

used

to

solve

these

complex

thermo-viscoelastic equations. Applying Laplace transformation

s f s   n



0

 n f  t   st e dt t n

(32)

to Eqs. (22)-(25), we can obtain:

 xxi  s  Ri  Rvi  

ui  sRvi Tw  T0  atii x

1

 u

(33)

 yyi   zzi  s  Rvi  Ri  i  sRvi Tw  T0  atii 2  x 

(34)

i  2u  sui   Ri  Rvi  2i  Rvi Tw  T0  ati i 1 x x

(35)

 2i u 1   ti s  2  bi1  s   qi s 2 i  bi 2 1   qi s i bi 3s 2 1   qi s  i x x

(36)

where

bi1  i ci vl / ki ,bi 2  bi wbi cbi l 2 / ki ,bi 3  ati K1RviT0vl / ki / Tw  T0  , Ri  t  

   4 i  1 K 1    A1 g  t   ,Rvi  t   i   A2 g  t   , g  t   .   3K1  s K1  s   ss   

Similarly, we apply Laplace transformation to Eqs. (26)-(29):

1  0, s  

0 s

; xx1  0, s   0

(37)

2  L, s   0;u2  L, s   0

(38)

1  h1 , s   2 (h1 , s);u1  h1 , s   u2 (h1 , s)

(39)

 xx1  h1 , s    xx 2 (h1 , s);q1(h1 , s)  q2 (h1 , s)

(40)

To investigate the effect of relaxation times on the interface condition, the continuous heat flux interface condition (Eq.(19)) is discussed in terms of the Fourier

10

and non-Fourier heat flux interface conditions. The continuous Fourier heat flux interface condition can be expressed as:

1

k1

x

 k2

x h

 2 x

(41) x h

The continuous non-Fourier heat flux interface condition with the DPL heat conduction equation (Eq.(10)) can be written as:

k11  s t1  1 k 1  s t 2   2  2 1  s q1 x x h 1  s q 2 x x h

(42)

Eliminating  i between Eqs. (35) and (36), we can obtain:

(

4 2  p  pi 2 )ui  0 i1 x 4 x 2

(43)

where pi1  i s / 1 /  R i  R vi   bi 3 s 2 ati Rvi Tw  T0  1   qi s  / p3   bi1s  bi 2  1   qi s  / 1   ti s  ; pi 2  i s  bi1s  bi 2  1   qi s  / 1 / pi 3 ; pi 3   R i  R vi  1   ti s  ;

Thus, the solution of Eq. (43) can be expressed as: 4

ui   uˆij e m x  j

i

(44)

j 1

where mij  i  1, 2; j  1, 2,3, 4  are the roots of the characteristic equations in tumor and normal tissues:

m 

j 4

i

 pi1  mij   pi 2  0 2

(45)

Accordingly, the temperature and stress in tumor and normal tissues can be expressed as: 4

i  ˆi j e m x j

(46)

i

j 1

4

 xxi   ˆ xxij e m x j

i

(47)

j 1

11

4

 yyi   zzi   ˆ yyij e m x j

(48)

i

j 1

j j j According to Eqs. (33)-(35), ˆi ,ˆ xxi ,ˆ yyi  i  1, 2; j  1, 2,3, 4 

can be

expressed as follows:

ˆi j  di j uˆij ,di j 

2 i s   Ri  Rvi   mij  1

Rvi Tw  T0  ati mij

ˆ xxij  bi j uˆij ,bi j  s  Ri  Rvi  mij  sRvi Tw  T0  ati di j  

 

1

j ˆ yyi  fi j uˆij , fi j  s  Rvi  Ri  mij  sRvi Tw  T0  ati di j  2

(49)

(50)

(51)

Substituting Eqs. (46)-(51) into the boundary and interface conditions (Eqs. (37) -(42)), we can obtain: 4

 d uˆ

j j 1 1

j 1

4

 b uˆ

j j 1 1

j 1

4

0 s

(52)

0

(53)



 d uˆ e j 1

4

j j j  m2 L 2 2

 uˆ e j 1

j j  m2 L 2

  d uˆ e 4

j j j  m1 h 1 1

j 1

  uˆ e 4

j 1

j j  m1 h 1

 b uˆ e 4

j 1

j j j  m1 h 1 1

  k m d uˆ e 4

j 1

1

j j j j  m1 h 1 1 1

0

(54)

0

(55)



 d 2j uˆ2j e m2 h  0 j

(56)



 uˆ2j e m2 h  0 j

(57)



 b2j uˆ2j e  m2 h  0 j

(58)



 k2 m2j d 2j uˆ2j e  m2 h  0 j

(59) 12

 k11  s t1  j j j  m j h k2 1  s t 2  j j j  m j h  m1 d1 uˆ1 e 1  m2 d 2 uˆ2 e 2   0     1  s  1  s  j 1  q1 q2  4

(60)

The coefficients uˆij  i  1, 2; j  1, 2,3, 4  can be obtained by solving Eqs. (52) -(60). Thus far, the solutions in Laplace domain are obtained. To obtain the solutions in time domain, an efficient numerical inversion of Laplace transform (NILT) algorithm [42] is adopted in this work, which is based on fast Fourier transformation.

4. Numerical results 4.1. Verification of numerical results Ramadan [43] obtained the temperature distribution in a multilayered media by using DPL heat conduction equation. To check the validity of the solutions in the present work, the results of Ramadan [43] are compared with the model of this study which ignores the coupled thermo-viscoelastic effect. The thermal parameters, geometries, dimensionless quantities, initial and boundary conditions are the same as those reported in [43]. Fig. 2 shows the temperature distribution of our DPL heat conduction model and Ramadan. Great agreements with each other are achieved, which indicate that the present numerical method is accurate and effective. And the dimensionless quantities of t ,x and  in Fig. 2 correspond to the parameters  , and  in [43].

13

Fig. 2 Comparisons with Ramadan [43] on the distribution of temperature

4.2. Influence of viscoelastic property on responses The material constants of tumor and normal tissues in the following calculation are shown in Table 1 [15, 33, 44]. The non-dimensional thickness of tumor tissue is h  0.2 , the non-dimensional length of the model is L  1.2 . Table 1 Material constants of biological tissue

Tumor

Normal Tissue Blood

1  2.481107 kg/m s 2

1  1.034 106 kg/m s2

1  1660kg/m3

c1  2540J/kgK

k1  0.778W/mK

wb1  0.009 / s

h  0.015m

at1  1104 /K

A1  0.106

2  8.27 107 kg/m s2

2  3.446 106 kg/m s2

2  1000kg/m3

c2  3720J/kgK

k2  0.642W/mK

wb 2  0.00018 / s

L  0.09m

at 2  1104 /K

A2  0.08

b  1060kg/m3

cb  3770J/kgK

Tb  T0  310K

Fig. 3-Fig. 6 show the results of dual phase lag thermo-elastic (DPL-TE) and dual phase lag thermo-viscoelastic (DPL-TVE) models at different instants to investigate the influence of viscoelastic parameter on the transient responses. The

14

classical Fourier heat flux interface condition (Eq. (41)) is used and  ti /  qi  0.5 in this case. From Fig. 3, we can see that the viscoelastic parameter has little effect on temperature distribution. The absolute values of displacement and stress in DPL-TVE model are smaller than DPL-TE model as time goes on. And there is a jump of stress at the tumor-normal tissue interface in Fig. 5 and Fig. 6.

Fig. 3 The influence of viscoelastic parameter on temperature

15

Fig. 4 The influence of viscoelastic parameter on displacement

Fig. 5 The influence of viscoelastic parameter on stress  xx

16

Fig. 6 The influence of viscoelastic parameter on stress  yy

4.3. Influence of relaxation times on responses Fig. 7-Fig. 10 show the influences of the ratio of  ti /  qi on the distributions of temperature, displacement and stress at the instant t  0.3 in the context of dual phase lag thermo-viscoelastic model. We can see that the larger ratio of  ti /  qi is, the farther the thermal and elastic waves arrive at, and the larger the absolute values of the maximum displacement and stress  yy  zz  are.

17

Fig. 7 Temperature distributions with different ratios of  ti /  qi

Fig. 8 Displacement distributions with different ratios of  ti /  qi 18

Fig. 9 Stress  xx distributions with different ratios of  ti /  qi

Fig. 10 Stress  yy distributions with different ratios of  ti /  qi 19

Tzou [8] has proposed that the ratio of  ti /  qi dominates the lagging behavior of heat transfer. In other words, the thermal responses are identical for a specified value of  ti /  qi . Comparisons of temperature variations under the same ratio



ti

/  qi  5.0  are performed for three different cases to check the suitability of this

assertion in these cases. From Fig. 11-Fig. 14, we can clearly see that the distributions of temperature, displacement and thermal stress are not exactly the same in these three cases, which implies that the lagging thermo-viscoelastic responses in biological tissue depend not only on the ratio of  ti /  qi but also on the absolute magnitudes of

 ti and  qi .

Fig. 11 The influence of relaxation times on temperature

20

Fig. 12 The influence of relaxation times on displacement

Fig. 13 The influence of relaxation times on stress  xx 21

Fig. 14 The influence of relaxation times on stress  yy

4.4. Influence of heat flux interface condition on responses Fig. 15-Fig. 18 show the influences of classical Fourier and non-Fourier heat flux interface conditions on the distributions of temperature, displacement and thermal stress at the instant t  0.3 in the context of dual phase lag thermo-viscoelastic model. It can be see that the absolute values of temperature, displacement and thermal stress in the Fourier heat flux interface condition are larger than that of in the non-Fourier heat flux interface condition when  ti /  qi  0.2 . But the opposite results can be observed when  ti /  qi  5.0 .

22

Fig. 15 The influence of Fourier and non-Fourier heat flux on temperature

Fig. 16 The influence of Fourier and non-Fourier heat flux on displacement 23

Fig. 17 The influence of Fourier and non-Fourier heat flux on stress  xx

Fig. 18 The influence of Fourier and non-Fourier heat flux on stress  yy 24

5. Discussion and conclusions In present work, the dual-phase-lag thermo-viscoelastic model is developed and the transient thermo-viscoelastic responses in the tumor and normal tissues are investigated. From the numerical results, one may conclude that: viscoelastic parameter has little effect on temperature, but the absolute values of displacement and stress in thermo-viscoelastic model are smaller than thermoelastic model, which indicates that under the same conditions of thickness and time, the thermal deformation damage of thermoelastic model is more serious than that of the thermo-viscoelastic

model

during

hyperthermia

treatment;

the

lagging

thermo-viscoelastic responses in biological tissue depend not only on the ratio of

 t /  q but also on the absolute values of  t and  q ; the absolute values of temperature, displacement and thermal stresses in the non-Fourier heat flux interface condition are larger than that of Fourier heat flux interface condition when

 ti /  qi  1.0 ; but the opposite results can be observed when  ti /  qi  1.0 . The limitation of the model in this study is assumed that the tumor and normal tissues have a linearly viscoelastic properties. This allows us to use a simpler linear dual phase lag thermo-viscoelastic model to describe the transient thermo-mechanical response. An obvious extension of the model is to develop the non-linear generalized thermo-viscoelastic model to approach the realistic viscoelastic properties of biological tissue as close as possible. But this work is an attempt to provide a thermo-viscoelastic theoretical framework to help researchers understand the complex thermo-mechanical phenomena that occur in thermal therapies.

Acknowledgements This work is supported by the National Science Foundation of China (11732007, 11572237) and the Fundamental Research Funds for the Central Universities and the State Scholarship Fund from the China Scholarship Council (CSC). 25

Reference [1] W. Agnelli, C. Padra, C.V. Turner, Shape optimization for tumor location, J. Comput. Math. Appl. 62 (11) (2011) 4068-4081. [2] H.H. Pennes, Analysis of tissue and arterial blood temperatures in the resting human forearm, J. Appl. Physiol. 1 (2) (1948) 93-122. [3] A.V. Luikov, Analytical Heat Diffusion Theory, Academic Press, New York, 1968. [4] A.M. Braznikov, V.A Karpychev, A.V Luikova, One engineering method of calculating heat conduction process, Inzh. Fiz. Zh. 28 (4) (1975) 677-680. [5] W. Kaminski, Hyperbolic heat conduction equation for material with a non-homogenous inner structure, ASME J. Heat. Transf. 112 (3) (1990) 555-560. [6] K. Mitra, S. Kumar, A. Vedavarz, et al., Experimental evidence of hyperbolic heat conduction in processed meat, ASME J. Heat Transf. 117 (3) (1995) 568-573. [7] Liu J., Ren Z.P., C.C. Wang, Interpretation of living tissue’s temperature oscillations by thermal wave theory, Chinese Sci. Bull. 40 (17) (1995) 1493-1495. [8] D.Y.

Tzou, Macro-to micro-scale heat transfer: the lagging behavior,

Washington: Taylor and Francis, 1997. [9] K.C. Liu, Y.N Wang, Y.S Chen, Investigation on the bio-heat transfer with the dual-phase-lag effect, Int. J. Therm. Sci. 58 (2012) 29-35. [10] H. Askarizadeh, H. Ahmadikia, Analytical study on the transient heating of a two-dimensional skin tissue using parabolic and hyperbolic bioheat transfer equations, Appl. Math. Model. 39 (13) (2015) 3704-3720. [11] C.P. Lau, Y.T. Tai, P.W. Lee, The effects of radio frequency ablation versus medical therapy on the quality-of-life and exercise capacity in patients with accessory pathway-mediated supraventricular tachycardia: a treatment comparison study, Pace. 18 (3) (1995) 424-32. [12] K.C. Liu, H.T. Chen, Investigation for the dual phase lag behavior of bio-heat transfer, Int. J. Therm. Sci. 49 (7) (2010) 1138-1146. 26

[13] L. Cao, Q.H. Qin, N. Zhao, An RBF–MFS model for analysing thermal behaviour of skin tissues, Int. J. Heat Mass Tran. 53 (7) (2010) 1298-1307. [14] D. Kumar, P. Kumar, K.N. Rai, A study of DPL model of heat transfer in bi-layer tissues during MFH treatment, Comput. Biol. Med. 75 (2016) 160-172. [15] K.C. Liu, H.T. Chen, Analysis of the thermal response and requirement for power dissipation in magnetic hyperthermia with the effect of blood temperature, Int. J. Heat. Mass Tran. 126 (2018) 1048-1056. [16] W. Shen, J. Zhang, F. Yang, Modeling and numerical simulation of bioheat transfer and biomechanics in soft tissue, Math. Comput. Model. 41 (11-12) (2005) 1251-1265. [17] F. Xu, T.J. Lu, K.A. Seffen, Biothermomechanical behavior of skin tissue, Acta Mech. Sin. 24 (1) (2008) 1-23. [18] F. Xu, K.A. Seffen, T.J. Lu, Non-Fourier analysis of skin biothermomechanics, Int. J. Heat Mass Tran. 51 (9-10) (2008) 2237-2259. [19] F. Xu, T.J. Lu, K.A. Seffen, Biothermomechanics of skin tissue. J. Mech. Phys. Solids 56 (2008) 1852-1884. [20] F. Xu, T.J. Lu, Introduction to skin biothermomechanics and thermal pain, New York: Science Press, 2011. [21] A. Iordana, S. Alexandru, Advanced thermo-mechanical analysis in the magnetic hyperthermia, J. Appl. Phys. 122 (16) (2017) 1-13. [22] M. Biot, Thermoelasticity and irreversible thermo-dynamics, J. Appl. Phys. 27 (3) (1956) 240-253. [23] V. Peshkov, Second sound in helium II, J. Phys. 8 (1944) 381-382. [24] H.W. Lord, Y. Shulman, A generalized dynamic theory of thermoelasticity, J. Mech. Phys. Solids 15 (5) (1967) 299-309. [25] A.E. Green, K. Lindsay, Thermoelasticity, J. Elast. 2 (1) (1972) 1-7. [26] A.E. Green, P.M. Naghdi, Thermoelasticity without energy dissipation, J. Elast. 31 (3) (1993) 189-208. 27

[27] M. Marin, A. Öchsner, The effect of a dipolar structure on the Holder stability in Green-Naghdi thermoelasticity, Continuum Mech. Therm. 29 (6) (2017) 1365-1374. [28] S. Chiriţă, M. Ciarletta, V. Tibullo, Qualitative properties of solutions in the time differential dual-phase-lag model of heat conduction, Appl. Math. Model. 50 (2017) 380-393. [29] M. Marin, Cesaro means in thermoelasticity of dipolar bodies, Acta Mech. 122(1-4) (1997) 155-168. [30] Z. Liu, L. Bilston, On the viscoelastic character of liver tissue: experiments and Mmodeling of the linear behavior, Biorheology 37 (3) (2000) 191-201. [31] S. Ocal, M.U. Ozcan, I. Basdogan, et al., Effect of preservation period on the viscoelastic material properties of soft tissue with implication for liver transplantation, J. Biomech. Eng.-T ASME 132 (2010) (1-7). [32] A. Majeed, A. Zeeshan, S.Z. Alamri, et al., Heat transfer analysis in ferromagnetic viscoelastic fluid flow over a stretching sheet with suction, Neural Comput. Appl. 30 (6) (2018) 1947-1955. [33] M.M. Attar, M. Haghpanahi, H. Shahverdi, et al., Thermo-mechanical analysis of soft tissue in local hyperthermia treatment, J. Mech. Sci. Technol. 30 (3) (2016) 1459-1469. [34] X.Y. Li, C.L. Li, Z.N Xue, et al., Analytical study of transient thermo-mechanical responses of dual-layer skin tissue with variable thermal material properties, Int. J. Therm. Sci. 124 (2018) 459-466. [35] X.Y. Li, C.L. Li, Z.N Xue, et al., Investigation of transient thermo-mechanical responses on the triple-layered skin tissue with temperature dependent blood perfusion rate, Int. J. Therm. Sci. 139 (2019) 339-349. [36] A. McBride, S. Baraman, D. Pond, et al., Thermoelastic modelling of the skin at finite derformations, J. Therm. Biol. 62 (2016) 201-209.

28

[37] J. Jiao, Z. Guo, Thermal interaction of short-pulsed laser focused beams with skin tissues, Phys. Med. Biol. 54 (13) (2009) 4225-4241. [38] A.Y. Sajjadi, K. Mitra, M. Grace, Ablation of subsurface tumors using an ultrashort pulse laser, Opt. Lasers Eng. 49 (3) (2011) 451-456. [39] M.A. Ezzat, A.S. El-Karamany, The uniqueness and reciprocity theorems for generalized thermo-viscoelasticity with two relaxation times, Int. J. Eng. Sci. 40 (2002) 1275-1284. [40] M.A. Ezzat, A.S. El-Karamany, A.A. El-Bary, Thermo-viscoelastic materials with fractional relaxation operators, Appl. Math. Model. 2015, 39 (23-24): 7499-7512. [41] M.A. Ezzat, A.S. El-Karamany, A.A. Samaan, Sate space formulation to generalized thermo-viscoelasticity with thermal relaxation, J. Thermal Stresses 24 (9) (2001) 823-846. [42] L. Brancik, Programs for fast numerical inversion of Laplace transforms in MATLAB language environment, Proceedings of the 7th conference MATLAB’ 99 Czech Republic Prague 1999, pp: 27-39. [43] K. Ramadan, Semi-analytical solutions for the dual phase lag heat conduction in multilayered media, Int. J. Therm. Sci. 48 (1) (2009) 14-25. [44] M. Mohajer, M.B. Ayani, H.B. Tabrizi, Numerical study of non-Fourier heat conduction in a biolayer spherical living tissue during hyperthermia, J. Therm. Biol. 62 (2016) 181-188.

29