Numerical model of biological tissue heating using the models of bio-heat transfer with delays

. The numerical model of thermal processes in domain of biological tissue subjected to an external heat source is discussed. The model presented is based on the second order dual-phase-lag equation (DPLE) in which the relaxation time and thermalization time (τ q and τ T ) are taken into account. In this paper the homogeneous, cylindrical skin tissue domain is considered. The most important aim of the research is to compare the results obtained using the classical model (the first-order DPLE) with the numerical solution resulting from the higher order form of this equation. At the stage of numerical computations the Finite Difference Method (FDM) is applied. In the final part of the paper the examples of computations are shown.


Introduction
The problem of thermal processes occurring in the domain of skin tissue subjected to an external heat source is discussed. Recently there is the view that, taking into account the specific internal structure of the tissue, the hyperbolic type equations better than parabolic ones reproduce the actual course of the thermal processes taking place in the domain considered, e.g. [1][2][3][4]. In this paper the heat transfer in the tissue is described by the single, second-order DPLE (a homogeneous domain). In the further stages of the research, the system of these equations (the multilayered skin tissue domain) will be considered. The energy equation contains the additional components (the source functions) related to the blood perfusion and metabolism [5][6][7] just like the classic Pennes equation [8]. The starting point for the considerations concerning the DPLE is (as one knows) the generalized form of the Fourier law in which the lag times (the relaxation and thermalization times) are introduced. The left and right hand sides of generalized Fourier law are developed into the Taylor series with accuracy to the first derivative and next, after using this development in the energy equation, the first-order DPLE equation can be obtained. In this paper the second-order equation is considered, in particular the Taylor series with the accuracy to the second derivative is taken into account (e.g. [9]). The external tissue heating is determined by the Neumann boundary condition (heat flux is given), the other boundary and initial conditions will be presented in the next chapter. The numerical algorithm presented in this paper is based on the FDM (the implicit scheme). Taking into account the form of assumed external heat flux, the axiallysymmetrical problem is considered. The original computer program that implements the numerical calculations has been prepared in such a way that it can be used for the first-and second-order equations as well as for the mixed variants. For example, in the several works -e.g. [10], the second-order Taylor expression of heat flux and the first-order Taylor expression of temperature gradient are applied to take into account the phase lagging behavior. The system of equations corresponding to the transition from t to t+∆t is solved using the iterative Gauss-type procedure. In the final part of the paper the examples illustrating the differences in the solutions resulting from the applications of the first-and the second-order models are presented.

Governing equations
The starting point for the formulation of the energy equation with delays is the generalized Fourier law. In particular, the relationship between heat flux q and the temperature gradient T  is given in the form where λ is a thermal conductivity, τ q and τ T are the relaxation time and thermalization time, respectively, X = {r, z} and t denote spatial co-ordinates (cylindrical coordinates) and time.
Using the Taylor series expansions, the following secondorder approximation of the formula (1) can be taken into account Depending on the assumed number of the components, the different forms of DPLE can be obtained. Using the known diffusion equation (in the case of the constant values of thermophysical parameters) one obtains where c is a volumetric specific heat, a=λ/c is a thermal diffusivity, Q is a capacity of internal heat sources, w q and w T which are equal to 0 or 1 determine the order of DPLE. In particular (0, 0) corresponds to the first order equation, (1, 1) corresponds to the second order equation, while (1, 0) or (0, 1) correspond to the mixed variants.
The internal heat source according to the Pennes theory [9] is a sum of two components where G B [m 3 blood/m 3 tissue/s] is a perfusion coefficient, c B is a volumetric specific heat of blood, T B is an arterial blood temperature, Q met is a metabolic heat source (treated here as a constant value). The basic assumption leading to the adopted form of the perfusion heat source is that the tissue domain is supplied by a large number of capillary blood vessels uniformly distributed in the domain considered. Taking into account the formula (4) one has The equation (5) can be written in the form T r z t T r z t T r z t r r r r z The energy equation (6) must be supplemented by the boundary and initial conditions. In the case of the problem considered, the boundary conditions take a form of the Dirichlet and Neumann ones. Thus, on the bottom of the cylinder (z=Z) the body core temperature T b is given. On the top of the cylinder (z=0) outside the domain of external heat source action and on the lateral surface (r=R) the no-flux condition is assumed. For r ≤ R 0 (radius of external heat flux action) the boundary heat flux is equal to q b (Figure 1). The condition discussed is of the form where ( , , ) T r z t  n denotes a normal derivative. One can see, that for the constant value of q b (r, z, t) and also for the no-flux condition the equation (8) where T p is an initial temperature, while u(r, z) and v(r, z) are the known functions. For the first-order DPLE the function v is not defined.

Numerical solution
The algorithm presented below is based on the implicit scheme of the finite difference method (FDM). Let , ( , , ) z j =jh (h is the mesh step in r and z directions) and f = 0, 1,…, F. Taking into account the initial conditions (9), under the assumption that u(r,z)=v(r,z)=0, one has the approximate form of equation (6) resulting from the introduction of the adequate differential quotients is as follows ,   2  1  2  ,  ,  ,  2   2  1  2  3  , , , where (c.f. formula (7))   , 1    In the model presented, on the external surface of the system the Neumann boundary condition (8) is taken into account. The FDM form of this condition is the following The problem is solved using the implicit scheme of the finite difference method under the assumption that the grid step is equal to h=0.0003 m (n=100, number of nodes 10201) and time step is equal to ∆t=0.01 s. Further compaction of the differential mesh did not change the results of numerical computations. In Figures 2, 3 and 4 the heating /cooling curves at the points (0, 0) and (0, 0006 m) are shown. The first-order DPLE solution is marked by a continuous line, while the second-order one is marked by a dashed line. In the case of the second-order solution, the higher maximum tissue temperatures are obtained, and these maxima are shifted to the right compared to the first-order DPLE solution. Differences in the solutions are not too large, but in the case of the tasks related to the heating of biological tissue they become important (e.g. in the case of the burns modeling or the simulation of hyperthermia treatment).        The computations above presented have been repeated for the other values of the lag times, namely τ q =2 s and τ T = 2 s (Figures 8, 9, 10). These significantly different from the previously assumed values of the delay times result from the considerations presented in the works [11,12] in which τ q and τ T are dependent mainly on the porosity of the tissue. Delays determined in this way are shorter than the most frequently reported in the literature (e.g. [5,13]) and are similar to each other. The differences between both solutions are visible, although all geometric and thermophysical parameters have been retained (except for delay times).

Conclusions
The basic goal of the paper was to compare the results of numerical solutions relating to the heating of biological tissue subjected to an external heat source. Two variants of dual-phase-lag equation have been taken into account, in particular the first-and the second-order DPLE have been considered. To achieve this aim, an algorithm based on the implicit scheme of the Finite Difference Method has been developed. The 3D axially-symmetrical tissue domain has been considered. A large external radius of the domain was assumed and then on the lateral surface of the cylinder the adiabatic condition could be accepted. The system of equations associated with the transition from time t to time t+∆t was solved using the Gaussian iteration method. Differences between solutions corresponding to the firstand the second-order DPLE are clearly visible. In particular, they have a place in the case of the significant external source capacity (more than 1 kW/m 2 ) and the rather short exposure time. Typical heating / cooling curves for the process under consideration and the small values of delays reach the maximum values on the heated surface (e.g. at the point (0, 0)) for t ≈ t exp . In the inner points (e.g. (0, 0.0005 m)) of the cylinder, the extreme value appears after a slightly longer time. The similar situation takes place for the larger values of delay times.
The assumption concerning the values of delay times is very important for the results of numerical simulations.
Comparison of the solutions obtained for τ q =15 s and τ T = 5 s and the results corresponding to τ q =2 s and τ T = 2 s confirms the qualitative consistency of the results, but the numerical values are strongly different for both the first-and second-order DPLE. It can also be seen that for similar values of delay times the differences between solutions (see: Figures 8 -10) are negligible. Summing up, it should be emphasized that the developed computer programs work correctly but their practical usefulness is determined by the adoption of the appropriate values of delay times. Further research will be carried out to take into account the layered structure of the skin tissue [12,14] and, as a consequence, the modeling of boundary conditions at the interface between the skin layers.