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

Authors: Hui Wang; Ming-Yue Han; Fang Yuan; Zhao-Ran Xiao
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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) [3–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 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 [15–17]. 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 Section2, 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(X1,X2). The heat equilibrium equation in the absence of internal heat generation is expressed as [18]
(1)∂q1∂X1+∂q2∂X2=0inΩ,
where Ω denotes the computing domain bounded by a simple closed curve Γ, and the quantity qi is the heat transfer rate along the ith direction governed by the Fourier’s law
(2)qi=-k∂T∂Xi.In (2), T is the temperature change and k is the thermal conductivity of the material. In the study, the thermal conductivity is assumed to be temperature-dependent, that is, k=k(T), and thus the governing equation (1) shows nonlinearity.Besides, the following boundary conditions:(3)T=T-onΓ1q≡qini=q-onΓ2
should be complemented to (1) to obtain a complete mathematic system. In (3), T- and q-, respectively, stand for the given distributions on the boundary Γ=Γ1+Γ2.


### 2.2. Kirchhoff Transformation

Usually, for the case that the material property is dependent on the temperature, that is,k=k(T), 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]
(4)dϕdT=k(T)orϕ=∫k(T)dT,
where ϕ is a new quantity named as Kirchhoff variable, which is related to the sought temperature field T.Then, substituting (4) into (2), one obtains
(5)qi=-k(T)∂T∂Xi=-k(T)∂T∂ϕ∂ϕ∂Xi=-∂ϕ∂Xi.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 ϕ:
(6)∂2ϕ∂X12+∂2ϕ∂X22=0
and the transformed boundary conditions
(7)ϕ=ϕ(T-)onΓ1ψ≡-∂ϕ∂Xini=q-onΓ2.In this paper, two types of material nonlinearities are considered.


#### 2.2.1. Power-Type Thermal Conductivity

Consider the following:(8)k=k0(α+βT)nwith constantsk0, α, β, and n.Integrating (8) by the Kirchhoff transformation (4) gives
(9)ϕ=k0β(n+1)(α+βT)n+1
from which the inverse manipulation produces
(10)T=[β(n+1)ϕ/k0]1/(n+1)-αβ.


#### 2.2.2. Exponent-Type Thermal Conductivity

Consider the following:(11)k=k0e(α+βT)Integrating (11) by the Kirchhoff transformation (4) gives
(12)ϕ=k0βe(α+βT)
from which the inverse manipulation produces
(13)T=ln(βϕ/k0)-αβ.


### 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 elemente, the weak integral functional can be written as follows [12, 13]:
(14)Πme=-12∫Ωe[(∂ϕ∂X1)2+(∂ϕ∂X2)2]dΩ+-∫Γ2eψ-ϕ~dΓ+∫Γeψ(ϕ~-ϕ)dΓ,
where ϕ and ϕ~ stand for independent element interior and boundary fields, respectively.Within the elemente, the assumed element interior field ϕ, also named as intraelement field, is usually approximated by the linear combination of the fundamental solution ϕ*(x,xs), that is,
(15)ϕ(x,xs)={N}{ce},
where
(16){N}={ϕ1*ϕ2*⋯ϕNs*},{ce}={c1c2⋯cNs}T,
and x and xs are field point and source field, respectively. Ns is the number of source points xs located outside the element domain, and ci 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(17)ϕ~(x)={N~}{de},
where {de} and {N~} 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:
(18)Πme=-12∫ΓeψϕdΓ-∫Γ2eψ-ϕ~dΓ+∫Γeψϕ~dΓ
in which only element boundary integrals are included.The substitution of (15) and (17) into (18) finally yields
(19)Πme=-12{ce}T[He]{ce}-{de}T{ge}+{ce}T[Ge]{de},
where
(20)[He]=∫Γe{Q}T{N}dΓ,[Ge]=∫Γe{Q}T{N~}dΓ,{ge}=∫Γ2e{N~}Tψ-dΓ
with
(21){Q}=-∂{N}∂xini=-{n1n2}[∂{N}∂x1∂{N}∂x2].The stationary conditions of the functionalΠme with respect to {ce} and {de} produce
(22)[Ke]{de}={ge},{ce}=[He]-1[Ge]{de}
in which the stiffness equation and the optional relationship between {ce} and {de} are well established.


## 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(X1,X2). The heat equilibrium equation in the absence of internal heat generation is expressed as [18]
(1)∂q1∂X1+∂q2∂X2=0inΩ,
where Ω denotes the computing domain bounded by a simple closed curve Γ, and the quantity qi is the heat transfer rate along the ith direction governed by the Fourier’s law
(2)qi=-k∂T∂Xi.In (2), T is the temperature change and k is the thermal conductivity of the material. In the study, the thermal conductivity is assumed to be temperature-dependent, that is, k=k(T), and thus the governing equation (1) shows nonlinearity.Besides, the following boundary conditions:(3)T=T-onΓ1q≡qini=q-onΓ2
should be complemented to (1) to obtain a complete mathematic system. In (3), T- and q-, respectively, stand for the given distributions on the boundary Γ=Γ1+Γ2.


## 2.2. Kirchhoff Transformation

Usually, for the case that the material property is dependent on the temperature, that is,k=k(T), 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]
(4)dϕdT=k(T)orϕ=∫k(T)dT,
where ϕ is a new quantity named as Kirchhoff variable, which is related to the sought temperature field T.Then, substituting (4) into (2), one obtains
(5)qi=-k(T)∂T∂Xi=-k(T)∂T∂ϕ∂ϕ∂Xi=-∂ϕ∂Xi.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 ϕ:
(6)∂2ϕ∂X12+∂2ϕ∂X22=0
and the transformed boundary conditions
(7)ϕ=ϕ(T-)onΓ1ψ≡-∂ϕ∂Xini=q-onΓ2.In this paper, two types of material nonlinearities are considered.


### 2.2.1. Power-Type Thermal Conductivity

Consider the following:(8)k=k0(α+βT)nwith constantsk0, α, β, and n.Integrating (8) by the Kirchhoff transformation (4) gives
(9)ϕ=k0β(n+1)(α+βT)n+1
from which the inverse manipulation produces
(10)T=[β(n+1)ϕ/k0]1/(n+1)-αβ.


### 2.2.2. Exponent-Type Thermal Conductivity

Consider the following:(11)k=k0e(α+βT)Integrating (11) by the Kirchhoff transformation (4) gives
(12)ϕ=k0βe(α+βT)
from which the inverse manipulation produces
(13)T=ln(βϕ/k0)-αβ.


## 2.2.1. Power-Type Thermal Conductivity

Consider the following:(8)k=k0(α+βT)nwith constantsk0, α, β, and n.Integrating (8) by the Kirchhoff transformation (4) gives
(9)ϕ=k0β(n+1)(α+βT)n+1
from which the inverse manipulation produces
(10)T=[β(n+1)ϕ/k0]1/(n+1)-αβ.


## 2.2.2. Exponent-Type Thermal Conductivity

Consider the following:(11)k=k0e(α+βT)Integrating (11) by the Kirchhoff transformation (4) gives
(12)ϕ=k0βe(α+βT)
from which the inverse manipulation produces
(13)T=ln(βϕ/k0)-αβ.


## 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 elemente, the weak integral functional can be written as follows [12, 13]:
(14)Πme=-12∫Ωe[(∂ϕ∂X1)2+(∂ϕ∂X2)2]dΩ+-∫Γ2eψ-ϕ~dΓ+∫Γeψ(ϕ~-ϕ)dΓ,
where ϕ and ϕ~ stand for independent element interior and boundary fields, respectively.Within the elemente, the assumed element interior field ϕ, also named as intraelement field, is usually approximated by the linear combination of the fundamental solution ϕ*(x,xs), that is,
(15)ϕ(x,xs)={N}{ce},
where
(16){N}={ϕ1*ϕ2*⋯ϕNs*},{ce}={c1c2⋯cNs}T,
and x and xs are field point and source field, respectively. Ns is the number of source points xs located outside the element domain, and ci 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(17)ϕ~(x)={N~}{de},
where {de} and {N~} 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:
(18)Πme=-12∫ΓeψϕdΓ-∫Γ2eψ-ϕ~dΓ+∫Γeψϕ~dΓ
in which only element boundary integrals are included.The substitution of (15) and (17) into (18) finally yields
(19)Πme=-12{ce}T[He]{ce}-{de}T{ge}+{ce}T[Ge]{de},
where
(20)[He]=∫Γe{Q}T{N}dΓ,[Ge]=∫Γe{Q}T{N~}dΓ,{ge}=∫Γ2e{N~}Tψ-dΓ
with
(21){Q}=-∂{N}∂xini=-{n1n2}[∂{N}∂x1∂{N}∂x2].The stationary conditions of the functionalΠme with respect to {ce} and {de} produce
(22)[Ke]{de}={ge},{ce}=[He]-1[Ge]{de}
in which the stiffness equation and the optional relationship between {ce} and {de} 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 Section3.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 Figure1.Figure 1
Sketch of material nonlinearity in unit square domain.For the case of power-type variation of the thermal conductivity with temperature, that is,k=k0(α+βT)n, an analytical solution of the form
(23)T=({(α+300β)n+12}{[(α+400β)n+1-(α+300β)n+1]X1+++(α+300β)n+1}1/(n+1)-α)(β)-1
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 k0=1, n=1, α=-2, and β=0.01, which corresponds to k=-2+0.01T and T=1003X1+1+200. Figure 3 displays the variation of the error defined by
(24)Arerr(u)=∑i=1N(uHFS-FEM-uexact)i2∑i=1N(uexact)i2
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 r 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 O(lnr), and O(r-2), respectively, in two dimensions so that inaccuracies in the evaluation of element boundary integrals (near-singular problems) would be caused when r 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 {He}, which is nearly singular due to the almost same entries [12]. Therefore, γ=5 is chosen in the following practical computation.Figure 2
Investigation of sensitivity to mesh distortion.Figure 3
Effect of location of source point to numerical accuracy.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 Table1. 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.Table 1
Sensitivity to distortion of sub-domain mesh.
Regular mesh
Irregular mesh
(|Tregular-Tirregular|/|Tregular|)×100%
A
324.9997
324.9990
0.000
B
374.9998
375.0006
0.000
C
339.1944
339.1934
0.000
D
385.4052
385.4057
0.000
E
358.1139
358.1139
0.000
Regular mesh
Irregular mesh
(
|
q
1regular
-
q
1irregular
|
/
|
q
1regular
|
)
×
100
%
A
−150.0053
−150.0102
0.003
B
−150.0053
−150.0034
0.001
C
−150.0053
−150.0033
0.001
D
−150.0053
−150.0103
0.003
E
−149.9287
−149.8496
0.053Next, the temperature and heat flux distributions with various power parametern 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 n 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.Figure 4
Variation of temperature along the lineX2=0.5 with different power order n.Figure 5
Variation of heat flux along the lineX2=0.5 with different power order n.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 casek=k0e(α+βT), it is easy to get the following analytical solutions:
(25)T=ln[(eα+400β-eα+300β)X1+eα+300β]-αβ,q1=-k0β[e(α+400β)-e(α+300β)].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 keepk0=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 q1 dramatically increase by almost fourteen orders of magnitude, that is, the value increases from −466.9718 to -3.1847×1017.Figure 6
Variation of temperature along the lineX2=0.5 with different coefficient β.Figure 7
Variation of heat flux along the lineX2=0.5 with different coefficient β.


### 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 ask=-2+0.01T.In the computation, two special elements, respectively, including 8 nodes and 16 nodes are compared (see Figure8) 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.Two special elements with different number of nodes.
(a)
8-node special circular hole element
(b)
16-node special circular hole elementFigure 9
Temperature distribution in the computing domain.


## 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 Figure1.Figure 1
Sketch of material nonlinearity in unit square domain.For the case of power-type variation of the thermal conductivity with temperature, that is,k=k0(α+βT)n, an analytical solution of the form
(23)T=({(α+300β)n+12}{[(α+400β)n+1-(α+300β)n+1]X1+++(α+300β)n+1}1/(n+1)-α)(β)-1
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 k0=1, n=1, α=-2, and β=0.01, which corresponds to k=-2+0.01T and T=1003X1+1+200. Figure 3 displays the variation of the error defined by
(24)Arerr(u)=∑i=1N(uHFS-FEM-uexact)i2∑i=1N(uexact)i2
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 r 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 O(lnr), and O(r-2), respectively, in two dimensions so that inaccuracies in the evaluation of element boundary integrals (near-singular problems) would be caused when r 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 {He}, which is nearly singular due to the almost same entries [12]. Therefore, γ=5 is chosen in the following practical computation.Figure 2
Investigation of sensitivity to mesh distortion.Figure 3
Effect of location of source point to numerical accuracy.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 Table1. 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.Table 1
Sensitivity to distortion of sub-domain mesh.
Regular mesh
Irregular mesh
(|Tregular-Tirregular|/|Tregular|)×100%
A
324.9997
324.9990
0.000
B
374.9998
375.0006
0.000
C
339.1944
339.1934
0.000
D
385.4052
385.4057
0.000
E
358.1139
358.1139
0.000
Regular mesh
Irregular mesh
(
|
q
1regular
-
q
1irregular
|
/
|
q
1regular
|
)
×
100
%
A
−150.0053
−150.0102
0.003
B
−150.0053
−150.0034
0.001
C
−150.0053
−150.0033
0.001
D
−150.0053
−150.0103
0.003
E
−149.9287
−149.8496
0.053Next, the temperature and heat flux distributions with various power parametern 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 n 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.Figure 4
Variation of temperature along the lineX2=0.5 with different power order n.Figure 5
Variation of heat flux along the lineX2=0.5 with different power order n.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 casek=k0e(α+βT), it is easy to get the following analytical solutions:
(25)T=ln[(eα+400β-eα+300β)X1+eα+300β]-αβ,q1=-k0β[e(α+400β)-e(α+300β)].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 keepk0=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 q1 dramatically increase by almost fourteen orders of magnitude, that is, the value increases from −466.9718 to -3.1847×1017.Figure 6
Variation of temperature along the lineX2=0.5 with different coefficient β.Figure 7
Variation of heat flux along the lineX2=0.5 with different coefficient β.


## 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 ask=-2+0.01T.In the computation, two special elements, respectively, including 8 nodes and 16 nodes are compared (see Figure8) 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.Two special elements with different number of nodes.
(a)
8-node special circular hole element
(b)
16-node special circular hole elementFigure 9
Temperature distribution 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.



---

*Source: 695457-2013-09-05.xml*