Meshless Method with Operator Splitting Technique for Transient Nonlinear Bioheat Transfer in Two-Dimensional Skin Tissues

A meshless numerical scheme combining the operator splitting method (OSM), the radial basis function (RBF) interpolation, and the method of fundamental solutions (MFS) is developed for solving transient nonlinear bioheat problems in two-dimensional (2D) skin tissues. In the numerical scheme, the nonlinearity caused by linear and exponential relationships of temperature-dependent blood perfusion rate (TDBPR) is taken into consideration. In the analysis, the OSM is used first to separate the Laplacian operator and the nonlinear source term, and then the second-order time-stepping schemes are employed for approximating two splitting operators to convert the original governing equation into a linear nonhomogeneous Helmholtz-type governing equation (NHGE) at each time step. Subsequently, the RBF interpolation and the MFS involving the fundamental solution of the Laplace equation are respectively employed to obtain approximated particular and homogeneous solutions of the nonhomogeneous Helmholtz-type governing equation. Finally, the full fields consisting of the particular and homogeneous solutions are enforced to fit the NHGE at interpolation points and the boundary conditions at boundary collocations for determining unknowns at each time step. The proposed method is verified by comparison of other methods. Furthermore, the sensitivity of the coefficients in the cases of a linear and an exponential relationship of TDBPR is investigated to reveal their bioheat effect on the skin tissue.


Introduction
It is widely recognized that accurate and fast prediction of temperature distribution in biological tissues is important in various practical diagnostics. For example, doctors would like to know the heat and temperature changes during surgery on a skin tumor or the human eye so that they can adjust the power of the laser therapy to avoid extra burning injury of healthy tissue [1,2]. Due to the complexity of the bioheat system itself, numerical methods have been overwhelmingly utilized for noninvasive diagnostics. Several linear or nonlinear steady-state bioheat models involving changed thermal conductivity and blood perfusion rate have been numerically solved to analyze the induced temperature distribution in biological tissues [3][4][5][6][7][8]. Besides, a non-Fourier heat conduction model in one-dimensional multilayered systems was analyzed by Laplace transform and the fast inversion technique [9][10][11]. In this paper, a model for describing transient nonlinear bioheat transfer model in two-dimensional (2D) skin tissue is developed.
In the context of transient nonlinear bioheat transfer, when there is a need to dynamically monitor the changes of temperature in time and space during the bioheat transfer process, some numerical models have been developed for various biological tissues. For instance, Arunn et al. [12] investigated the variation of transient temperature of a 2D human eye computational model using the finite volume method. Trakic et al. [13] predicted the transient temperature rise in a nonlinear heat transfer model of tumor and healthy tissue of mouse by the commercial finite element software FEMLAB. Feng et al. [14] applied finite element technology coupled with a nested-blocked optimization algorithm to predict the temperature distribution in a prostate during a nanoshell-mediated laser surgery.
Boundary-type methods, which are different from the domain-type methods above, have also been presented to solve such problems. Among them, the dual reciprocity boundary element formulation has been proposed for analysis of the transient nonlinear 2D bioheat transfer system subjected to a sinusoidal heat flux on the skin surface [15]. An axisymmetric boundary element formulation using the time-dependent fundamental solution was derived by Majchrzak for the analysis of freezing and thawing processes in biological tissues [16]. Alternatively, Cao et al. [17] developed a mixed meshless method RBF-MFS by coupling the radial basis function (RBF) and the method of fundamental solutions (MFS) to analyze the transient linear thermal behavior in skin tumor tissue. It is noted that in these numerical methods, the Laplace transform method or finite difference technology with respect to time were applied to handle the time variable in the bioheat transfer governing equation. However, the Laplace transform method is usually limited to linear transient problems [18]. For other cases, the finite difference scheme needs carefully consideration of the time step length to obtain accurate, stable, and convergent results [19]. These difficulties have motivated researchers to develop other methods for effectively handling the time derivative term and for approximating the nonlinear source terms in the governing equation. For example, the operator splitting method for uniformly handling the transient term and nonlinear source term explicitly using two-level higher-order time step schemes has received considerable attention [20,21].
In this work, we aim to develop a mixed meshless method for analyzing transient nonlinear bioheat transfer in 2D skin tissue by way of the operator splitting method. In the proposed solution procedure, the operator splitting method is employed first to isolate the transient and nonlinear terms in the original Penney bioheat governing equation by the explicit second-order Adams-Bashforth time marching for the half time step and the second-order Adams-Moulton scheme for the next time step. Then, the new equation in the form of a modified Helmholtz equation can be derived and solved at each time step. Next, the mesh-free dual reciprocity method implemented by the RBF interpolation and the mesh-free MFS in terms of fundamental solution kernels are respectively utilized to determine the particular solutions and the homogeneous solutions of the modified Helmholtz problem at each time step by simple internal and boundary collocations. The advantages of this meshless scheme are that only the boundary integrals are included in the calculation process. Thus the computational time is less than the corresponding finite element method. Secondly there will be no singular integrals generated in the process, because the source points and field points are placed outside the solution domain respectively. There is no complicated and time-consuming mesh generation required during the process. However, the disadvantage of the proposed meshless scheme is that it is difficult to deal with multi-material problems, compared to the conventional finite element method based on feasible element material definition.
The paper is organized as follows: In Section 2, the numerical results from the proposed numerical method are verified by comparison with those from ANSYS software. Sensitivity analysis of the effect of coefficients on temperature distribution in the case of temperature-dependent blood perfusion rate is performed in the same section. In Section 3, the mathematical bioheat model of 2D skin tissue is described and in Section 4, the solution procedure including the operator splitting method, the dual reciprocity method, and the MFS is presented. Finally, some conclusions are drawn in Section 5.

Verification of the Proposed Method
In this section, the efficiency and accuracy of the proposed method for analyzing transient nonlinear bioheat transfer in 2D skin tissue are validated by the finite element software ANSYS through a benchmark example. In this work, we consider the nonlinear term induced by the blood perfusion rate, which is a function of the tissue temperature. The thermal parameters of the 2D skin tissue model used in the calculation are given in Table 1 [17,22,23]. The ANSYS transient thermal toolbox is employed to simulate bioheat transfer in the biological material. The mesh generated by ANSYS is shown in Figure 1, in which 647 elements and 733 nodes are generated for finite element analysis.
For the purpose of comparison, we consider that blood perfusion rate is a linear function of tissue temperature: 1 2 ( ) b T a a T ω = + , where 1 0.0005 a = and 2 0.0001 a = . 63 interpolation points and 32 boundary collocations (see Figure 2) are used to calculate the transient temperature distribution. Numerical results along the x-axis at three time instants ∆t = 50, 80, and 100 s are presented in Figure 3 to show the accuracy and stability of the second-order Adams-Bashforth and Adams-Moulton schemes. From Figure 3, it can be seen that the results from the proposed algorithm with fewer collocation points are in good agreement with the results from the ANSYS Transient thermal toolbox. The relative error of the results from the proposed method with respect to those from the Transient thermal toolbox of ANSYS is less than 0.5%.   Next, the exponential case of the temperature-dependent blood perfusion rate  Further, Figure 5 presents the temperature variation from t = 0 s to t = 2560 s at the point (1.875 mm, 0) on the x-axis for the case of a linear blood perfusion rate. It can be seen from Figure 5 that the variation of temperature with time from the proposed meshless method is almost identical to that obtained from ANSYS, although much less unknowns are used in the proposed method.
Thus, the convergence and accuracy of the present meshless method with the higher-order AB and AM time-stepping schemes is validated for transient nonlinear bioheat analysis in the rectangular model of skin tissue. More numerical results are now presented to illustrate temperature distribution in the solution domain caused by different temperature-dependent blood perfusion rates. In Figures 6 and 7, the temperature distribution in the skin tissue along the x-axis at different times is presented. It is found that the steady state in Figure 6 is reached much earlier (the linear case, at about 1600 s) than that in Figure 7 (the exponential case, at about 8000 s). It is also noted that the slope of the steady-state temperature curve along the x-axis increases and then decreases from the left side to the right side for both linear and exponential cases. However, the slope for the linear case appears greater than that for the exponential case in the region close to the left surface, which has a lower environmental temperature, whereas the slope for the linear case becomes less than that for the exponential case in the region close to the right surface, which has a higher body core temperature. Moreover, the exponential-form blood perfusion rate produces a higher interior temperature in the region close to x = 18.75 mm than that for the linear-form rate. The main reason is that the exponential-form blood perfusion rate generally has a lower value of the blood perfusion rate than the linear-form with the coefficients given above. In the region close to the left surface, where the skin tissue temperature is evidently lower than the blood temperature, the greater blood perfusion rate means that more heat flows from blood to skin tissue, causing a rapid increase of the tissue temperature. Thus there is greater temperature gradient in this region for the linear case than the exponential case. When the tissue temperature exceeds the blood temperature, a greater blood perfusion rate causes more heat to flow from tissue to blood and causes the tissue temperature to decrease.

Sensitivity of Temperature to Variation of Constants ai in Linear Case of Blood Perfusion Rate
In this section, the linear case of temperature-dependent blood perfusion rate 1 2 ( ) b T a a T ω = + is considered for the sensitivity analysis of temperature to the constant coefficients 1 a and 2 a . Firstly, the coefficient a2 is set to be constant 0.0002 and the constant 1 a is assumed to be 0.005, 0.0005, and 0.00005, respectively. As we can see in Figure 8, the steady-state temperatures are quite close to each other when 1 0.00005 a = and 1 0.0005 a = . But the skin temperature curve has a relative larger gap with the two curves mentioned above when 1 0.005 a = . In addition, the three temperature curves intersect at the point (12.65 mm, 0), at which the skin temperature is approximately 37 °C. Therefore, from the left surface of the skin tissue to the approximate location point (12.65 mm, 0), the greater blood perfusion rate indicates that more heat flow transfer occurs between the blood and skin tissue. If the blood temperature is higher than the skin tissue, more heat flow transferred from blood to skin tissue causes a rapid increase in skin temperature. If the blood temperature is lower than the skin tissue, more heat flow transferred from skin tissue to blood causes a rapid decrease in skin temperature. It can be seen that blood perfusion protects the skin tissue from extreme temperature increases or decreases caused by the environment. To study the effect of a2 on skin temperature, we assume the first constant 1 a to be 0.0005 and the second constant 2 a is set as 0.00002, 0.0002, and 0.002. From Figure 9, it can be seen that variation of the constant 2 a causes a more rapid change in the steady-state temperature curve than that due to variation of the constant 1 a . In particular, when the constant 2 a = 0.002, the curve of the tissue temperature is steeper than the other two curves. The highest value of the skin temperature appears at the approximate location (7.5 mm, 0), which is closer to the left hand side boundary of the skin tissue than the other two curves. From the location (15 mm, 0) to (26.25 mm, 0), the skin tissue temperature is stable at a certain level when the constant 2 a = 0.002.

Sensitivity of Temperature to Variation of ai in the Exponential Case of Blood Perfusion Rate
In this section, the sensitivity analysis of temperature to the constant coefficients i a is investigated by considering the exponential case of temperature-dependent blood perfusion rate When the constant 2 a is assumed to be 0.01, the constant 1 a is respectively tested at 0.005, 0.0005, and 0.00005.
Compared with the linear case shown in Figure 8, the difference or gap between each skin temperature curve is relative larger, as shown in Figure 10. Similarly, the three temperature curves with different values of constant 1 a intersect at almost the same point (the distance from the left hand side boundary being roughly 13.125 mm). This finding means that, at the location (13.125 mm, 0), the skin temperature has almost the same value of 37.75 °C for different values of the constant 1 a . Figure 10 illustrates the stronger regulatory and protective effect of the exponential-form blood perfusion rate than that in the linear case (see Figure 8).  Again, we assume constant a1 to be 0.0005, while constant 2 a is set to be 0.03, 0.01, and 0.003. As we can see from Figure 11, when constant 2 a = 0.03, the temperature of the skin tissue becomes increasingly steep before the point (11.25 mm, 0), but the curve is flatter than temperature curves with smaller values of 2 a . Compared with the effect of the different values of 1 a in Figure 10, the increase in the value of 2 a causes a larger reduction of the peak value of the skin tissue temperature and the temperature becomes more stable from the location (11.25 mm, 0) to (26.25 mm, 0). In summary, an increase in the value of constant 2 a has a higher sensitivity to the temperature of skin tissue than an increase in the value of constant 1 a . Simultaneously, it is found that an increase in the blood perfusion rate causes the temperature of the skin tissue to reach its final steady state more quickly and reduces the peak value of the tissue temperature. That means that if the skin tissue absorbs a large quantity ofbiological heat from its environment, the blood perfusion effect causes the temperature to reach a certain value quickly and reduces the risk of burning of the skin tissue. Figure 11. Sensitivity to constant a2 in the exponential case of blood perfusion rate.

Mathematical Bioheat Transfer Model in 2D Skin Tissue
The transient heat transfer in biological tissue is usually governed by the well-known Pennes's bioheat transfer equation [7,24]: where T is the temperature, ρ the tissue density, c the tissue specific heat, k the tissue thermal conductivity, b ρ the blood density, b c the blood specific heat, b ω the blood perfusion rate, b T the arterial temperature, r Q the spatial heat sources, m Q the metabolic heat generation rate, t the time, and 2 ∇ the standard Laplacian operator.
In practice, blood flow accelerates with the increase in temperature of the environment. Thus, blood perfusion rate can be viewed as a function of tissue temperature. In this case, the governing Equation (1) can be written in the form of nonlinear equation as follows: In hyperthermia treatment, the blood perfusion is usually assumed to vary linearly [15,22,25]: or exponentially [5,22,26]: with the tissue temperature T. 1 a and 2 a are positive constants.
For the sake of convenience, we introduce a new temperature variable θ: then, the nonlinear governing Equation (2) can be rewritten in terms of the new variable as: represents the generalized interior heat source term including the metabolic heat of the tissue and the spatial heat source caused by laser heating. Further, Equation (7) can be expressed in the general unsteady Poisson equation form as: Besides the governing Equation (9), boundary conditions of the problem that describes bioheat transfer in a rectangular skin domain as shown in Figure 12 include [4,15,27]: (1) the temperature condition at the right boundary 1 Γ , that is assumed to be the body core temperature c θ : (2) the thermal insulation conditions at the upper boundary 2 Γ and the lower boundary 3 Γ : 2 3 0 at boundaries and T k n ∂ − = Γ Γ ∂ (12) (3) the temperature condition at the left boundary 4 Γ , that is assumed to be constant: 4 at boundary Moreover, the initial condition of the problem is given by: The governing Equation (9), the boundary conditions (11)- (13), and the initial condition (14) consist of a complete partial differential equation system, which will be solved by the meshless method developed in the next section.

The Operator Splitting Method
Making use of the concept of operator splitting [20], the time-dependent governing Equation (9) can be expressed as a sum of two operators 1 L and 2 L : From Equation (15), a solution in time by a two-level time-stepping scheme is used in this work. Typically, the second-order Adams-Bashforth scheme [20]: and the second-order Adams-Moulton scheme [20]: are, respectively, employed to model the nonlinear operator 2 L and the Laplacian operator 1 L .
In Equations (18) and (19), If a new variable * θ defined by: is introduced, Equation (21) can be transformed to: that can be rearranged in the form: Equation (24) is a type of modified Helmholtz equation and * θ is a generalized function to be determined at each time step. The right nonhomogeneous term in Equation (24)  can be obtained through Equation (22). Unlike the backward time-stepping scheme, this scheme needs the function values at step (n) and the previous step (n-1). Therefore, it cannot start by itself. The function value at the first time step can be evaluated by the extrapolated explicit forward Euler scheme presented below [20] 1 0 t t Then we have: If we set the assumed initial guess 0 θ , the value 1 θ for the first time step can be calculated according to Equation (26). Furthermore, the time iteration can be commenced from Equation (24). For the sake of simplicity, Equation (24) is rewritten as: and ( ) ( ) As well, the boundary conditions (11)-(13) should be modified for the time iteration so that a complete partial differential equation system can be formed with the conjunction of the governing Equation (27) and the modified boundary conditions below.

Solution of the Modified Helmholtz System
In this subsection, the dual reciprocity method (DRM) using RBFs and the MFS using fundamental solutions are used to solve the modified Helmholtz equation system (27)- (30). Both methods are based on boundary or internal collocation and have been successfully applied to similar nonhomogeneous problems [17,[28][29][30]. The solution method is described in detail below.
First, the DRM is introduced by simply setting: and then Equation (27) can be expressed as the following nonhomogeneous Laplace equation: According to the linear feature of the Laplacian operator, the solution to Equation (32) can be expressed as [28,29]: where N is the number of boundary collocation points, j β are source intensity and ( ) ( , ) is the fundamental solution to the linear Laplacian operator [35]: and has the form: Finally, with the obtained particular and homogeneous approximations, the full solution can be written in the form: The normal derivative of the full solution can then be given by: For the purpose of simplicity, Equations (43) and (44) are written in the matrix form: where ( 1,2,3,4) i N i = are respectively the number of collocation points on the four edges of the rectangular domain (see Figure 2) and 1 The unknown coefficient vector c can be determined from linear equation system (50), and then the temperature variable * θ at each time step can be calculated from Equation (43) Figure 2 shows an illustration of the 32 collocations, 32 source points and 63 interpolation points for the half rectangular domain in our calculation.

Conclusions
In this paper, an operator splitting technique coupled with the dual reciprocity method and the method of fundamental solutions is presented to develop a mesh-free algorithm for solving the transient nonlinear bioheat transfer in a 2D model of skin tissue with a temperature-dependent blood perfusion rate. Use of the operator splitting technique, including two-level second-order time-stepping schemes, makes it possible to establish an accurate and convergent solution procedure for transient and nonlinear cases, and then the dual reciprocity method and the method of fundamental solutions are respectively employed to solve the obtained modified Helmholtz equation system at each time step. This meshless method is dependent on the internal interpolating points and boundary collocation points of the domain only and thus is really meshless and dimension-independent. The numerical results demonstrate the accuracy and efficiency of the meshless method in the analysis of the transient nonlinear bioheat transfer problem under consideration, with very few interpolation and collocation points. Moreover, the sensitivity analysis of temperature sensitivity to the constants in the linear and exponential expressions of the blood perfusion rate demonstrates the increase in the constant 2 a in the linear case. It is found that the exponential case has a more significant influence on the tissue temperature distribution than the constant 1 a , and an increase in its value results in a relatively fast increase in the tissue temperature in the region close to the outer surface and simultaneously, the peak temperature value decreases. This reflects the regulation and protection effect of the blood perfusion rate in biological tissue.