Fundamental-Solution-Based Hybrid Element Model for Nonlinear Heat Conduction Problems with Temperature-Dependent Material Properties

The boundary-type hybrid finite element formulation coupling the Kirchhoff transformation is proposed for the two-dimensional nonlinear heat conduction problems in solids with or without circular holes, and the thermal conductivity of material is assumed to be in terms of temperature change. The Kirchhoff transformation is firstly used to convert the nonlinear partial differential governing equation into a linear one by introducing the Kirchhoff variable, and then the new linear system is solved by the present hybrid finite element model, in which the proper fundamental solutions associated with some field points are used to approximate the element interior fields and the conventional shape functions are employed to approximate the element frame fields. The weak integral functional is developed to link these two fields and establish the stiffness equation with sparse and symmetric coefficient matrix. Finally, the algorithm is verified on several examples involving various expressions of thermal conductivity and existence of circular hole, and numerical results show good accuracy and stability.


Introduction
Compared to conventional materials, materials with temperature-dependent thermal conductivities usually serve in the high-temperature environment, for example, refractory materials for blast furnace, and thus the thermal analysis is important for such materials.However, the temperaturedependent feature of material properties leads to the nonlinearity of the governing equation, and thus the difficulties of the derivation of analytical solutions increase.
Different to analytical solutions, numerical results can be easily obtained by numerical methods for such nonlinear heat transfer problems, for example, the finite element method (FEM) [1,2], the boundary element method (BEM) [3][4][5], or the dual-reciprocity BEM (DRBEM) [6,7], the method of fundamental solution (MFS) [8,9], and the meshless element free Galerkin method [10,11].As an alternative to numerical approaches mentioned above, the fundamentalsolution-based hybrid finite element method, named as HFS-FEM, was initially developed for linear heat transfer problem [12] and then was extended to complicated thermal analysis in composite structures [13] and biological tissues [14].More recently, the elastic analysis was also performed by means of the present HFS-FEM and several types of special elements were developed for problems associated with circular or elliptical holes, nonhomogeneous materials, and discontinuous loads [15][16][17].Generally, as one of the domaintype methods based on element division in the domain, the HFS-FEM has some advantages over the conventional FEM.For instance, arbitrary-shape elements including element boundary integrals only can be constructed in the HFS-FEM, and thus special-purposed elements can be developed for problems with local defects like holes and inclusions [13,15,17] to achieve the effort of significant mesh reduction.Besides, in contrast to the boundary-type BEM or DRBEM, the HFS-FEM based on element discretization of the domain of interest can be easily applied to multidomain problems and each element can have itself material definition during the analysis.However, in the formulation of BEM, the treatment of interface continuity conditions between adjacent subdomains can highly weaken the symmetry of the coefficient matrix of the final solving equation in the formulation of BEM.
In this study, the hybrid finite element formulation with fundamental solution kernels (HFS-FEM) will be extended to the nonlinear heat conduction problems with temperaturedependent thermal conductivity, by coupling the noniteration technique, Kirchhoff transformation [18].During the computation, the Kirchhoff transformation is firstly employed to remove the temperature dependence of thermal conductivity, and then a new linear system in terms of the introduced Kirchhoff variable can be obtained.Subsequently, the fundamental-solution-based hybrid finite element formulation is presented to solve this new linear system.Once the Kirchhoff variable distribution in the domain under consideration is determined, the inverse transformation can be available to derive the desired temperature distribution.
A brief outline of the paper is arranged as follows: the mathematical models including the basic equations of heat conduction, the Kirchhoff transformation, and the hybrid finite element formulation are described in Section 2, and numerical results are presented and discussed in Section 3. Finally, conclusions are summarized in Section 4.

Mathematical Models
2.1.Basic Equations.In the paper, the two-dimensional steady-state heat conduction in isotropic media is taken into consideration with reference to the Cartesian coordinate system ( 1 ,  2 ).The heat equilibrium equation in the absence of internal heat generation is expressed as [18] where Ω denotes the computing domain bounded by a simple closed curve Γ, and the quantity   is the heat transfer rate along the th direction governed by the Fourier's law In (2),  is the temperature change and  is the thermal conductivity of the material.In the study, the thermal conductivity is assumed to be temperature-dependent, that is,  = (), and thus the governing equation (1) shows nonlinearity.
Besides, the following boundary conditions: should be complemented to (1) to obtain a complete mathematic system.In (3),  and , respectively, stand for the given distributions on the boundary

Kirchhoff Transformation.
Usually, for the case that the material property is dependent on the temperature, that is,  = (), necessary iteration procedures, for instance, Newton iteration, Nash iteration, Picard iteration, and so on, can be employed to perform linearization for such nonlinear problems.In the present work, a noniteration technology is employed by introducing the Kirchhoff transformation such as [18] d d =  () or  = ∫  () d, where  is a new quantity named as Kirchhoff variable, which is related to the sought temperature field .
Then, substituting (4) into (2), one obtains As a result, the nonlinear governing equation ( 1) and boundary conditions (3) reduce to the following linear system consisting of the linear partial differential equation associated with the introduced Kirchhoff variable : and the transformed boundary conditions In this paper, two types of material nonlinearities are considered.

Regular mesh
Irregular mesh
Integrating (8) by the Kirchhoff transformation (4) gives from which the inverse manipulation produces

Exponent-Type Thermal Conductivity.
Consider the following: Integrating ( 11) by the Kirchhoff transformation (4) gives from which the inverse manipulation produces 2.3.Hybrid Finite Element Formulation.Currently, the hybrid finite element with fundamental solution as trial function has been successfully formulated to solve the linear heat transfer problems using general or special elements.In this section, the HFS-FEM is applied to solve the new linear system consisting of ( 6) and ( 7) to determine the induced Kirchhoff variable.
The overall computing domain is firstly discretized with some elements.For a typical element , the weak integral functional can be written as follows [12,13]: where  and φ stand for independent element interior and boundary fields, respectively.Within the element , the assumed element interior field , also named as intraelement field, is usually approximated by the linear combination of the fundamental solution where and x and x  are field point and source field, respectively.  is the number of source points x  located outside the element domain, and   presents unknown source intensity.
In the paper, two different fundamental solutions are involved.One is the general fundamental solution associated with the problem without hole, and another is the special fundamental solution associated with the problem with circular hole.For completeness, we present, in the Appendix, the involved expressions of fundamental solutions.According to the physical definitions of fundamental solutions given in the Appendix, we can easily find that the linear combination of fundamental solutions with different source points can exactly satisfy the linear partial differential equation ( 6) and the specified interfacial condition along the circular hole.This feature is important to simplify the hybrid functional ( 14) by removing the domain integral in it.
Besides, the independent frame field defined over the element boundary is constructed by where {d  } and { Ñ} represent the element nodal degree of freedom (DOF) and conventional interpolating shape functions, respectively.Subsequently, applying the Gaussian divergence theorem to (14), we have the following formula: in which only element boundary integrals are included.The substitution of ( 15) and ( 17) into (18) finally yields where with The stationary conditions of the functional Π  with respect to {c  } and {d  } produce in which the stiffness equation and the optional relationship between {c  } and {d  } are well established.

Numerical Examples
In this section, some numerical examples are presented to illustrate the validity of the proposed hybrid finite element formulation with fundamental solution kernels.There are few nonlinear problems that have the known exact solutions.Among them is the first example in Section 3.1.The numerical results computed by the present method are compared with the exact solutions.The second example in Section 3.2 is a problem with a centered circular cut, and the numerical results by the present algorithm are compared with those by the MATLAB PDE Toolbox.

Nonlinear Heat Transfer in a Unit
Square.As a first example, let us consider the numerical solutions for plane heat transfer over a unit square domain.The temperatures on the left and right edges of the square remain 300 K and 400 K, respectively, while the remainders are assumed to be insulated, as shown in Figure 1.For the case of power-type variation of the thermal conductivity with temperature, that is,  =  0 ( + )  , an analytical solution of the form is obtained to verify the numerical solutions.In the practical computation, two different mesh schemes (regular mesh and irregular mesh) are employed to model the square domain and validate the stability of the proposed approach (see Figure 2).For this purpose, we first set  0 = 1,  = 1,  = −2, and  = 0.01, which corresponds to  = −2 + 0.01 and  = 100√3 1 + 1 + 200. Figure 3 displays the variation of the error defined by with the parameter , which controls the location of source point used for intraelement approximation.It can be seen that there are large stable range to choose the parameter , for either regular or irregular meshes.At the same time, one also observes that too small values of , meaning that the source points lie so closely to the element boundary that the distance  between the source points and element nodes approaches to zero, bring bad results.The main reason is that the kernel functions of the fundamental solution vary according to (ln ), and ( −2 ), respectively, in two dimensions so that inaccuracies in the evaluation of element boundary integrals (near-singular problems) would be caused when  approaches to zero.On the other hand, too large value of  corresponding to the large distance of the source points and element nodes will affect the inverse manipulation of the element matrix {H  }, which is nearly singular due to the almost same entries [12].Therefore,  = 5 is chosen in the following practical computation.
To investigate the sensitivity to mesh distortion of the presented formulation, numerical results at five test points for uniform and distorted meshes are obtained by the formulation presented in the paper and displayed in Table 1.The coordinates of the five test points are, respectively, A(0.1875, 0.2500), B(0.6875, 0.2500), C(0.3125, 0.7500), D(0.8125, 0.7500), and E(0.5000, 0.5000).As expected, a marked insensitivity to the mesh distortion is observed in the table.Besides, we find that the results inside the subdomains, that is, at points A, B, C, and D, are less affected by the irregularity of the mesh than those at the subdomains boundary point E.
Next, the temperature and heat flux distributions with various power parameter  are, respectively, given in Figures 4 and 5 to illustrate the effect of nonlinear material property.It can be seen that the temperature profile becomes steeper, and the temperature gradient has larger value, as  increases.Besides, the good agreement between numerical results and exact solutions means that the proposed approach can capture the nonlinear effect, even with coarse mesh.
In contrast to the power-type variation of thermal conductivity, another form of thermal conductivity with exponent expression is also investigated in this example.For the case  =  0  (+) , it is easy to get the following analytical solutions: During computation, the same 2 by 2 regular mesh as the one used in the previous test is used to model the computing domain.Simultaneously, we keep  0 = 1 and  = −2 invariant.The results obtained by the HFS-FEM are displayed in Figures 6 and 7, from which we notice that the results of HFS-FEM and exact solutions agree well.At the same time, it is found that the temperature curve shows stronger nonlinearity, as the parameter  becomes larger.Meanwhile, the average values of heat flux component  1 dramatically increase by almost fourteen orders of magnitude, that is, the value increases from −466.9718 to −3.1847 × 10 17 .

Nonlinear Heat Transfer in a Unit Square with a Circular
Cut.To illustrate the advantage of the present approach over the conventional FEM in the aspect of mesh reduction, let us consider a unit square with a centered circular cut.The diameter of the circular hole is taken to be 0.5.The boundary conditions along the outer edges of the square are the same as those in the previous example, and the rim of the hole is assumed to be insulated.The thermal conductivity of the material is assumed as  = −2 + 0.01.
In the computation, two special elements, respectively, including 8 nodes and 16 nodes are compared (see Figure 8) to investigate the change of accuracy of the present algorithm.Figure 9 displays the temperature isoline maps corresponding to the two different elements shown in Figure 8.In the figure, the numerical results of conventional FEM implemented by MATLAB PDE toolbox are also provided to make comparison.Total 4448 triangle finite elements with 2339 nodes are employed to discretize the computing domain in the MATLAB PDE toolbox.It can be observed that the numerical accuracy of the present special elements increases, as the number of nodes of the special element becomes large.Moreover, there is a good agreement between the numerical results of the present special element with 16 nodes and that of MATLAB PDE toolbox.More interestingly, a great effect on mesh reduction is achieved by use of the proposed special elements.Thus, the present hybrid finite element model with special elements can obtain better efficiency than the conventional FEM, when the circular hole exists in the computing domain.

Conclusions
The present study proposes the fundamental-solution-based hybrid finite element formulation to solve the nonlinear heat role.Based on the mathematical definition of them, they are employed to construct the element interior field, which exactly satisfies the governing equation, and simultaneously convert the domain integral into the boundary integral in the hybrid functional.
General Fundamental Solution.For the general case without holes in the 2D infinite plane, the fundamental solution  * (x, x  ) of the problem is required to satisfy the following partial differential equation: )  * (x, x  ) +  (x, x  ) = 0 ∀ x, x  ∈ R 2

Figure 1 :
Figure 1: Sketch of material nonlinearity in unit square domain.

Figure 4 :
Figure 4: Variation of temperature along the line  2 = 0.5 with different power order .

Figure 5 :
Figure 5: Variation of heat flux along the line  2 = 0.5 with different power order .

Figure 6 :
Figure 6: Variation of temperature along the line  2 = 0.5 with different coefficient .

Figure 7 :
Figure 7: Variation of heat flux along the line  2 = 0.5 with different coefficient .

3 )
* (x, x  ) = − 1 2Re [ln ( −   )] .(A.2)In (A.1) and (A.2),  represents the Dirac delta function, Re denotes the real part of the bracketed expression, and  =  1 +  2  and   =   1 +   2  are the field point and source point expressed in the complex space, respectively.The derivative of (A.2) is given by Special Fundamental Solution.Specially, when there is a central circular hole with radius  in the 2D infinite plane, the fundamental solution satisfies both partial differential equation (A.1) and the specific boundary condition along the rim of circular hole.Here, the insulated boundary condition  * / = 0 is assumed, and the corresponding special fundamental solutions can be written as * (,   ) = − 1 2 Re [ln ( −   ) ( 2 −   )   2 −  2    ( −   ) ( 2 −   )] .(A.5)

Table 1 :
Sensitivity to distortion of sub-domain mesh.