Analysis of 2D heat conduction in nonlinear functionally graded materials using a local semi-analytical meshless method

1.   Introduction

The functionally graded material is a composite with continuously varying microstructure and mechanical properties and is often applied to reduce the thermal stress or residual stress in dissimilar material joints. With this regard, it is very important to analyze the heat transfer characteristics of functionally graded materials. Generally, it is not easy work to derive an exact expression of temperature distribution in nonlinear functionally graded materials. Therefore, the numerical simulation has become one of the main methods for addressing such problems.

Over the past few decades, various numerical methods have been developed for the analysis of heat transfer in functionally graded materials. The main methods include the finite element method (FEM) ^[1,2,3,4], the boundary element method (BEM) ^[5,6,7], and the finite difference method (FDM) ^[8,9,10]. Among them, the FEM occupies dominant position due to its mature theory system and good stability. These traditional mesh-type methods have their advantages, but at the same time they have many deficiencies, such as the complexity of preprocessing, especially for complex structures.

In order to reduce the time-consuming task of grid partition, numerous meshless/meshfree methods and coupled finite element-meshfree methods ^[11,12] have been proposed in recent years. At present, the common meshless methods include the meshless local Petrov-Galerkin method (MLPG) ^[13,14,15], the generalized finite difference method (GFDM) ^{[16,17,18,19]}, the singular boundary method (SBM) ^{[20,21,22,23,24]}, the boundary knot method (BKM) ^{[25,26,27,28]} and the method of fundamental solution (MFS) ^{[29,30,31,32]}, etc. The successful application of these approaches shows their development prospect. Meanwhile, these methods also have some limitations. For example, as the semi-analytical and boundary-type meshless schemes, the BKM and MFS encounter the challenge of simulating large-scale and/or high-dimensional problems. With these reasons in mind, the localized method of fundamental solution (LMFS) ^[33] and the localized boundary knot method (LBKM) ^[34] have been proposed and successfully applied to large-scale acoustics ^[35,36], elastic mechanics ^[37], inverse problems ^[38], steady-state heat transfers ^[39,40] and other mechanics problems ^[41,42,43].

The LBKM is a domain-type meshfree approach using the non-singular general solution as the basis function. In fact, it can be regarded as a localized version of traditional boundary-type BKM. When this method was first proposed, the artificial boundaries and boundary nodes on it are required for every local subdomain. Later, the LBKM is modified into the local knot method (LKM) without designing the artificial boundaries. Like the local radial basis function collocation method and the element-free Galerkin method ^[44,45,46], the LKM is a typical localized domain-type meshless method. Unlike them, the LKM is a semi-analytical and strong-form approach, in which the nonsingular general solution is taken as the basis function. Compared with the local radial basis function collocation method, the LKM avoids the issue of selecting shape parameters, and has a higher numerical accuracy. Recently, the LKM has been successfully applied to solving convection-diffusion-reaction equations ^[47] and acoustic problems ^[48], and has demonstrated the advantages of truly meshless, high-accuracy, and large-scale calculation.

In view of its above merits, this study makes a first attempt to employ the LKM in conjunction with the Kirchhoff transformation to solve the heat conduction problems of two-dimensional non-linear functionally graded materials. In the calculation, the Kirchhoff transformation is applied to transform the nonlinear heat conduction equation into a linear equation. The transformed equation is a modified Helmholtz equation whose non-singular general solution is available. The LKM can be directly used to approximate the solution of the modified Helmholtz equation. The solution of the original problem can be acquired by using the inverse Kirchhoff transformation.

The general outline of this paper is as follows. The heat conduction problem of functionally graded nonlinear material is introduced in Section 2. Section 3 describes the procedures of solving the anisotropic modified Helmholtz equation via the LKM. In Section 4, three benchmark numerical examples are presented to demonstrate the feasibility and accuracy of the proposed methodology. Finally, some conclusions are drawn in Section 5.

2.   Problem statement

Consider the steady-state heat conduction problem in a two-dimensional nonlinear functionally graded material with the domain $\mathit{\Omega } \subset {R^\mathit{2}}$ and boundary $\mathit{\Gamma }$, the governing equation can be represented as

with the following boundary conditions:

where $T({\boldsymbol{x}})$ represents the temperature at the point ${\boldsymbol{x}}$, $K = {\left\{ {{K_{ij}}({\boldsymbol{x}}, T)} \right\}_{1 \leqslant i, j \leqslant 2}}$ the heat conduction coefficient matrix satisfying positive definite $({\Delta _K} = \det (K) = {K_{12}}{K_{21}} - {K_{11}}{K_{22}} < 0)$ and symmetry $({K_{21}} = {K_{12}})$ conditions, ${{\boldsymbol{n}}_{\boldsymbol{x}}}$ the unit outer normal vector at the boundary point ${\boldsymbol{x}}$, ${h_e}$ the conductivity factor, ${T_f}$ the ambient temperature.${\mathit{\Gamma }_D}$, ${\mathit{\Gamma }_N}$ and ${\mathit{\Gamma }_R}$ denote boundary segments according to Dirichlet, Neumann and Robin boundaries, and $\mathit{\Gamma } = {\mathit{\Gamma }_D} + {\mathit{\Gamma }_N} + {\mathit{\Gamma }_R}$.

In this study, we focus on the exponentially functionally graded material whose heat conductivity matrix could be represented as follows

where $a(T) > 0$, ${\beta _1}$ and ${\beta _2}$ are the material parameters, ${\boldsymbol{\bar K}} = {\left\{ {{{\bar K}_{ij}}} \right\}_{1 \leqslant i, j \leqslant 2}}$ is a symmetric positive definite matrix.

By using the Kirchhoff transformation

Equations (1)–(4) can be converted to

where ${\mathit{\Phi }_T}({\boldsymbol{x}}) = \varphi (T({\boldsymbol{x}}))$.

Let ${\mathit{\Phi }_T} = \mathit{\Psi }{e^{ - \sum\nolimits_{i = 1}^2 {{\beta _i}({x_i} + {s_i})} }}$, Eq (7) can be simplified to the following formula

where $\lambda = \sqrt {\sum\nolimits_{i = 1}^2 {\sum\nolimits_{j = 1}^2 {{\beta _i}{{\bar K}_{ij}}{\beta _j}} } }$, $s = ({s_1}, {s_2})$ is source point. Since ${e^{\sum\nolimits_{i = 1}^2 {{\beta _i}({x_i} + {s_i})} }} > 0$, the above formula can degenerate into an anisotropic modified Helmholtz equation. Employing the following transformation

with ${\Delta _{\bar K}} = \det (\bar K) = {\bar K_{11}}{\bar K_{22}} - {\bar K_{12}}{\bar K_{21}} > 0$, an isotropic modified Helmholtz equation will be obtained, namely,

The non-singular general solution of Eq (13) is available ^[49], and thus the non-singular general solution of Eq (11) can be derived by using the inverse transformation of Eq (12),

where $R = \sqrt {\sum\nolimits_{i = 1}^2 {\sum\nolimits_{j = 1}^2 {{r_i}\bar K_{ij}^{ - 1}{r_j}} } }$, ${r_1} = {x_1} - {s_1}$, ${r_2} = {x_2} - {s_2}$, ${\boldsymbol{x}} = ({x_1}, {x_2})$ and $s = ({s_1}, {s_2})$ are collocation node and source point, and ${I_0}$ indicates the zero-order modified Bessel function of first kind. Using the transform ${\mathit{\Phi }_T} = \mathit{\Psi }{e^{ - \sum\nolimits_{i = 1}^2 {{\beta _i}({x_i} + {s_i})} }}$ again, the non-singular general solution of Eq (7) is yielded, and written as

3.   Local knot method

This section will introduce the implementation of the LKM in solving Eq (7), with the help of the non-singular general solution (15). According to the basic theory of the LKM, $\left({{{\boldsymbol{x}}_i}} \right)_{i = 1}^N$ are discretized inside the domain $\mathit{\Omega }$ and along its boundary $\mathit{\Gamma }$, here $N = {n_p} + {n_{b1}} + {n_{b2}} + {n_{b3}}$, the numbers of boundary nodes satisfying the Dirichlet, Neumann and Robin boundary conditions are ${n_{b1}}$, ${n_{b2}}$ and ${n_{b3}}$, respectively. ${n_p}$ represents the number of interior nodes. Figure 1 shows the schematic diagram of the LKM, which indicates the distribution of the central node ${{\boldsymbol{x}}^{(0)}}$ and its m supporting points (${{\boldsymbol{x}}^{(p)}}, p = 1, 2, ..., m$) in a corresponding supporting domain ${\mathit{\Omega }_\mathit{s}}$. For the $m + 1$ local nodes inside the supporting domain ${\mathit{\Omega }_\mathit{s}}$, we can express the unknown at the center node using a linear combination of the non-singular general solutions

or for brevity

where $({\alpha _q})_{q = 0}^m$ are the unknown coefficients, $G({R_{pq}})$ is the non-singular general solution given in Eq (15).

Based on the essential ideas of moving least squares (MLS) approximation, in each local subdomain, the residual function can be defined as follows

where ${\omega ^{(p)}}$ is the weight function, many kinds of weight functions can be adopted. We use the spline weight function ^[50] in the present study.

where ${d_p}$ is the distance between the central node ${\rm{ }}{{\boldsymbol{x}}^{(0)}}$ and the pth supporting node ${\rm{ }}{{\boldsymbol{x}}^{(p)}}$, ${d_{{\rm{max}}}}$ is denotes the maximum value of distances between ${\rm{ }}{{\boldsymbol{x}}^{(0)}}$ and m supporting nodes, i.e., ${d_{\max }} = \mathop {\max }\limits_{p = 1, 2, ..., m} ({d_p})$.

On the basis of the MLS theory, the undetermined coefficients ${\boldsymbol{\alpha }} = {({\alpha _0}, {\alpha _1}, ..., {\alpha _m})^T}$ could be decided by minimizing the functional $B(\mathit{\Phi })$, namely,

From Eq (20), the following linear system can be consisted

where

The unknown vector ${\boldsymbol{b}}$ in Eq (22) could be recast as

In accordance with Eqs (21)–(23), the coefficient vector ${\boldsymbol{\alpha }} = {({\alpha _0}, {\alpha _1}, ..., {\alpha _m})^T}$ can be repressed as

Substituting Eq (24) into Eq (17) as $p = 0,$ the temperature at central node ${\rm{ }}{{\boldsymbol{x}}^{(0)}}$ is expressed as

or

in which $({f^{(q)}})_{q = 0}^m = {{\boldsymbol{G}}^{(0)}}{{\boldsymbol{A}}^{ - 1}}{\boldsymbol{B}}.$ Follow the above process, the unknowns $\mathit{\Phi }({{\boldsymbol{x}}_i})$ at all nodes $\left({{{\boldsymbol{x}}_i}} \right)_{i = 1}^N$ can be expressed in the form of Eq (26).

In addition, the normal heat flux at the boundary node can be calculated by

in which

or for brevity

where $({e^{(l)}})_{l = 0}^m = Q({{\boldsymbol{x}}^{(0)}}, {{\boldsymbol{x}}^{(q)}}){{\boldsymbol{A}}^{ - 1}}{\boldsymbol{B}}.$ If the central node ${{\boldsymbol{x}}^{(0)}}$ is a boundary node with Robin boundary condition, we have

Substituting Eq (29) into Eq (30) yields

For the interior nodes, the temperature distribution should satisfy Eq (26), namely,

where the subscript p of $f_p^{(q)}$ is a coefficient used to differentiate the different internal nodes. For the boundary nodes satisfying Dirichlet boundary conditions, the following equation should be satisfied

For the boundary nodes with Neumann conditions, we have

For the boundary nodes with Robin conditions, we have

In Eqs (34) and (35), the subscript "p" (indicates the node number) is used to distinguish the different boundary nodes.

By using given boundary data and combining Eqs (32)–(35), the following sparse system of linear equations is obtained

in which ${{\boldsymbol{A}}_{N \times N}}$ is the matrix of coefficient, ${{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_T} = {\left[{{\mathit{\Phi }^{(1)}}, \; {\mathit{\Phi }^{(2)}}, ..., \; {\mathit{\Phi }^{(N)}}} \right]^T}$ represents the unascertained vector of variables at all nodes, and ${\mathit{\boldsymbol{b}}_{N \times 1}}$ denotes the known vector. After obtaining the solution ${{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_T}$ of Eq (36), the solution of the original problem (Eq (1)) can be acquired by using the inverse Kirchhoff transformation

To summarize and to make clearer the procedure of the developed method for solving the heat conduction problem in nonlinear functionally graded material, a computational flow chart is given in Figure 2.

4.   Numerical examples

In this section, three typical numerical examples are provided to verify the applicability and accuracy of the proposed method for solving the heat conduction problem of two-dimensional nonlinear functionally graded materials. In order to estimate the accuracy, the measured errors are defined as

where W is the number of all test points, $u_n^k$ and $u_e^k$ are numerical and analytical solutions at test points, respectively. Remarkably, all the calculations were done on a laptop (Intel® CoreTMi5-9300H CPU at 2.40 GHz, 8G RAM) by using the MATLAB program.

4.1.   Example 1

The typical heat transfer problem of exponential heterogeneous gradient materials is considered. In a high-temperature environment, the coefficient of heat conduction is defined as $a(T) = {e^T}$. We can acquire ${\mathit{\Phi }_\mathit{T}}\mathit{ = }{\mathit{e}^\mathit{T}}$, $T = {\varphi ^{ - 1}}({\mathit{\Phi }_T}) = \ln ({\mathit{\Phi }_T})$ by the Kirchhoff transformation. Here, we consider a square domain $\mathit{\Omega = }\left[{- 1, 1} \right] \times \left[{- 1, 1} \right]$ of orthotropic material with $\bar K = \left(\begin{gathered} 2{\rm{ \; \; 0}} \\ 0{\rm{ \; \; 1}} \\ \end{gathered} \right)$ and ${\beta _1}{\rm{ = }}0$, ${\beta _2}{\rm{ = 1}}$. The lower and right boundaries of the structure respectively satisfy the Dirichlet and Neumann boundary conditions, and the other boundaries satisfy the Robin boundary condition. The analytical solution is

where $Tx = {{{x_1}} / {\sqrt 2 }} - 1$, $Ty = {x_2}$, and $Tr = \sqrt {T{x^2} + T{y^2}}$.

For the sake of investigating the influence of the total number of nodes on the calculation results, 2300 regular and irregular nodes are used in the calculation (as shown in Figure 3). Irregular nodes are derived by jiggling the regular nodes, i.e., assigning a perturbation on regular nodes along the x and y directions. The temperature distribution in the domain can be calculated with different numbers of supporting node ($10 \leqslant m \leqslant 40$). Figure 4 shows the error curves as the number of supporting node increases, under regular/irregular nodal distributions. As can be seen from Figure 4, although the errors show fluctuation for some specific numbers of supporting node, it can still be seen that the global errors and the maximum errors generally have a decreasing trend with the increase of supporting nodes. In this example, $m = 35$ is a relatively optimal value to achieve the highest accuracy. This indicates that the number of supporting nodes has a certain influence on the calculation results, but a relatively ideal numerical result can be obtained in a larger range.

In fact, the number of supporting nodes ($m$) is related to the node spacing ($\Delta h = \mathop {\max }\limits_{1 \leqslant i \leqslant N} \mathop {\min }\limits_{1 \leqslant j \leqslant N} \left| {{x_i} - {x_j}} \right|$). In order to numerically analyze the relationship between $m$ and $\Delta h$, we investigate the global errors under different values of $m$ and $\Delta h$. Figure 5 displays the 3D error surface and the 2D error plane. The red dots in Figure 5(b) indicate the positions of ${\rm{Global \;error}} = {10^ - }^3$. Using the curve fitting method, we can obtain an estimation formula that ensures the global error is almost equal to 10^-3, $m = \left\lceil {0.02297 \cdot \Delta {h^{({\rm{ - }}1.885)}}} \right\rceil$, where $\left\lceil \cdot \right\rceil$ denotes the round up operation. Furthermore, the positions corresponding to the relatively optimal value of $m$ are plotted by white squares. In the same way, the relatively optimal value of $m$ can be formulated by

According to Eq (41), we can easily estimate the number of supported nodes according to the node spacing. This example uses 2300 nodes, and thus $\Delta h$ is about 0.038. It can be observed from Figure 5(b) that ${\rm{Global \;error}} \leqslant {10^ - }^3$ when $10 \leqslant m \leqslant 40$.

Figure 6 shows the variations of errors with the increase of the total number of nodes under the regular and irregular nodes, where $m = 30$. We can observe that both the maximum errors and global errors are less than ${\rm{5}}{\rm{.621}} \times {10^{ - 2}}$, and both decrease with the augment of the total number of nodes, indicating that the LKM has the good accuracy and convergence in dealing with the heat transfer problem of nonlinear functionally graded materials. It can also be obtained from Figure 6 that the calculation precision with regular node distribution is higher than that with irregular node distribution, which verifies that the proposed method, as a local meshless method, can simulate the heat conduction problem with arbitrary node distribution, and the calculation accuracy with regular node distribution is slightly better.

In the case of regular nodes, Figure 7 displays the distributions of the exact solution and absolute error in the calculation domain. From Figure 7, we can clearly observe that the LKM can reach a higher precision, and the absolute error is less than $1.104 \times {10^{ - 3}}$.

To validate the accuracy and condition number of the developed method, we compared the LKM with the BKM. Table 1 lists the errors and the condition numbers of the above two methods with increasing number of boundary nodes. It is worth noting that the BKM is a boundary meshless method, while the LKM is a local meshless method. In the BKM, the number of boundary nodes (NB) is set to the same value as the LKM to ensure the fairness of comparison. From Table 1, it is observed that the numerical veracity of the LKM is marginally better than that of the BKM for various values of NB. Moreover, the condition number of the LKM is obviously less than that of the BKM. It could be carefully concluded that the presented LKM is accurate and stable for solving two-dimensional nonlinear heat conduction problems.

4.2.   Example 2

The second example considers a heat transfer problem on a square functionally graded material plate with four circular holes of the same size, as shown in Figure 8(a). The temperatures are assumed to be known on the boundaries a and b, the heat fluxes are given on the boundaries c and d, and the boundaries of four circular holes are adiabatic. The relevant parameters are set as: $L = 1.0$, ${L_1} = 0.3$, ${L_2} = 0.2$, $r = 0.1$. In fact, the thermal conductivity can be expressed as a linear relationship with temperature.

The exact solution is indicated as

Using the Kirchhoff transformation, we can get

where $\mu = {1 / 4},$ $Tx = \frac{{{x_1}\sqrt {{\Delta _{\bar K}}} }}{{{{\bar K}_{11}}}}, {\rm{ }}Ty = - \frac{{{x_1}{{\bar K}_{12}}}}{{{{\bar K}_{11}}}} + {x_2},$ $\tau = \sqrt {{{\bar K}_{11}}{{(\frac{{\sqrt {{\Delta _{\bar K}}} - {{\bar K}_{12}}}}{{{{\bar K}_{11}}}})}^2} + 2{{\bar K}_{12}}(\frac{{\sqrt {{\Delta _{\bar K}}} - {{\bar K}_{12}}}}{{{{\bar K}_{11}}}}) + {{\bar K}_{22}}}.$$\bar K = \left(\begin{gathered} 1{\rm{ 0}}{\rm{.25}} \\ {\rm{0}}{\rm{.25 3}} \\ \end{gathered} \right),$ ${\beta _1} = 0.1,$ ${\beta _2} = 0.8$.

Under $N = 1226$ (as shown in Figure 8(b)) and $m = 60$, Figure 9 illustrates the comparison of numerical results on the computational domain. It can be found from the figure that the numerical results obtained by the LKM are extremely in agreement with the analytical solutions. Furthermore, the absolute errors are less than ${\rm{8}}{\rm{.306}} \times {10^{ - 5}}$, and the maximum error appears at the boundary d. It can be seen from Table 2, both global errors and maximum errors decrease with increasing number of supporting node, showing a convergence trend. As can be expected, the CPU times gradually increase with increasing number of supporting nodes, but not too much even for the relatively larger value $m = 60$.

Finally, we set $m = 60$, Figure 10 depicts the absolute errors of the LKM under different numbers ($N = 960, {\rm{ }}1890, {\rm{ }}2620$) of total nodes. Noted that the error becomes smaller and smaller as the number of total nodes increases, indicating a convergence trend. The above numerical experiment with complicated geometries and mixed boundary conditions confirms the capacity, accuracy, and convergence of the developed methodology in solving the 2D heat conduction in nonlinear functionally graded material.

4.3.   Example 3

In the last example, the LKM is applied to a nonlinear heat conduction problem in an irregular domain. Figure 11 shows the geometry model of the problem. In this case, we use the same analytical solution and parameters as for example 2.

To investigate the influence of the number of supporting nodes on the computational accuracy, we set N to be 1228 (irregular nodes) and plot the relative errors of the temperatures at all points in Figure 12, when m is equal to 30, 40, and 50. It is noticed that the error gradually decreases with increasing number of supporting nodes.

To investigate the convergence of the LKM for nonlinear heat conduction problem in an irregular structure, Figure 13 shows the variations of global errors and maximum errors as the increase of node number. Numerical results in the figure demonstrate the good convergence properties of the proposed scheme. Furthermore, even if the number of nodes exceeds 100000, the numerical simulation work is still performed on a regular laptop. A great number of numerical tests indicate that the proposed methodology in this paper is accurate, stable and convergent, and is suitable for large-scale simulations in an arbitrary domain.

5.   Conclusions

This paper firstly presented a novel local semi-analytical meshless method, the local knot method (LKM) in conjunction with the Kirchhoff transformation, to numerically simulate heat conduction problems of nonlinear functionally graded materials. As a local semi-analytical meshless algorithm, the LKM uses the non-singular general solution as the basis function. Compared with the traditional BKM, this method not only retains the simplicity and high accuracy, but also avoids an ill-conditioned system of equations, and is more appropriate for large-scale simulations associated with complex structures.

Numerical experiments including simply/multi-connected and regular/irregular domains shown that the proposed method is accurate, stable, and convergent. Furthermore, the resultant matrix is well-conditioned and has a smaller condition number than the traditional BKM. Compared with the existing methods for solving the heat conduction problem of nonlinear functionally graded materials, the present scheme could be regarded as a highly competitive one.

Acknowledgments

The work described in this paper was supported by the National Natural Science Foundation of China (No.11802151), the Natural Science Foundation of Shandong Province of China (No. ZR2019BA008), and the China Postdoctoral Science Foundation (No. 2019M652315).

Conflict of interest

The authors declare that they have no conflicts of interest to report regarding the present study.