Finite Element Analysis of Nonlinear Bioheat Model in Skin Tissue Due to External Thermal Sources

: In this work, numerical estimations of a nonlinear hyperbolic bioheat equation under various boundary conditions for medicinal treatments of tumor cells are constructed. The heating source components in a nonlinear hyperbolic bioheat transfer model, such as the rate of blood perfusions and the metabolic heating generations, are considered experimentally temperature-dependent functions. Due to the nonlinearity of the governing relations, the ﬁnite element method is adopted to solve such a problem. The results for temperature are presented graphically. Parametric analysis is then performed to identify an appropriate procedure to select signiﬁcant design variables in order to yield further accuracy to achieve efﬁcient thermal power in hyperthermia treatments.


Introduction
Recent research developments have shown that the problem of heating transfers in skin tissues is an incredibly complicated problem. There are many controversial parameters found in this field. To describe the current status [1], at temperatures above 60 • C, the time needed to check irreversible damages decreases rapidly, while, at temperatures around 40-45 • C, irreversible tissue damage happens after only prolonged exposures. A human body behaves differently under different environmental conditions (e.g., different air temperatures, humidity levels, and wind velocities). Thermotherapy operations, such as laser tissues welding [2], hyperthermy [3], and lasers operations [4], have been openly used in modern medicine. Since the temperature distribution in the living tissues depends on complex phenomena, like metabolic heating generations and blood circulation, investigators have expanded some formulations. In clinical therapy, various contemporary thermotherapeutic techniques have been widely used for microwave, laser, ultrasound, and radiofrequency technologies. First, using an objective thermal lens, the laser focuses on the tumor. One of the main challenges in the treatment process is to provide adequate thermal power to diseased tissue without affecting healthy tissue. Therefore, it is important to consider the impact of temperature and stress fields on heat treatment kinetics. Van and Gybels [5] demonstrated that deformation due to heating and cooling can also contribute to a feeling of discomfort. Therefore, exact predictions of the Sun, mechanical reaction, and thermal damage in organic tissues are required for the treatment planning and development of new clinical heating systems. In 1948, Pennes [6] investigated the thermal behavior of skin temperature in the forearm. The investigation included many phenomenological mechanizations, such as metabolic heat generation, radiation, blood perfusion, thermal conduction, and phases change. In biological tissues, the phase changes occur in wide ranges. The adjusted Penne's bioheat models are presented by various techniques of numerical approach available in the literature: the homotopy perturbation method [7], the Legendre wavelet Galerkin approach, the finite element approaches, the Galerkin approach, and the approach of variations and iterations [8,9], as well as the finite-decomposition approach [10]. Hobiny et al. [11] presented analytical estimations of temperature in biological tissue under TPL bioheat model with experimental verifications. Dillenseger and Esneault [12] studied temperature improvement over time by using the approach of finite differences in the case where the body temperature is abnormally low. Hobiny and Abbas [13] studied the temperature increment and variations of thermal damages in living tissues due to mobile heating sources. Zhu et al. [14] considered the rate process models for the resulting thermal injury and the sedimentations of lighting energy in tissues by using the diffusion theory. Diaz et al. [15] applied the finite element scheme to solve the thermos-diffusions model type in the biological tissues to determine the thermal damages. Alzahrani and Abbas [16] investigated the analytical solution of thermal damage in living tissue due to laser irradiations. In considering the study of real phenomenon regarding thermal transfers in finite mediums, the linear/nonlinear models of heat transfers have been expanded, and their numerical or analytical solutions have been solved by several authors [17][18][19][20][21][22][23][24][25][26][27][28][29].
This work explores the numerical solutions of the thermal injury of biological tissue using a nonlinear bioheat theory in the living tissue. Due to the nonlinearity of the governing relations, the finite element method is used in order to determine solutions for this problem. The numerical outcomes for temperature behaviors are provided graphically. Furthermore, a comparison is made with the results obtained in the presence and absence of delay times. Additionally, comparisons between the nonlinear and linear influences on bioheat transfers are made when the nonlinear parameter are zeroes.

Mathematical Model
We consider that the medium is a semi-infinite living tissue under thermal insulation. Based on Cattaneo [30], with the delay time τ o , the bioheat relation in skin tissues can be given by the following equation [31,32]: where τ o is the thermal delay time, ρ is the mass density of tissues, T b is the temperature of blood, t is the time, T is the temperature of tissue, c is the tissue-specific heat, k is the tissue thermal conductivity, and Q ext represents the thermal generated per unit volume of tissue due to the absorption of electromagnetic radiations, which is expressed by the following equation [7]: where P is the transmitted power, S and η are the antenna constants, x p is the probe radius, x = L − x is the distance of tissues from external surface, and Q b refers to the blood perfusion heat source, which is expressed as follows: where ω b (T) is the blood perfusion rate dependent on temperature, ρ b is the mass density of blood, and c b is the specific heat of blood. Depending on the anatomical structure of the biological living tissue containing blood vessels, this point will extend with the increasing temperature of the local tissues. The perfusion rate of blood with the temperature of local tissues is presumed as in [33] as follows: where ω bo is the reference blood perfusions rate, α is the constant associated with blood perfusions, T o is the initial temperature of normal tissues, and Q m is the heat generated by metabolics process. Mitchell et al. [33] noticed that the metabolic heat generations is a function of the temperature of local tissues and is expressed as follows: where β is the associated metabolic constant and Q mo is the metabolic heat source reference.

Application
The appropriate initial and boundary conditions are defined by the following: For its adaptations, the non-dimensional variables can be expressed by the following: In terms of this nondimensional parameters (8), the formulation (1) as well the initial (6) and the boundary (7) conditions can be inserted as (the scripts have been abandoned for its suitability) as follows:

Numerical Scheme
In this section, the highly nonlinear formulation of bioheat transfer based upon thermal delay time τ o in biological tissues is summarized using the finite element method (FEM). This technique is a powerful approach that is primarily used to obtain the numerical solution of complex problem and the finite element method is the selected method for nonlinear systems in several fields. The finite element scheme is applied here to obtain the solution of nonlinear formulation (9) under the initial condition (10) and the boundary (11) condition. Abbas and his colleagues [34][35][36][37][38][39][40][41][42] presented solutions for various problems under deference generalized thermoelasticity models. The finite element formulation of nonlinear bioheat transfer can be obtained by using the standard procedure. The non-dimension weak formulations of bioheat model are derived. The set of independent test functions resulting in the temperature δT is specified. The governing formulation is multiplied by independent weighting functions and then the boundary conditions are used to integrate over the spatial domain. The applications of integration by parts and the use of the divergence theorem to decrease the order of the spatial derivative paves the way for the application of the boundary condition. Using the Galerkin procedure, the unknown temperature T and the corresponding test function are approximated by the same shape function.
where N refers to the shape function and n refers to the node numbers per element. The quadratic element with three nodes is used. In this problem, a one-dimensional quadratic element is used, and we suppose that the local coordinates in the range [−1, 1] are used for the master element, which can be defined as follows: Furthermore, the time derivatives of the unknown variables should be determined by the implicit approaches. Thus, the weak formulations for the finite element scheme corresponding to the basic Equation (9) can be written as follows: Symbolically, the discredited equation of Equation (15) can be given as follows: where (16) can be expressed as follows:

Numerical Results
The temperature distribution in a biological tissue under nonlinear bioheat transfer models was investigated. For numerical computation, exemplary values of thermal properties for biological tissue are given below (as in [43]): The calculations were made by the MATLAB software (R2018a), and the outcomes are graphically presented. The impacts of the external heating sources on the skin's surface were integrated. In this study, it was considered that T b = 37 • C is the normal temperature and a slab of tissues is 0.05 m thick. The grid sizes were refined and, consequently, the values of various parameters were stabilized. Further refinement of mesh sizes over 18,000 elements did not change the values considerably, therefore the mesh sizes were accepted as the grid sizes for computing purposes. In Figures 1 and 2 Figures 1 and 2. Figures 3 and 4 display the influences of thermal delay time on the temperature distribution. As expected, the thermal delay time has a major impact on the temperature distributions. Figures 1 and 3 show the temperature variations with respect to the distances x at t = 1.5 min. It is observed from the graphs that the temperature increases from the normal temperature T b = 37 • C to the maximum values on the tissue surface. A timeline of the surface temperature is presented in Figures 2 and 4. It was observed that after starting from the normal temperature T b , temperature increases with time until the highest temperature values are obtained. Finally, Figures 5 and 6 display the effects of the rate of blood perfusion ω b under bioheat model with one relaxation time on the temperature variation. The higher the rate of blood perfusion, the greater the convective heat loss due to faster blood flow. 021, 9, x FOR PEER REVIEW 5 of 9 time on the temperature variation. The higher the rate of blood perfusion, the greater the convective heat loss due to faster blood flow.   time on the temperature variation. The higher the rate of blood perfusion, the greater the convective heat loss due to faster blood flow.   x FOR PEER REVIEW 6 of 9   x FOR PEER REVIEW 6 of 9

Conclusions
Based on a highly nonlinear bioheat model of skin tissues, temperature behaviors were studied when applying thermotherapy. Numerical solutions were derived for the nonlinear bioheat with one relaxation time under an external heat source. The finite element method was applied to solve the second-order nonlinear differential equation by using a quadratic element. The theoretical outcomes can be used as confirmation tools in order to examine practical operations like laser treatment.

Conclusions
Based on a highly nonlinear bioheat model of skin tissues, temperature behaviors were studied when applying thermotherapy. Numerical solutions were derived for the nonlinear bioheat with one relaxation time under an external heat source. The finite element method was applied to solve the second-order nonlinear differential equation by using a quadratic element. The theoretical outcomes can be used as confirmation tools in order to examine practical operations like laser treatment.

Conclusions
Based on a highly nonlinear bioheat model of skin tissues, temperature behaviors were studied when applying thermotherapy. Numerical solutions were derived for the nonlinear bioheat with one relaxation time under an external heat source. The finite element method was applied to solve the second-order nonlinear differential equation by using a quadratic element. The theoretical outcomes can be used as confirmation tools in order to examine practical operations like laser treatment. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.