NUMERICAL MODELING OF HEAT TRANSFER IN BIOLOGICAL TISSUE DOMAIN USING THE INTERVAL FINITE DIFFERENCE METHOD

In the paper the numerical analysis of heat transfer process proceeding in the domain of a biological tissue is presented. In particular, the two-dimensional problem is considered, in which the thermophysical parameters (volumetric specific heat and thermal conductivity) are given as intervals. The problem discussed has been solved using the interval finite difference method using the rules of directed interval arithmetic. In the final part of the paper the results of numerical computation are shown.


INTRODUCTION
Thermophysical parameters of biological tissue (thermal conductivity, volumetric specific heat, perfusion coefficient etc.) can change in the wide range.For instance, the value of dermis thermal conductivity (the 'second layer' of skin tissue) can be assumed as a number from 0.37 W/mK to 0.52 W/mK.It results from the fact that the tissue parameters depend on the numerous individual traits such as the age, sex, occupation etc.
So, the paper concerns imprecisely defined transient bio-heat transfer problems, when in the mathematical description the uncertain parameters are defined and treated as directed interval numbers (e.g.[1]).The base of mathematical model is given by the Pennes interval set of equations supplemented by the boundary-initial conditions.
The first part of the paper is devoted to a short presentation of governing equations describing the heat transfer processes proceeding in biological tissue domain (continuous models [2] are taken into account).
Next, the generalization of FDM algorithm in the case of directed interval values of tissue parameters are presented, at the same time the approach close to the control volume algorithm (proposed by Mochnacki [3]) is applied.
In the final part of the paper the examples of numerical simulations are presented.In particular the solutions concerning the thermal processes proceeding in the skin tissue domain subjected to external heat source is shown.

GOVERNING EQUATIONS
Thermal processes proceeding in the heterogeneous skin tissue domain (Figure 1) can be described by the following system of interval energy equations where e = 1, 2, 3 corresponds to the successive layers of skin (epidermis, dermis, subcutaneous region -as in Figure 1  In Figure 2 the tissue domain oriented in cylindrical co-ordinate system is shown.
where [ , ] Be Be G G − + is the interval perfusion coefficient, B c is the volumetric specific heat of blood, B T is the arterial blood temperature, [ , ] me me is the interval metabolic heat source.
Interval equations (1) should be supplemented by the boundary and initial conditions.So, the skin surface ( 0 r R ≤ ) is subjected to an external heat source, for 0 r R > the boundary condition of the 3 rd type is assumed, while for the others parts of the boundary the no-flux conditions are taken into account , where b q is the given interval boundary heat flux, α is the heat transfer coefficient, a T is the Between the successive sub-domains the continuity condition can be taken into account The initial condition is also given The equations ( 1) -( 5) create the mathematical model of the process discussed.
The problem formulated has been solved by means of interval finite difference method using the rules of directed interval arithmetic [1].In this arithmetic the set of proper intervals is extended by improper intervals: it is possible to obtain the number zero by subtraction of two identical intervals and the number one as the result of the division of two identical intervalswhich was impossible applying classical interval arithmetic.The directed interval arithmetic allows one to obtain the temperature intervals much narrower than the classical interval arithmetic does and the intervals width does not increase in time.

NUMERICAL MODEL
Numerical model of thermal processes proceeding in domain of heating tissue bases on the finite difference method in the version presented in [4,5].At first, the time grid is introduced with a constant step ∆t .The geometrical mesh is shown in Figure 3.One can see that the 'boundary' nodes are located at the distance 0.5h or 0.5k with respect to the real boundary (h, k are the steps of regular mesh in directions r and z), respectively.This approach gives the better approximation of the Neumann and Robin boundary conditions [4].The operator ( ) for the axially-symmetrical problem at the internal node (i, j) and time f-1 (the explicit differential scheme is considered) is of the form ( ) The FDM approximation of expression (7) at the node (i, j) (see: Figure 4) is constructed in the similar way as in [4].We introduce the nodes marked by 'empty circles' and the approximation of differential operators using the mean quotient is applied Interval thermal conductivities 1 , 0.5 are assumed in the form of harmonic means of values corresponding to the basic nodes and then in formulas (8), (9) the classically defined thermal resistances 1 , 1 − between the central node (i, j) and nodes (i, j+1), (i, j−1) appear .5 0.5 0.5 , The nodes (i, j), (i, j+1), (i, j−1) can belong to the different sub-domains [4,5].It should be pointed out that the similar formulas can be obtained in the case of 'boundary' nodes.Let us assume that the 'boundary' node is located at the distance of 0.5h from the external boundary on which the Robin condition is given.Then the interval thermal resistance 1 , 1 Additionally in place of 1 , 1 + the ambient temperature should be introduced.For a very small value of heat transfer coefficient, e.g.
) the no-flux condition can be taken into account.The problem of expression (7) approximation in the case of Neumann boundary conditions is discussed, in details, in [4. 5].
Using again the mean differential quotient one obtains The same approach can be used for the second term of formula (7).In particular and next .5 0.5 0.5 , , According to the rules of explicit differential scheme for transition t f−1 → t f one obtains the interval FDM equation in the form From which the temperature at the point (i, j) for time level f can be found.The stability conditions [4] must be fulfilled, of course.All these interval values must be calculated according to the rules of the directed interval arithmetic [1]. (e = 1, 2, 3), the time of external heat source exposition has been assumed as 5 s.In Figure 5 the heating and cooling curves at the selected nodes are presented.
Figure 5. Heating and cooling curves at nodes 1 ( ) , L R and 4 ( ) Figure 7 illustrates a comparison between the heating and cooling curves obtained using interval FDM and the results obtained using classical FDM for thermophysical parameters defined without intervals (dashed line).

FINAL REMARKS
The application of interval FDM gives the solution for which the heating/cooling curves at the set of points selected from the tissue domain are obtained in a 'fuzzy' form and the differences between border courses are visible.A such information seems to be very interesting and allows one to look at the process more closely.One can see that the solution obtained using the classical approach fits in the interior of interval and the decrease of intervals causes the decrease of the distance between border curves.
the capacity of interval internal heat sources, ( , ) e T x t , x and t denote the temperature, spatial co-ordinates and time.