Abstract

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.

1. 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 temperature-dependent 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) [35], 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 fundamental-solution-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 [1517]. Generally, as one of the domain-type 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 temperature-dependent 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.

2. 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 . 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 .

2.2. 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] 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.

2.2.1. Power-Type Thermal Conductivity

Consider the following:

with constants , , , and .

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

2.2.2. 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 , that is, where and and are field point and source field, respectively. is the number of source points 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 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 and produce in which the stiffness equation and the optional relationship between and are well established.

3. 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.

3.1. 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, , 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 , , , and , which corresponds to and . 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 , and , 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 , which is nearly singular due to the almost same entries [12]. Therefore, 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 , 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 and 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 dramatically increase by almost fourteen orders of magnitude, that is, the value increases from −466.9718 to .

3.2. 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 .

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.

4. Conclusions

The present study proposes the fundamental-solution-based hybrid finite element formulation to solve the nonlinear heat conduction problems with temperature-dependent thermal conductivity. It is a direct extension of the HFS-FEM for linear heat conduction problems by combining the Kirchhoff transformation and it is easy to be numerically implemented. The numerical results show that the present hybrid finite element formulation is very effective and convenient to solve the nonlinear heat conduction problem. It is shown that, for the analysis of nonlinear heat conduction problems, especially those related to the circular hole, the improvements of both the efficiency of mesh discretization and the numerical accuracy of the present hybrid finite element model are significant, for the sake of the use of general or special fundamental solutions.

Appendix

Fundamental Solutions

In the formulation of the present hybrid finite element, the fundamental solution of the problem plays an important 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 of the problem is required to satisfy the following partial differential equation: which produces

In (A.1) and (A.2), represents the Dirac delta function, Re denotes the real part of the bracketed expression, and and 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 is assumed, and the corresponding special fundamental solutions can be written as whose derivative is

Acknowledgments

This research work was partially supported by the Natural Science Foundation of China under Grant nos. 11102059 and 51178165. This financial support is gratefully acknowledged. Also, the authors would like to thank the support of the basic research programs from Henan Province Office of Education and Henan University of Technology.